# Coupled S waves in inhomogeneous weakly anisotropic media using first-order ray tracing

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

► To cite this version:

Véronique Farra, Ivan Pšenčík. Coupled S waves in inhomogeneous weakly anisotropic media using first-order ray tracing. Geophysical Journal International, 2010, 180, pp.405-417. 10.1111/j.1365-246X.2009.04423.x. insu-03605350

HAL Id: insu-03605350

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

Submitted on 11 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
# Coupled *S* waves in inhomogeneous weakly anisotropic media using first-order ray tracing

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

Accepted 2009 October 19. Received 2009 October 19; in original form 2009 April 3

## SUMMARY

We present a modification of our recently proposed approximate procedure for computing coupled *S* waves in inhomogeneous weakly anisotropic media. The new procedure can be used to compute *S* waves propagating in smooth inhomogeneous isotropic or anisotropic media. In isotropic media, it reduces to standard ray theory procedure for *S* waves. In anisotropic media, it can be used to study coupled as well as decoupled *S* waves. As the previous procedure, the new one is also based on the approximately computed common *S*-wave ray. First-order ray tracing and dynamic ray tracing, originally developed for computations of *P*-wave fields, is used to compute common *S*-wave rays and the dynamic ray tracing along them. The principal difference between the previous and new procedure consists in implicit incorporation of the second-order common *S*-wave traveltime correction and more accurate estimate of traveltime difference in the modified coupling equations. This leads to a substantial increase of accuracy of the coupling equations, which are solved along the common ray to evaluate *S*-wave amplitudes. The new coupling equations provide, first of all, more accurate traveltimes, but they also allow modelling of decoupled *S* waves, which could hardly be done with the original coupling equations. There is no need for a choice of a reference medium. The reference medium is determined uniquely from the actual medium.

The new procedure has all the advantages of the previous procedure. Among the basic advantages is that it can describe the coupling of *S* waves. The procedure eliminates problems with ray tracing in the vicinity of singularities; the common *S*-wave ray tracing is as stable as *P*-wave ray tracing. Due to the use of perturbation formulae, the ray tracing, dynamic ray tracing and coupling equations are much simpler and more transparent than in the exact case. As a byproduct of both coupling procedures, we get formulae for approximate evaluation of traveltimes of separate *S* waves. These formulae can find applications in migration and traveltime tomography.

The accuracy of the previous and modified coupling procedures is studied on several models of varying strength of anisotropy. First, we investigate the accuracy of perturbation formulae in homogeneous models, in which coupling does not exist. Then we study both coupling and accuracy of perturbation formulae in inhomogeneous models. We compare the results obtained by the coupling procedures with the results of the quasi-isotropic approach and standard ray theory.

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

## 1 INTRODUCTION

*S* waves propagating in inhomogeneous, weakly anisotropic media are usually coupled. In the far field, coupling can be described, for example, by the ‘connection’ formulae (Thomson *et al.* 1992) or by the coupling ray theory (Coates & Chapman 1990; Bulant & Klimeš 2002). The coupling ray theory represents an extension of the zero-order ray theory. There are various versions of the coupling

ray theory, which differ in the choice of the trajectory—common *S*-wave ray—along which coupling effects are evaluated, or by various approximations of the coupling equations. For an extensive review of such approximations, see Klimeš. & Bulant (2004).

One of the approximations of the coupling ray theory is the quasi-isotropic approach. It was proposed by Kravtsov (1968), see also Kravtsov & Orlov (1990). Later it was applied by Pšenčík (1998) to elastic media, see also Pšenčík & Dellinger (2001). In the

<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>
quasi-isotropic approach, the common S-wave ray is traced in a reference isotropic medium approximating the actual anisotropic medium. This reduces the accuracy of the corresponding computations. In order to increase their accuracy, Bakker (2002) proposed the use of the common S-wave ray, which is traced in the actual anisotropic medium, see also Klimeš (2006). As the polarization plane of S waves, Bakker (2002) used the plane perpendicular to the direction of the exact P-wave polarization vector. Farra & Pšenčík (2008), hereinafter referred to as Paper I, used Bakker's (2002) approach and combined it with the first-order ray tracing (FORT) and first-order dynamic ray tracing (FODRT) concept, which they used before for computing P waves in inhomogeneous weakly anisotropic media (Pšenčík & Farra 2005, 2007). In contrast to Bakker (2002), the S-wave polarization plane used in Paper I was only an approximation of the actual one. It was defined as a plane perpendicular to the first-order slowness vector.

As in the papers on FORT and FODRT for P waves, the computation of common S-wave rays and dynamic ray tracing along them in Paper I is also based on the perturbation theory, in which anisotropy is considered to be weak. The trajectory of the common S-wave ray corresponds to the Hamiltonian obtained from the average of the first-order eigenvalues of the Christoffel matrix, corresponding to the two S-waves propagating in anisotropic media. The common S-wave ray is computed in the actual medium and, therefore, it does not require the specification of a reference medium as in the quasi-isotropic approach. The reference medium in Paper I is only necessary if we wish to compute the second-order common S-wave traveltime corrections. The great advantage of such common S-wave rays is that their computation does not collapse in the vicinity of S-wave singularities as is typical for standard S-wave ray tracing in anisotropic media, see, for example, Vavryčuk (2003). The use of the common S-wave rays computed in the actual anisotropic medium is advantageous even in homogeneous media, where all common rays between the source and the receiver coincide with exact S-wave rays. The quantities calculated along the common ray (for example, the slowness vector) in the actual anisotropic medium approximate the actual quantities better.

The approximate coupling equations are derived in Paper I under the assumption that the deviations of anisotropy from isotropy are, simply speaking, of the order $O(\omega^{-1})$, where $\omega$ is the circular frequency. The coupling equations consist of two coupled, frequency-dependent, linear, ordinary differential equations for the S-wave amplitude coefficients and are solved along common S-wave rays. The coupling equations do not require definition of a reference medium.

In this paper, we propose a simple modification of the coupling equations proposed in Paper I, which considerably increases their accuracy. The modification consists in replacing two mutually perpendicular unit vectors, which define the zero-order polarization plane of the common S wave in Paper I, by vectors, which define the polarization plane more accurately. The expressions for the new vectors are very simple. Their use, however, leads not only to an improved specification of polarization plane, but, most importantly, to a considerable increase of the accuracy of traveltimes. Although the expressions for the new vectors require a reference medium, the reference medium enters the equations implicitly, and does not need to be specified.

There is a considerable difference between the formulation of the coupling ray theory studied by Klimeš. & Bulant (2004) or Bulant & Klimeš (2008) and the approach presented in this paper. In our approach, the common ray computation, dynamic ray tracing and coupling equations along it are approximate. The traveltime correc-

tions are incorporated in the coupling equations. There is no need for additional numerical quadratures along the common S-wave ray. Despite the differences in formulation, general conclusions are very similar. The proposed scheme can describe S-wave propagation in smooth inhomogeneous isotropic media (for this situation it reduces to the standard ray tracing for isotropic media), and it describes, with sufficient accuracy, the coupling of S waves in smooth inhomogeneous, weakly anisotropic media and, as the presented synthetic tests illustrate, can even properly describe well separated S waves.

Interesting and useful byproducts of the coupling equations are approximate formulae for computing traveltimes and polarizations of separate S waves propagating in inhomogeneous anisotropic media. The traveltime formulae are simple and transparent. We can distinguish in them the terms responsible for S1- and S2-wave traveltime separation and terms responsible for corrections of the traveltime along the common S-wave ray.

