# First-order ray computations of coupled S waves in inhomogeneous weakly anisotropic media

Véronique Farra, Ivan Pšenčík

► To cite this version:

Véronique Farra, Ivan Pšenčík. First-order ray computations of coupled S waves in inhomogeneous weakly anisotropic media. Geophysical Journal International, 2008, 173, pp.979-989. 10.1111/j.1365-246X.2008.03778.x. insu-03603728

HAL Id: insu-03603728

https://insu.hal.science/insu-03603728

Submitted on 10 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# First-order ray computations of coupled S waves in inhomogeneous weakly anisotropic media

Véronique Farra<fnref n="1" /> and Ivan Pšenčík<fnref n="2" />

Accepted 2008 March 3. Received 2008 March 3; in original form 2007 October 3

## SUMMARY

We propose an approximate procedure for computing coupled S waves in inhomogeneous weakly anisotropic media. The procedure is based on the first-order ray tracing (FORT) and dynamic ray tracing (FODRT), which was originally developed for P waves. We use the so-called common ray tracing concept to derive approximate ray tracing and dynamic ray tracing equations, and an approximate solution of the transport equation for coupled S waves propagating in laterally varying, weakly anisotropic media. In our common ray tracing, ray equations are governed by the first-order Hamiltonian formed by the average of first-order eigenvalues of the Christoffel matrix, corresponding to the two S-wave modes propagating in anisotropic media. The solution of the transport equation for the coupled S waves leads to a system of two coupled frequency-dependent, linear ordinary differential equations for amplitude coefficients, which is evaluated along the S-wave common ray. For derivation of the FORT and FODRT equations, we use the perturbation theory in which deviations of anisotropy from isotropy are considered to be first-order perturbations. To derive the coupled differential equations for S-wave amplitudes, we assume that the first-order perturbations are of order $O(\omega^{-1})$, where $\omega$ is the circular frequency. This makes it possible to express the amplitude coefficients in the coupled differential equations in terms of geometrical spreading and other quantities related to the common ray. The proposed procedure removes problems of most currently available ray tracers, which yield distorted results or even collapse when shear waves propagating in weakly anisotropic media are computed. The first-order approximation leads to simpler ray tracing and dynamic ray tracing equations than in the exact case. For anisotropic media of higher-symmetry than monoclinic, all equations involved differ only slightly from the corresponding equations for isotropic media. If the anisotropy vanishes, the equations reduce to standard, exact ray tracing and dynamic ray tracing equations for S waves propagating in isotropic media. The proposed ray tracing and dynamic ray tracing equations, corresponding traveltimes and geometrical spreading are all given to the first order. The accuracy of the traveltimes along S-wave first-order common rays can be increased by calculating a second-order traveltime correction.

**Key words:** Body waves; Seismic anisotropy; Wave propagation.

## 1 INTRODUCTION

In our two recent papers, Pšenčík & Farra (2005, 2007), we proposed the so-called first-order ray tracing (FORT) and first-order dynamic ray tracing (FODRT) for P waves propagating in smoothly varying weakly anisotropic media. We used the perturbation theory in which deviations of anisotropy from isotropy were considered to be first-order perturbations. In this paper, we propose an extension of the above procedures to coupled S waves propagating in smooth, laterally varying, weakly anisotropic structures. It removes problems of standard ray tracers, which deal with S waves separately, and yield distorted results or even collapse. The main application of the proposed procedure will be in modelling S-wavefields in inhomogeneous weakly anisotropic media. Further applications are expected in multicomponent seismic imaging, in Kirchhoff migrations including S or PS waves, in tomography. In all these applications, the proposed procedure can simply substitute currently used isotropic ray tracers.

Farra (2005) proposed FORT and FODRT equations applicable separately to each of the two S waves propagating in anisotropic media. In inhomogeneous weakly anisotropic media, however, S waves are usually coupled. Coupling can be described by various techniques, for

<footnote n="1">¹Institut de Physique du Globe de Paris, 4 Place Jussieu, 75252 Paris Cedex 05, France</footnote>

<footnote n="2">²Institute of Geophysics, Acad. Sci. of Czech Republic, Boční II, 141 31 Praha 4, Czech Republic. E-mail: ip@ig.cas.cz</footnote>
example, by the ‘connection’ formulae (Thomson et al. 1992) or by the coupling ray theory (Coates & Chapman 1990; Bulant & Klimeš 2002). In the coupling ray theory, the S waves are calculated along a common ray. There are various simplified versions of the coupling ray theory, one of them being the quasi-isotropic (QI) approximation (Kravtsov 1968; Pšenčík 1998). The basic idea of the QI approximation is to seek the asymptotic solution of the elastodynamic equation as an expansion with respect to two small parameters of the same order: a small parameter $\epsilon_1$ ($\epsilon_1 \sim \omega^{-1}$, where $\omega$ is the circular frequency) used in the standard ray series expansion (Červený 2001) and a parameter $\epsilon_2$ characterizing differences of tensors of elastic parameters of the weakly anisotropic medium and of the nearby ‘reference’ isotropic medium. The common ray for both S waves in the QI approximation is constructed in a reference isotropic medium. Along this S-wave common ray, coupled amplitude coefficients are calculated by integrating two coupled frequency-dependent ordinary differential equations. The QI approach has several limitations, leading to its lower accuracy. The most important of these limitations are the common ray traced in the isotropic reference medium (such a ray may differ considerably from the trajectory along which the energy of the two coupled S waves propagates in the anisotropic medium), the inaccuracy of the calculated traveltime and the geometrical spreading calculated along the reference ray in the reference isotropic medium. To reduce these limitations, Bakker (2002) proposed the use of an artificial Hamiltonian obtained as the average of the actual Hamiltonians of the two S waves propagating in an anisotropic medium. This approach led to an S-wave common ray, which was traced in the studied anisotropic medium. In contrast to the rays of individual S waves (obtained from their actual Hamiltonians), tracing the S-wave common ray in an anisotropic medium does not encounter problems in S-wave singular regions. In a way, tracing an S-wave common ray is like tracing a ray of a P wave. A detailed review and analysis of S-wave common ray tracing and dynamic ray tracing can be found in Klimeš (2006b).

In this paper, we use Bakker’s (2002) approach to extend applicability of the QI approximation. As in the QI approximation, we consider the two small parameters $\epsilon_1$ and $\epsilon_2$. In the first step, we construct the common S-wave ray as Bakker (2002) did, by using averaged S-wave eigenvalues of the Christoffel matrix. Instead of exact eigenvalues used by Bakker (2002), we, however, use their first-order approximations as in Pšenčík & Farra (2005, 2007). We use the perturbation theory, in which the deviation of anisotropy from isotropy is considered to be a first-order perturbation. In this way, we derive the FORT and FODRT equations along the common rays of coupled S waves. In the second step, we derive coupled differential equations for the amplitude coefficients along the common rays. For their derivation, we assume that both small parameters $\epsilon_1$ and $\epsilon_2$ are of the same order, $O(\omega^{-1})$. Described two-step procedure leads to considerably more accurate approach than the QI approximation, applicable to weakly anisotropic media of arbitrary symmetry. As in the case of P waves, the accuracy of the first-order traveltime calculated along a first-order common ray can be increased by calculation of a simple second-order correction.

In Section 2, we review the basic equations of the ray theory used in the derivations in the following sections. In Section 3, an explicit expression for the mean value of the first-order S-wave eigenvalues of the Christoffel matrix is given. The mean value is used in deriving the first-order common-ray tracing for S waves. The formula for the second-order traveltime correction, which can be evaluated by quadratures along the first-order common ray, is also given in Section 3. Section 4 is devoted to FODRT along the S-wave first-order common ray. It is shown how the results of dynamic ray tracing can be used to determine the first-order relative geometrical spreading. Section 5 gives the solution of the transport equation for a specific ansatz containing amplitude coefficients of the coupled S waves. The ansatz is very similar to ansatzs used by Pšenčík (1998) and Bakker (2002). In Section 5, the basic formula of the paper, the system of two coupled ordinary differential equations for computing the amplitude coefficients of the coupled S waves along the first-order common ray, is given. The paper ends with a short concluding section. Appendix A contains expressions for weak-anisotropy (WA) parameters and for elements of matrix $\tilde{B}$, which controls wave propagation in weakly anisotropic media. Appendix B contains expressions for averaged S-wave eigenvalues and their derivatives for orthorhombic media and for transversely isotropic media with vertical axes of symmetry.

The lower-case indices i, j, k, l, ... take the values of 1,2,3, the upper-case indices I, J, K, L, ... take the values of 1, 2. The Einstein summation convention over repeated indices is used. The index following the comma in the subscript denotes the partial derivative with respect to the relevant Cartesian coordinate. The upper index [M] is used to denote quantities related to the S-wave common ray.

## 2 BASIC EQUATIONS

In the frequency domain, the homogeneous elastodynamic equation (without a source term) for an inhomogeneous anisotropic medium reads
$$ (c_{ijkl}u_{k,l})_j + \rho\omega^2 u_i = 0. \quad (1) $$

In eq. (1), $c_{ijkl} = c_{ijkl}(x_m)$ are elements of the tensor of the elastic moduli, $\rho = \rho(x_m)$ is the density and $\omega$ is the circular frequency. Symbol $u_i = u_i(x_m, \omega)$ denotes the components of the sought displacement vector $\mathbf{u}$. The solution of the elastodynamic equation (1) is sought in the form of the zero-order term of the ray series:

$$ u_i(x_m, \omega) = U_i(x_m)\exp[i\omega\tau(x_m)]. \quad (2) $$

Here i is the imaginary unit, $\tau(x_m)$ is the eikonal, which also serves as the traveltime, and $U_i(x_m)$ are components of the zero-order vectorial amplitude coefficient $\mathbf{U}(x_m)$. By inserting (2) into (1), we arrive at

