On occurrence of point singularities in orthorhombic media

Degeneracies of the slowness surfaces of shear (and compressional) waves in low‐symmetry anisotropic media (such as orthorhombic), known as point singularities, pose difficulties during modelling and inversion, but can be potentially used in the latter as model parameter constraints. I analyse the quantity and spatial arrangement of point singularities in orthorhombic media, as well as their relation to the overall strength of velocity anisotropy. A classification scheme based on the number and spatial distribution of singularity directions is proposed. In normal orthorhombic models (where the principal shear moduli are smaller than the principal compressional moduli), point singularities can only be arranged in three distinct patterns, and media with the theoretical minimum (0) and maximum (16) number of singularities are not possible. In orthorhombic models resulting from embedding vertical fractures in transversely isotropic background, only two singularity distributions are possible, in contrast to what was previously thought. Although the total number of singularities is independent of the overall anisotropy strength, for general (non‐normal) orthorhombic models, different spatial distributions of singularities become more probable with increasing magnitude of anisotropy.


I N T R O D U C T I O N
Seismic wave propagation in anisotropic media is characterized by many interesting and challenging features, among which are so-called singularity directions (also degeneracies, acoustic axes)-points in the phase space where the slowness surfaces of two or more wave modes come in contact. Three types of singularities exist in anisotropic media of different symmetry: tangential ('kiss') and line (intersection) singularities are mostly encountered in hexagonal (transversely isotropic) systems, and conical (or point) singularities are often present in lower symmetry (orthorhombic, monoclinic, and triclinic) systems (Crampin 1991). Although all types of singularities pose certain challenges in seismic modelling (Vavryčuk 2001;Grechka 2017), only conical points have significant influence on the geometrical structure of the * E-mail: yus.ivanov@gmail.com wavefronts and on the polarization in their vicinities (Musgrave 1985;Vavryčuk 2003a;Grechka 2015).
Due to increasing interest in microseismic monitoring, in which a recorded wavefield from small-scale earthquakes contains direct shear waves of high energy, and the medium symmetry is often lower than hexagonal (e.g. orthorhombic or monoclinic as that of fractured shale), the conical singularities have to be treated with care in both modelling and inversion (Grechka and Yaskevich 2014;Grechka 2015). On the other hand, they can potentially be used to constrain model parameters in inversion (Vavryčuk 2013). The goals of the present work are to study (a) the conditions of occurrence of conical singularities in orthorhombic (ORT) media, (b) their spatial distribution, (c) potential dependence of their number and arrangement on the anisotropy strength and (d) to provide relevant numerical examples that illustrate the theoretical findings. This paper is an extended and corrected version of Ivanov and Stovas (2018 This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
In the text below, I alternate between common four-index notation and two-index Voigt notation for the components of the stiffness tensor. I use lower case indices to describe the former and capital indices for the latter, and no summation over repeated indices is assumed.

T H E O R Y
In orthorhombic media, the singularity directions form a pattern that is symmetric about the coordinate planes, and, taking into account the symmetry about the origin, the maximum allowed number of real-valued singularity directions is 16 (Musgrave 1985;Ivanov and Stovas 2017). In fact, that is also the maximum number for any symmetry lower than orthorhombic (Khatkevich 1977;Vavryčuk 2005;Grechka 2015). Theoretically, there also exist orthorhombic models with no singularity directions (Alshits and Lothe 1979). Physical experiments, however, have not identified any real materials with such properties.
Studying the condition of multiple eigenvalues of the Christoffel tensor, Boulanger and Hayes (1998) present equations needed to calculate all singularity directions in an orthorhombic medium. I rewrite their equations using the index notation. In order to calculate the singularity directions n = (n 1 , n 2 , n 3 ) within the [x i , x j ] plane, given n k = 0, one has to solve the following equation for either n i , n j or n i /n j : where i = j = k = i. Assuming that all stiffness coefficients are different, equation (1) results in 0, 2 or 4 real-valued singularity directions. Equal values of any two principal stiffness coefficients (c LL = c MM in Voigt notation) result in multiple roots of equation (1) meaning that the singularity occurs along one of the coordinate axes.
Equation (1) is bi-quadratic in the components of n, and the sign of its discriminant ( ) determines the nature of its roots. If > 0, all four roots are either real valued or complex valued. Writing this condition for the three symmetry planes, one arrives at the following system of inequalities (in Voigt notation): (2) It is evident upon closer inspection of equations (2) that unless the shear stiffness moduli (c LL , L = 4, 5, 6) are larger than the compressional moduli (c LL , L = 1, 2, 3), the system of equations has no solutions. Hence, in orthorhombic materials with the allowed number of singularity directions within the symmetry planes is between 2 and 12. I shall call any orthorhombic medium satisfying condition 3 'normal'. Let us now examine the aforementioned condition from the geometrical standpoint. I consider an arbitrary orthorhombic medium whose symmetry planes are aligned with the coordinate frame. Taking advantage of the mirror symmetry, I can focus on one of the octants only, say, defined by x i > 0. In each of the three quadrants of this octant, there can be zero, one or two intersections (i.e. singularity directions) of the phase surfaces of three wave modes as has been discussed above. Assuming that all principal stiffness moduli are different, and that the medium is normal, there are six possible relations among the principal shear stiffness moduli: c 55 ≶ c 66 along x 1 , c 44 ≶ c 66 along x 2 and c 44 ≶ c 55 along x 3 (two impossible cases, e.g. c 55 > c 44 , c 44 > c 66 , c 66 > c 55 , are excluded). The intersection of the coordinate planes with the phase surfaces of the shear waves produces two smooth curves in each plane. One of the curves corresponds to a shear wave that has the same velocity along the coordinate axes (e.g. √ c 55 along x 1 and x 2 ) and is commonly referred to as an SV-wave in that symmetry plane following the limited equivalence between orthorhombic and transversely isotropic (TI) media. The other curve corresponds to the shear wave that has different velocities along the coordinate axes, also known as the SH-wave in TI media. It is then impossible for these six line segments (two smooth curves in each of the three planes) to have an even number of intersections in the considered octant ( Fig. 1). The out-of-plane singularity directions are obtained explicitly as where and the factor λ is determined from the condition |n| = 1. Equation (4) results in either four or zero real singularity directions (one in each octant in the upper half-space x 3 > 0). To investigate the dependence of the distribution of the singularity directions upon the anisotropy strength, I propose a classification of orthorhombic media based on the number of singularity directions, as follows from equations (1). A similar, although more general and slightly less intuitive classification is proposed in Musgrave (1985). I define 10 classes (Table 1) based on the number of the singularity directions within the symmetry planes. Numbering of the symmetry planes is arbitrary due to their equivalence, and possible permutations of the symmetry planes are indicated in the table. Additionally, each class can have zero or four out-of-plane singularities. Based on the aforementioned observations, for normal orthorhombic media, one of the symmetry planes must contain just one singularity direction, thus allowing only three classes, namely, II, V and IX, to be present.
To quantify the overall elastic anisotropy strength, I use a recently introduced parameter known as the universal elastic anisotropy index (Ranganathan and Ostoja-Starzewski 2008). It accounts for both shear and bulk (compressional) contributions, is uniquely defined and is invariant to rotational transformations: where G and K are shear and bulk moduli calculated using the Voigt (V) and Reuss (R) averages, The indices i, j and k in the second column indicate different symmetry planes. The number of possible permutations for a given distribution is specified in the third column. Each class can have zero or four out-of-plane singularity directions, as specified in the parentheses.
s = c −1 is the compliance matrix in Voigt notation. In the isotropic limit, A U = 0. (1) and 'kiss' singularities.

Multiple roots of equation
If velocities of any two waves along any symmetry axis become equal, the degeneracy known as the 'kiss' singularity (Crampin 1981;Vavryčuk 1999) occurs (same as the on-axis singularity in TI media). At this point, the slowness sheets touch tangentially and do not exhibit the anomalous curvature observed in the vicinity of a conical point. Let us investigate the case when the shear stiffness moduli are equal along x i -axis, assuming that the orthorhombic medium is normal.
Then c ijij = c ikik , and equation (1) written for the two symmetry planes that contain axis x i simplifies to Equation (8) has a double root that corresponds to the on-axis singularity (n i = 1, n j = n k = 0). Solving it for n j , one arrives at Since the denominator on the right-hand side of equation (9) is positive in a normal medium, the condition that the real singularity directions exist in the symmetry planes adjacent to x i -axis requires a positive numerator. Employing an idea somewhat similar to that of the anisotropic notation of Thomsen (1986) and Tsvankin (1997), one can rewrite the condition above as where v P i = √ c iiii and v S i = √ c ijij are the on-axis velocities of the P-wave and S-wave polarized in the x j -direction, respectively, and δ j , j , γ j are the anisotropy parameters defined in the [x i , x j ] plane with respect to the x i axis, Further inspecting inequality 10, one can notice that the parameters inside the brackets are responsible for the shear-wave normal-moveout (NMO) velocities in TI media and in the symmetry planes of orthorhombic media. Hence, the inequality simplifies to where v SV j and v S H j are the NMO velocities of the shear waves in the [x i , x j ] plane defined for a reflector orthogonal to the Here, following the limited equivalence between the symmetry planes of ORT media and TI media, SV and SH notations are used. Since the NMO velocity v n is related to the curvature of the slowness surface, k (k ∝ −v 2 n ), condition 12 can be written as sgn(γ j )k SV j < sgn(γ j )k S H j , where k SV j and k S H j are the curvatures of the slowness surfaces of the SV-and SH-waves calculated at x j = x k = 0 in the x j -direction. The geometric interpretation of this condition then becomes clear and intuitive: in order for the ellipse (SH-wave) and a quartic curve (SV-wave) to have an intersection outside the x i -axis, their curvatures have to be related accordingly. This point is illustrated in Figure 2. Following the classification in Table 1, and counting the on-axis singularity once in each symmetry plane adjacent to x i , the allowed classes for normal media are I-V and VIII. If all on-axis shear-wave moduli are equal, equation (1) has no solutions other than n = x, and no in-plane singularities can occur. It is also possible to encounter a 'kiss' singularity outside the symmetry axes (Schoenberg and Helbig 1997). The condition for this situation inside the symmetry planes follows from equation (1) when the latter has a double root for n 2 j /n 2 i and is given by the following expressions: Note that equation (4) specifies a direction along which the phase velocities of two (or, potentially, three, see, e.g. Musgrave (1981)) modes coincide, but does not distinguish between the singularity types. It is possible (but highly unlikely) to encounter a 'kiss' singularity outside the symmetry planes in orthorhombic media (Crampin 1991), and the necessary condition (for arbitrary symmetry) in terms of the polarization vectors is given in Shuvalov (1998) and Vavryčuk (2003b).

Randomly generated orthorhombic models
Following Vavryčuk (2005), I conduct four numerical experiments to study the occurrence of the shear-wave singularities in orthorhombic media and its relation to the anisotropy strength. In the first experiment, the orthorhombic stiffness tensor is obtained by the following perturbation: where c ISO is the isotropic stiffness tensor with the Lamé parameters λ = μ = 1 GPa, ε determines the anisotropy strength andc ORT is a randomly generated orthorhombic stiffness tensor, such as the probability density function of each stiffness coefficient corresponding to a continuous uniform distribution in the interval [−3, 3] GPa, (c ORT ) i jkl ∼ U(−3, 3). Generated stiffness tensors are tested for the thermodynamic stability. Using equations (1) and (4), all singularity directions for the obtained model are calculated and their types (conical or tangential) are determined. Based on the number of singularities, a class is assigned to the model in accordance with Table 1, and it is also noted if out-of-plane singularities are present. The total of 10 5 thermodynamically stable models are generated and analysed. The resulting distribution of the orthorhombic models based on the number of singularities as a function of the anisotropy strength, ε and A U , is shown in Figure 3. The frequency of each class is normalized by the number of all possible permutations depending on which symmetry plane contains the singularities (the multiplicity in Table 1). For practically realizable models that correspond to ε < 1/3 and result in normal orthorhombic media, there exist only three classes (Figs. 3a,b), as predicted in Section 2, and it is almost equally probable to encounter a model with or without out-of-plane singularities. Classes II and IX have a multiplicity of 3, whereas class V has a multiplicity of 6. It is difficult to relate ε and A U since the latter has a large variance for a fixed value of ε.
The ε − A U cross-plot is shown in Figure 4, where the maximum threshold for A U is indicated in the graph. Importantly, all models with ε < 1/3 have singularities associated with shear waves only. Once ε becomes larger than 1/3, the shear stiffnesses (c 44 , c 55 , c 66 in Voigt notation) can become larger than the compressional stiffnesses (c 11 , c 22 , c 33 ). Then all classes become present in the generated models, including I (the maximum number of singularities) and X (no in-plane singularities) classes. Interestingly, there are no models with only out-of-plane singularities. For classes I, III, IV, VI-VIII  Table 1. The horizontal axis is the class number (total number of the in-plane singularities is indicated in parentheses), the vertical axis is the relative frequency in percent, purple and pink colours correspond to models with only in-plane singularities and models with in-plane and out-of-plane singularities, respectively. The total relative count of these two groups is shown inside the legend boxes. The frequencies are normalized by the class multiplicity. The average values of A U are displayed in the individual panels (a-d). and X, at least one singularity is associated with the P-wave. In all tested models, no 'kiss' singularities are encountered due to the very low probability of satisfying the necessary conditions.

Variable strength of anisotropy
In the second numerical experiment, one of the randomly generated orthorhombic stiffness matricesc ORT withc 11 = 2.88,c 22 = 2.8,c 33 = 2.09,c 44 = 1.9,c 55 = 1.6,c 66 = 0.26, c 23 = 1.38,c 13 = −0.72,c 12 = −1.47 km 2 s −2 is fixed and the orthorhombic model is obtained by varying the strength of anisotropy in equation (15), with ε ∈ [−5000, 5000]. For each obtained model (both stable and unstable), the out-of-plane singularity directions are calculated and plotted on the unit sphere (Fig. 5). A similar experiment for generally anisotropic media was performed by Vavryčuk (2005), who obtained singularity directions that followed closed and selfintersecting curves (comprised of the singularity directions in stable and unstable media). In the simpler orthorhombic case presented here, these curves form relatively complicated self-intersecting patterns and pass through the directions associated with the in-plane singularities. The complexity of the pattern is model dependent, and it is possible to encounter non-intersecting curves where the singularity directions that correspond to stable and unstable media are completely separated.

Anisotropy due to aligned fractures
In the third numerical experiment, the effective anisotropy due to aligned fractures is considered. I use the linear-slip theory (Schoenberg 1980;Schoenberg and Helbig 1997) to  Figure 3. Since A U has a large variance, its maximum value is limited by a threshold indicated in the plot. model orthorhombic media by embedding a single set (or multiple orthogonal sets) of parallel fractures in a transversely isotropic background or multiple orthogonal fracture sets in an isotropic background (Bakulin, Grechka and Tsvankin 2000). The occurrence of point singularities for the constructed media is then studied. The effective compliance tensor is obtained by adding the excess-compliance tensor due to the fractures to the compliance tensor of the background medium. The effective stiffness matrix of the orthorhombic medium obtained by putting two vertical orthogonal sets of parallel fractures into a transversely isotropic background has the form: with and where c i jb are the stiffness elements of the transverselyisotropic background medium with a vertical symmetry axis as many models with out-of-plane singularities, and most of the generated media fall into class V (six in-plane singularities), even after correcting for the multiplicity due to mirror symmetry. When only one fracture set in the VTI background is present (Fig. 6b) or both fracture sets in the isotropic background are rotationally invariant (Fig. 6d), there are no models that belong to class II. Schoenberg and Helbig (1997) provide a detailed discussion on the number and distribution of singularity directions in orthorhombic models with a single set of fractures in VTI background. The authors state that there is a maximum number of 14 degenerate points (class II in the present classification). The numerical observations discussed above, however, do not support this conclusion. Interestingly, the probability to encounter an out-of-plane singularity in class IX in the case of two sets of circular fractures in isotropic background is much lower than in any other case or class. As in the experiment in Section 3.1, no 'kiss' singularities are encountered in any of the tested models.

Samples from materials database
In the final numerical example, I apply the developed classification methodology to materials reported in Jong et al. (2015).   The database contains the elastic tensors of inorganic compounds taken from the Materials Project (an online database of inorganic compounds; see Jain et al. 2013). Out of the 1181 materials, I select 193 samples with orthorhombic symmetry. The resulting distribution based on the number of singularities is shown in Figure 7. Only one material, germanium selenide (GeSe), falls into class I, and two materials in class III are found, yttrium aluminium (YAl) and tellurium oxide (TeO 2 ). Although this database is of little interest in exploration seismology, it is nevertheless important to understand what type of singularity pattern one might expect.

Examples of phase and group surfaces
As a final remark, I would like to illustrate the peculiar behaviour and complexity of the slowness and group-velocity surfaces in some of the orthorhombic media encountered in numerical experiments. Many scholars note the complicated shape of these surfaces in the vicinity of singularity points in orthorhombic and lower-symmetry media (Crampin 1991;Brown et al. 1993;Grechka and Obolentseva 1993;Vavryčuk 2005;Grechka 2015Grechka , 2017Ivanov and Stovas 2017). Conical points and anomalous curvatures of the slowness sheets associated with them cause multi-valued wave fronts, caustics and strongly non-linear and rapidly varying particle motion (polarization) directions. Two realizations of orthorhombic media from the first experiment that belong to classes I and X and contain the maximum (16) and minimum (0) number of singularity directions, respectively, are chosen somewhat arbitrarily. The phase and group surfaces of all wave modes in these media are displayed in Figures 8 and 9, respectively. One can see that the corresponding surfaces are highly complex. However, I should note that neither of these two cases has been observed.

C O N C L U S I O N S
I analysed the conditions that govern the occurrence of singularity directions in orthorhombic media. The relationship between the anisotropy strength and the number and distribution of conical singularities in orthorhombic media is studied. It is demonstrated theoretically and confirmed numerically that for models in which the principal compressional stiffnesses are larger than the shear stiffnesses (referred to as normal orthorhombic media), the singularities can be distributed in only three different patterns (classes II, V and IX), and all singularities are associated with the shear waves only. The distribution of singularity directions between the classes in nor-mal orthorhombic media does not depend on the anisotropy strength due to a limited number of possible relations between the stiffness moduli. If both normal and non-normal orthorhombic models are considered, all possible singularity patterns are present, and the stronger the anisotropy, the more uniform is the distribution of models among the classes. However, the singularity directions are not restricted to the shear phase surfaces only. No relation is observed between the anisotropy strength and the total number of the singularity directions. The probability to encounter out-of-plane singularities is close to 50% regardless of the anisotropy strength. Tangential ('kiss') singularities can only be encountered in special cases, and neither in-plane nor out-of-plane 'kiss' singularities have been identified in the analysed models. By varying the anisotropy strength, it was shown that the out-of-plane singularities follow closed, potentially intersecting trajectories and pass through the symmetry planes (merging with the in-plane singularities).
In normal orthorhombic media obtained using the linearslip theory, all possible (II, V, IX) classes of singularity distributions are present only if two orthogonal sets of parallel fractures are embedded in a VTI or isotropic background. If either one set of fractures in a VTI medium or two sets of circular fractures in a purely isotropic background are considered, orthorhombic models of only two classes are observed (in contrast to previously published conclusions of Schoenberg and Helbig 1997).
The above results can potentially be used to constrain model parameters in the inversion of laboratory, microseismic and surface seismic data, provided that shear-wave singularities can be observed.

A C K N O W L E D G E M E N T S
Petromaks2 project is acknowledged for the financial support. The author is grateful to Vladimir Grechka for fruitful discussions related to the topic and his valuable comments on the initial version of the manuscript. Editor Ilya Tsvankin and an anonymous reviewer are acknowledged. Petr Bulant is thanked for help with the software used to triangulate group-velocity surfaces, openly available from Seismic Waves in Complex 3D Structures consortium research project.