In Section 2, we review the main results of Paper I. In Section 3, we derive modified coupling equations, which represent the basic result of this paper. Section 4 is devoted to the illustration of the performance of both types of approximate coupling equations for several models of varying strength of anisotropy. For homogeneous models, in which coupling does not exist and for inhomogeneous models, in which it does exist, we compare the results obtained by the coupling equations of Paper I and by the modified coupling equations with the results of the standard ray theory or the quasi-isotropic approach. Advantages and limitations and future plans are briefly discussed in Section 5. Finally, in Appendix A, we study the behaviour of coupling equations, derived in Paper I and herein, in special cases. We describe how equations reduce in isotropic media or in media characterized by stronger anisotropy and/or weaker inhomogeneity and/or for waves with higher frequencies. Approximate formulae for computing separate S-wave traveltimes are also given there.

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 upper index [*M*] is used to denote quantities related to the S-wave common ray.

## 2 BASIC EQUATIONS

In the frequency domain, the zero-order ray approximation of the displacement vector **u** of an arbitrary wave propagating in an inhomogeneous anisotropic medium can be expressed as

$$ \mathbf{u}(x_m, \omega) = \mathbf{U}(x_m)\exp[i\omega\tau(x_m)]. \quad (1) $$

Here *i* is the imaginary unit, $\omega$ is the circular frequency, $\tau(x_m)$ is the eikonal, which also serves as the traveltime, and $\mathbf{U}(x_m)$ is the zero-order vectorial ray amplitude coefficient. Expression (1) is a good approximation of the displacement vector if the variations of the model and wave parameters within a wavelength are small, in other words if characteristic length *L* (the distance, on which the parameters change by an amount comparable with their size) is large with respect to the wavelength. This condition can be expressed as the requirement that parameter $\epsilon_1 \sim c/(\omega L)$ is small (Pšenčík 1998). This requirement is often simplified to $\epsilon_1 \sim \omega^{-1}$, see, for example, Paper I. Here we use the more general former definition of $\epsilon_1$, that is, $\epsilon_1 \sim c/(\omega L)$.

As in Paper I, we concentrate on the coupled S waves computed along a common S-wave ray and assume that in an inhomogeneous
weakly anisotropic medium the vectorial amplitude coefficient $U(x_m) = U[x_m(\tau)]$ has the form

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

In (2), $\mathcal{A}$ and $\mathcal{B}$ are S-wave amplitude coefficients and vectors $\mathbf{e}^{[1]}$ and $\mathbf{e}^{[2]}$ are mutually perpendicular unit vectors, to which the amplitude coefficients are related. Parameter $\tau$ is the traveltime. All the mentioned quantities are computed along the S-wave common ray, obtained by solving the first-order ray tracing (FORT) equations

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

Here $x_i$ and $p_i$ are the Cartesian coordinates of the S-wave common first-order ray and the components of the corresponding first-order slowness vectors, respectively. Parameter $\tau = \tau^{[\mathcal{M}]}(x_m)$ is the first-order traveltime. Symbol $G^{[\mathcal{M}]}$ denotes the first-order S-wave mean eigenvalue

$$\begin{aligned} G^{[\mathcal{M}]}(x_m, p_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]}). \end{aligned} \quad (4)$$

In (4), $G^{[1]}$ and $G^{[2]}$ are the first-order approximations of two smaller eigenvalues of the generalized Christoffel matrix $\Gamma$ (we call it generalized because it contains components of slowness vector $\mathbf{p}$ instead of the unit vector in the direction of $\mathbf{p}$ used in the standard Christoffel matrix):

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

Symbols $a_{ijkl}$ denote density-normalized elastic moduli,

$$a_{ijkl} = c_{ijkl}/\rho, \quad (6)$$

$c_{ijkl}$ being elements of the fourth-order tensor of elastic moduli and $\rho$ the density. The first-order S-wave eigenvalues $G^{[1]}$ and $G^{[2]}$ can also be expressed in terms of elements $B_{11}$, $B_{12}$ and $B_{22}$ 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]}. \quad (7)$$

Symbols $e_i^{[j]}$ in (4) and (7) denote the components of unit vectors $\mathbf{e}^{[j]}$. Vectors $\mathbf{e}^{[1]}$ and $\mathbf{e}^{[2]}$, see (2), are perpendicular to the third vector $\mathbf{e}^{[3]}$ chosen so that $\mathbf{e}^{[3]} = \mathbf{n}$. Here $\mathbf{n}$ is a unit vector specifying the direction of the first-order slowness vector $\mathbf{p}$. Vectors $\mathbf{e}^{[K]}$ can be chosen arbitrarily in the plane perpendicular to $\mathbf{n}$. Vector $\mathbf{e}^{[3]}$ can be determined from the second set of FORT eqs (3). Along the common S-wave ray, vectors $\mathbf{e}^{[K]}$ can be computed as the vectorial basis of the ‘wavefront orthonormal coordinate system’, see, for example, Červený (2001):

$$\frac{\mathrm{d}e_i^{[K]}}{\mathrm{d}\tau} = -(c^{[\mathcal{M}]})^2 \left( e_k^{[K]} \frac{\mathrm{d}p_k}{\mathrm{d}\tau} \right) p_i. \quad (8)$$

Here, $c^{[\mathcal{M}]} = c^{[\mathcal{M}]}(x_m, n_m)$ is the first-order common S-wave phase velocity corresponding to $G^{[\mathcal{M}]}$, $(c^{[\mathcal{M}]})^2 = G^{[\mathcal{M}]}(x_m, n_m)$.

For more details about FORT and FODRT along common S-wave rays see Paper I.

It was shown in Paper I that the accuracy of traveltime computations can be enhanced by calculating a correction along the first-order common S-wave ray. Although it is not strictly second-order traveltime correction, we refer to it as such, for the sake of simplicity. It reads

$$\Delta\tau^{[\mathcal{M}]} = \frac{1}{4} \int_{\tau_0}^{\tau} [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(x_m) - V_S^2(x_m)} d\tau. \quad (9)$$

Quantities $B_{13}$ and $B_{23}$ in eq. (9) are elements of the symmetric matrix $\mathbf{B}(x_m, p_m)$, see (7). Let us mention that the traveltime correction (9) does not depend on the choice of basis vectors $\mathbf{e}^{[K]}$, but it does depend on the choice of $V_P$ and $V_S$. Symbols $V_S$ and $V_P$ denote the S- and P-wave velocities of the reference isotropic medium, closely approximating the weakly anisotropic medium along the considered ray.

The zero-order vectorial ray amplitude coefficient $U(\tau)$ of the common S wave given in eq. (2), calculated along the S-wave common ray, reads

$$U(\tau) = \frac{A_0(\tau)\mathbf{e}^{[1]}(\tau) + B_0(\tau)\mathbf{e}^{[2]}(\tau)}{[\rho(\tau)c^{[\mathcal{M}]}(\tau)]^{1/2} \mathcal{L}^{[\mathcal{M}]}(\tau)}. \quad (10)$$

