G J e Ray tracing of multiple transmitted/reﬂected/converted waves in 2-D/3-D layered anisotropic TTI media and application to crosswell traveltime tomography

To overcome the deﬁciency of some current grid-/cell-based ray tracing algorithms, which are only able to handle ﬁrst arrivals or primary reﬂections (or conversions) in anisotropic media, we have extended the functionality of the multistage irregular shortest-path method to 2-D/3-D tilted transversely isotropic (TTI) media. The new approach is able to track multiple transmitted/reﬂected/converted arrivals composed of any kind of combinations of transmis-sions, reﬂections and mode conversions. The basic principle is that the seven parameters (ﬁve elastic parameters plus two polar angles deﬁning the tilt of the symmetry axis) of the TTI media are sampled at primary nodes, and the group velocity values at secondary nodes are obtained by tri-linear interpolation of the primary nodes across each cell, from which the group velocities of the three wave modes (q P , q SV and qSH ) are calculated. Finally, we conduct grid-/cell-based wave front expansion to trace multiple transmitted/reﬂected/converted arrivals from one region to the next. The results of calculations in uniform anisotropic media indicate that the numerical results agree with the analytical solutions except in directions of SV -wave triplications, at which only the lowest velocity value is selected at the singularity points by the multistage irregular shortest-path anisotropic ray tracing method. This veriﬁes the accuracy of the methodology. Several simulation results show that the new method is able to efﬁciently and accurately approximate situations involving continuous velocity variations and undulating discontinuities, and that it is suitable for any combination of multiple transmit-ted/reﬂected/converted arrival tracking in TTI media of arbitrary strength and tilt. Crosshole synthetic traveltime tomographic tests have been performed, which highlight the importance of using such code when the medium is distinctly anisotropic.