$$ (i\omega)^2 N_i(U_n) + i\omega M_i(U_n) + L_i(U_n) = 0. \quad (3) $$

The differential vectorial operators $\mathbf{N}$, $\mathbf{M}$ and $\mathbf{L}$ are given by the relations:

$$ N_i(U_n) = a_{ijkl} U_k \tau_{l,j} - U_i, \quad (4) $$
$$M_i(U_n) = \rho^{-1}(\rho a_{ijkl} U_k \tau_{,l})_j + a_{ijkl} U_{k,l} \tau_{,j}, \qquad (5)$$

$$L_i(U_n) = \rho^{-1}(\rho a_{ijkl} U_{k,l})_j \qquad (6)$$

(see Červený 2001). In eqs (4)–(6), $a_{ijkl}$ are density normalized elastic moduli,

$$a_{ijkl} = c_{ijkl} / \rho. \qquad (7)$$

Since eq. (3) should be satisfied for arbitrarily high frequency $\omega$, operators $\mathbf{N}$, $\mathbf{M}$ and $\mathbf{L}$ should be zero. From $\mathbf{N} = \mathbf{0}$, we get the Christoffel equation,

$$(\Gamma_{ik} - \delta_{ik})U_k = 0, \qquad (8)$$

where

$$\Gamma_{ik}(x_m, p_m) = a_{ijkl}(x_m) p_j p_l. \qquad (9)$$

In eq. (9), $\Gamma_{ik}$ are the elements of the generalized Christoffel matrix $\Gamma$. We call it generalized because it contains components $p_j = \tau_{,j} = n_j/c$ of slowness vector $\mathbf{p}$ instead of components $n_j$ of unit vector $\mathbf{n}$, used in the standard Christoffel matrix; $c$ is the phase velocity. For the sake of brevity, we call $\Gamma(x_m, p_m)$ the Christoffel matrix in the following. Each of the three eigenvalues $G(x_m, p_m)$ of the Christoffel matrix (eq. 9) can be used in the eikonal equation $G = 1$ of the corresponding wave. The eikonal equation can be solved using the ray tracing equations, from which the traveltimes of the considered wave can be determined. Equation $\mathbf{M} = \mathbf{0}$, with $\mathbf{M}$ given in (5), yields the transport equation, from which ray amplitudes along the corresponding ray can be calculated. Equations $\mathbf{N} = \mathbf{0}$ and $\mathbf{M} = \mathbf{0}$ are sufficient for the complete determination of the displacement vector $\mathbf{u}$ in eq. (2). The magnitude of the operator $\mathbf{L}$ could be used to check precision of the approximate solution (2).

Below, we introduce the FORT and FODRT equations with which we can determine the S-wave common rays, the first- or second-order common traveltime $\tau^{[\mathcal{M}]}$ or $\tau^{[\mathcal{M}]} + \Delta\tau^{[\mathcal{M}]}$ and the corresponding first-order relative geometrical spreading $\mathcal{L}^{[\mathcal{M}]}$, a quantity important for determining the ray amplitudes.

### 3 FIRST-ORDER COMMON RAY TRACING FOR S WAVES

To calculate the S-wave first-order common rays and ray amplitudes, we use a procedure similar to that used by Pšenčík & Farra (2005, 2007), who obtained the FORT and FODRT equations for P waves. Instead of using the largest of the three first-order eigenvalues of the Christoffel matrix (9), $G^{[3]} = G^{[3]}(x_m, p_m)$, which refers to P waves, we use the mean value, $G^{[\mathcal{M}]} = G^{[\mathcal{M}]}(x_m, p_m)$, of the smaller first-order eigenvalues $G^{[1]}$ and $G^{[2]}$ of the Christoffel matrix $\Gamma$, which refer to S waves. Farra & Pšenčík (2003) and Farra (2005) give expressions for the first-order S-wave eigenvalues $G^{[1]}$, $G^{[2]}$ and P-wave eigenvalue $G^{[3]}$:

$$G^{[1]} = \frac{1}{2}(B_{11} + B_{22} + \sqrt{D}), \qquad G^{[2]} = \frac{1}{2}(B_{11} + B_{22} - \sqrt{D}), \qquad G^{[3]} = B_{33}, \qquad (10)$$

where

$$D = (B_{11} - B_{22})^2 + 4B_{12}^2. \qquad (11)$$

In eqs (10) and (11), symbols $B_{jl}$ are elements of the symmetric matrix $\mathbf{B}(x_m, p_m)$:

$$B_{jl}(x_m, p_m) = \Gamma_{ik}(x_m, p_m) e_i^{[j]} e_k^{[l]}, \qquad (12)$$

(see Farra & Pšenčík 2003). Symbols $e_i^{[j]}$ denote the components of three mutually perpendicular unit vectors $\mathbf{e}^{[1]}$, $\mathbf{e}^{[2]}$ and $\mathbf{e}^{[3]}$. Vector $\mathbf{e}^{[3]}$ is chosen so that $\mathbf{e}^{[3]} = \mathbf{n}$, where $\mathbf{n}$ is a unit vector specifying direction of the slowness vector $\mathbf{p}$. Vectors $\mathbf{e}^{[1]}$ and $\mathbf{e}^{[2]}$ can be chosen arbitrarily in the plane perpendicular to $\mathbf{n}$. Matrix $\mathbf{B}(x_m, p_m)$ [or $\mathbf{B}(x_m, n_m)$] can be obtained by simple rotation,

$$\mathbf{B} = \mathbf{R}^T \bar{\mathbf{B}} \mathbf{R}, \qquad (13)$$

from matrix $\bar{\mathbf{B}}(x_m, p_m)$ [or $\bar{\mathbf{B}}(x_m, n_m)$], specified explicitly in Appendix A, eqs (A3) and (A4). Rotation matrix $\mathbf{R}$ has the following form:

$$\mathbf{R} = \begin{pmatrix} \cos \Phi & -\sin \Phi & 0 \\ \sin \Phi & \cos \Phi & 0 \\ 0 & 0 & 1 \end{pmatrix}. \qquad (14)$$

In eq. (14), angle $\Phi$ is the acute angle between the basis vectors $\mathbf{e}^{[K]}$ introduced above and the vectors $\bar{\mathbf{e}}^{[K]}$ used in Appendix A (see eq. A5). Matrices $\mathbf{B}$ and $\bar{\mathbf{B}}$ are independent of the choice of the reference velocities $\alpha$ and $\beta$ (see Appendix A). The reference velocities $\alpha$ and $\beta$ appear in the matrices $\mathbf{B}$ and $\bar{\mathbf{B}}$ as a consequence of the use of $\alpha$ and $\beta$ in definitions of WA parameters (see Appendix A, eq. A1). The reference velocities can thus be chosen as arbitrary non-zero quantities. From eqs (12) to (14) it is easy to prove that terms $B_{11} + B_{22}$ and $B_{13}^2 + B_{23}^2$ are independent of the choice of vectors $\mathbf{e}^{[1]}$ and $\mathbf{e}^{[2]}$.

Using eqs (10) and (12), we can express the S-wave first-order mean value $G^{[\mathcal{M}]}$ in the following way:

$$G^{[\mathcal{M}]} = \frac{1}{2} (G^{[1]} + G^{[2]}) = \frac{1}{2} (B_{11} + B_{22}) = \frac{1}{2} \Gamma_{ik} (e_i^{[1]} e_k^{[1]} + e_i^{[2]} e_k^{[2]}) = \frac{1}{2} \Gamma_{ik} (\delta_{ik} - e_i^{[3]} e_k^{[3]}) = \frac{1}{2} [\mathrm{Tr}(\Gamma) - G^{[3]}]. \qquad (15)$$
The FORT equations for S-wave common rays can then be expressed as follows:

$$ \frac{dx_i}{d\tau} = \frac{1}{2} \frac{\partial G^{[\mathcal{M}]}[x_i]}{\partial p_i}, \qquad \frac{dp_i}{d\tau} = -\frac{1}{2} \frac{\partial G^{[\mathcal{M}]}[x_i]}{\partial x_i}. \qquad (16) $$

Here $G^{[\mathcal{M}]}$ is given by eq. (15). In eq. (16), $x_i = x_i(\tau)$ are the Cartesian coordinates of the trajectory of the first-order common ray of an S wave. Symbols $p_i = p_i(\tau)$ denote components of the first-order slowness vectors, $p_i = n_i/c^{[\mathcal{M}]}$. Here $n_i$ are again components of the unit vector $\mathbf{n}$ specifying direction of $\mathbf{p}$, and $c^{[\mathcal{M}]} = c^{[\mathcal{M}]}(x_m, n_m)$ is the S-wave first-order common phase velocity corresponding to $G^{[\mathcal{M}]}$; $(c^{[\mathcal{M}]})^2 = G^{[\mathcal{M}]}(x_m, n_m)$. The parameter along the common ray, $\tau$, is the first-order common traveltime $\tau^{[\mathcal{M}]}$. At each point along the ray, the first-order slowness vector $\mathbf{p}$ determined from eq. (16) should satisfy the first-order eikonal equation

$$ G^{[\mathcal{M}]}(x_m, p_m) = [c^{[\mathcal{M}]}(x_m, n_m)]^{-2} G^{[\mathcal{M}]}(x_m, n_m) = 1. \qquad (17) $$

The explicit expression for $G^{[\mathcal{M}]}$ could be obtained from eq. (15), using, for example, expressions for $\tilde{B}_{11}$ and $\tilde{B}_{22}$ given in Appendix A, and exploiting the fact that $B_{11}(x_m, p_m) + B_{22}(x_m, p_m) = \tilde{B}_{11}(x_m, p_m) + \tilde{B}_{22}(x_m, p_m)$. Instead, we use the last term on the right-hand side of eq. (15) and express the S-wave first-order mean value $G^{[\mathcal{M}]}(x_m, p_m)$ in terms of Tr($\Gamma$) and $G^{[3]}$. The explicit formula for the trace of the Christoffel matrix Tr($\Gamma$) reads