Parameter $\tau = \tau^{[\mathcal{M}]}(x_m)$ is again the first-order traveltime along the S-wave common ray, obtained by solving FORT eqs (3). The symbols $\rho$, $c^{[\mathcal{M}]}$ and $\mathcal{L}^{[\mathcal{M}]}$ denote the density and first-order common S-wave phase velocity and geometrical spreading. The spreading $\mathcal{L}^{[\mathcal{M}]}$ represents the first-order approximation of the exact common S-wave spreading computed along an exact common S-wave ray. Vectors $\mathbf{e}^{[1]}$ and $\mathbf{e}^{[2]}$ are the basis vectors of the wavefront orthonormal coordinate system (8). Vectors $\mathbf{e}^{[K]}$ define the zero-order polarization plane of the common S wave. Thus the vectorial amplitude coefficient $U(\tau)$ is situated in the zero-order polarization plane. Frequency-dependent amplitude coefficients $A_0$ and $B_0$ can then be obtained by solving the system of coupled differential equations

$$\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} A_0 \\ B_0 \end{pmatrix}. \quad (11)$$

In (11), $B_{jL}$ are again the elements of the symmetric matrix $\mathbf{B}(x_m, p_m)$, see (7). Matrix $\mathbf{B}(x_m, p_m)$ can be obtained by simple rotation from matrix $\tilde{\mathbf{B}}(x_m, p_m)$, specified explicitly in Paper I. Because of the relation of matrix $\mathbf{B}$ to the first-order eigenvalues of the Christoffel matrix (Farra & Pšenčík 2003), we call eqs (11) ‘first-order coupling equations’.

If the wavefield is generated by point force $\mathbf{F}$ acting at time $\tau = \tau_0$, the initial conditions for coupling equations (11) read

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

Here $c_0^{[\mathcal{M}]}$ 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.

In Paper I, the system of coupling equations (11) was derived under the assumption that another small parameter $\epsilon_2 \sim ||\Delta a_{ijkl}||/||a_{ijkl}|| \sim \Delta c/c (||.||$ denotes the tensor norm), characterizing weak anisotropy, is of the same order as the earlier introduced small parameter $\epsilon_1$. Here, $\Delta a_{ijkl}$ are perturbations of the density-normalized elastic parameters from the reference isotropic medium, $\Delta c$ is the maximum difference of the phase velocities of the two S waves propagating in the weakly anisotropic medium.

## 3 SECOND-ORDER COUPLING EQUATIONS

For greater accuracy of the results, it was proposed in Paper I to replace vectors $\mathbf{e}^{[K]}$, which define the zero-order polarization plane of the common S wave, by their first-order counterparts $\mathbf{f}^{[K]}$

$$\mathbf{f}^{[K]} = \mathbf{e}^{[K]} - (c^{[\mathcal{M}]})^2 \frac{B_{K3}(x_m, p_m)}{V_P^2 - V_S^2} \mathbf{e}^{[3]}. \quad (13)$$
Vectors $\mathbf{f}^{[K]}$ are situated in the plane perpendicular to the first-order eigenvector $\mathbf{f}^{[3]}$ corresponding to the largest eigenvalue of the generalized Christoffel matrix (5)

$$ \mathbf{f}^{[3]} = (c^{[\mathcal{M}]})^2 \frac{B_{13}(x_m, p_m)}{V_P^2 - V_S^2} \mathbf{e}^{[1]} + (c^{[\mathcal{M}]})^2 \frac{B_{23}(x_m, p_m)}{V_P^2 - V_S^2} \mathbf{e}^{[2]} + \mathbf{e}^{[3]}. \quad (14) $$

Quantities $B_{13}$, $B_{23}$ are again elements of matrix $\mathbf{B}(x_m, p_m)$, see (7). Symbols $V_P$ and $V_S$ denote the P- and S-wave velocities of the reference isotropic medium, also used in eq. (9).

In this paper, increase of accuracy of the specification of the S-wave polarization plane is not the only reason for replacing vectors $\mathbf{e}^{[K]}$ by $\mathbf{f}^{[K]}$. More important is the use of the vectors $\mathbf{f}^{[K]}$ in the rederivation of coupling equations, which leads to their significantly increased accuracy. After replacement of vectors $\mathbf{e}^{[K]}$ by $\mathbf{f}^{[K]}$ in (10), the zero-order vectorial ray amplitude coefficient $\mathbf{U}(\tau)$ reads

$$ \mathbf{U}(\tau) = \frac{A_1(\tau)\mathbf{f}^{[1]}(\tau) + B_1(\tau)\mathbf{f}^{[2]}(\tau)}{[\rho(\tau)c^{[\mathcal{M}]}(\tau)]^{1/2} \mathcal{L}^{[\mathcal{M}]}(\tau)}. \quad (15) $$

In (15), $A_1$ and $B_1$ are the amplitude coefficients, which we are seeking instead of coefficients $A_0$ and $B_0$ used in (10). We can now proceed as in Paper I: insert eq. (15) into eq. (1), and the result into the elastodynamic equation, and neglect the terms of order $O(1)$, $O(\omega^{-1})$ and less. Within this approximation, the replacement of vectors $\mathbf{e}^{[K]}$ by $\mathbf{f}^{[K]}$ only affects the coefficient of $(i\omega)^2$ in eq. (31) of Paper I. The effects on the coefficient of $i\omega$ in eq. (31) can be neglected. Thus, if we take into account the definition of matrix $\mathbf{B}(x_m, p_m)$ in (7), eq. (32) of Paper I now reads

$$ (\Gamma_{ik} - \delta_{ik}) C f_i^{[K]} f_k^{[L]} = C [B_{KL} - \delta_{KL} - 2(c^{[\mathcal{M}]})^2 \frac{B_{K3} B_{L3}}{V_P^2 - V_S^2} + (c^{[\mathcal{M}]})^4 (B_{33} - 1) \frac{B_{K3} B_{L3}}{(V_P^2 - V_S^2)^2}]. \quad (16) $$

For a specific choice of velocities $V_S$ and $V_P$

$$ V_S^2(x_m) = (c^{[\mathcal{M}]})^2, \quad V_P^2(x_m) = (c^{[\mathcal{M}]})^2 B_{33}(x_m, p_m), \quad (17) $$

eq. (16) can be rewritten to read

$$ (\Gamma_{ik} - \delta_{ik}) C f_i^{[K]} f_k^{[L]} = C(M_{KL} - \delta_{KL}). \quad (18) $$

The elements of the $2 \times 2$ matrix $\mathbf{M}$ in (18) read

$$ M_{KL}(x_m, p_m) = B_{KL}(x_m, p_m) - (c^{[\mathcal{M}]})^2 \frac{B_{K3}(x_m, p_m) B_{L3}(x_m, p_m)}{V_P^2 - V_S^2}. \quad (19) $$

Matrix $\mathbf{M}$ was introduced in Farra (2001) and Farra & Pšenčík (2003) for calculating the second-order S-wave eigenvalues of the Christoffel matrix and the corresponding first-order S-wave eigenvectors, which represent first-order polarization vectors. Eq. (17) yields $B_{33} - 1 = (c^{[\mathcal{M}]})^{-2}(V_P^2 - V_S^2)$. Eq. (19) can thus be rewritten to the alternative form

$$ M_{KL}(x_m, p_m) = B_{KL}(x_m, p_m) - \frac{B_{K3}(x_m, p_m) B_{L3}(x_m, p_m)}{B_{33}(x_m, p_m) - 1}. \quad (20) $$

Note that the element $B_{33}$ of the matrix $\mathbf{B}$ represents the first-order approximation of the largest eigenvalue of the Christoffel matrix (5). The coupled system of differential equations (11) then transforms into the system