pseudospectral method (e.g. Crampin 1984) and the spectral element method (e.g. Komatitsch et al. 2000) to the governing partial differential wave equations, subject to appropriate source and boundary conditions. This yields the full waveform synthetic seismograms in which individual modes or arrivals are not explicitly identified, but rather the full elastodynamic response is computed. The other approach is the ray tracing method, which tracks the traveltime and trajectory of wave energy-flux for specific modes, subject to the high-frequency assumption (e.g.Červený 2001). Compared to the full wave solution, ray tracing is far less computationally expensive, especially for large 3-D simulations. Although it does not yield the full wave response, it still provides an effective way to approximately model the wavefield kinematics and dynamics, so important in the processing and interpretation of seismic data. Examples where ray tracing is frequently used in seismic exploration include traveltime tomography (e.g. , seismic migration (e.g. Alkhalifah & Larner 1994) and synthetic seismograms computations (e.g. Guest & Kendall 1993). Two different classes of ray tracing can be distinguished. One uses a fixed grid and computes the traveltime field for all nodes in space as time elapses, and is referred as to grid-/cell-based ray tracing method. The other uses the high-frequency approximation to the wave equation to trace rays by following traveltime fields as they evolve. Such a high-frequency ray tracing method can be a shooting or a raybending method. There are several examples of this approach to calculating traveltimes, their corresponding ray paths and even amplitudes in anisotropic media. Červený (1972) applied the method of characteristics and the numerical Runge-Kutta method to the first-order differential equations. Unfortunately, the equations for ray tracing in a general anisotropic medium are much more complicated than those for an isotropic medium, or a transversely isotropic medium involving a vertical axis of symmetry (VTI media). For a weakly anisotropic medium (Červený & Jech 1982) demonstrated a linearized approach to traveltime computation without any ray tracing performed. Gajewski & Pšenčik (1987) presented a 3-D ray tracing algorithm for the computation of high-frequency synthetic seismograms in a layered anisotropic medium. Shearer & Chapman (1988) gave the solutions of the ray path and traveltime in a constant gradient anisotropic medium. In some cases, such shooting or bending methods encounter difficulties. The former may miss the target even if different initial directions are applied. The latter may get trapped in a local minimum, even if many initial ray paths are perturbed, when the complexity of the anisotropic medium increases.
For the VTI case, involving a vertical symmetry axis, the mathematical treatment for the phase and group velocities (only five elastic parameters and with decoupled qSV and qSH wave modes), and shear wave singularity directions (Vavryčuk 2001) is greatly simplified. For example, Faria & Stoffa (1994) presented an iterative approach to the traveltime computation for a 2-D/3-D VTI medium. Rüger & Alkhalifah (1996) applied simple versions of the eigenvectors in a VTI medium and developed an efficient 2-D ray tracing method. Ettrich & Gajewski (1998) suggested a perturbation theory to compute traveltimes in VTI media, starting with a background velocity model that exhibits elliptical anisotropy. Cardarelli & Cerreto (2002) presented a ray tracing method in elliptical anisotropic media using linear traveltime interpolation.
Although Lagrangian algorithms can be used to compute reflections and/or refractions, and even multivalued rays in an anisotropic medium (Červený 2001), their applications will be impractical when the wavefield becomes too complicated. In such cases, it is impossible to distinguish whether the computed ray is a first or later arrival.
Alternatively, for tracing first arrivals of some seismic phases, the Eulerian grid-/cell-based wave front expanding algorithms, such as finite-difference eikonal equation solver (Vidale 1988(Vidale , 19902-D/3-D) and the shortest-path method (2-D, Moser 1991;3-D, Klimeš & Kvasnička 1994), are a favourable choice to conduct this kind of ray tracing in complex isotropic media and are also valid in anisotropic media. For example, Qian & Symes (2002) developed a finite-difference method to calculate the traveltime of qP waves in anisotropic media. Alkhalifah (2002) demonstrated another traveltime computation scheme by solving a linearized eikonal equation for 2-and 3-D VTI media. However, in general anisotropic media, the ray paths cannot immediately be retrieved from the computed traveltimes, because the phase-slowness vector (the gradient of the traveltime) may differ from the group-velocity vector (direction of the ray path). Such ray paths are required for the purpose of tomography.
It should be appreciated that the subsurface structure may be a composite assembly of tilted transversely isotropic (TTI) media with variable orientations of the symmetry axis due to folding, dislocation and other tectonic movements subsequent to deposition. This means that the direction of the symmetry axis is not constant (e.g. vertical or horizontal) in a complex anisotropic medium. Therefore, general anisotropic TI media should be used to define any symmetry axis of the anisotropic media. Zhou & Greenhalgh (2004, 2006 investigated first arrival and primary reflected arrival tracking in the most general 2-D/3-D anisotropic media using the shortest-path method. Later they developed a non-linear traveltime tomography algorithm for 3-D strongly anisotropic media (Zhou & Greenhalgh 2008). Following this research direction, we extend in this contribution the functionality of the multistage irregular shortest-path method (multistage ISPM, Bai et al. 2009Bai et al. , 2010, previously used for 2-D/3-D complex layered isotropic media, to 2-D/3-D general anisotropic TTI media and trace multiple transmitted/reflected/converted arrivals. The main advantages of the multistage ISPM method are its simplicity, robustness and ability to deal with a general anisotropic medium. Here, for reasons of simplicity, we restrict ourselves to general TTI media. The structure of the paper is as follows. First, we give the formulae for calculating the phase and group velocities of the three different wave modes (qP, qSV and qSH), and briefly summarize how the multistage ISPM works in 2-D/3-D anisotropic TI media. Then we use synthetic examples having known analytical solutions to demonstrate the accuracy of the method, the reduced central processing unit (CPU) time consumption and the overall performance of the method. Finally, we give a crosshole synthetic data example of multiple phase tomographic inversion in TTI media, and compare with the results of commonly applied isotropic inversion of such anisotropic data.

P H A S E A N D G RO U P V E L O C I T Y I N T I LT E D T I M E D I A
To obtain the phase and group velocities for a general anisotropic TTI medium, we first consider the relatively simple case of a VTI medium (only five elastic parameters). Subsequently, we will rotate the symmetry axis into an arbitrary orientation defined by the inclination angle (θ 0 ) and azimuth angle (ϕ 0 ) to obtain the more general result for a TTI medium. Here the inclination angle θ 0 is measured from the z-axis and the azimuth angle ϕ 0 is measured from the x-axis (see Fig. 1). In this coordinate system the wave front normal direction (or slowness vector) is given by the unit vector n = (sin θ cos ϕ, sin θ sin ϕ, cos θ ). This is the direction of phase velocity propagation. The group velocity vector direction, or ray path direction, is given by the unit vectorr = (sin θ 1 cos ϕ 1 , sin θ 1 sin ϕ 1 , cos θ 1 ).
For a VTI medium there is no azimuthal dependence of the phase velocity and so we can set ϕ = 0 in the above expression. The phase velocity for each mode is then simply a function of the direction of the slowness vector (or wave front normal direction) as well as the elastic constants, and given by (Daley & Hron 1977): where c 1 , c 2 and c 3 are the phase velocities of the qP, qSV and qSH waves, respectively. The quantity υ is the inclination angle of the slowness vector [(i.e. the angle between the wave front normal and the symmetry axis (a 11 , a 13 , a 33 , a 44 , a 66 )] are the five elastic moduli parameters, which are functions of spatial position. In eq. (1) the parameters P and Q are defined as follows: Q 1 = a 44 + (a 11 − a 44 ) sin 2 υ Q 2 = a 33 + (a 44 − a 33 ) sin 2 υ Q 3 = 1 4 (a 13 + a 44 ) 2 sin 2 2υ . ( Applying Crampin's (1981) formulae (see also Berryman 1979;Thomsen 1986) we can obtain the group velocities for the three different modes, The subscripts m = 1, 2, 3 stand for qP, qSV and qSH, respectively. It is possible to express the group velocity vector in terms of its horizontal (superscript h) and vertical (superscript z) components (Berryman 1979), In eq. (5) the phase velocity derivative is given by, where ∂ P ∂υ = 1 2 (a 11 − a 33 ) sin 2υ, − 1 2 (a 44 + a 13 ) 2 sin 4υ.
To make eqs (1)-(7) applicable to a general TTI medium (having arbitrary tilt of the symmetry axis), we can use the standard cosine expression for the dot product between two unit vectors to obtain the angle υ subtended between the arbitrary symmetry axis direction and the slowness vector direction: cos υ = sin θ 0 sin θ cos(ϕ 0 − ϕ) + cos θ 0 cos θ.
Here the two pairs of angles (θ 0 , ϕ 0 ) and (θ, ϕ) specify the orientation of the symmetry axis of the TTI medium and the phase-slowness (wave front normal) direction, given by n = (sin θ cos ϕ, sin θ sin ϕ, cos θ ), respectively. Substituting eq. (8) into eqs (1)-(7) and using the standard trigonometric expressions for projecting the vector components from the natural rock frame into the geographic (recording) frame yields the group velocities for a general TI medium.
Here the superscripts z and h denote values in directions parallel and perpendicular to the symmetry axis (given by eq. 5) whereas superscripts 1, 2 and 3 denote values in the geographic frame (or x, y and z). In the 2-D case, ϕ = ϕ 0 = 0, and eq. (9) simplifies to, From the above discussion it is clear that a 3-D general TTI medium may be described by the model vector m(x) given as where all the components of the model vector vary with the spatial coordinates x = (x, y, z). In the 2-D case it is assumed that the symmetry axis lies entirely in the x-z plane and so there is only one angle θ 0 to be considered (φ 0 = 0). Substituting (11) into (1)-(7), one obtains three group velocities {U m (x, n), m = 1, 2, 3} for the 2-D case, which correspond to the three wave modes (qP, qSV and qSH-see Zhou & Greenhalgh 2004). In addition to the slowness direction n, the group velocities are functions of the spatial position coordinates x. The mapping from the elastic moduli parameters to the group-velocity is fundamental for the ray tracing method described in the next section.

Model parametrization
To account for irregular interfaces and velocity discontinuities (such as low-or high-velocity layers), we can utilize irregular cells near the interfaces, and use regular rectangular cells elsewhere. Fig. 2(a) is a schematic display of the model parametrization for the multistage ISPM in 2-D layered media having undulating interfaces (including the surface topography) and an embedded low-velocity layer, in which triangular (Fig. 2c) and trapezoidal (Fig. 2d) cells are exploited near the irregular interfaces and rectangular (Fig. 2b) cells are used elsewhere. In the 3-D case, the trapezoidal mesh can be divided into rectangular prismatic (or cubic) cells and a tetrahedral mesh (see Fig. 3, referred to Bai et al. 2010).
In both the 2-D and 3-D cases, the nodes at the cell corners are referred to as primary nodes and the nodes along the cell boundaries (or cell surfaces in 3-D) are called secondary nodes. The benefit of introducing the secondary nodes is that one can use a relatively large cell size in the model parametrization and still maintain a good ray coverage (more secondary nodes are applied) and hence a high computational accuracy in the ray tracing process (Bai et al. 2007). The velocity field is in general sampled at the primary nodes and a specified velocity function (a bilinear function for two dimensions and a tri-linear Lagrangian interpolation function for three dimensions) is defined across the cell, which links the primary and secondary nodes (including the source and receiver positions in a cell). For secondary node velocity interpolation, the reader is referred to Bai et al. (2010), but here the only difference is that the group velocity and the ray direction in an anisotropic medium are used. For example, to obtain the group velocity at a fixed ray direction of the secondary node position, the group velocities for the same ray direction at the primary node positions within the same cell are used to conduct the bi-linear or tri-linear interpolation.

Ray tracing
There are well-established algorithms for ray tracing in anisotropic media, for example, it is easy to trace rays when the group velocity is obtained for a given slowness vector direction (i.e.Červený 2001). In the context of the shortest-path method (i.e. a grid-/cell-based ray tracing algorithm), to calculate the minimum traveltimes and locate the associated ray paths for all grid nodes, one can gradually expand the volume of the computed nodes by continually adding the undetermined neighbouring nodes of the cells to the computed Figure 2. A schematic explanation for computing the transmitted, reflected and mode-converted arrivals using the multistage ISPM in 2-D layered anisotropic media with undulating interfaces (diagram a). Diagrams (b-d) are grid cells used in the model parametrization (black and grey circles indicate primary and secondary nodes, respectively). In (a) the white line indicates the ray path for the reflected and mode-converted phase SV 1 P 1 P 1 P 2 P 2 SV 2 SV 3 P 3 P 2 P 1 in the six consecutive stages of the ray tracing process.  nodes. In this process one should start with the node that has a minimum traveltime in the subset N j (whereN j is the total set of computed nodes in the current computed wave front) to keep track of the first arrival times for the undetermined nodes. By using an interval sorting method (Klimeš & Kvasnička 1994), it is possible to delete the larger traveltimes, while recording the minimum traveltime and corresponding ray path.
The minimum traveltime from a source to an undetermined node j in a cell is expressed by the following equation: where U m = U m (x, θ 0 , r i j ) for two dimensions, and U m = U m (x, θ 0 , φ 0 , r i j ) for three dimensions. r i j and R(x i , x j ) are the ray direction and distance from node x i to node x j , respectively. Furthermore, in this wave front expanding process the order number of the incident (or previous) node i * (treated as a new source at the wave front) giving the minimum traveltime to node j is recorded for establishing the coordinates of the related ray path. In such a fashion, the minimum traveltimes and the corresponding ray paths within the gridded model are calculated simultaneously. Note that at present the ISPM algorithm cannot deal with multivalued qSV group velocities, which can arise in certain slowness directions. With the understanding that only the lowest (minimum) group velocity is selected in such cases, the method acts as a usual first arrival modelling algorithm. Since triplications of the qSV group velocity sheets show precursors (see Figs 4 and 7), it is questionable whether the ISPM algorithm generates first qSV arrivals in cases of strong anisotropy. If a wave front triplication occurs then an eigenvalue or eigenvector method can be used to deal with the phase or group velocity at the singularity points (for details, see Zhou & Greenhalgh 2004).  top and bottom interfaces of the middle layer, and then parametrize each layer by means of regular and irregular cells. Secondary nodes are added according to the specified secondary node density or spacing. For simplicity of wave identification, in our phase convention: P or SV or SH represents a qP or qSV or qSH wave, respectively, and the number in the subscript or superscript indicates a downward or upward propagating seismic wave in the different layers, given by the number, respectively. For example, the ray path for the multiple reflected and converted phase indicated in Fig. 2 is denoted as SV 1 P 1 P 1 P 2 P 2 SV 2 SV 3 P 3 P 2 P 1 in this phase naming convention.

Multistage computational scheme
In short, if we consider the ray that transmits, reflects and/or mode-converts at the upper interface, then a simulated downward wave front (P 1 phase) is propagated through the upper layer until it impinges on all sampled nodes of the first interface. At this stage the independent computational domain is halted at the upper layer and we are left with traveltime values defined along the sampled upper interface. From here, a downward propagation of a pure transmitted qP-wave (P 1 P 2 ) or a transmitted and mode-converted phase (P 1 SV 2 ) can be simulated by reinitializing it, starting at the sampled node position with the minimum traveltime on the first interface. From Huygens Principle, the minimum time node is treated as a new source point on the wave front to propagate into the middle layer. Meanwhile, an upward-propagating wave front consisting of a pure reflected phase (P 1 P 1 ) or a mode-converted reflected phase (P 1 SV 1 ) can also be predicted by reinitializing the wave front and starting at the sampled node position with the minimum traveltime, from the upper interface back into the upper layer (incident  In summary, the multiple transmitted/reflected/converted arrivals are the different combinations or conjugations, via the velocity discontinuities (i.e. interfaces) of the incident, transmitted, reflected (or refracted) and converted phases which obey Snell's Law, Fermat's Principle and Huygen's Principle. Note that in some cases, it may happen that the reflected/transmitted wave hits the interface at some location earlier than the recorded arrival time of the incoming wave (e.g. a head wave is generated). In such circumstances, we can define the interface not as a sampled line (or surface) but as a sampled region for two dimensions (or volume for three dimensions) of reasonably narrow width so that one can account for head wave ray tracing (see Table 1, Pn or Sn, Huang et al. 2013).

B E N C H M A R K T E S T S A N D C O M P U TAT I O N A L E F F I C I E N C Y
To check the algorithm accuracy in 2-D/3-D anisotropic TI media, we first show a comparison between the numerical results and the analytical solutions in a uniform 2-D TI medium, and then provide benchmark test results against the analytical solutions in two uniform 3-D TI media and finally discuss its performance in terms of computational accuracy and efficiency.  (a) Primary reflections from interface 1 (P 1 P 1 ), from interface 2 (P 1 P 2 P 2 P 1 ) and from model base (P 1 P 2 P 3 P 3 P 2 P 1 ); (b) primary converted reflections from interface 1 (P 1 SV 1 ), from interface 2 (P 1 P 2 SV 2 SV 1 ) and from model base (P 1 P 2 P 3 SV 3 SV 2 SV 1 ); (c) three-times pure P-wave reflections between interfaces 1 and 2 (P 1 P 2 P 2 P 2 P 2 P 1 ) and between the interfaces 1 and model base (P 1 P 2 P 3 P 3 P 2 P 2 P 3 P 3 P 2 P 1 ) and (d) multiple transmitted and mode converted reflections from interface 2 (SV 1 P 2 SV 2 P 2 SV 2 P 1 ) and from the model base (SV 1 P 2 SV 3 P 3 SV 2 P 2 SV 3 P 3 SV 2 P 1 ). × 0.2 km and the secondary node spacing is 20 m. The source is located at the centre of the model surface (2.5 km, 0.0 km). The first arrivals of P, SV and SH for the three cases are calculated and the corresponding contour lines of the traveltime field are depicted in Fig. 4, in which the white lines are the analytical solutions whereas the black lines are the numerical results. The analytical solutions are obtained from the formulae given in Section 2. From Fig. 4 it is clear that the numerical results coincide with the analytical solutions, except for the singular directions (triplications) in the SV wave. The traveltime contours change with the orientation angle of the symmetry axis. The results for the SV wave show that the ISPM method fails to calculate the cusps (triplications of traveltime), because we are unable to use all three values of the group velocity in a given ray direction. The shortest path method is limited to just one value, the minimum. The analytical solutions circumvent the problem because they can express group velocity in terms of the wave front normal direction (rather than the ray direction) which is a single-valued function. Inability to handle SV triplications is a disadvantage of the ISPM numerical method, but it works perfectly well for the P and SH waves. Note that the traveltime contours for SH wave are concentric circles due to the fact that for this model a 44 = a 66 = 1.0 (i.e. the velocity value for the SH wave is constant).
To visualize the spatial distribution in three dimensions of the wave front patterns for all three wave modes we have selected four different models. The elastic moduli parameters and the corresponding Thomsen parameters for the four different models are tabulated in Tables 1 and 2, respectively. Fig. 6 is a plot of the analytic wave fronts for the three modes (from left to right is P, SV and SH) at t = 1 s in the two VTI models (upper panels for Model 1 and lower panel for Model 2). From the figure the significant anisotropy is evident by the departure in shape from a simple sphere. The three wave modes Figure 10. Ray paths for different kinds of multiple transmitted, mode-converted reflections in the undulating layered TTI medium. (a) Pure primary P wave reflections (P 1 P 1 , P 1 P 2 P 2 P 1 and P 1 P 2 P 3 P 3 P 2 P 1 ) from three different interfaces; (b) twice converted reflections (SV 1 P 1 P 1 , SV 1 P 1 P 2 P 2 P 1 and SV 1 P 1 P 2 P 3 P 3 P 2 P 1 ) from three different interfaces; and (c, d) transmitted, converted and reflected arrivals from interface 3, corresponding to the incident P wave for (c) (phase P 1 SV 2 P 3 SV 3 P 2 SV 1 ) and incident SV for (d) (phase SV 1 P 2 SV 3 P 3 SV 2 P 1 ). have quite different wave front shapes or group velocity patterns. In particular, the two quasi-shear waves (SV and SH) split, propagating in different wave fronts. The SV wave has cusps (triplications) in the direction of the symmetry axis. The wave velocities of the three modes (for each VTI model) are constant in the plane perpendicular to the symmetry axes, but change with increasing departure angle from this plane. Fig. 7 is a plot of the equal analytic group velocity traveltime surfaces for the three different wave modes (from left to right is P, SV and SH) at t = 1 s for the two TTI models (upper panels for Model 3 and lower panels for Model 4). Comparing Figs 6 and 7, it is apparent that the wave front shapes are the same except for the rotation with respect to inclination angle θ 0 = 45 0 and azimuth angle ϕ 0 = 315 0 . The combination of elastic parameters underlying Figs 6 and 7 yields negative values for the Tsvankin & Thomsen (1994) parameter [σ = ( α 0 β 0 ) 2 (ε − δ)], which gives rise to the triplications along the symmetry axis in each case.

Computational efficiency in the 3-D case
In previous studies (Bai et al. 2007(Bai et al. , 2009(Bai et al. , 2011 we concluded that the efficiency of the ISPM depends on both the cell size and the number of secondary nodes. For the 2-D situation the computational efficiency is not an important issue, but for 3-D problems it can be important and so it is instructive to analyse the CPU time of ray tracing in 3-D anisotropic media. A uniform TTI model of size 10 km × 10 km × 10 km having the same elastic parameter values as underlying Fig. 5 (see previous section) and symmetry axis angles θ 0 = 45 0 and ϕ 0 = 0 0 was used for the test. We selected four different cell sizes (500 m, 250 m, 200 m and 125 m) in the model parametrization. For each different cell-sized model five different numbers of secondary nodes (3, 5, 7, 9 and 11) were used to progressively increase the ray coverage. The source is located at the centre of model base and 101 × 101 receivers were arranged over the top surface of the model (uniformly distributed). Results for all 20 different combinations of cell size and secondary node density were obtained for direct P arrivals. The traveltime errors Figure 11. Wave front simulation of multiple transmitted, mode-converted reflections (phase: P 1 P 2 P 3 P 3 P 2 P 2 P 3 P 3 P 3 P 3 P 2 P 1 ) in successive fashion in undulating layered TTI media (wave front interval is 0.2 s). The white lines are the two buried interfaces.  Figure 12. (a) 3-D ray paths of primary reflected phase (P 1 P 1 , blue line) from subsurface interface 1 and (P 1 P 2 P 2 P 1 , green line) from the model base and (b) 3-D ray paths of twice converted reflection phase (P 1 SV 1 P 1 , blue line) from the subsurface interface and (P 1 SV 1 P 2 SV 2 P 1 , green line) and from the model base.
(measured against the analytical solutions) along with the CPU time for each run are graphed in Fig. 8 (the PC computer used is a Lenovo Pentium(R) D-2.39 GHz). The average relative error drops significantly as the number of secondary nodes increases (or as the secondary node spacing decreases) for a fixed cell size. (Fig. 8a), but slightly increases as the cell size increases for a fixed secondary node density (Fig. 8b). That is, the computing error is dependent on the secondary node spacing if a suitable cell size is used in the model parametrization. We also observe from Figs 8(b) and (d) that the CPU time increases dramatically with an increase in the number of secondary nodes and as the cell size decreases. There is a clear trade-off between accuracy and efficiency. It is recommended that a relatively large cell size be used in the model parametrization but more secondary nodes be applied to increase the ray coverage. In this way it is possible to obtain the required level of accuracy without undue extra computational effort.

Multiple transmitted/reflected/converted arrivals in 2-D TTI media
For tracking multiple transmitted/reflected/converted arrivals in anisotropic media, we selected a four-layered TTI model having two horizontal interfaces (upper and basal) and one tilted interface (see Fig. 9). The elastic moduli and the inclination angle of the symmetry axis for the three distinct layers above the basal interface are listed in Table 3. The source is located at the centre of the model surface (4.0 km, 0.0 km), and 33 receivers were uniformly arranged along the top model surface. Fig. 9 shows the ray paths for various transmitted, reflected and converted waves from the three different interfaces-diagram (a): ray paths for pure primary P reflections (phases: P 1 P 1 , P 1 P 2 P 2 P 1 and P 1 P 2 P 3 P 3 P 2 P 1 ) from the three interfaces; diagram (b): ray paths for primary converted (P→SV) reflections (phases: P 1 SV 1 , P 1 P 2 SV 2 SV 1 and P 1 P 2 P 3 SV 3 SV 2 SV 1 ) from the three interfaces; diagram (c): ray paths for three-times reflected pure P waves between interfaces 1 and 2 (P 1 P 2 P 2 P 2 P 2 P 1 ) and between interfaces 1 and 3 at the base of the model (P 1 P 2 P 3 P 3 P 2 P 2 P 3 P 3 P 2 P 1 ); and diagram (d): ray paths for multiple transmitted and converted reflections from interface 2 (SV 1 P 2 SV 2 P 2 SV 2 P 1 ) and from interface 3 (SV 1 P 2 SV 3 P 3 SV 2 P 2 SV 3 P 3 SV 2 P 1 ). The multistage ISPM scheme (Bai et al. 2010) is capable of computing traveltimes in the presence of undulating interfaces with lateral variations in the model parameters. To illustrate its performance in such cases for anisotropic media, we have selected the same scale length model and elastic parameters as for Fig. 9, but now impose undulating surface topography and layer interfaces (see Fig. 10). The tilt angles for the anisotropy symmetry axes are 0 • , 45 • and 90 • for layers 1, 2 and 3, respectively. The source is located at the far left boundary at a depth of 1.0 km and 21 receivers are uniformly arranged along the model top surface. Fig. 10(a) shows the ray paths for pure primary P reflections (P 1 P 1 , P 1 P 2 P 2 P 1 and P 1 P 2 P 3 P 3 P 2 P 1 ) from the three layer boundaries, whereas Fig. 10(b) shows the ray paths for the two-times mode-converted reflections (SV 1 P 1 P 1 , SV 1 P 1 P 2 P 2 P 1 and SV 1 P 1 P 2 P 3 P 3 P 2 P 1 ) from the same interfaces. Figs 10(c) and (d) depict the same ray paths for transmitted, converted and reflected arrivals from the interface 3, assuming for incident P in diagram (c) (phase P 1 SV 2 P 3 SV 3 P 2 SV 1 ) and incident SV in diagram (d) (phase SV 1 P 2 SV 3 P 3 SV 2 P 1 ).
To visualize wave propagation in the anisotropic medium, we selected the same sized model and the undulating interfaces as in the Fig. 10, but with different elastic moduli parameters and symmetry axes for each layer (see Table 4), and shown in Fig. 11 are the wave fronts for the pure transmitted and reflected phase P 1 P 2 P 3 P 3 P 2 P 2 P 3 P 3 P 3 P 3 P 2 P 1 . The simulated wave emanates from the source, propagates downwards and transmits through interfaces 1, 2 and reaches the bottom interface (diagram a). It then reflects back from the bottom interface and transmits upward through interface 2 and arrives at interface 1 (diagram b). From here it reflects downwards from interface 1 and transmits through interface 2 and reaches the basal interface (diagram c). Next it reflects back upwards again from this interface and propagates to interface 2 (diagram d). It then reflects back downwards from the interface 2 and reaches the bottom interface again (diagram e), from where it reflects back upwards and propagates all the way to the top surface, passing through interfaces 2 and 1 (diagram f).

Multiple transmitted/reflected/converted arrivals in 3-D TTI media
In this section we consider a 3-D TTI model having dimensions of 50 km × 50 km × 50 km and comprising two distinct layers above a half-space. The surface topography and the bottom interface are undulating whereas the top interface separating layers 1 and 2 is planar and dipping. The elastic parameters and symmetry axis angles for the two layers are as listed in Table 5. For the model parametrization, we selected a cell size (2.0 km × 2.0 km × 2.0 km) and used seven secondary nodes (at a spacing of 0.25 km) in each direction of the cell surface. The source is located at a depth of 10 km on the centre line of the model. The 31 receivers are arranged along the back edge and two side boundaries of the top topographic surface. Fig. 12(a) shows the ray paths computed for the pure P-wave primary reflections from the first interface (P 1 P 1 ) and the second interface (P 1 P 2 P 2 P 1 ). The twice mode-converted reflections between the model surface and the first interface (P 1 SV 1 P 1 ), and between the model surface and the second interface (P 1 SV 1 P 2 SV 2 P 1 ) are computed and the corresponding ray paths shown in Fig. 12(b). For clarity of presentation we have shown only these primary and twice mode-converted reflections but more complicated multiple transmitted/reflected/converted arrivals can easily be modelled and incorporated with other subroutines, such as an inversion solver, for tomographic purposes. Note that for simplicity of presentation purposes, we only show some examples for media with homogeneous layers in the 2-D and 3-D cases, but lateral variations in the elastic parameters are also possible. This is one distinct advantage of the ISPM algorithm; it is capable of finding the ray paths in a medium with high contrast elastic parameters.   Zhou & Greenhalgh (2008) developed two kinds of non-linear anisotropic inversion algorithm to directly reconstruct the elastic moduli or the Thomsen parameters for various anisotropic TTI models by using the first arrivals of the P, SV and SH waves. One is referred to as the Thomsen parameter scheme and the other is called the elastic modulus scheme. In this section we used the elastic Figure 16. Results of the elastic model parameter inversion using direct qSV arrivals and pure reflected qSVqSV arrivals. modulus scheme and extended its functionality by exploiting more arrivals (direct waves + reflected waves) to image the structure of the five elastic constant model parameters. Fig. 13 shows a faulted VTI model and its elastic parameters. The corresponding Thomsen parameters for background medium and the faulted layer are listed in Table 6. The crosswell recording geometry for this synthetic tomography experiment entails 21 sources located in the left borehole (left edge of the model) at an equal spacing over the 600 m depth range and 21 receivers uniformly arranged in the opposite borehole at the right edge of the model. The model base is taken to be a Figure 17. Results of the elastic model parameter inversion using direct qSH arrivals and pure reflected qSHqSH arrivals. reflecting interface. Traveltimes for both direct P, SV and SH waves and pure reflected PP, SVSV and SHSH waves from the model base were computed and used as the input data for tomographic purposes. To appreciate the ray coverage, we show, in Fig. 14, the ray paths for three kinds of direct arrivals (grey lines) and pure reflections (black lines) from the model base for the P-wave (left panel), the SV-wave (middle panel) and the SH-wave (right panel). From Fig. 14, it is clear that the ray coverage is good in the central bottom part of the model, but not so good in the central top part and in the neighbourhood of the two embedded anomalies.

A P P L I C AT I O N T O C RO S S H O L E T O M O G R A P H Y
To simulate the real case, we added ±50 ms and ±100 ms maximum random noise levels to the traveltimes of the direct waves (qP, qSV and qSH) and pure reflected waves (qPqP, qSVqSV and qSHqSH), respectively. First, we use the background velocitiesV P = √ a 33 , V SV = √ a 44 , V S H = √ a 66 ) as the initial model parameters and undertake a non-linear isotropic inversion (see Bai & Greenhalgh 2005) to reconstruct the P, SV and SH velocity structures (results not shown for reasons of economy of presentation). It is clear from these results that the faulted layer is not recovered, and even the background velocity structure is incorrect. Inverting anisotropic data with an isotropic algorithm totally fails despite exploiting both direct and reflected arrivals. By contrast, it is possible to recover the five elastic moduli parameters for each layer (background plus faulted bed) by using an anisotropic inversion code (see Zhou & Greenhalgh 2008) with the combined arrival time data. By using the background model (i.e. a 11 = 15.1, a 13 = 1.6, a 33 = 10.8, a 44 = 3.1 and a 66 = 4.3) as the starting model, Fig. 15 shows the four elastic model parameters (a 11 , a 13 , a 33 and a 44 ) reconstructed by combining the direct qP and reflected qPqP traveltimes. From the figure it is clear that the four elastic moduli parameters are well recovered both in terms of anomaly shape and magnitude, except for the parameter a 33 , which is distorted to some degree. This is partially due to the high contrast in the a 33 values between the anomaly and the background. Note that the above elastic moduli inversion scheme is able to recover the structure of the 'shear wave-related' parameter (a 44 ) from the P-wave data. Fig. 16 displays the same four elastic moduli parameters reconstructed by combining direct the qSV and reflected qSVqSV data. Only the elastic modulus parameter a 44 is recovered in the right position and of the correct strength, whereas parameter a 13 is roughly in the right position but not of the correct magnitude. Fig. 17 depicts the tomograms obtained for the other two elastic moduli parameters (a 44 and a 66 ) by combining the direct qSH and reflected qSHqSH arrivals times. From Fig. 17 it can be seen that the parameter a 66 is recovered in the right position and magnitude, but for a 44 only the anomaly strength is reliably recovered. Furthermore, there exist two image artefacts (see left panel of the Fig. 17).
From Figs 15-17 it is evident that all elastic model parameters are recovered in the right positions and approximately the right magnitudes . Furthermore, the reconstructed parameter a 44 for the faulted layer is in the correct position, regardless of whether combined qP + qPqP, qSV + qSVqSV, or qSH + qSHqSH arrival times are used. From this synthetic data inversion test it is obvious that the tomographic solution is strongly dependent on using the appropriate Figure 18. Comparisons between the non-linear VTI and TTI inversions using TTI data of direct qP arrivals and pure reflected qPqP arrivals.
inversion algorithm if anisotropy is present. The traditional (isotropic) inversion method cannot be expected to work in such situations. It was assumed in the above anisotropic inversion that the orientation angles for the axis of symmetry (VTI case) is known in advance.
Next we show the difference between the VTI and TTI inversions when the data conform to a tilted symmetry axis. The five elastic moduli parameters are the same as mentioned above (Fig. 13). The corresponding Thomsen parameters are equivalent for the model (see Table 6), but the symmetry axis is now rotated through an angle of 45 • (TTI case). Figs 18-20 show the tomograms for both the VTI and TTI inversions using three different combinations of data (Fig. 18, results for qP + qPqP; Fig. 19, results for qSV + qSVqSV and Fig. 20, results for qSH + qSHqSH). The VTI inversion is unable to properly image the structure, regardless of which combination of arrival time data qP + qPqP, or qSV + qSVqSV, or qSH + qSHqSH, is used. Both the position and magnitude of the anomalous faulted bed are incorrect. In fact, the faulted bed gets rotated into the Figure 19. Comparisons between the non-linear VTI and TTI inversions using TTI data of direct qSV arrivals and pure reflected qSVqSV arrivals. direction of the TTI symmetry axis (see two left panels of Figs 18 and 19, and the top panels of Fig. 20). Note that parameter a 11 is not sensitive to the symmetry axis, but only the shape of the anomaly is recovered. By contrast, TTI inversion of the data (assuming the correct tilt angle) recovers most elastic moduli parameters (see the two right panels of Figs 18 and 19, and the bottom panels of Fig. 20). It is therefore important to know in advance the rough model structure (i.e. whether the medium is anisotropic or not and if so, the approximate tilt angle). Of course, the tilt angle could be an additional unknown to solve for in TTI inversion but as shown by Zhou & Greenhalgh (2008) the large differences in sensitivities between the angles and the elastic parameters make this problematic. Figure 20. Comparisons between the non-linear VTI and TTI inversions using TTI data of direct qSH arrivals and pure qSHqSH arrivals.

C O N C L U S I O N S
Multiple transmitted, mode-converted and reflected (or critically refracted, as a special case) wave fronts of seismic waves are caused by discontinuous variations in the wave speed and/or ray directions in anisotropic media. The additional phases yield more subsurface information than the direct (first arrival) waves alone. However, in traditional earthquake location, traveltime tomography and migration imaging, only the information from the first (primary) arrivals is normally used. Nevertheless, the primary reflections or other later arriving phases are important for improving the hypocentre location accuracy and also the tomographic resolution because improved angular coverage of the target and more constraint parameters are provided in seeking a physical solution. We have therefore extended in this study the multistage ISPM scheme to anisotropic media and developed the multiple transmitted/reflected/converted arrival tracking algorithm for general TTI media (including as special cases VTI and HTI media). The numerical tests, in terms of accuracy and computational economy, suggest that it is a practical and efficient algorithm to trace multiple transmitted/reflected/converted arrivals in general 2-D/3-D TTI media. With the irregular cell interpolation technique and the multistage computational scheme, it is a promising way to predict the multiple transmitted, converted and reflected waves in general 2-D/3-D layered TTI media involving curved interfaces. Numerical (synthetic) experiments of crosshole tomography data clearly show the significant influence of anisotropy and the need to incorporate it into the inversion. Note that resolving the orientation of the symmetry axis in the inversion remains difficult and it is not attempted in this work.