$$ \begin{aligned} \operatorname{Tr}(\Gamma) ={}& \alpha^2 p_i p_i + 2\alpha^2 [\epsilon_x p_1^2 + \epsilon_y p_2^2 + \epsilon_z p_3^2 \\ & + (\epsilon_{16} + \epsilon_{26})p_1 p_2 + (\epsilon_{15} + \epsilon_{35})p_1 p_3 + (\epsilon_{24} + \epsilon_{34})p_2 p_3] \\ & + 2\beta^2 [p_k p_k + \epsilon_{45} p_1 p_2 + \epsilon_{46} p_1 p_3 + \epsilon_{56} p_2 p_3 \\ & + \gamma_x (p_2^2 + p_3^2) + \gamma_y (p_1^2 + p_3^2) + \gamma_z (p_1^2 + p_2^2)]. \end{aligned} \qquad (18) $$

The explicit formula for the first-order P-wave eigenvalue $G^{[3]} = B_{33}$ of the Christoffel matrix (9) can be found from Appendix A. For $G^{[\mathcal{M}]}$, eq. (15) then yields

$$ G^{[\mathcal{M}]} = \beta^2 [p_i p_i + \epsilon_{45} p_1 p_2 + \epsilon_{46} p_1 p_3 + \epsilon_{56} p_2 p_3 + \gamma_y (p_1^2 + p_3^2) + \gamma_x (p_2^2 + p_3^2) + \gamma_z (p_1^2 + p_2^2)] - \alpha^2 (p_i p_i)^{-1} [\eta_x p_2^2 p_3^2 + \eta_y p_1^2 p_3^2 \\ & + \eta_z p_1^2 p_2^2 + (\chi_x - \epsilon_{16})p_1 p_2 (p_2^2 + p_3^2 - p_1^2) + (\chi_y - \epsilon_{26})p_1 p_2 (p_1^2 + p_3^2 - p_2^2) + (\chi_z - \epsilon_{24})p_2 p_3 (p_1^2 + p_3^2 - p_2^2) \\ & + (\chi_y - \epsilon_{15})p_1 p_3 (p_2^2 + p_3^2 - p_1^2) + (\chi_x - \epsilon_{34})p_2 p_3 (p_1^2 + p_2^2 - p_3^2) + (\chi_y - \epsilon_{35})p_1 p_3 (p_1^2 + p_2^2 - p_3^2)]. \qquad (19) $$

Eq. (19) holds for weak anisotropy of any symmetry. The coefficients of eq. (19) contain linear combinations of the WA parameters. For the definition of the WA parameters and parameters $\eta_x, \eta_y$ and $\eta_z$, see Appendix A. The S-wave first-order mean value $G^{[\mathcal{M}]}(x_m, p_m)$ and its derivatives are independent of reference velocities $\alpha$ and $\beta$. Another interesting feature of eq. (19) is that, as in the P-wave case (see Pšenčík & Farra 2005), $G^{[\mathcal{M}]}$ is expressed in terms of 15 independent coefficients, which are formed by 21 WA or elastic parameters. The S-wave coefficients differ from the P-wave coefficients. As shown by Farra & Pšenčík (2003), the whole set of P- and S-wave coefficients can be used for estimation of all 21 WA parameters.

The initial conditions for the ray tracing equations (16) for $\tau = 0$ read:

$$ x_i(0) = x_i^0, \qquad p_i(0) = p_i^0. \qquad (20) $$

Here, $x_i^0$ are the coordinates of source point $\mathbf{x}^0$, and $p_i^0 = n_i^0/c_0^{[\mathcal{M}]}$ are the components of the first-order slowness vector $\mathbf{p}^0$ at the source. Symbol $c_0^{[\mathcal{M}]}$ denotes the first-order approximation of the S-wave common phase velocity in direction $\mathbf{n}^0$ at source point $\mathbf{x}^0$. Velocity $c_0^{[\mathcal{M}]}$ is given by the square root of $G^{[\mathcal{M}]}(x_m^0, n_m^0)$ (see eq. 17). Vector $\mathbf{n}^0$ can be specified by two ray parameters, $\gamma^{(J)}$, chosen as two take-off angles, $\phi_0$ and $\delta_0$, so that

$$ n_1^0 = \cos \phi_0 \cos \delta_0, \quad n_2^0 = \sin \phi_0 \cos \delta_0, \quad n_3^0 = \sin \delta_0. \qquad (21) $$

The explicit expressions for $G^{[\mathcal{M}]}$ and their derivatives, required in equation (16), for orthorhombic and VTI media, can be found in Appendix B. If the anisotropy vanishes, equations (16) reduce to the exact ray tracing equations for an S wave in an inhomogeneous isotropic medium.

By integrating FORT system in eq. (16), we obtain the first-order common traveltime $\tau^{[\mathcal{M}]}$. The accuracy of the traveltime can be simply enhanced by calculating the second-order traveltime correction $\Delta\tau^{[\mathcal{M}]}$ along the first-order common ray (16):

$$ \Delta\tau^{[\mathcal{M}]} = -\frac{1}{2} \int_{\Omega_0} \Delta G^{[\mathcal{M}]}(x_m, p_m) d\tau. \qquad (22) $$

Here $\Omega_0$ denotes the first-order common ray of an S wave, $x_m$ are its coordinates and $p_m$ are the components of the first-order slowness vectors along $\Omega_0$. Symbol $\Delta G^{[\mathcal{M}]}$ denotes the difference between the exact mean value $G_{ex}^{[\mathcal{M}]}(x_m, p_m)$ and its first-order approximation $G^{[\mathcal{M}]}(x_m, p_m)$ along $\Omega_0$. To estimate $G_{ex}^{[\mathcal{M}]}$, we use the mean of the second-order expressions for the S-wave eigenvalues. The expressions for the second-order S-wave eigenvalues given by Farra & Pšenčík (2003) lead to

$$ \Delta G^{[\mathcal{M}]}(x_m, p_m) = -\frac{1}{2} [c^{[\mathcal{M}]}(x_m, n_m)]^2 \frac{B_{13}^2(x_m, p_m) + B_{23}^2(x_m, p_m)}{V_p^2 - V_S^2}. \qquad (23) $$

Inserting (23) into (22), we arrive at the final expression for the second-order traveltime correction:

$$ \Delta\tau^{[\mathcal{M}]} = \frac{1}{4} \int_{\Omega_0} [c^{[\mathcal{M}]}(x_m, n_m)]^2 \frac{B_{13}^2(x_m, p_m) + B_{23}^2(x_m, p_m)}{V_p^2 - V_S^2(x_m)} d\tau. \qquad (24) $$
Quantities $B_{13}$ and $B_{23}$ in eq. (24) are elements of the symmetric matrix **B** (see eq. 12). Because the term $B_{13}^2 + B_{23}^2$ is independent of the choice of vectors **e**[1] and **e**[2], we can use directly $\tilde{B}_{13}^2 + \tilde{B}_{23}^2$, where $\tilde{B}_{13}$ and $\tilde{B}_{23}$ are elements of matrix $\tilde{\mathbf{B}}(x_m, p_m)$ given in eq. (A3).

Symbols $V_P$ and $V_S$ in (24) denote the P- and S-wave velocities corresponding to the reference isotropic medium, closely approximating the studied weakly anisotropic medium along ray $\Omega_0$. Similarly to the case of P waves, we can use $V_S^2 = (p_k p_k)^{-1}$ and $V_P^2 = 3V_S^2$, where $p_k$ are the components of the first-order slowness vector obtained in integrating the FORT system (eq. 16). The values of $V_P$ and $V_S$ determined in this way vary along the ray $\Omega_0$. The above specification yielded the most accurate results in the case of P waves and it is expected to do the same in the case of S waves.

Note the difference in signs between eq. (24) and the corresponding second-order correction for the P-wave traveltimes (Pšenčík & Farra 2005, 2007). Obviously, the S-wave common traveltime correction is always positive. It agrees with the observation of Farra & Pšenčík (2003) that the first-order S-wave eigenvalues are always larger than or equal to the exact ones.

## 4 FIRST-ORDER DYNAMIC RAY TRACING ALONG THE S-WAVE COMMON RAY

To calculate the amplitudes of the coupled S waves (see Section 5), we need to know the first-order relative geometrical spreading $\mathcal{L}^{[\mathcal{M}]}(R, S)$ along the first-order common ray from source S to receiver R. For the sake of brevity, we shall call $\mathcal{L}^{[\mathcal{M}]}(R, S)$ the geometrical spreading below. It is defined as

$$ \mathcal{L}^{[\mathcal{M}]}(R, S) = |\mathbf{X}^{(1)} \times \mathbf{X}^{(2)}|^{1/2}. \quad (25) $$

Vectors $\mathbf{X}^{(1)}$ in eq. (25) and other important vectors $\mathbf{Y}^{(1)}$, defined below, are determined from the FODRT system with especially chosen initial conditions, given below (see also Pšenčík & Teles 1996). If we denote the ray parameters specifying the S-wave common ray by $\gamma^{(1)}$, where $\gamma^{(1)}$ are the initial angles, $\phi_0$ and $\delta_0$, specifying the orientation of vector $\mathbf{n}^0$ (see eq. 21), the components of vectors $\mathbf{X}^{(1)}$ and $\mathbf{Y}^{(1)}$ are:

$$ X_i^{(1)} = \left[ \frac{\partial x_i}{\partial \gamma^{(1)}} \right]_{\tau=\text{const}}, \quad Y_i^{(1)} = \left[ \frac{\partial p_i}{\partial \gamma^{(1)}} \right]_{\tau=\text{const}}. \quad (26) $$