$$ \begin{pmatrix} \mathrm{d}A_1/\mathrm{d}\tau \\ \mathrm{d}B_1/\mathrm{d}\tau \end{pmatrix} = -\frac{i\omega}{2} \begin{pmatrix} M_{11} - 1 & M_{12} \\ M_{12} & M_{22} - 1 \end{pmatrix} \begin{pmatrix} A_1 \\ B_1 \end{pmatrix}. \quad (21) $$

We call eqs (21) ‘second-order coupling equations’. They are obtained from the first-order coupling equations (11) just by replacing the elements $B_{11}$, $B_{12}$ and $B_{22}$ of matrix $\mathbf{B}$ by the elements $M_{11}$, $M_{12}$ and $M_{22}$ of matrix $\mathbf{M}$. In Appendix A, we show that the replacement of eqs (11) by (21) leads to the increase of accuracy of the computed traveltimes along the S-wave common ray. From the analysis made in Appendix A, we can see that the traveltimes computed implicitly by eqs (21) contain the second-order traveltime correction (9). We illustrate this fact also on numerical examples.

The initial conditions have a form similar to those in (12). Only vectors $\mathbf{e}^{[K]}$ are replaced by vectors $\mathbf{f}^{[K]}$

$$ A_1(\tau_0) = \frac{f_k^{[1]}(\tau_0) F_k}{4\pi (\rho_0 c_0^{[\mathcal{M}]})^{1/2}}, \quad B_1(\tau_0) = \frac{f_k^{[2]}(\tau_0) F_k}{4\pi (\rho_0 c_0^{[\mathcal{M}]})^{1/2}}. \quad (22) $$

## 4 NUMERICAL EXAMPLES

In order to illustrate the accuracy of the first-order and second-order coupling equations (11) and (21), we consider the VSP configuration, which we used in our previous studies of FORT and FODRT for P waves, see Pšenčík & Farra (2005, 2007). The source and the borehole are situated in a vertical plane $(x, z)$. The borehole is parallel to the $z$ axis, the vertical single-force source is located on the surface at $z = 0$ km, at a distance of 1 km from the borehole. The source-time function is a windowed symmetric Gabor wavelet, $\exp[-(2\pi f/\gamma)^2 t^2] \cos(2\pi f t)$, with the dominant frequency $f = 50$ Hz and $\gamma = 4$. There are 29 three-component receivers in the borehole, distributed with a uniform step of 0.02 km, with receiver depths ranging from 0.01 to 0.57 km. The receivers record the vertical (positive downwards), transverse and radial (along the line connecting the source and the top of the borehole; positive away from the source) components of the wavefield. The recording system is right-handed. All calculated seismograms are shown with no differential scaling between components and traces, so that true relative amplitudes can be seen.

We consider three models, QI, QI2 and QI4, used by Klimeš & Bulant (2004) and Bulant & Klimeš (2008). Model QI coincides with the WA model of Pšenčík & Dellinger (2001). The models are vertically inhomogeneous HTI media with constant vertical gradients of the elastic moduli. The axis of symmetry is rotated everywhere in the horizontal plane by $45^\circ$ from the vertical plane $(x, z)$. The S-wave anisotropy in the $(x, z)$ plane defined as $(c_{S1} - c_{S2})/c_{\text{average}} \times 100$ per cent ranges, from the horizontal to vertical direction, from 1 to 4 per cent, from 4 to 7 per cent and from 11 to 13 per cent for the QI, QI2 and QI4 models, respectively. The matrices of the density-normalized elastic moduli can be found in the above references. The variations of the S-wave phase velocities in the $(x, z)$ plane for all three models are shown in Fig. 1. The left-hand plots correspond to $z = 0$ km, the right-hand plots to $z = 1$ km. Model QI is shown in the top, QI2 in the middle and QI4 in the bottom plot. Velocities are shown as functions of the angle of propagation. They vary from 0° (horizontal propagation) to 90° (vertical propagation). Although the coupling method based on the FORT can accommodate arbitrary lateral variations of the elastic moduli, the models used exhibit only vertical variations.

In addition to inhomogeneous QI, QI2 and QI4 models, we also consider their homogeneous counterparts specified by the elastic parameters of the above models at $z = 0$ km. The homogeneous models are thus represented by the plots in the left-hand column of Fig. 1. We call these models QI HOM, QI2 HOM and QI4 HOM. Since there is no coupling in homogeneous media (the two S waves
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The S-wave phase-velocity sections in the (x, z) plane for models QI (top panel), QI2 (middle panel) and QI4 (bottom panel). The velocities vary from the horizontal (0°) to vertical (90°) direction of the wave normal. Left-hand plots correspond to *z* = 0 km, right-hand plots to *z* = 1 km.</figcaption></figure>

are decoupled there), homogeneous models allow us to investigate
the accuracy of the traveltime and spreading approximations.

Fig. 2 shows the relative traveltime and geometrical-spreading
differences for the homogeneous QI2 HOM model. The differences
are calculated from quantities computed along common S-wave
rays and from exact quantities computed along rays of S1 and S2
waves generated by the ray tracer for anisotropic media—a modified
program package ANRAY (Gajewski & Pšenčík 1990), which is
based on the standard ray theory. The upper plot shows relative
traveltime differences of the first-order traveltime τ[<sup>M</sup>] (blue – 1)
obtained by solving FORT eqs (3) along the common S-wave ray and
second-order traveltime τ[<sup>M</sup>] + Δτ[<sup>M</sup>] (red – 2) from the average
of exact S1- and S2-wave traveltimes. For Δτ[<sup>M</sup>] see eq. (9). The
V<sub>S</sub> and V<sub>P</sub> velocities in eq. (9) are determined from eq. (17). In the
bottom plot of Fig. 2, we show relative differences of the first-order
common S-wave geometrical spreading L[<sup>M</sup>] from the spreading
of exact S1 (blue), S2 (red) wave and average of S1 and S2 waves
(black).

We can see that the relative difference of the first-order travel-
time τ[<sup>M</sup>] from average of exact traveltimes of S1 and S2 waves is
less than 1 per cent. This difference is substantially reduced if the
second-order traveltime correction (9) is used. Both relative differ-
ences remain nearly constant for the studied range of depths. The
relative difference of the first-order geometrical spreading L[<sup>M</sup>]

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Model QI2 HOM. Top panel: the relative traveltime differences of the first-order (blue – 1) and the second-order (red – 2) common S-wave traveltimes from average of exact S1- and S2-wave traveltimes. Bottom panel: the relative differences of the first-order common S-wave geometrical spreading from the exact S1-wave (blue), S2-wave (red) and average spreading (black).</figcaption></figure>

from the average of exact spreadings of S1 and S2 waves is approx-
imately –3 per cent. The first-order geometrical spreading L[<sup>M</sup>],
computed along the common S-wave ray, is smaller than the exact
spreading of each of the two S waves (the relative differences, red
and blue, are negative).

For generation of approximate synthetic seismograms in the re-
maining figures, we used solely first- and second-order coupling
equations (11) and (21). We did use neither eq. (9) nor equations
from Appendix A.

In Fig. 3, we compare the synthetic seismograms, computed from
the first-order coupling equations (11) (red) and the standard ray
theory seismograms (black) in the QI2 HOM model. The red seis-
mograms are plotted over the black ones in this and the following
comparisons. We can see that the faster S1 wave (observable in the
vertical component and, for deeper receivers, also on the transverse
component) is relatively well approximated by eqs (11). The approx-
imation of the S2 wave is worse. We can see a very poor fit of the
approximate seismograms in the radial and transverse components,
especially for deeper receivers. This is the consequence of the trav-
eltime approximation incorporated in eqs (11), see eqs (A5)–(A6)
in Appendix A. The discrepancies between the approximate and
standard ray theory seismograms (black) are removed, if the
<figure><img src="image_4.png" /></figure>

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

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

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

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

