## Abstract

The wave quantities needed in constructing wave fields propagating in anisotropic elastic media are usually calculated as a function of the slowness vector, or of its direction called the wave normal. In some applications, however, it is desirable to calculate the wave quantities as a function of the ray direction. In this paper, a method of calculating the slowness vector for a specified ray direction is proposed. The method is applicable to general anisotropy of arbitrary strength with arbitrary complex wave surface. The slowness vector is determined by numerically solving a system of multivariate polynomial equations of the sixth order. By solving the equations, we obtain a complete set of slowness vectors corresponding to all wave types and to all branches of the wave surface including the slowness vectors along the acoustic axes. The wave surface can be folded to any degree. The system of equations is further specified for rays shot in the symmetry plane of an orthorhombic medium and for a transversely isotropic medium. The system is decoupled into two polynomial equations of the fourth order for the *P**–SV* waves, and into equations for the *SH* wave, which yield an explicit closed-form solution. The presented approach is particularly advantageous in constructing ray fields, ray-theoretical Green functions, wavefronts and wave fields in strong anisotropy.

## 1. Introduction

The phase and group velocities of elastic waves are, generally, different in anisotropic media. The phase velocity (or its inverse called the slowness) describes the propagation of plane waves and is directed along the wave normal. The group velocity determines the signal propagation and energy transport and is directed along a ray (Červený 2001). The wave quantities required to construct wave fields are usually calculated as a function of the wave normal. This way is simple and satisfactory in many cases. Nevertheless, in some applications, calculating the wave quantities for a specified ray vector is more appropriate and desirable. Such approach, however, has not yet been fully developed, because it conceals complications. Whereas the slowness surface (defined by the directional variation of the slowness) is of the sixth order, the wave surface (defined by the directional variation of the group velocity) can be of up to the 150th order (Fedorov 1968; Musgrave 1970; Helbig 1994). The slowness surface is formed by three sheets corresponding to three types of waves (*P*, *S*1 and *S*2), and each wave quantity is unique for a specified slowness vector and a wave type. However, the wave surface can be folded (Musgrave 1970; Payton 1983; Every 1986; Wolfe 1995, 1998; Every *et al*. 1998; Vavryčuk 2003*a*,*b*) and several different values of a wave quantity (including the slowness vector) can correspond to one specified ray and wave type (see figure 1).

The problem of calculating the phase and group velocities for a specified ray has been treated by several authors. First, approximate formulae valid under various symmetries of weak anisotropy have been suggested by Vavryčuk (1997), Song & Every (2000), Pšenčík & Vavryčuk (2002) and Farra (2004). The weak anisotropy condition is employed by the authors in order to linearize the problem and make it easier to solve. Second, formulae valid under two-dimensional anisotropy and derived using the Stroh formalism were suggested by Wang (1995). Third, the group velocity along a general direction in three-dimensional strong triclinic anisotropy was calculated by Sharma (2002), who applied Newton's method of solving nonlinear simultaneous equations. His approach, however, encounters difficulties if the wave surface is folded.

In this paper, the problem is solved with the use of a system of three multivariate polynomial equations of the sixth order. The system allows calculation of all slowness vectors corresponding to a specified ray direction. Subsequently, all the other wave quantities can be easily computed. The approach can find applications in tracing rays, in constructing wavefronts, and in calculating ray-theoretical Green functions and wave fields in anisotropic media of arbitrary strength and symmetry (Every & Kim 1994; Wolfe 1998; Červený 2001).

## 2. Slowness surface

The propagation of waves in elastic anisotropic media is described by the Christoffel equation which reads (Červený 2001, eqn 2.2.35–37)(2.1)Tensor *Γ*_{jk} is the Christoffel tensor (also called the Christoffel matrix), *δ*_{jk} is the Kronecker delta, and *G* is the eigenvalue of *Γ*_{jk}. The Christoffel tensor *Γ*_{jk} can be defined either in terms of the slowness direction * n*,(2.2)or in terms of slowness vector

*=*