Quantities $X_i^{(1)} = X_i^{(1)}(\tau)$ and $Y_i^{(1)} = Y_i^{(1)}(\tau)$ describe the variations along the wave front of coordinates $x_i$ of the first-order common ray and of components $p_i$ of the first-order slowness vector due to the variations of parameters $\gamma^{(1)}$. Values of $X_i^{(1)}$ and $Y_i^{(1)}$ can be found by solving a system of linear differential equations obtained by differentiating the FORT equations (16) with respect to $\gamma^{(1)}$. The resulting system of differential equations is the FODRT system along the common ray of an S wave. It can be developed in the same way as for P waves (see Pšenčík & Farra 2007):

$$ \begin{aligned} \frac{dX_i^{(1)}}{d\tau} &= \frac{1}{2} \left( \frac{\partial^2 G^{[\mathcal{M}]}(x_m, p_m)}{\partial p_i \partial x_j} X_j^{(1)} + \frac{\partial^2 G^{[\mathcal{M}]}(x_m, p_m)}{\partial p_i \partial p_j} Y_j^{(1)} \right), \\ \frac{dY_i^{(1)}}{d\tau} &= -\frac{1}{2} \left( \frac{\partial^2 G^{[\mathcal{M}]}(x_m, p_m)}{\partial x_i \partial x_j} X_j^{(1)} + \frac{\partial^2 G^{[\mathcal{M}]}(x_m, p_m)}{\partial x_i \partial p_j} Y_j^{(1)} \right). \end{aligned} \quad (27) $$

Here, $G^{[\mathcal{M}]}(x_m, p_m)$ is given in eqs (15) or (19).

The explicit expressions for $G^{[\mathcal{M}]}$ and its derivatives, required in eq. (27), for orthorhombic and VTI media can be found in Appendix B. If the anisotropy vanishes, equations (27) reduce to the corresponding exact dynamic ray tracing equations for an S wave propagating in an isotropic medium.

We obtain the quantities $X_i^{(J)}$ required to calculate the geometrical spreading (see eq. 25), by specifying the initial conditions for the FODRT equations (27) for $\tau = 0$ in the following way:

$$ X_i^{(J)}(0) = 0, \quad Y_i^{(J)}(0) = Z_{i1} - p_i^0 v_{0k}^{[\mathcal{M}]} Z_{k1}, \quad (28) $$

where

$$ \begin{aligned} Z_{11} &= -\sin \phi_0, & Z_{21} &= \cos \phi_0, & Z_{31} &= 0, \\ Z_{12} &= -\cos \phi_0 \sin \delta_0, & Z_{22} &= -\sin \phi_0 \sin \delta_0, & Z_{32} &= \cos \delta_0 \end{aligned} \quad (29) $$

(see eq. (A.5) of Pšenčík & Teles 1996). In eq. (28), $v_0^{[\mathcal{M}]}$ denotes the components of the first-order ray-velocity vector of the common S wave, $v_i^{[\mathcal{M}]} = dx_i/d\tau$, at point S, at which $\tau = 0$. Symbols $\phi_0$ and $\delta_0$ in eq. (29) again denote the take-off angles introduced in eq. (21). Initial conditions of the FODRT equations for P waves, discussed by Pšenčík & Farra (2007), should be specified in the same way as the above initial conditions for common S waves (see eqs 28 and 29).

As in the case of P waves, the solution of the FODRT equations (27) can be effectively used for two-point ray tracing. Quantities $X_i^{(J)}$ can be used to estimate ray parameters $\gamma^{(J)}$ (take-off angles $\phi_0$ and $\delta_0$) from the distance between the termination points of rays shot at a receiver surface (surface containing receivers) and the receiver, at which the two-point ray should terminate.

## 5 SOLUTION OF THE TRANSPORT EQUATION ALONG THE COMMON RAY

In this section, we make use of an ansatz similar to ansatzs used by Pšenčík (1998) or Bakker (2002), in which both S waves are considered as one wave computed along the S-wave common ray. As in the mentioned papers, we require that the ansatz satisfies, asymptotically, the
elastodynamic equation (1). Throughout this section, we assume that the small parameter $\epsilon_2$, representing weak anisotropy perturbation, $\epsilon_2 \sim |\Delta a_{ijkl}|/|a_{ijkl}|$, is of the same order as the small parameter $\epsilon_1$, that is, of the order $O(\omega^{-1})$. This has important consequences. The diagonal elements of the matrix **B** are of the order $O(1)$, the off-diagonal elements are of the order $O(\omega^{-1})$. The difference $B_{11} - B_{22}$ is also of the order $O(\omega^{-1})$, which implies that the deviations of the eigenvalues $G^{[1]}$ and $G^{[2]}$ of the $S1$ and $S2$ waves from their mean value $G^{[M]}$, corresponding to the $S$-wave common ray, are of the order $O(\omega^{-1})$. Along the $S$-wave common ray, we seek the solution of the elastodynamic equation (1) in the form of eq. (2) with

$$ \mathbf{U}(x_m) = \mathcal{A}(x_m)\mathbf{e}^{[1]}(x_m, p_m) + \mathcal{B}(x_m)\mathbf{e}^{[2]}(x_m, p_m), \quad \tau(x_m) = \tau^{[M]}(x_m). \quad (30) $$

In eq. (30), $\mathcal{A}$ and $\mathcal{B}$ are scalar amplitude factors, $\tau^{[M]}(x_m)$ is the first-order traveltime along the $S$-wave common ray, obtained by solving FORT equations (16), $x_m$ are the coordinates along the common ray and $p_m$ are the components of the first-order slowness vector **p** along it. Vectors $\mathbf{e}^{[1]}$, $\mathbf{e}^{[2]}$ are the two mutually perpendicular unit vectors introduced after eq. (12). Our task is to find $\mathcal{A}$ and $\mathcal{B}$ so that (2) with (30) satisfies, asymptotically, elastodynamic equation (1).

Inserting (30) into (2), the result into (1) and neglecting the terms of order $O(1)$, $O(\omega^{-1})$ and less, we arrive at:

$$ (i\omega)^2(\mathbf{\Gamma} - \mathbf{I})(\mathcal{A}\mathbf{e}^{[1]} + \mathcal{B}\mathbf{e}^{[2]}) + i\omega\mathbf{M}(\mathcal{A}\mathbf{e}^{[1]} + \mathcal{B}\mathbf{e}^{[2]}) = 0. \quad (31) $$

In eq. (31), symbol **M** denotes the vectorial differential operator introduced in eq. (5), $\mathbf{\Gamma}$ denotes the Christoffel matrix, $\mathbf{\Gamma} = \mathbf{\Gamma}(x_m, p_m)$, where $p_m$ are the components of the first-order slowness vector.

To satisfy eq. (31), we require that the scalar product of its left-hand side with vectors $\mathbf{e}^{[m]}$ is zero for any of $m = 1, 2, 3$. We concentrate here only on equations resulting from scalar multiplication by $\mathbf{e}^{[1]}$ and $\mathbf{e}^{[2]}$. Multiplication by $\mathbf{e}^{[3]}$ yields the equation for the so-called additional component (Pšenčík 1998), which we do not study here. Using $\mathcal{C}\mathbf{e}^{[K]}$ instead of $\mathcal{A}\mathbf{e}^{[1]} + \mathcal{B}\mathbf{e}^{[2]}$ in the term with $(i\omega)^2$ in eq. (31), and multiplying the result by $\mathbf{e}^{[L]}$, we get

$$ (\Gamma_{ik} - \delta_{ik})\mathcal{C}e_i^{[K]}e_k^{[L]} = \mathcal{C}(B_{KL} - \delta_{KL}). \quad (32) $$

Symbol $B_{KL} = B_{KL}(x_m, p_m)$ denotes again elements of matrix **B** (see eq. 12). Doing the same with the term with $i\omega$ in (31) and taking into account eq. (5), we get, after some algebra

$$ M_i(C\mathbf{e}^{[K]})e_i^{[L]} = C\left[\rho^{-1}(\rho a_{ijkl}e_k^{[K]}e_i^{[L]}p_l)_{,j} + a_{ijkl}(e_{k,l}^{[K]}e_i^{[L]}p_j - e_k^{[K]}e_{i,j}^{[L]}p_l)\right] + a_{ijkl}e_k^{[K]}e_i^{[L]}(C_{,j}p_l + C_{,l}p_j). \quad (33) $$

As terms of order $O(1)$ and less are neglected in eq. (31), the term (32) can be calculated to $O(\omega^{-1})$ and the term (33) to $O(1)$. This means that eq. (32) can remain as it is. Eq. (33) can be, however, simplified. With accuracy to $O(1)$, we can rewrite the terms $a_{ijkl}e_k^{[K]}e_i^{[L]}p_l$ in the following way

$$ a_{ijkl}e_k^{[K]}e_i^{[L]}p_l = \delta_{KL}v_j^{[M]}. \quad (34) $$

The remaining term in (33) can be simplified if we take into account that for the required accuracy it is sufficient to substitute $a_{ijkl}$ by its counterpart specifying isotropic reference medium:

$$ a_{ijkl} \left( e_{k,l}^{[K]} e_i^{[L]} p_j - e_k^{[K]} e_{i,j}^{[L]} p_l \right) = 2(1 - \delta_{KL}) e_i^{[L]} \left( \frac{\mathrm{d}e_i^{[K]}}{\mathrm{d}\tau} \right). \quad (35) $$

Multiplication of eq. (31) successively by $\mathbf{e}^{[1]}$ and $\mathbf{e}^{[2]}$ and consideration of eqs (32)–(35) result in the following system of coupled linear differential equations:

$$ \begin{pmatrix} \mathrm{d}A/\mathrm{d}\tau \\ \mathrm{d}B/\mathrm{d}\tau \end{pmatrix} = -\frac{i\omega}{2} \begin{pmatrix} B_{11}-1 & B_{12} \\ B_{12} & B_{22}-1 \end{pmatrix} \begin{pmatrix} \mathcal{A} \\ \mathcal{B} \end{pmatrix} - \frac{1}{2} \begin{pmatrix} \rho^{-1} \mathrm{div}(\rho \mathbf{v}^{[M]}) & 2\mathbf{e}^{[1]} \cdot \frac{\mathrm{d}\mathbf{e}^{[2]}}{\mathrm{d}\tau} \\ 2\mathbf{e}^{[2]} \cdot \frac{\mathrm{d}\mathbf{e}^{[1]}}{\mathrm{d}\tau} & \rho^{-1} \mathrm{div}(\rho \mathbf{v}^{[M]}) \end{pmatrix} \begin{pmatrix} \mathcal{A} \\ \mathcal{B} \end{pmatrix}. \quad (36) $$

The first matrix term on the right-hand side of eq. (36) is responsible for frequency-dependent coupling. If this term is small, the studied $S$ waves tend to decouple. This can occur either for low frequencies or for the case of a medium close to isotropic, for which $B_{11}$ and $B_{22}$ approach 1 and $B_{12}$ approaches 0. Tendency to decoupling also occurs for high frequencies and/or stronger anisotropy. Coupling is strongest for intermediate frequencies. See more detailed discussions of these phenomena in Pšenčík (1998) and Bakker (2002).

The off-diagonal elements $\mathbf{e}^{[2]} \cdot \frac{\mathrm{d}\mathbf{e}^{[1]}}{\mathrm{d}\tau}$ of the second matrix term on the right-hand side of eq. (36) depend on the choice of vectors $\mathbf{e}^{[K]}$. We show later that the vectors $\mathbf{e}^{[K]}$ can be chosen so that $\mathbf{e}^{[2]} \cdot \frac{\mathrm{d}\mathbf{e}^{[1]}}{\mathrm{d}\tau} = 0$. The diagonal elements $\rho^{-1} \mathrm{div}(\rho \mathbf{v}^{[M]})$ are related to the geometrical spreading $\mathcal{L}^{[M]}$ given in (25). For details of this relation see sections 2.4.3. and 3.10.6 in Červený (2001). If we express the amplitude factors $\mathcal{A}$ and $\mathcal{B}$ in the form common in the ray method (see also (Pšenčík 1998, eq. 37)):

$$ \mathcal{A}(\tau) = \frac{A_0(\tau)}{[\rho(\tau)c^{[M]}(\tau)]^{1/2} \mathcal{L}^{[M]}(\tau)}, \quad \mathcal{B}(\tau) = \frac{B_0(\tau)}{[\rho(\tau)c^{[M]}(\tau)]^{1/2} \mathcal{L}^{[M]}(\tau)}, \quad (37) $$

the coupled equations (36) simplify to

$$ \begin{pmatrix} \mathrm{d}A_0/\mathrm{d}\tau \\ \mathrm{d}B_0/\mathrm{d}\tau \end{pmatrix} = -\frac{i\omega}{2} \begin{pmatrix} B_{11}-1 & B_{12} \\ B_{12} & B_{22}-1 \end{pmatrix} \begin{pmatrix} \mathcal{A}_0 \\ \mathcal{B}_0 \end{pmatrix} - \begin{pmatrix} 0 & \mathbf{e}^{[1]} \cdot \frac{\mathrm{d}\mathbf{e}^{[2]}}{\mathrm{d}\tau} \\ \mathbf{e}^{[2]} \cdot \frac{\mathrm{d}\mathbf{e}^{[1]}}{\mathrm{d}\tau} & 0 \end{pmatrix} \begin{pmatrix} \mathcal{A}_0 \\ \mathcal{B}_0 \end{pmatrix}. \quad (38) $$

In (37), $c^{[M]}$ and $\mathcal{L}^{[M]}$ are the phase velocity (see eq. 17) and the geometrical spreading (see eq. 25), associated with the first-order mean value $G^{[M]}$ along the common ray of the $S$ wave. We can see that coefficients $A_0$ and $B_0$ in (38) are frequency-dependent. To solve the coupled system (38), it is necessary to specify its initial conditions and how to calculate vectors $\mathbf{e}^{[K]}$ along the ray.
The initial conditions can be specified in the very same way as in Pšenčík (1998, eq. 39). Specifically, if the wavefield is generated by a point force F acting at time $\tau = 0$, the comparison of eq. (37) with the corresponding expressions of Pšenčík & Teles (1996) yields the initial conditions for eq. (38):

$$A_0(0) = \frac{e_k^{[1]}(0)F_k}{4\pi(\rho_0 c_0^{[M]})^{1/2}}, \quad B_0(0) = \frac{e_k^{[2]}(0)F_k}{4\pi(\rho_0 c_0^{[M]})^{1/2}} \quad (39)$$

Here $c_0^{[M]}$ again denotes the first-order approximation of the S-wave common phase velocity at the source, and $\rho_0$ denotes the density at the same point.

If we specify the point force at point S as a unit vector oriented successively along all coordinate axes $x_n$, that is, $F_k = \delta_{kn}$, solve eq. (38) with the initial conditions (39) along the common ray from point S to point R, and insert the results into eq. (37), we get frequency-dependent coefficients $A_n(\tau) = A(\tau)$ and $B_n(\tau) = B(\tau)$. These coefficients can be used in the expression for the first-order common S-wave Green's function:

$$G_{in}^{[M]}(R, S, \omega) = \left[ A_n(\tau)e_i^{[1]}(\tau) + B_n(\tau)e_i^{[2]}(\tau) \right] \exp \left[ iT^G(R, S) + i\omega\tau^{[M]}(R, S) \right]. \quad (40)$$

Here $T^G(R, S)$ denotes the complete phase shift due to the caustics along the S-wave common ray from S to R. All other quantities have been defined above. To increase accuracy of the traveltime $\tau^{[M]}$, we can substitute it by $\tau^{[M]} + \Delta\tau^{[M]}$ (see eq. 24).

There are several possibilities how to calculate vectors $e^{[K]}$; vector $e^{[3]} = n$ is calculated automatically along the common ray. One possibility how to calculate $e^{[K]}$ is to specify them as the vectorial basis of the so-called wave front orthonormal coordinate system (see Červený 2001, section 4.2.2; Klimeš 2006a). Červený (2001, eq. 4.2.17) gives differential equations, which in our notation read

$$\frac{de_i^{[K]}}{d\tau} = -(c^{[M]})^2 \left( e_k^{[K]} \frac{dp_k}{d\tau} \right) p_i. \quad (41)$$

Solution of eq. (41) along a selected ray yields vectors $e^{[K]}$ at any point of the ray. Note that vectors $e^{[K]}$ are perpendicular to the first-order slowness vector $p$, and not to the ray as in the QI approximation of Pšenčík (1998). Note also that vectors $e^{[K]}$ are calculated along an S-wave common ray in the studied weakly anisotropic medium and not in a reference isotropic medium as in the QI approximation.

The use of the basis vectors of the wave front orthonormal coordinate system (41) yields $e^{[K]} \cdot de^{[L]}/d\tau = 0$ for $K \neq L$ and it thus leads to the coupled system of differential equations whose form reminds form of the QI coupled equations of Pšenčík (1998, eq. 38):

$$\begin{pmatrix} dA_0/d\tau \\ dB_0/d\tau \end{pmatrix} = -\frac{i\omega}{2} \begin{pmatrix} B_{11} - 1 & B_{12} \\ B_{12} & B_{22} - 1 \end{pmatrix} \begin{pmatrix} A_0 \\ B_0 \end{pmatrix}. \quad (42)$$

Matrix $B(x_m, p_m)$ used in (42) can be obtained by rotation (13) from matrix $\tilde{B}(x_m, p_m)$ given in Appendix A.

Vectors $e^{[K]}$ define zero-order polarization plane of the common S wave. For greater accuracy of the formula (40), we can substitute them by their first-order counterparts $f^{[K]}$,

$$f^{[K]} = e^{[K]} - \frac{B_{K3}(x_m, n_m)}{V_P^2 - V_S^2} e^{[3]}. \quad (43)$$

Vectors $f^{[K]}$ are situated in the plane perpendicular to the first-order P-wave eigenvector $f^{[3]}$. Vector $f^{[3]}$ is given by (see Farra & Pšenčík 2003):

$$f^{[3]} = \frac{B_{13}(x_m, n_m)}{V_P^2 - V_S^2} e^{[1]} + \frac{B_{23}(x_m, n_m)}{V_P^2 - V_S^2} e^{[2]} + e^{[3]}. \quad (44)$$

Symbols $V_P$ and $V_S$ in eqs (43) and (44) have exactly the same meaning as in eq. (24). Elements $B_{13}$ and $B_{23}$ of the matrix $B(x_m, n_m)$ can be obtained by rotation (13) from $\tilde{B}_{13}$ and $\tilde{B}_{23}$ given in (A4). Since $B_{K3}$ are $O(\omega^{-1})$, the differences between vectors $f^{[K]}$ and vectors $e^{[K]}$ are $O(\omega^{-1})$ too.

# 6 DISCUSSION AND CONCLUSIONS

The preceding sections describe an approximate procedure for computing S waves in inhomogeneous weakly anisotropic media. In such media, the two S waves are coupled and must, therefore, be calculated together. We have used common ray tracing and dynamic ray tracing (Bakker 2002; Klimeš 2006a) with the common ray situated in the studied weakly anisotropic medium. We have combined the S-wave common ray tracing with the perturbation approach, which we applied recently to P waves (Pšenčík & Farra 2005, 2007) and obtained an approximate procedure, which has similar advantages as FORT and FODRT for P waves and differs significantly from the standard QI approach (Kravtsov 1968; Pšenčík 1998). The most important differences from the QI approach are:

(1) The common ray is traced in the studied inhomogeneous, weakly anisotropic medium, not in a reference isotropic medium.