<figure><img src="image_9.png" /><figcaption>**Figure 3.** Comparison of the seismograms computed with the first-order coupling equations (11) (red) and the standard ray theory seismograms (black) for the vertical single-force source in the QI2 HOM model.</figcaption></figure>

<caption>**Figure 4.** Comparison of the seismograms computed with the second-order coupling equations (21) (red) and the standard ray theory seismograms (black) for the vertical single-force source in the QI2 HOM model.</caption>

approximate seismograms (red) are computed from the second-order coupling equations (21), see Fig. 4. This is the consequence of the use of a better traveltime approximation incorporated in eqs (21), see eqs (A15)–(A16) in Appendix A. We can see that the second-order coupling equations (21) (red) yield a good fit with the standard ray theory seismograms (black) even for model QI4 HOM with the S-wave anisotropy between 11 and 13 per cent, shown in Fig. 5. The S1 wave computed by eqs (21) is slightly faster and has a slightly overestimated amplitude with respect to the S1 wave computed by the standard ray method. The relative difference of the first-order geometrical spreading $\mathcal{L}^{[M]}$ from the S1-wave spreading

is about –11 per cent in this case. Let us emphasize that the two well separated S waves shown in red in Fig. 5 are computed along a ‘single’ common S-wave ray.

Let us now consider QI models with vertical inhomogeneity, in which we can observe coupling effects. We do not show the plots of the relative differences of the first- and second-order traveltimes and of the geometrical spreading. For model QI2, they differ only a little from the plots shown in Fig. 2.

In Fig. 6, we compare the seismograms computed from the first-order coupling equations (11) (black) and from the second-order coupling equations (21) (red) for the QI model. For the faster S1
<figure><img src="image_10.png" /></figure>

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

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

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

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

<figure><img src="image_15.png" /><figcaption>**Figure 5.** Comparison of the seismograms computed with the second-order coupling equations (21) (red) and the standard ray theory seismograms (black) for the vertical single-force source in the QI4 HOM model.</figcaption></figure>

<caption>**Figure 6.** Comparison of the seismograms computed with the first-order coupling equations (11) (black) and the second-order coupling equations (21) (red) for the vertical single-force source in the QI model.</caption>

wave (mostly observable in the vertical and radial components), we have a nearly perfect fit of both approximations. For the S2 wave (mostly observable in the transverse component), we can observe certain differences. The S2 wave computed by the first-order coupling equations is slightly faster than the S2 wave computed by the second-order coupling equations.

As mentioned above, model QI is identical with model WA studied by Pšenčík & Dellinger (2001). They used the quasi-isotropic approximation with a common S-wave ray traced in a reference isotropic medium and compared the results with the results computed by the reflectivity method. In Fig. 7, we compare the seismo-

grams computed by the second-order coupling equations (21) (red) with the seismograms computed by the quasi-isotropic approach (black) for the QI model. No amplitude normalization (which had to be used by Pšenčík & Dellinger 2001) is used. Except for the slightly higher amplitudes of both S waves and slightly faster S2 wave computed by the quasi-isotropic approach, Fig. 7 has a similar character as Fig. 6. It means that the first-order coupling equations and quasi-isotropic approach generate results of comparable accuracy.

Since the anisotropy of the QI model is relatively weak, the results of the quasi-isotropic approach and the first-order and second-order coupling equations presented in this paper are comparable. Let us
<figure><img src="image_16.png" /></figure>

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

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

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

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

<figure><img src="image_21.png" /><figcaption>**Figure 7.** Comparison of the seismograms computed with the second-order coupling equations (21) (red) and the quasi-isotropic approach (black) for the vertical single-force source in the QI model.</figcaption></figure>

<caption>**Figure 8.** Comparison of the seismograms computed with the second-order coupling equations (21) (red) and the standard ray theory seismograms (black) for the vertical single-force source in the QI model.</caption>

now compare the seismograms computed by the second-order coupling equations (21) (red) with the standard ray theory seismograms computed by the package ANRAY for the QI model. Fig. 8 shows this comparison. We can see the nearly perfect fit of the S1 wave, observable again mostly in the vertical and radial components, except for the slightly higher amplitudes computed by ANRAY at deeper receivers, but a strong misfit of the S2 wave, mostly observable in the transverse component. The most pronounced differences, mainly due to the phase shift, can be seen at the shallow receivers. This misfit indicates the failure of the standard ray theory (used

in the ANRAY package) to describe properly the phenomenon of coupling.

Let us now consider the stronger anisotropy, specifically model QI2. Similarly as in Fig. 6, we compare the results of the computations based on the first-order coupling equations (11) (black) with those based on the second-order coupling equations (21) (red) in Fig. 9. While in Fig. 6 both equations yielded comparable results, in the medium with stronger anisotropy their results differ. The most dramatic difference (both in phase shift and amplitude) can be observed for the S2 wave. The differences are so pronounced that they can be observed in all three components, mostly in the transverse.
<figure><img src="image_22.png" /><figcaption>**Figure 9.** Comparison of the seismograms computed with the first-order coupling equations (11) (black) and the second-order coupling equations (21) (red) for the vertical single-force source in the QI2 model.</figcaption></figure>

As superior, we consider, of course, the seismograms computed using eqs (21). In contrast to the first-order coupling equations, eqs (21) describe, specifically for the deeper receivers in the transverse component, the separation of the S1 and S2 waves.

As mentioned above, coupling equations (11) and (21) yield significantly different results for the slower S2 wave while results for the faster S1 wave are nearly identical. Explanation is related to the accuracy of the approximation of polarization vectors. In the studied models, faster S1 wave is polarized perpendicularly to the plane defined by its slowness vector and the axis of symmetry. From the analysis of polarization and traveltime formulae given in

Appendix A, we can find that vectors **e**[1] and **f**[1] are identical and approximate rather well the exact polarization vector. Analysis of eqs (A6) and (A16) then shows that the terms $\Delta\tau_{S1}$ in these equations represent identical traveltime corrections (therefore, the S1 wave is described equally well by both systems of coupling equations). The terms $\Delta\tau_{S2}$ obtained from (A6) and (A16), however, differ. The difference equals twice the second-order traveltime correction (9). This indicates that the second-order traveltime difference may play equally important role in coupling equations (21) as the second-order traveltime correction, which alone is not sufficient for good traveltime approximation.

If the wave with the polarization perpendicular to the plane defined by the slowness vector and the axis of symmetry were the slower one, coupling equations (11) and (21) would yield nearly identical results for the S2 wave. For more general anisotropy, coupling equations (11) and (21) would yield different results for each of the S waves.

The separation of the S1 and S2 waves is visible more clearly in Fig. 10, which shows the same comparison as Fig. 8. Specifically, Fig. 10 compares the seismograms, computed by the second-order coupling equations (21) (red), with the standard ray theory seismograms computed by the ANRAY package for the QI2 model. We can see that the fit is now much better than in Fig. 8. This is because the anisotropy is now sufficiently strong so that the coupling effects are not as pronounced (the coupling still affects wave field at shallow receivers). In the transverse component, we can observe very clear separation of the S1 and S2 waves.

