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

We present a modiﬁcation 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 ﬁelds, 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 modiﬁed 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, ﬁrst 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 ﬁnd applications in migration and traveltime tomography


406
V. Farra and I. Pšenčík 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 Swave 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(ω −1 ), where ω 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 Swave 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 quasiisotropic 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.

B A S I C E Q UAT I O N S
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 Here i is the imaginary unit, ω is the circular frequency, τ (x m ) is the eikonal, which also serves as the traveltime, and 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 1 ∼ c/(ωL) is small (Pšenčík 1998). This requirement is often simplified to 1 ∼ ω −1 , see, for example, Paper I. Here we use the more general former definition of 1 , that is, 1 ∼ c/(ω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 Coupled S waves in weak anisotropy 407 weakly anisotropic medium the vectorial amplitude coefficient (2) In (2), A and B are S-wave amplitude coefficients and vectors e [1] and e [2] are mutually perpendicular unit vectors, to which the amplitude coefficients are related. Parameter τ 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 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 τ = τ [M] (x m ) is the firstorder traveltime. Symbol G [M] denotes the first-order S-wave mean eigenvalue In (4), G [1] and G [2] are the first-order approximations of two smaller eigenvalues of the generalized Christoffel matrix (we call it generalized because it contains components of slowness vector p instead of the unit vector in the direction of p used in the standard Christoffel matrix): Symbols a ijkl denote density-normalized elastic moduli, c ijkl being elements of the fourth-order tensor of elastic moduli and ρ 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 B(x m , p m ) Symbols e [j] i in (4) and (7) denote the components of unit vectors e [j] . Vectors e [1] and e [2] , see (2), are perpendicular to the third vector e [3] chosen so that e [3] = n. Here n is a unit vector specifying the direction of the first-order slowness vector p. Vectors e [K] can be chosen arbitrarily in the plane perpendicular to n. Vector e [3] can be determined from the second set of FORT eqs (3). Along the common S-wave ray, vectors e [K] can be computed as the vectorial basis of the 'wavefront orthonormal coordinate system', see, for example,Červený (2001): 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 secondorder traveltime correction, we refer to it as such, for the sake of simplicity. It reads Quantities B 13 and B 23 in eq. (9) are elements of the symmetric matrix B(x m , p m ), see (7). Let us mention that the traveltime correction (9) does not depend on the choice of basis vectors 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(τ ) of the common S wave given in eq. (2), calculated along the S-wave common ray, reads Parameter τ = τ [M] (x m ) is again the first-order traveltime along the S-wave common ray, obtained by solving FORT eqs (3). The symbols ρ, c [M] and L [M] denote the density and first-order common S-wave phase velocity and geometrical spreading. The spreading L [M] represents the first-order approximation of the exact common S-wave spreading computed along an exact common S-wave ray. Vectors e [1] and e [2] are the basis vectors of the wavefront orthonormal coordinate system (8). Vectors e [K] define the zero-order polarization plane of the common S wave. Thus the vectorial amplitude coefficient U(τ ) 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 In (11), B JL are again the elements of the symmetric matrix B(x m , p m ), see (7). Matrix B(x m , p m ) can be obtained by simple rotation from matrixB(x m , p m ), specified explicitly in Paper I. Because of the relation of matrix 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 F acting at time τ = τ 0 , the initial conditions for coupling equations (11) read Here c [M] 0 denotes the first-order approximation of the S-wave common phase velocity at the source, and ρ 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 2 ∼ || a ijkl ||/ ||a ijkl || ∼ c/c (||.|| denotes the tensor norm), characterizing weak anisotropy, is of the same order as the earlier introduced small parameter 1 . Here, a ijkl are perturbations of the density-normalized elastic parameters from the reference isotropic medium, c is the maximum difference of the phase velocities of the two S waves propagating in the weakly anisotropic medium.

S E C O N D -O R D E R C O U P L I N G E Q UAT I O N S
For greater accuracy of the results, it was proposed in Paper I to replace vectors e [K] , which define the zero-order polarization plane of the common S wave, by their first-order counterparts f [K]

V. Farra and I. Pšenčík
Vectors f [K] are situated in the plane perpendicular to the firstorder eigenvector f [3] corresponding to the largest eigenvalue of the generalized Christoffel matrix (5) Quantities B 13 , B 23 are again elements of matrix 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 e [K] by f [K] . More important is the use of the vectors f [K] in the rederivation of coupling equations, which leads to their significantly increased accuracy. After replacement of vectors e [K] by f [K] in (10), the zero-order vectorial ray amplitude coefficient U(τ ) reads 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(ω −1 ) and less. Within this approximation, the replacement of vectors e [K] by f [K] only affects the coefficient of (iω) 2 in eq. (31) of Paper I. The effects on the coefficient of iω in eq. (31) can be neglected. Thus, if we take into account the definition of matrix B(x m , p m ) in (7), eq. (32) of Paper I now reads For a specific choice of velocities V S and V P eq. (16) can be rewritten to read The elements of the 2 × 2 matrix M in (18) read Matrix 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 (19) can thus be rewritten to the alternative form Note that the element B 33 of the matrix 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 We call eqs (21) 'second-order coupling equations'. They are obtained from the first-order coupling equations (11) (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 e [K] are replaced by vectors f [K]

N U M E R I C A L E X A M P L E S
In order to illustrate the accuracy of the first-order and secondorder 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 singleforce 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 with the dominant frequency f = 50 Hz and γ = 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 • from the vertical plane (x, z). The S-wave anisotropy in the (x, z) plane defined as (c S1 − c S2 )/c average × 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 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 τ [M] (blue -1) obtained by solving FORT eqs (3) along the common S-wave ray and second-order traveltime τ [M] + τ [M] (red -2) from the average of exact S1-and S2-wave traveltimes. For τ [M] see eq. (9). The V S and V P 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 [M] 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 traveltime τ [M] 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 differences remain nearly constant for the studied range of depths. The relative difference of the first-order geometrical spreading L [M] 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). from the average of exact spreadings of S1 and S2 waves is approximately −3 per cent. The first-order geometrical spreading L [M] , 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 remaining 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 seismograms 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 approximation 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 traveltime approximation incorporated in eqs (11)  approximate seismograms (red) are computed from the secondorder 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 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 firstorder coupling equations (11) (black) and from the second-order coupling equations (21)  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.  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 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. 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 τ S1 in these equations represent identical traveltime corrections (therefore, the S1 wave is described equally well by both systems of coupling equations). The terms τ 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 secondorder 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.

D I S C U S S I O N A N D C O N C L U S I O N S
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  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 Figure 11. Comparison of the seismograms computed with the secondorder coupling equations (21) (red) and the standard ray theory seismograms (black) for the vertical single-force source in the QI4 model. 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.

A C K N O W L E D G M E N T S
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. 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 2 = 0, eqs (11) yield A 0 = const, B 0 = 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 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 ( c large), higher frequencies (ω large) or weaker inhomogeneity, including homogeneity (characteristic length L large or infinite), that is, for the case 1 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) Here A 0 and B 0 are the amplitude factors, τ are deviations of the traveltimes from the first-order traveltime τ [M] calculated along the common S-wave ray. Inserting (A1) into eqs (11) and considering 1 2 , we arrive at dA 0 /dτ ∼ 0, dB 0 /dτ ∼ 0 and 2 d τ dτ Eq. (A2) can be expressed in the form of an eigenvalue problem Eq. (A3) represents a system of two equations for two eigenvalues 1 − 2d τ /dτ 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 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 τ S1 and τ S2 of the S waves, between points corresponding to τ 0 and τ on the common S-wave ray τ S1,S2 (τ, τ 0 ) = τ [M] (τ, τ 0 ) + τ S1,S2 (τ, τ 0 ).
The term τ [M] (τ, τ 0 ) is the first-order traveltime calculated between points corresponding to τ 0 and τ on the common S-wave ray. If the common S-wave ray does not pass through a singularity, the term τ S1,S2 , which represents deviations of the traveltimes of S1 and S2 waves from τ [M] (τ, τ 0 ), has the following meaning τ S1 = − V. Farra and I. Pšenčík 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 1 2 τ τ 0 √ (B 11 −B 22 ) 2 +4B 2 12 dτ . Since the difference is related to the firstorder traveltime calculated along the common ray, we call it the 'first-order traveltime difference'. Formulae (A6) are independent of the choice of e [K] .
The eigenvectors of eq. (A3) specify the directions of the zeroorder polarization vectors of the two decoupled waves. We can introduce the polarization vectors e [K] as mutually orthogonal unit vectors situated in the plane specified by vectors e [K] , see eq. (8). At any point of the common S-wave ray, the vectors e [K] Angle 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 0 changes the matrix B into B , in which B 12 = 0. From this condition we get (Pšenčík 1998;Farra 2001;Farra & Pšenčík 2003) tan 2 0 = 2B 12 Angle 0 is chosen so that vector e [1] corresponds to the S1 wave. The polarization vector e [2] corresponds to the S2 wave.
The factors A 0 , B 0 can be determined from the initial conditions, the angle 0 (τ 0 ) from eq. (A8) applied at τ = τ 0 . The traveltime deviations τ S1 , τ S2 are given in (A6), the polarization vectors 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 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 D 0 and the polarization vector 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 τ 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 secondorder coupling equations (21) in the form Angle 1 is chosen so that vector f [1] corresponds to the S1 wave.
The polarization vector f [2] corresponds to the S2 wave.