(2) The wave front orthonormal basis related to the common ray in the studied inhomogeneous, weakly anisotropic medium is used instead of the basis related to the common ray in the reference isotropic medium.

(3) The accuracy of traveltimes can be increased by considering the second-order traveltime correction.
(4) The common geometrical spreading is calculated along the common ray in the studied inhomogeneous, weakly anisotropic medium, not along the common ray in the reference isotropic medium.

The most interesting features and advantages of the first-order common ray and dynamic ray tracing of coupled S waves in inhomogeneous weakly anisotropic media are:

(1) The FORT and FODRT equations are expressed in terms of the WA parameters, which represent a more natural description of a weakly anisotropic medium than standard elastic moduli.

(2) The proposed S-wave common ray tracing procedure is regular everywhere. It avoids problems in S-wave singular regions, known from tracing separate S waves in anisotropic media.

(3) The FORT and FODRT equations have a simpler structure than the exact ones, which leads to a reduction in the number of algebraic operations involved.

(4) In the most general case, the FORT and FODRT equations for S waves depend only on 15 independent coefficients (which differ from the coefficients of P waves).

(5) The FORT and FODRT equations allow S waves to be treated separately from P waves, as in isotropic media.

(6) The accuracy of the traveltimes calculated along S-wave common rays can be easily increased by introducing the second-order traveltime correction (eq. 24).

(7) The procedure is applicable to inhomogeneous, weakly anisotropic media of arbitrary symmetry.

(8) The FORT and FODRT equations derived in this paper are applicable to inhomogeneous isotropic as well as anisotropic media. In isotropic media they yield exact rays and exact dynamic ray tracing results for S waves; in anisotropic media, first-order approximations along S-wave common rays.

(9) The procedure described above can be simply generalized for laterally varying, layered, weakly anisotropic structures.

(10) FORT and FODRT for inhomogeneous weakly anisotropic media of higher symmetry (orthorhombic, TI) represent a natural substitute for ray tracing and dynamic ray tracing procedures in routine processing codes.

## ACKNOWLEDGMENTS

The authors are grateful to Joe Dellinger for his review and stimulating comments. Comments by the editor, J-M.Kendall, an anonymous referee and our colleagues E.Iversen and P.Jílek are also appreciated. Special thanks belong to Peter Bakker for his contribution to this work. A substantial part of this work was done during IP's stay at the IPG Paris at the invitation of the IPGP. We are grateful to the Consortium Project 'Seismic waves in complex 3-D structures' (SW3D) and Research Projects 205/05/2182 and 205/08/0332 of the Grant Agency of the Czech Republic for support.

## REFERENCES

Bakker, P., 2002. Coupled anisotropic shear-wave ray tracing in situations where associated slowness sheets are almost tangent, *Pure appl. Geophys.*, **159**, 1403–1417.

Bulant, P. & Klimeš, L., 2002. Numerical algorithm of the coupling ray theory in weakly anisotropic media, *Pure appl. Geophys.*, **159**, 1419–1435.

Červený, V., 2001. *Seismic Ray Theory*, Cambridge University Press, Cambridge.

Coates, R.T. & Chapman, C.H., 1990. Quasi-shear wave coupling in weakly anisotropic 3-D media, *Geophys. J. Int.*, **103**, 301–320.

Farra, V., 2005. First-order ray tracing for qS waves in inhomogeneous weakly anisotropic media, *Geophys. J. Int.*, **161**, 309–324.

Farra, V. & Pšenčík, I., 2003. Properties of the zero-, first- and higher-order approximations of attributes of elastic waves in weakly anisotropic media, *J. acoust. Soc. Am.*, **114**, 1366–1378.

Klimeš, L., 2006a. Ray-centred coordinate systems in anisotropic media, *Stud. Geophys. Geod.*, **50**, 431–447.

Klimeš, L., 2006b. Common-ray tracing and dynamic ray tracing for S waves in a smooth elastic anisotropic medium, *Stud. Geophys. Geod.*, **50**, 449–461.

Kravtsov, Yu. A., 1968. 'Quasiisotropic' approximation to geometrical optics, *Dokl. AN SSSR*, **183**(1), 74–77 (in Russian).

Pšenčík, I., 1998. Green's functions for inhomogeneous weakly anisotropic media, *Geophys. J. Int.*, **135**, 279–288.

Pšenčík, I. & Farra, V., 2005. First-order ray tracing for qP waves in inhomogeneous weakly anisotropic media, *Geophysics*, **70**, D65–D75.

Pšenčík, I. & Farra, V., 2007. First-order P-wave ray synthetic seismograms in inhomogeneous weakly anisotropic media. *Geophys. J. Int.*, **170**, 1243–1252.

Pšenčík, I. & Teles, T.N., 1996. Point source radiation in inhomogeneous anisotropic structures, *Pure appl. Geophys.*, **148**, 591–623.

Thomson, C.J., Kendall, J.-M. & Guest, W.S., 1992. Geometrical theory of shear wave splitting: corrections to ray theory for interference in isotropic/anisotropic transitions. *Geophys. J. Int.*, **108**, 339–363.

## APPENDIX A: WA PARAMETERS AND ELEMENTS OF MATRIX B

If we denote by $\alpha$ and $\beta$ the P- and S-wave velocities of a reference isotropic medium, we can introduce 21 WA parameters in a way slightly different from that of Farra & Pšenčík (2003):

$$ \epsilon_x = \frac{A_{11} - \alpha^2}{2\alpha^2}, \quad \epsilon_y = \frac{A_{22} - \alpha^2}{2\alpha^2}, \quad \epsilon_z = \frac{A_{33} - \alpha^2}{2\alpha^2}, $$

$$ \delta_x = \frac{A_{23} + 2A_{44} - \alpha^2}{\alpha^2}, \quad \delta_y = \frac{A_{13} + 2A_{55} - \alpha^2}{\alpha^2}, \quad \delta_z = \frac{A_{12} + 2A_{66} - \alpha^2}{\alpha^2}, $$
$$
\begin{align*}
\chi_x &= \frac{A_{14} + 2A_{56}}{\alpha^2}, & \chi_y &= \frac{A_{25} + 2A_{46}}{\alpha^2}, & \chi_z &= \frac{A_{36} + 2A_{45}}{\alpha^2}, \\
\epsilon_{15} &= \frac{A_{15}}{\alpha^2}, & \epsilon_{16} &= \frac{A_{16}}{\alpha^2}, & \epsilon_{24} &= \frac{A_{24}}{\alpha^2}, & \epsilon_{26} &= \frac{A_{26}}{\alpha^2}, & \epsilon_{34} &= \frac{A_{34}}{\alpha^2}, & \epsilon_{35} &= \frac{A_{35}}{\alpha^2}, \\
\epsilon_{46} &= \frac{A_{46}}{\beta^2}, & \epsilon_{56} &= \frac{A_{56}}{\beta^2}, & \epsilon_{45} &= \frac{A_{45}}{\beta^2}, & \gamma_x &= \frac{A_{44} - \beta^2}{2\beta^2}, & \gamma_y &= \frac{A_{55} - \beta^2}{2\beta^2}, & \gamma_z &= \frac{A_{66} - \beta^2}{2\beta^2}. \tag{A1}
\end{align*}
$$

Please, note that parameters $\delta_x, \delta_y$ and $\gamma_x, \gamma_y$ are defined here differently from Farra & Pšenčik (2003) and Pšenčik & Farra (2005, 2007). In the text, the following variables are also used.

$$ \eta_x = \delta_x - \epsilon_y - \epsilon_z, \quad \eta_y = \delta_y - \epsilon_x - \epsilon_z, \quad \eta_z = \delta_z - \epsilon_x - \epsilon_y. \tag{A2} $$

The matrix $\bar{\mathbf{B}}(x_m, p_m)$ used in eq. (13) is related to the matrix $\bar{\mathbf{B}}(x_m, n_m)$ by the relation

$$ \bar{\mathbf{B}}(x_m, p_m) = (c^{[\mathcal{M}]})^{-2} \bar{\mathbf{B}}(x_m, n_m). \tag{A3} $$

The elements of $\bar{\mathbf{B}}(x_m, n_m)$ read:

$$
\begin{align*}
\bar{B}_{11} &= \beta^2 + 2\alpha^2 D^{-2} \left\{ n_3^3 \left[ (\epsilon_{15} - \epsilon_{35})n_1^3 + (\epsilon_{24} - \epsilon_{34})n_2^3 + (\chi_x - \epsilon_{34})n_1^2 n_2 + (\chi_y - \epsilon_{35})n_1 n_2^2 \right] \right. \\
&\qquad \left. + n_3^2 \left[ (\eta_z - \eta_x - \eta_y)n_1^2 n_2^2 - \eta_y n_1^4 - \eta_x n_2^4 + 2(\epsilon_{16} - \chi_z)n_1^3 n_2 + 2(\epsilon_{26} - \chi_z)n_1 n_2^3 \right] \right. \\
&\qquad \left. + n_3 \left[ (\epsilon_{35} - \epsilon_{15})n_1^5 + (\epsilon_{34} - \epsilon_{24})n_2^5 + (\epsilon_{34} - \chi_x)n_1^4 n_2 + (\epsilon_{35} - \chi_y)n_1 n_2^4 \right. \right. \\
&\qquad \left. \left. + (2\epsilon_{35} - \chi_y - \epsilon_{15})n_1^3 n_2^2 + (2\epsilon_{34} - \chi_x - \epsilon_{24})n_1^2 n_2^3 \right] \right\} + 2\beta^2 D^{-2} (\gamma_x n_2^2 + \gamma_y n_1^2 + \epsilon_{45} n_1 n_2),
\end{align*}
$$