Fig. 11 shows the same as Fig. 10, but for the QI4 model. We can now observe clear separation (of about 0.06 s; approximately three wave periods) of the two S waves in all components. The fit is not as good as in Fig. 10. The S1 wave computed from eqs (21) is slightly faster and has a slightly overestimated amplitude with respect to the S1 wave computed by the standard ray method, see a similar observation in Fig. 5. Let us emphasize again that, while the two well-separated S waves shown in black in Fig. 11 are calculated each along a different S-wave ray, the S waves shown in red in Fig. 11 are computed along a 'single' common S-wave ray.

## 5 DISCUSSION AND CONCLUSIONS

We have modified approximate coupling equations for computing the coupled S-waves propagating in smooth inhomogeneous, weakly anisotropic media, proposed in Paper I. The modification consists in replacing two mutually perpendicular unit vectors, which define the zero-order polarization plane of the common S wave in Paper I by vectors, which define the polarization plane more accurately. The modification leads to a substantial increase in the accuracy of the computed S-wavefield as illustrated by numerical examples.

The proposed procedure is applicable to the S-waves propagating in inhomogeneous, isotropic or weakly anisotropic media of arbitrary symmetry. In isotropic media, it reduces to exact ray tracing and dynamic ray tracing. Eqs (8) reduce to the equations for the vector basis of the ray-centred coordinate system, in which the vectors $e^{[K]}$ specify the polarization vectors of the S wave propagating in an inhomogeneous isotropic medium. In anisotropic media, the proposed procedure can deal with coupled as well as decoupled S waves. The common S-wave ray tracing is regular everywhere including singular regions.

The second-order coupling equations have a simple and transparent form. They differ from the first-order coupling equations only by the substitution of elements $B_{11}$, $B_{12}$ and $B_{22}$ of matrix **B** by elements $M_{11}$, $M_{12}$ and $M_{22}$ of matrix **M**. Computation of elements of
<figure><img src="image_23.png" /></figure>

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

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

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

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

<figure><img src="image_28.png" /><figcaption>**Figure 10.** Comparison of the seismograms computed with the second-order coupling equations (21) (red) and the standard ray theory seismograms (black) for the vertical single-force source in the QI2 model.</figcaption></figure>

<caption>**Figure 11.** Comparison of the seismograms computed with the second-order coupling equations (21) (red) and the standard ray theory seismograms (black) for the vertical single-force source in the QI4 model.</caption>

matrix **M** is simple and does not require much extra computational effort. The second-order coupling equations contain the corrections of the common *S*-wave traveltime implicitly. Simple, approximate formulae for traveltimes of separate *S* waves can be obtained as a byproduct of the coupling procedure. These formulae may find applications in migration and traveltime tomography based on *S* waves. The coupling procedure described in this paper is designed for smooth media without interfaces, but it can be simply generalized for laterally varying, layered, weakly anisotropic structures.

The computations of common *S*-wave rays and of dynamic ray tracing along them are based on FORT and FODRT (Pšenčík & Farra

2005, 2007). The coupling equations are also approximate, based on perturbation formulae. Despite this, the proposed scheme generates quite accurate results; compare the results presented in this paper with the results of Bulant & Klimeš (2008). The scheme proposed in this paper also provides satisfactory results for decoupled *S* waves as shown in the comparisons with the results of the standard ray theory for anisotropic media.

As the next step, we plan to study the behaviour of the proposed procedure in more complicated models. We plan to concentrate specifically on *S*-wave computations close to singularities. We also plan to test the traveltime formulae given in Appendix A.
ACKNOWLEDGMENTS

A substantial part of this work was done during IP's stay at the IPG Paris at the invitation of the IPGP. We appreciate the useful discussions with Einar Iversen and his comments. Comments of Peter Bakker helped to improve the text significantly. We also thank Václav Vavryčuk for his comments. We are grateful to the Consortium Project 'Seismic waves in complex 3-D structures' (SW3D) and Research Project 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.

Bulant, P. & Klimeš, L., 2008. Numerical comparison of the isotropic-common-ray and anisotropic-common-ray approximations of the coupling ray theory, *Geophys. J. Int.*, **175**, 357–374.

Červený, V., 2001. *Seismic Ray Theory*, Cambridge Univ. 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., 2001. High order expressions of the phase velocity and polarization of qP and qS waves in anisotropic media, *Geophys. J. Int.*, **147**, 93–105.

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.

Farra, V. & Pšenčík, I., 2008. First-order ray computations of coupled S waves in inhomogeneous weakly anisotropic media, *Geophys. J. Int.*, **173**, 979–989.

Gajewski, D. & Pšenčík, I., 1990. Vertical seismic profile synthetics by dynamic ray tracing in laterally varying layered anisotropic structures, *J. geophys. Res.*, **95**, 11301–11315.

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

Klimeš, L. & Bulant, P., 2004. Errors due to the common ray approximations of the coupling ray theory, *Stud. Geophys. Geod.*, **48**, 117–142.

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

Kravtsov, Yu.A. & Orlov, Yu.I., 1990. *Geometrical Optics of Inhomogeneous Media*, Springer–Verlag, Heidelberg.

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

Pšenčík, I. & Dellinger, J., 2001. Quasi-shear waves in inhomogeneous weakly anisotropic media by the quasi-isotropic approach: a model study, *Geophysics*, **66**, 308–319.

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.

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.

Vavryčuk, V., 2003. Behavior of rays near singularities in anisotropic media, *Phys. Rev.*, **B67**, 054 105-1–054 105-8.

APPENDIX A: SPECIAL SITUATIONS

Modified coupling equations were derived for inhomogeneous weakly anisotropic media under the assumption that the two small parameters $\epsilon_1 \sim c/(\omega L)$ and $\epsilon_2 \sim |\Delta c|/c$ are of the same order, $\epsilon_1 \sim \epsilon_2$. Here we analyse how these equations reduce in various

special situations like, for example, for S waves propagating in inhomogeneous isotropic media or for weakly coupled waves in anisotropic media including homogeneous anisotropic media. Study of weakly coupled S waves offers us possibility to analyze structure of traveltime corrections, polarization, etc. See also Coates & Chapman (1990), Pšenčík (1998).

In isotropic media, that is, for $\epsilon_2 = 0$, eqs (11) yield $A_0 = \text{const}$, $B_0 = \text{const}$. If we insert this into eq. (10), and the result into eq. (1), we obtain the well-known zero-order ray expression for the displacement vector $\mathbf{u}$ of an S wave propagating in an inhomogeneous isotropic medium. The vectors $e^{[K]}$ computed from eq.(8) represent the S-wave polarization vectors.

For stronger anisotropy ($\Delta c$ large), higher frequencies ($\omega$ large) or weaker inhomogeneity, including homogeneity (characteristic length $L$ large or infinite), that is, for the case $\epsilon_1 \ll \epsilon_2$, the coupling is expected to be weaker. In the following, we first seek solutions of coupling equations (11) and then (21) under the assumption of weak coupling.

We can seek solutions of eqs (11) in the following form (Pšenčík 1998)

$$A_0(\tau) = A'_0 \exp(i\omega\Delta\tau), \quad B_0(\tau) = B'_0 \exp(i\omega\Delta\tau). \quad (A1)$$

Here $A'_0$ and $B'_0$ are the amplitude factors, $\Delta\tau$ are deviations of the traveltimes from the first-order traveltime $\tau^{[M]}$ calculated along the common S-wave ray. Inserting (A1) into eqs (11) and considering $\epsilon_1 \ll \epsilon_2$, we arrive at $dA'_0/d\tau \sim 0$, $dB'_0/d\tau \sim 0$ and