**p***/*

**n***c*,(2.3)where

*a*

_{ijkl}=

*c*

_{ijkl}/

*ρ*are the density-normalized elastic parameters,

*c*

_{ijkl}are the elastic parameters,

*ρ*is the density of the medium, and

*c*is the phase velocity of the wave. Tensor

*Γ*

_{jk}is positive-definite. It has three eigenvalues

*G*, which are real and positive, and three unit eigenvectors

*, called polarization vectors. Eigenvalues*

**g***G*(

*) and*

**n***G*(

*) read(2.4)(2.5)*

**p**The slowness surface is defined by the set of all slowness vectors * p*(

*) and can be determined by equation (2.5). The degree of the slowness surface is 6, thus any straight line intersects the surface at six points at the most. Since, the Christoffel tensor has three eigenvalues, the slowness surface consists of three sheets corresponding to three different waves (*

**n***P*,

*S*1 and

*S*2 waves). The inner sheet (corresponding to the

*P*wave) must be wholly convex (see Musgrave 1970, p. 92), but the other sheets (corresponding to the

*S*1 and

*S*2 waves) can also be locally concave or saddle shaped. The lines which separate convex, concave and saddle-shaped areas on the slowness surface are called the parabolic lines (Every & Kim 1994; Shuvalov & Every 1996, 1997; Vavryčuk 2003

*b*).

The sheets of the *P*, *S*1 and *S*2 waves can be fully detached, but they can also touch or intersect each other. Points common to sheets of different waves are called acoustic axes or singularities (see Khatkevich 1963, 1977; Alshits & Lothe 1979*a*, 2004; Musgrave 1985; Darinskij 1994; Boulanger & Hayes 1998; Vavryčuk 2001, 2003*c*). For the singularities, two eigenvalues of the Christoffel tensor coincide,(2.6)and the tensor is degenerate. Exceptionally, all three eigenvalues can coincide but this type of singularity is very rare and will not be considered here. In general, the behaviour of waves is anomalous at the singularity and its vicinity (Alshits & Lothe 1979*b*; Shuvalov 1998) and thus the singularities frequently complicate modelling of waves.

## 3. Wave surface

The group velocity vector * v*=

*(*

**v***) satisfies the equation (Červený 2001, eqn 2.2.65)(3.1)The direction*

**p***of the group velocity, , is called the ray direction or the ray vector. If we differentiate equation (3.1) using the theorem on implicit functions, we obtain (Červený 1972, eqn 15; Červený 2001, eqn 3.6.15)(3.2)where*

**N***D*

_{jk}is the matrix of cofactors of ,(3.3)where

*Γ*

_{jk}=

*Γ*

_{jk}(

*),*

**p***G*=

*G*(

*)=1 and*

**p***D*=

*D*(

*) is defined as(3.4)The group velocity can be calculated by equation (3.2) only for the so-called regular directions, which are defined by*

**p***D*≠0. This inequality is equivalent to the condition of the non-degenerate Christoffel tensor. On the contrary, if

*D*=0, the Christoffel tensor is degenerate, and the slowness vector lies along a singularity. At the singularity, the group velocity is calculated as (Červený 2001, eqn 3.6.10)(3.5)where

*=*

**g***(*

**g***) is the polarization vector of the degenerate wave. The polarization vector is not unique at the singularity, but can be of arbitrary direction in the plane perpendicular to the polarization vector of the non-degenerate wave. Polarization vectors near the singularity have anomalous properties, which depend on the type of singularity.*

**p**The wave surface (also called the group velocity surface) is defined by the set of all group velocity vectors * v*. They can be parameterized by slowness vector

*,*

**p***=*

**v***(*

**v***), by slowness direction*

**p***,*

**n***=*

**v***(*

**v***), or by ray direction*

**n***,*

**N***=*

**v***(*

**v***). Compared with the slowness surface, the wave surface is much more complicated. The degree of the wave surface may be significantly higher than that of the slowness surface (but must be less than 150; see Musgrave 1970, p. 92). If the slowness sheet of the wave in question contains saddle-shaped or concave areas, function*

**N***=*

**v***(*

**v***) is multivalued. This means that many group-velocity vectors can correspond to one specified ray direction. This effect is known as the ‘folding’ or the ‘triplication’ of the wave surface. The existence of folding of the wave surface is conditioned by the existence of parabolic lines on the slowness surface. The parabolic lines on the slowness surface are mapped on to caustics on the wave surface (also called cusps, cuspidal lines or cuspidal edges). Triplications and caustics are observed, particularly, near conical and wedge singularities (see Vavryčuk 2003*

**N***b*).

## 4. The slowness vector from a ray vector

From equations (3.2) and (3.5), we can calculate ray vector * N* as a function of slowness vector

*,*

**p***=*

**N***(*

**N***). In this section, we attempt to find the inverse relation,*

**p***=*

**p***(*

**p***). We will derive the equations valid for general triclinic anisotropy and discuss their applicability.*

**N**### (a) System of equations

Equation (3.2) can be altered to read(4.1)Taking into account that(4.2)and consequently,(4.3)we get(4.4)This equation can be viewed as a system of three linear equations in unknowns *N*_{1}, *N*_{2} and *N*_{3} provided that elastic parameters *a*_{ijkl} and components of the slowness vector *p*_{1}, *p*_{2} and *p*_{3} are known. However, the equation can also be viewed as a system of three coupled polynomial equations of the sixth order in three unknowns *p*_{1}, *p*_{2} and *p*_{3}, provided that elastic parameters *a*_{ijkl} and components of the ray direction *N*_{1}, *N*_{2} and *N*_{3} are known. Since the equation is invariant to inserting * N* or −

*, it should yield solutions for a ray direction with both orientations ±*

**N***.*

**N**### (b) Spurious solutions

Equation (4.4) yields not only true slowness vectors for a specified ray direction, but also some spurious solutions, which must be rejected. First, we reject complex-valued solutions that correspond to inhomogeneous plane waves. Second, we have to reject some of the solutions with *D*=0. These solutions were incorporated during the derivation of (4.4), when equation (3.2) was multiplied by *D*. Even though all true slowness vectors for a specified ray lie in regular directions (*D*≠0), the complete set of solutions of (4.4) will also comprise some spurious solutions characterized by *D*=0. However, we cannot simply scrap all slowness vectors yielding *D*=0, because some of them might be true and correspond to a singularity. In order to decide whether the solution with *D*=0 is spurious or a true singular solution, we have to proceed in the following way: first, we calculate the group velocity vector from equation (3.5) and check, whether equation (4.2) is satisfied. Second, we check whether the calculated group velocity vector is parallel to the specified ray direction. If equation (4.2) is satisfied and the group velocity is parallel to the specified ray, the slowness vector represents a true solution, which is pointed at the singularity. In the other case, the solution is spurious. Finally, we have also to reject all slowness vectors, which yield an orientation opposite to the specified ray direction.

## 5. Anisotropy of higher symmetry

### (a) Orthorhombic anisotropy

Let us specify the previously derived equations for a symmetry plane of orthorhombic anisotropy (ORT). The elastic parameters of ORT are defined in two-index notation as follows (see Musgrave 1970, p. 117): *a*_{11}, *a*_{22}, *a*_{33}, *a*_{44}, *a*_{55}, *a*_{66}, *a*_{12}, *a*_{13}, *a*_{23}. The other elastic parameters are zero. The equations will be studied for the symmetry plane *x*_{1}−*x*_{3}. Since *p*_{2}=0, the system (4.4) of three equations in three unknowns *p*_{1}, *p*_{2} and *p*_{3} reduces to two equations in two unknowns *p*_{1} and *p*_{3}. The equations read(5.1)where(5.2)The system of equations can be further simplified. Since two waves are always polarized in the symmetry plane (the *P* and *SV* waves), and the third wave is always transverse (the *SH* wave), the system can be decoupled into two systems of equations: equations for the *P**–SV* waves, and separately for the *SH* wave. The equations for the *SH* wave can be solved explicitly. If we apply the formulae for the group velocities derived by Musgrave (1970, eqn 9.3.6–9.3.8), we obtain(5.3)for the *P**–SV* waves, and(5.4)for the *SH* wave.

### (b) Transverse isotropy

Let us assume that the medium is transversely isotropic (TI) with a vertical axis of symmetry. The elastic parameters are specified in two-index notation as follows: *a*_{11}=*a*_{22}, *a*_{33}, *a*_{44}=*a*_{55}, *a*_{66}, *a*_{13}=*a*_{23}, *a*_{12}=*a*_{11}−2*a*_{66}, the other parameters being zero. Since the medium is rotationally symmetric about the vertical axis, it is sufficient to study all quantities and equations in the *x*_{1}*–x*_{3} plane only. Similarly to ORT, the system (4.4) of three equations in three unknowns *p*_{1}, *p*_{2} and *p*_{3} reduces to two equations in two unknowns *p*_{1} and *p*_{3}. The form of equations is identical to (5.1), but coefficients (5.2) slightly differ:(5.5)

Also, for TI, the system of equations can be further simplified. As for orthorhombic media in a symmetry plane, the system can be decoupled into two systems of equations: equations for the *P**–SV* waves, and separately for the *SH* wave. If we apply the formulae for the group velocities derived by Musgrave (1970, eqn 8.2.2*a*, 8.2.8, 8.2.9), we obtain(5.6)for the *P**–SV* waves, and(5.7)for the *SH* wave. Similar formulae can be derived also for tetragonal and cubic symmetries.

## 6. Numerical procedure

Equations (4.4), (5.3) and (5.6) represent systems of multivariate polynomial equations, which can be solved numerically using, for example, the MATHEMATICA software package. We obtain a set of real- and complex-valued slowness vectors * p*. We reject all complex-valued solutions. Then we calculate the group velocity vectors

*corresponding to all real-valued slowness vectors*

**v***using equation (3.5). We select only those solutions for which the following conditions are satisfied:(6.1)where is the group velocity. Obviously, the conditions in (6.1) can only be satisfied approximately, i.e. within the precision of the procedure for calculating the roots of the multivariate polynomial equations. Using equation (6.1), we also skip the solutions corresponding to the opposite orientation of the ray. Finally, we have to skip all multiplicities in the solutions. Once the set of true slowness vectors is known, we can trace rays or calculate all other wave quantities necessary for constructing wave fields such as phase velocities, polarization vectors, wave metric tensor (Červený 2002; Vavryčuk 2003*

**p***b*) or the Gaussian curvatures of the slowness and wave surfaces (Klimeš 2002).

Equation (4.4) also has some limitations. First, it cannot be applied to calculating slowness vectors at conical points of the wave surface, when an infinite set of slowness vectors corresponds to a single ray direction. Second, the equation is not suitable for calculating the slowness vectors under extremely weak anisotropy. In this case, the medium is closely degenerate (*D*→0), and the solution can be distorted by numerical errors. Instead, a modified approach should be applied such as iterations, or the linearized formulae obtained from the perturbation theory.

## 7. Example

The efficiency of the procedure is exemplified on tracing rays generated by a point source situated in a vertically inhomogeneous TI medium with a vertical axis of symmetry. The TI medium is used just for simplicity and clarity of the interpretation of the numerical results. Even though the example is very simple and confined to two-dimensions only, the anisotropy displays complexities such as singularities on the slowness surface and triplications on the wave surface, hence the example is complex enough to demonstrate the efficiency of the proposed approach. Obviously, the approach can be applied, not only to two-dimensional, but also to truly three-dimensional media, but the interpretation and presentation of the results would then be more involved.

In the presented example, the approach is applied in order to compare ray fields calculated by specifying two different initial conditions of rays at the source: the rays specified by initial slowness directions **n**_{0} and by initial ray directions **N**_{0}. The comparison will show that the standard ray tracing procedure, which utilizes the specification of the initial slowness directions, can produce ray fields with a rather non-uniform distribution of rays, while the specification of the initial ray directions produces a distribution of rays which is more uniform and thus more suitable for a wavefront construction.

As the medium, we assume the Bazhenov shale (Vernik & Liu 1997, appendix A, 12.507 ft) with the following elastic parameters (in km^{2} s^{−2}): *a*_{11}=*a*_{22}=23.52, *a*_{33}=10.89, *a*_{44}=*a*_{55}=5.29, *a*_{66}=9.42, *a*_{12}=4.69, *a*_{13}=*a*_{23}=9.46. These values characterize the medium at the source, which is situated at a depth of 1 km. The elastic parameters at other depths are calculated as(7.1)where *ϵ*=0.3 km^{−1}. The medium is strongly anisotropic with anisotropy strength of 38.0, 26.1 and 28.6% for *P*, *SV* and *SH* waves, respectively. The anisotropy strength is calculated as(7.2)where *c*^{max} and *c*^{min} are the maximum and minimum phase velocities of the wave. The medium displays a kiss singularity along the vertical axis. Since standard ray-tracing algorithms may fail for rays near the singularity, we use the ray-tracing algorithm proposed by Vavryčuk (2001, 2003*c*). This algorithm is suitable for tracing rays in strongly as well as weakly anisotropic media and in all kinds of singularities and their vicinities.

We will study the rays of the *SV* wave. The *SV*-wave surface displays two triplications: one is near the symmetry axis, and the other is near the symmetry plane (see figure 2). The triplication near the symmetry axis is more noticeable. The rays are shot from the source in the *x*_{1}*–x*_{3} plane in the angular interval 〈0°, 360°〉 in steps of 15°. Figure 3*a* shows rays shot with an equidistant step of initial directions of the ray, figure 3*b* shows rays with an equidistant step of initial directions of the slowness vector. In both cases the angular step is 15°. The figures show a very irregular distribution of rays when the equidistant angular step in slowness directions is imposed (see figure 3*b*). The rays concentrate in directions close to the vertical and horizontal axes. The other directions are covered very sparsely. If we perform the ray tracing and impose an equidistant step in initial directions of rays, we obtain a more regular and dense ray field (see figure 3*a*). This is even clearer in figure 3*c*,*d*, which depict rays in the vicinity of the source. Figure 3*c* also shows that for some initial directions, we trace not one but three different rays. This happens for rays along the vertical axis and for rays deviated from the vertical axis by ±15°. For these directions, the rays hit the triplicate wave surface, and correspond to its three different branches. A similar situation occurs for rays shot in the horizontal plane, but the triplication in the horizontal plane is almost indistinguishable and the three rays coincide.

Figure 3 also shows that the algorithm is capable of tracing the singular rays. As mentioned, the studied TI medium has the kiss singularity along the vertical axis. Figure 3*a* shows six rays shot from the source which have the vertical initial direction. Three rays go up and three rays go down. From the six rays, four are regular and two are singular. The regular rays start to deviate gradually from the vertical axis, the singular rays are vertical straight lines. Figure 4 provides a detailed view of the slowness and wave surfaces near the upper kiss singularity, and the correspondence between the vertical rays and the slowness vectors. The singular ray corresponds to the vertical slowness vector. The other two vertical rays are regular and have slowness vectors, which do not point at the singularity. Since the singular ray is defined by the vertical ray vector and the vertical slowness vector, the vertical gradient of elastic parameters cannot change their direction. Consequently, the singular ray is the vertical straight line. The initial directions of the regular rays are defined by the vertical ray vectors and by inclined slowness vectors, and the vertical gradient of elastic parameters causes the rays to bend gradually (see figure 5).

## 8. Conclusions

The calculation of a slowness vector from a specified ray direction in general anisotropy of arbitrary strength leads to solving a system of three coupled equations of the sixth order in three unknowns. Such system can be solved numerically using, for example, the MATHEMATICA software package. The system of equations is simplified for a higher symmetry of anisotropy. In transverse isotropy and in the symmetry plane of orthorhombic, tetragonal or cubic anisotropy, the system can be further decoupled into equations for the *P**–SV* waves and into equations for the *SH* wave. The system of equations for the *P**–SV* waves consists of two equations of the fourth order in two unknowns. For the *SH* wave, an explicit closed-form solution can be obtained.

The basic advantage of the proposed approach is that we obtain a full set of slowness vectors corresponding to a specified ray even in the case of complex wave surface with the folding and multiple cusp edges. From the slowness vectors, we can calculate all wave quantities needed for modelling waves in anisotropic media. The approach works properly in regular directions as well as in all kinds of singularities (kiss, intersection, wedge or conical singularity). The approach is not applicable to conical points on the wave surface when an infinite set of slowness vectors correspond to a single ray direction. Also, for extremely weak anisotropy, the approach can yield distorted results due to numerical errors. In this case, it is more appropriate to apply iterations or the linearized formulae obtained from the perturbation theory.

## Acknowledgements

I thank Klaus Helbig and one anonymous referee for their comments. The work was supported by the Consortium Project SW3D ‘Seismic waves in complex 3-D structures’, and by the Grant Agency of the Academy of Sciences of the Czech Republic, grant A3012309.

## Footnotes

- Received May 9, 2005.
- Accepted October 31, 2005.

- © 2006 The Royal Society