$$
\begin{align*}
\bar{B}_{12} &= \alpha^2 D^{-2} \left\{ n_3^2 \left[ (\chi_x - \epsilon_{34})n_1^3 - (\chi_y - \epsilon_{35})n_2^3 + (\epsilon_{35} + \chi_y - 2\epsilon_{15})n_1^2 n_2 + (2\epsilon_{24} - \chi_x - \epsilon_{34})n_1 n_2^2 \right] \right. \\
&\qquad \left. + n_3 \left[ (\epsilon_{16} - \chi_x)n_1^4 - (\epsilon_{26} - \chi_z)n_2^4 + (\eta_y - \eta_x + \eta_z)n_1^3 n_2 + (\eta_y - \eta_x - \eta_z)n_1 n_2^3 + 3(\epsilon_{26} - \epsilon_{16})n_1^2 n_2^2 \right] \right. \\
&\qquad \left. + (\epsilon_{15} - \chi_y)n_1^4 n_2 + (\chi_x - \epsilon_{24})n_1 n_2^4 + (\epsilon_{15} - \chi_y)n_1^2 n_2^3 + (\chi_x - \epsilon_{24})n_1^3 n_2^2 \right\} \\
&\qquad + \beta^2 D^{-2} \left\{ n_3 \left[ 2(\gamma_x - \gamma_y)n_1 n_2 + \epsilon_{45} n_1^2 - \epsilon_{45} n_2^2 \right] + \epsilon_{46} n_1^2 - \epsilon_{56} n_1^3 + \epsilon_{46} n_1^2 n_2 - \epsilon_{56} n_1 n_2^2 \right\},
\end{align*}
$$

$$
\begin{align*}
\bar{B}_{22} &= \beta^2 + 2\beta^2 D^{-2} \left[ n_3^2 (\gamma_x n_1^2 + \gamma_y n_2^2 - \epsilon_{45} n_1 n_2) + n_3 (\epsilon_{46} n_1^3 + \epsilon_{56} n_2^3 + \epsilon_{56} n_1^2 n_2 + \epsilon_{46} n_1 n_2^2) \right] \\
&\quad + 2\alpha^2 D^{-2} \left\{ n_3 \left[ (\epsilon_{24} - \chi_x)n_1^2 n_2 + (\epsilon_{15} - \chi_y)n_1 n_2^2 \right] - \eta_z n_1^2 n_2^2 + (\epsilon_{26} - \epsilon_{16})n_1^3 n_2 - (\epsilon_{26} - \epsilon_{16})n_1 n_2^3 \right\} + 2\beta^2 \gamma_z (n_1^2 + n_2^2),
\end{align*}
$$

$$
\begin{align*}
\bar{B}_{13} &= \alpha^2 D^{-1} \left\{ n_3^3 (\epsilon_{34} n_2 + \epsilon_{35} n_1) + n_3^3 (\eta_y n_1^2 + \eta_x n_2^2 + 2\chi_z n_1 n_2) \right. \\
&\qquad + n_3^2 \left[ (4\chi_x - 3\epsilon_{34})n_1^2 n_2 + (4\chi_y - 3\epsilon_{35})n_1 n_2^2 + (4\epsilon_{15} - 3\epsilon_{35})n_1^3 + (4\epsilon_{24} - 3\epsilon_{34})n_2^3 \right] \\
&\qquad + n_3 \left[ (2\eta_z - \eta_x - \eta_y)n_1^2 n_2^2 + 2(2\epsilon_{16} - \chi_z)n_1^3 n_2 + 2(2\epsilon_{26} - \chi_z)n_1 n_2^3 \right. \\
&\qquad \left. - \eta_y n_1^4 - \eta_x n_2^4 + (\epsilon_x - \epsilon_z)n_1^2 + (\epsilon_y - \epsilon_z)n_2^2 \right] - \chi_x n_1^2 n_2 - \chi_y n_1 n_2^2 - \epsilon_{15} n_1^3 - \epsilon_{24} n_2^3 \},
\end{align*}
$$

$$
\begin{align*}
\bar{B}_{23} &= \alpha^2 D^{-1} \left\{ n_3^3 (\epsilon_{34} n_1 - \epsilon_{35} n_2) + n_3^2 \left[ (\eta_x - \eta_y)n_1 n_2 + \chi_z n_1^2 - \chi_z n_2^2 \right] + n_3 \left[ (2\chi_y - 3\epsilon_{15})n_1^2 n_2 - (2\chi_x - 3\epsilon_{24})n_1 n_2^2 + \chi_x n_1^3 - \chi_y n_1^3 \right] \right. \\
&\qquad \left. + \eta_z n_1^3 n_2 - \eta_z n_1 n_2^3 + 3(\epsilon_{26} - \epsilon_{16})n_1^2 n_2^2 + \epsilon_{16} n_1^4 - \epsilon_{26} n_1^4 + (\epsilon_y - \epsilon_x)n_1 n_2 \right\},
\end{align*}
$$

$$
\begin{align}
\bar{B}_{33} &= \alpha^2 + 2\alpha^2 \left[ 2n_3^3(\epsilon_{34}n_2 + \epsilon_{35}n_1) + n_3^2(\eta_y n_1^2 + \eta_x n_2^2 + 2\chi_z n_1 n_2 + \epsilon_z) \right. \nonumber \\
&\qquad \left. + 2n_3(\chi_x n_1^2 n_2 + \chi_y n_1 n_2^2 + \epsilon_{15} n_1^3 + \epsilon_{24} n_2^3) + \epsilon_x n_1^2 + \epsilon_y n_2^2 + \eta_z n_1^2 n_2^2 + 2\epsilon_{16} n_1^3 n_2 + 2\epsilon_{26} n_1 n_2^3 \right]. \tag{A4}
\end{align}
$$

The expressions in eq. (A4) correspond to the following choice of vectors $\mathbf{e}^{[i]}$:

$$ \bar{\mathbf{e}}^{[1]} \equiv D^{-1}(n_1 n_3, n_2 n_3, n_3^2 - 1), \quad \bar{\mathbf{e}}^{[2]} \equiv D^{-1}(-n_2, n_1, 0), \quad \bar{\mathbf{e}}^{[3]} \equiv \mathbf{n} \equiv (n_1, n_2, n_3), \tag{A5} $$

where

$$ D = (n_1^2 + n_2^2)^{1/2}, \qquad n_1^2 + n_2^2 + n_3^2 = 1. \tag{A6} $$

Symbols $n_i$ denote the components of unit vector $\mathbf{n}$ specifying direction of the slowness vector.

## APPENDIX B: EXPRESSIONS FOR $G^{[\mathcal{M}]}$ AND ITS FIRST- AND SECOND-ORDER DERIVATIVES IN ORTHORHOMBIC AND VTI MEDIA

The formulae for the ray tracing and dynamic ray tracing along an S-wave common ray simplify considerably in media with higher symmetry anisotropy. Below we present formulae for the S-wave first-order mean value $G^{[\mathcal{M}]}$ and its first- and second-order derivatives, which are the basic elements on the right-hand sides of the ray tracing (16) and dynamic ray tracing (27) equations. We present them for the two most frequently used types of media, media of orthorhombic and VTI symmetry. We use the following notations,

$$ p_{12} = p_1^2 + p_2^2, \quad p_{13} = p_1^2 + p_3^2, \quad p_{23} = p_2^2 + p_3^2. \tag{B1} $$
## B.1 Orthorhombic symmetry

In a medium of orthorhombic symmetry that has planes of symmetry coinciding with the coordinate planes, we have

$$ \chi_x = \chi_y = \chi_z = \epsilon_{15} = \epsilon_{16} = \epsilon_{24} = \epsilon_{26} = \epsilon_{34} = \epsilon_{35} = \epsilon_{45} = \epsilon_{46} = \epsilon_{56} = 0. \qquad (B2) $$

Formula (19) for the first-order eigenvalue $G^{[M]}$ then reduces to

$$ G^{[M]} = \beta^2(p_i p_i + \gamma_x p_{23} + \gamma_y p_{13} + \gamma_z p_{12}) - \alpha^2(p_k p_k)^{-1}(\eta_x p_2^2 p_3^2 + \eta_y p_1^2 p_3^2 + \eta_z p_1^2 p_2^2). \qquad (B3) $$

The formulae for the first derivatives of $G^{[M]}$, which appear on the right-hand sides of eq. (16), read:

$$ \frac{\partial G^{[M]}}{\partial p_1} = 2p_1 \{ \beta^2(1 + \gamma_y + \gamma_z) - \alpha^2(p_k p_k)^{-2} [p_{23}(\eta_y p_3^2 + \eta_z p_2^2) - \eta_x p_2^2 p_3^2] \}, $$

$$ \frac{\partial G^{[M]}}{\partial p_2} = 2p_2 \{ \beta^2(1 + \gamma_x + \gamma_z) - \alpha^2(p_k p_k)^{-2} [p_{13}(\eta_x p_3^2 + \eta_z p_1^2) - \eta_y p_1^2 p_3^2] \}, $$

$$ \frac{\partial G^{[M]}}{\partial p_3} = 2p_3 \{ \beta^2(1 + \gamma_x + \gamma_y) - \alpha^2(p_k p_k)^{-2} [p_{12}(\eta_x p_2^2 + \eta_y p_1^2) - \eta_z p_1^2 p_2^2] \}, $$

$$ \frac{\partial G^{[M]}}{\partial x_i} = \beta^2 [(\gamma_{x,i} + \gamma_{z,i})p_2^2 + (\gamma_{y,i} + \gamma_{z,i})p_1^2 + (\gamma_{x,i} + \gamma_{y,i})p_3^2] - \alpha^2(p_k p_k)^{-1}(\eta_{x,i} p_2^2 p_3^2 + \eta_{y,i} p_1^2 p_3^2 + \eta_{z,i} p_1^2 p_2^2). \qquad (B4) $$

The formulae for the second derivatives of $G^{[M]}$, which appear on the right-hand sides of eq. (27), read:

$$ \frac{\partial^2 G^{[M]}}{\partial p_1^2} = 2\beta^2(1 + \gamma_y + \gamma_z) - 2\alpha^2(p_k p_k)^{-3}[(\eta_y p_3^2 + \eta_z p_2^2)p_{23} - \eta_x p_2^2 p_3^2](p_{23} - 3p_1^2), $$

$$ \frac{\partial^2 G^{[M]}}{\partial p_1 \partial p_2} = -4\alpha^2 p_1 p_2 (p_k p_k)^{-3} [\eta_z (p_{23} p_{13} + p_1^2 p_2^2) + \eta_x p_3^2 (p_2^2 - p_{13}) + \eta_y p_3^2 (p_1^2 - p_{23})], $$

$$ \frac{\partial^2 G^{[M]}}{\partial p_1 \partial p_3} = -4\alpha^2 p_1 p_3 (p_k p_k)^{-3} [\eta_y (p_{23} p_{12} + p_1^2 p_3^2) + \eta_x p_2^2 (p_3^2 - p_{12}) + \eta_z p_2^2 (p_1^2 - p_{23})], $$

$$ \frac{\partial^2 G^{[M]}}{\partial p_2^2} = 2\beta^2(1 + \gamma_x + \gamma_z) - 2\alpha^2(p_k p_k)^{-3}[(\eta_x p_3^2 + \eta_z p_1^2)p_{13} - \eta_y p_1^2 p_3^2](p_{13} - 3p_2^2), $$

$$ \frac{\partial^2 G^{[M]}}{\partial p_2 \partial p_3} = -4\alpha^2 p_2 p_3 (p_k p_k)^{-3} [\eta_x (p_{13} p_{12} + p_2^2 p_3^2) + \eta_y p_1^2 (p_3^2 - p_{12}) + \eta_z p_1^2 (p_2^2 - p_{13})], $$

$$ \frac{\partial^2 G^{[M]}}{\partial p_3^2} = 2\beta^2(1 + \gamma_x + \gamma_y) - 2\alpha^2(p_k p_k)^{-3}[(\eta_x p_2^2 + \eta_y p_1^2)p_{12} - \eta_z p_1^2 p_2^2](p_{12} - 3p_3^2), $$

$$ \frac{\partial^2 G^{[M]}}{\partial x_i \partial p_1} = 2p_1 \{ \beta^2(\gamma_{y,i} + \gamma_{z,i}) - \alpha^2(p_k p_k)^{-2} [(\eta_{y,i} p_3^2 + \eta_{z,i} p_2^2)p_{23} - \eta_{x,i} p_2^2 p_3^2] \}, $$

$$ \frac{\partial^2 G^{[M]}}{\partial x_i \partial p_2} = 2p_2 \{ \beta^2(\gamma_{x,i} + \gamma_{z,i}) - \alpha^2(p_k p_k)^{-2} [(\eta_{x,i} p_3^2 + \eta_{z,i} p_1^2)p_{13} - \eta_{y,i} p_1^2 p_3^2] \}, $$

$$ \frac{\partial^2 G^{[M]}}{\partial x_i \partial p_3} = 2p_3 \{ \beta^2(\gamma_{x,i} + \gamma_{y,i}) - \alpha^2(p_k p_k)^{-2} [(\eta_{y,i} p_1^2 + \eta_{x,i} p_2^2)p_{12} - \eta_{z,i} p_1^2 p_2^2] \}, $$

$$ \frac{\partial^2 G^{[M]}}{\partial x_i \partial x_j} = \beta^2(\gamma_{x,ij} p_{23} + \gamma_{y,ij} p_{13} + \gamma_{z,ij} p_{12}) - \alpha^2(p_k p_k)^{-1}(\eta_{x,ij} p_2^2 p_3^2 + \eta_{y,ij} p_1^2 p_3^2 + \eta_{z,ij} p_1^2 p_2^2). \qquad (B5) $$

The initial conditions for ray-tracing equations (16) and dynamic ray-tracing equations (27) are given by eqs (20) and (28). The square of the first-order S-wave common phase velocity $c^{[M]}$ appearing in the initial conditions reads:

$$ (c^{[M]})^2 = \beta^2 [1 + \gamma_x (n_2^2 + n_3^2) + \gamma_y (n_1^2 + n_3^2) + \gamma_z (n_1^2 + n_2^2)] - \alpha^2 (\eta_x n_2^2 n_3^2 + \eta_y n_1^2 n_3^2 + \eta_z n_1^2 n_2^2). \qquad (B6) $$

## B.2 VTI symmetry

In a medium of a hexagonal symmetry with a vertical axis of symmetry, we have, in addition to conditions (B1),

$$ \epsilon_x = \epsilon_y = \frac{1}{2}\delta_z, \quad \delta_x = \delta_y, \quad \eta_x = \eta_y, \quad \eta_z = 0, \quad \gamma_x = \gamma_y. \qquad (B7) $$

Formula (19) for the first-order eigenvalue $G^{[M]}$ reads in this case:

$$ G^{[M]} = \beta^2 [p_i p_i + (\gamma_x + \gamma_z)p_{12} + 2\gamma_x p_3^2] - \alpha^2 (p_k p_k)^{-1} \eta_x p_{12} p_3^2. \qquad (B8) $$
The formulae for the first derivatives of $G^{[\mathcal{M}]}$, which appear on the right-hand sides of eq. (16), simplify to:

$$ 
\begin{aligned}
\frac{\partial G^{[\mathcal{M}]} }{\partial p_1} &= 2p_1 [\beta^2(1 + \gamma_x + \gamma_z) - \alpha^2(p_k p_k)^{-2} \eta_x p_3^4], \\
\frac{\partial G^{[\mathcal{M}]} }{\partial p_2} &= 2p_2 [\beta^2(1 + \gamma_x + \gamma_z) - \alpha^2(p_k p_k)^{-2} \eta_x p_3^4], \\
\frac{\partial G^{[\mathcal{M}]} }{\partial p_3} &= 2p_3 [\beta^2(1 + 2\gamma_x) - \alpha^2(p_k p_k)^{-2} \eta_x p_{12}^2], \\
\frac{\partial G^{[\mathcal{M}]} }{\partial x_i} &= \beta^2 [(\gamma_{x,i} + \gamma_{z,i})p_{12} + 2\gamma_{x,i} p_3^2] - \alpha^2(p_k p_k)^{-1} \eta_{x,i} p_{12} p_3^2.
\end{aligned}
\qquad (B9) $$

The formulae for the second derivatives of $G^{[\mathcal{M}]}$, which appear on the right-hand sides of eq. (27), read:

$$ 
\begin{aligned}
\frac{\partial^2 G^{[\mathcal{M}]} }{\partial p_1^2} &= 2\beta^2(1 + \gamma_x + \gamma_z) - 2\alpha^2(p_k p_k)^{-3} \eta_x p_3^4 (p_{23} - 3p_1^2), \\
\frac{\partial^2 G^{[\mathcal{M}]} }{\partial p_1 \partial p_2} &= 8\alpha^2 \eta_x p_1 p_2 p_3^4 (p_k p_k)^{-3}, \qquad \frac{\partial^2 G^{[\mathcal{M}]} }{\partial p_1 \partial p_3} = -8\alpha^2 \eta_x p_1 p_3^3 (p_k p_k)^{-3} p_{12}, \\
\frac{\partial^2 G^{[\mathcal{M}]} }{\partial p_2^2} &= 2\beta^2(1 + \gamma_x + \gamma_z) - 2\alpha^2(p_k p_k)^{-3} \eta_x p_3^4 (p_{13} - 3p_2^2), \\
\frac{\partial^2 G^{[\mathcal{M}]} }{\partial p_2 \partial p_3} &= -8\alpha^2 \eta_x p_2 p_3^3 (p_k p_k)^{-3} p_{12}, \\
\frac{\partial^2 G^{[\mathcal{M}]} }{\partial p_3^2} &= 2\beta^2(1 + 2\gamma_x) - 2\alpha^2(p_k p_k)^{-3} \eta_x p_{12}^2 (p_{12} - 3p_3^2), \\
\frac{\partial^2 G^{[\mathcal{M}]} }{\partial x_i \partial p_1} &= 2p_1 [\beta^2(\gamma_{x,i} + \gamma_{z,i}) - \alpha^2(p_k p_k)^{-2} \eta_{x,i} p_3^4], \\
\frac{\partial^2 G^{[\mathcal{M}]} }{\partial x_i \partial p_2} &= 2p_2 [\beta^2(\gamma_{x,i} + \gamma_{z,i}) - \alpha^2(p_k p_k)^{-2} \eta_{x,i} p_3^4], \\
\frac{\partial^2 G^{[\mathcal{M}]} }{\partial x_i \partial p_3} &= 2p_3 [2\beta^2 \gamma_{x,i} - \alpha^2(p_k p_k)^{-2} \eta_{x,i} p_{12}^2], \\
\frac{\partial^2 G^{[\mathcal{M}]} }{\partial x_i \partial x_j} &= \beta^2 [(\gamma_{x,ij} + \gamma_{z,ij})p_{12} + 2\gamma_{x,ij} p_3^2] - \alpha^2(p_k p_k)^{-1} \eta_{x,ij} p_{12} p_3^2.
\end{aligned}
\qquad (B10) $$

The initial conditions for ray-tracing equations (16) and dynamic ray-tracing equations (27) are again given by eqs (20) and (28). The square of the first-order S-wave common phase velocity $c^{[\mathcal{M}]}$, appearing in the initial conditions, now has the following simple form:

$$ (c^{[\mathcal{M}]} )^2 = \beta^2 [1 + (\gamma_x + \gamma_z)(n_1^2 + n_2^2) + 2\gamma_x n_3^2] - \alpha^2 \eta_x (n_1^2 + n_2^2) n_3^2. \qquad (B11) $$