$$2 \frac{d\Delta\tau}{d\tau} \begin{pmatrix} A'_{0} \\ B'_{0} \end{pmatrix} + \begin{pmatrix} B_{11} - 1 & B_{12} \\ B_{12} & B_{22} - 1 \end{pmatrix} \begin{pmatrix} A'_{0} \\ B'_{0} \end{pmatrix} = 0. \quad (A2)$$

Eq. (A2) can be expressed in the form of an eigenvalue problem

$$\begin{bmatrix} B_{11} - (1 - 2d\Delta\tau/d\tau) & B_{12} \\ B_{12} & B_{22} - (1 - 2d\Delta\tau/d\tau) \end{bmatrix} \begin{pmatrix} A'_{0} \\ B'_{0} \end{pmatrix} = 0. \quad (A3)$$

Eq. (A3) represents a system of two equations for two eigenvalues $1 - 2d\Delta\tau/d\tau$ and two corresponding eigenvectors situated in the plane specified by vectors $e^{[K]}$. The two eigenvalues and corresponding eigenvectors correspond to two independent S waves that we call S1 and S2 in the following. From the eigenvalues of eq. (A3), we get

$$\frac{d\Delta\tau}{d\tau} = \frac{1}{2} - \frac{1}{4} \left[ (B_{11} + B_{22}) \pm \sqrt{(B_{11} - B_{22})^2 + 4B_{12}^2} \right]. \quad (A4)$$

Taking into account that $B_{11}(x_m, p_m) + B_{22}(x_m, p_m) = 2$, see eqs (4) and the first-order eikonal equation along the common S-wave ray, $G^{[M]}(x_m, p_m) = 1$, we can use (A4) in the expression for the approximate traveltimes $\tau_{S1}$ and $\tau_{S2}$ of the S waves, between points corresponding to $\tau_0$ and $\tau$ on the common S-wave ray

$$\tau_{S1,S2}(\tau, \tau_0) = \tau^{[M]}(\tau, \tau_0) + \Delta\tau_{S1,S2}(\tau, \tau_0). \quad (A5)$$

The term $\tau^{[M]}(\tau, \tau_0)$ is the first-order traveltime calculated between points corresponding to $\tau_0$ and $\tau$ on the common S-wave ray. If the common S-wave ray does not pass through a singularity, the term $\Delta\tau_{S1,S2}$, which represents deviations of the traveltimes of S1 and S2 waves from $\tau^{[M]}(\tau, \tau_0)$, has the following meaning

$$\begin{align} \Delta\tau_{S1} &= -\frac{1}{4} \int_{\tau_0}^{\tau} \sqrt{(B_{11} - B_{22})^2 + 4B_{12}^2} d\tau, \nonumber \\ \Delta\tau_{S2} &= \frac{1}{4} \int_{\tau_0}^{\tau} \sqrt{(B_{11} - B_{22})^2 + 4B_{12}^2} d\tau. \tag{A6} \end{align}$$
In this way, the S1 wave is faster than the S2 wave. The zero value of the argument of the integrals in eq. (A6) indicates either a singularity or S-wave propagation in an isotropic medium. The approximate traveltime difference of the S1 and S2 waves is thus $\frac{1}{2} \int_{\tau_0}^{\tau} \sqrt{(B_{11}-B_{22})^2+4B_{12}^2} d\tau$. Since the difference is related to the first-order traveltime calculated along the common ray, we call it the ‘first-order traveltime difference’. Formulae (A6) are independent of the choice of $\mathbf{e}^{[K]}$.

The eigenvectors of eq. (A3) specify the directions of the zero-order polarization vectors of the two decoupled waves. We can introduce the polarization vectors $\mathbf{e}^{[K]}$ as mutually orthogonal unit vectors situated in the plane specified by vectors $\mathbf{e}^{[K]}$, see eq. (8). At any point of the common S-wave ray, the vectors $\mathbf{e}^{[K]}$ are given by the expressions

$$ 
\begin{aligned}
\mathbf{e}^{[1]} &= \mathbf{e}^{[1]} \cos \Phi_0 + \mathbf{e}^{[2]} \sin \Phi_0, \\
\mathbf{e}^{[2]} &= -\mathbf{e}^{[1]} \sin \Phi_0 + \mathbf{e}^{[2]} \cos \Phi_0.
\end{aligned}
\quad (A7) $$

Angle $\Phi_0$ is the angle, which leads to the diagonalization of the matrix on the left-hand side of eq. (A3). In other words, rotation through angle $\Phi_0$ changes the matrix $\mathbf{B}$ into $\mathbf{B}'$, in which $B_{12}' = 0$. From this condition we get (Pšenčík 1998; Farra 2001; Farra & Pšenčík 2003)

$$ \tan 2\Phi_0 = \frac{2B_{12}}{B_{11} - B_{22}}. \quad (A8) $$

Angle $\Phi_0$ is chosen so that vector $\mathbf{e}^{[1]}$ corresponds to the S1 wave. The polarization vector $\mathbf{e}^{[2]}$ corresponds to the S2 wave.

We can now rewrite eq. (10) in the following ‘uncoupled’ form

$$ \mathbf{U}(\tau) = \frac{\mathcal{D}'_0 \mathbf{e}^{[1]}(\tau) \exp(i\omega\Delta\tau_{S1}) + \mathcal{F}'_0 \mathbf{e}^{[2]}(\tau) \exp(i\omega\Delta\tau_{S2})}{[\rho(\tau)c^{[M]}(\tau)]^{1/2} \mathcal{L}^{[M]}(\tau)}. \quad (A9) $$

In (A9), $\mathcal{D}'_0$ and $\mathcal{F}'_0$ are factors constant along the ray, which can be expressed in terms of the factors $\mathcal{A}'_0$ and $\mathcal{B}'_0$, given in (A1), in the following way

$$ 
\begin{aligned}
\mathcal{D}'_0 &= \mathcal{A}'_0 \cos \Phi_0(\tau_0) + \mathcal{B}'_0 \sin \Phi_0(\tau_0), \\
\mathcal{F}'_0 &= -\mathcal{A}'_0 \sin \Phi_0(\tau_0) + \mathcal{B}'_0 \cos \Phi_0(\tau_0).
\end{aligned}
\quad (A10) $$

The factors $\mathcal{A}'_0$, $\mathcal{B}'_0$ can be determined from the initial conditions, the angle $\Phi_0(\tau_0)$ from eq. (A8) applied at $\tau = \tau_0$. The traveltime deviations $\Delta\tau_{S1}$, $\Delta\tau_{S2}$ are given in (A6), the polarization vectors $\mathbf{e}^{[K]}$ in (A7). After inserting eq. (A9) into eq. (1), eq. (1) describes two decoupled S waves. Each of them is given in the form of the zero-order ray expression for the displacement vector $\mathbf{u}$ of an S wave propagating in an inhomogeneous anisotropic medium. Both waves share the same spreading factor, but they differ by their polarizations and traveltimes, see eq. (A9). The wave specified by the factor $\mathcal{D}'_0$ and the polarization vector $\mathbf{e}^{[1]}$ is the S1 wave, the other is the S2 wave. Note that for the evaluation of (A9), it is only necessary to solve FORT and FODRT equations and to evaluate the deviations $\Delta\tau_{S1,S2}$ from eqs (A6) and polarizations from eqs (A7) and (A8). There is no need to solve the coupling equations (11).

Similarly as in eq. (A1), we can seek the solution of the second-order coupling equations (21) in the form

$$ \mathcal{A}_1(\tau) = \mathcal{A}'_1 \exp(i\omega\Delta\tau), \quad \mathcal{B}_1(\tau) = \mathcal{B}'_1 \exp(i\omega\Delta\tau). \quad (A11) $$

Inserting (A11) into eqs (21) and considering $\epsilon_1 \ll \epsilon_2$, we arrive at $d\mathcal{A}'_1/d\tau \sim 0$, $d\mathcal{B}'_1/d\tau \sim 0$ and

$$ 2 \frac{d\Delta\tau}{d\tau} \begin{pmatrix} \mathcal{A}'_1 \\ \mathcal{B}'_1 \end{pmatrix} + \begin{pmatrix} M_{11} - 1 & M_{12} \\ M_{12} & M_{22} - 1 \end{pmatrix} \begin{pmatrix} \mathcal{A}'_1 \\ \mathcal{B}'_1 \end{pmatrix} = 0. \quad (A12) $$

As in case of eq. (A2), we can consider (A12) as an eigenvalue problem for eigenvalues $1 - 2d\Delta\tau/d\tau$ and the corresponding eigenvectors situated in the plane specified by vectors $\mathbf{f}^{[K]}$. From the eigenvalues we can find

$$ \frac{d\Delta\tau}{d\tau} = \frac{1}{2} - \frac{1}{4}[(M_{11} + M_{22}) \pm \sqrt{(M_{11} - M_{22})^2 + 4M_{12}^2}]. \quad (A13) $$

Using eq. (19) and again taking into account that $B_{11}(x_m, p_m) + B_{22}(x_m, p_m) = 2$, we get

$$ M_{11} + M_{22} = 2 - (c^{[M]})^2 \frac{B_{13}^2 + B_{23}^2}{V_P^2 - V_S^2}. \quad (A14) $$

Using eq. (A14), and assuming that the common S-wave ray does not pass through a singularity, we can use eq. (A13) in the expressions for the approximate traveltimes $\tau_{S1}$ and $\tau_{S2}$ of faster S1 and slower S2 waves between points corresponding to $\tau_0$ and $\tau$ on the common S-wave ray

$$ \tau_{S1,S2}(\tau, \tau_0) = \tau^{[M]}(\tau, \tau_0) + \Delta\tau_{S1,S2}(\tau, \tau_0). \quad (A15) $$

The term $\tau^{[M]}(\tau, \tau_0)$ is the first-order traveltime calculated between points corresponding to $\tau_0$ and $\tau$ on the common S-wave ray. The term $\Delta\tau_{S1,S2}$ has now the following meaning

$$ 
\begin{aligned}
\Delta\tau_{S1} &= \Delta\tau^{[M]} - \frac{1}{4} \int_{\tau_0}^{\tau} \sqrt{(M_{11} - M_{22})^2 + 4M_{12}^2} d\tau, \\
\Delta\tau_{S2} &= \Delta\tau^{[M]} + \frac{1}{4} \int_{\tau_0}^{\tau} \sqrt{(M_{11} - M_{22})^2 + 4M_{12}^2} d\tau.
\end{aligned}
\quad (A16) $$

The term $\Delta\tau^{[M]}$ is the second-order traveltime correction given in eq. (9). We can thus see that the second-order traveltime correction is implicitly incorporated in the coupling equations (21). The second terms on the right-hand sides of eqs (A16) control separation of S waves. The approximate traveltime difference of the S1 and S2 waves is now $\frac{1}{2} \int_{\tau_0}^{\tau} \sqrt{(M_{11}-M_{22})^2+4M_{12}^2} d\tau$. Since it is related to the expressions, which contain second-order traveltime corrections calculated along the common ray, we call it the ‘second-order traveltime difference’.

Comparing eqs (A15) and (A16) with (A5) and (A6), we can see that the increased accuracy of (A15) is caused by two effects. First, by the second-order traveltime correction $\Delta\tau^{[M]}$ and second, by the more accurate estimate of traveltime difference, see the second term on the right-hand sides of eqs (A16). Formulae (A16) are independent of the choice of vectors $\mathbf{f}^{[K]}$.

The eigenvectors of eq. (A12) specify the direction of the first-order polarization vectors $\mathbf{f}^{[K]}$ of the two decoupled waves. They are situated in the plane specified by vectors $\mathbf{f}^{[K]}$, see eq. (13); at any point of the common S-wave ray, the vectors $\mathbf{f}^{[K]}$ can be expressed as follows:

$$ 
\begin{aligned}
\mathbf{f}^{[1]} &= \mathbf{f}^{[1]} \cos \Phi_1 + \mathbf{f}^{[2]} \sin \Phi_1, \\
\mathbf{f}^{[2]} &= -\mathbf{f}^{[1]} \sin \Phi_1 + \mathbf{f}^{[2]} \cos \Phi_1.
\end{aligned}
\quad (A17) $$

The polarization vectors $\mathbf{f}^{[K]}$ are not unit since vectors $\mathbf{f}^{[K]}$ are also not unit. Angle $\Phi_1$ is specified so that the element $M'_{12}$ of the rotated matrix $\mathbf{M}$ equals zero, $M'_{12} = 0$. This condition yields (Farra 2001; Farra & Pšenčík 2003)

$$ \tan 2\Phi_1 = \frac{2M_{12}}{M_{11} - M_{22}}. \quad (A18) $$
Angle $\Phi_1$ is chosen so that vector $\mathbf{f}^{[1]}$ corresponds to the S1 wave. The polarization vector $\mathbf{f}^{[2]}$ corresponds to the S2 wave.

We can now rewrite eq. (15) in the following ‘uncoupled’ form:

$$ \mathbf{U}(\tau) = \frac{\mathcal{D}'_1 \mathbf{f}^{[1]}(\tau)\exp(i\omega\Delta\tau_{S1}) + \mathcal{F}'_1 \mathbf{f}^{[2]}(\tau)\exp(i\omega\Delta\tau_{S2})}{[\rho(\tau)c^{[M]}(\tau)]^{1/2} \mathcal{L}^{[M]}(\tau)} \quad (A19) $$

In (A19), $\mathcal{D}'_1$ and $\mathcal{F}'_1$ are factors constant along the ray, which can be expressed in terms of the factors $\mathcal{A}'_1$ and $\mathcal{B}'_1$ given in (A11) in the following way:

$$ \begin{aligned} \mathcal{D}'_1 &= \mathcal{A}'_1 \cos \Phi_1(\tau_0) + \mathcal{B}'_1 \sin \Phi_1(\tau_0), \\ \mathcal{F}'_1 &= -\mathcal{A}'_1 \sin \Phi_1(\tau_0) + \mathcal{B}'_1 \cos \Phi_1(\tau_0). \end{aligned} \quad (A20) $$

The factors $\mathcal{A}'_1$, $\mathcal{B}'_1$ can be determined from the initial conditions, the angle $\Phi_1(\tau_0)$ from eq. (A18) applied at $\tau = \tau_0$. The terms $\Delta\tau_{S1}$, $\Delta\tau_{S2}$ are given in (A16), the polarization vectors $\mathbf{f}^{[k]}$ in (A17). We can see from eq. (A19) that, for $\epsilon_1 \ll \epsilon_2$, also the second-order coupling equations yield two separate S waves with the same spreading factor, but different polarizations and traveltimes.