Ray tracing of turning wave in elliptically anisotropic media with an irregular surface

Seismic ray tracing in anisotropic media with irregular surface is crucial for the exploration of the fine crustal structure. Elliptically anisotropic medium is the type of anisotropic media with only four independent elastic parameters. Usually, this medium can be described by only the vertical phase velocity and the horizontal phase velocity for seismic wave propagation. Model parameterization in this study is described by flexible triangular grids, which is beneficial for the description of irregular surface with high degree of approximation. Both the vertical and horizontal phase velocities are defined in the triangular grids, respectively, which are used for the description of phase velocity distribution everywhere in the model by linear interpolation. We develop a shooting ray tracing method of turning wave in the elliptically anisotropic media with irregular surface. Runge-Kutta method is applied to solve the partial differential equation of seismic ray in elliptically anisotropic media. Linearly modified method is used for adjusting emergent phase angles in the shooting scheme. Numerical tests demonstrate that ray paths coincide well with analytical trajectories in transversely homogeneous elliptically anisotropic media. Seismic ray tracing results in transversely inhomogeneous elliptically anisotropic media demonstrate that our method is effective for further first-arrival tomography in elliptically anisotropic media with an irregular surface.

The methods to solve the propagation of seismic waves mainly include wave-field simulation (Carcione et al. 1992;Lan and Zhang 2011;Liu et al. 2014a, b), seismic ray tracing (Cerveny 2001;Cardarelli and Cerreto 2002;Xu et al. 2006Xu et al. , 2008Xu et al. , 2010Xu et al. , 2014 and traveltime calculation using eikonal equation (Faria and Stoffa 1994;Lan and Zhang 2013a, b). Compared with the method of traveltime calculation using eikonal equation, ray tracing methods not only can obtain the traveltime of the seismic wave field, but also can get the ray trajectories in the ground (Shearer and Chapman 1989;Cerveny 2001). Traditional ray tracing methods include shooting methods (Langan et al. 1985;Virieux and Farra 1991;Xu et al. 2004Xu et al. , 2007Xu et al. , 2008 and bending methods (Julian and Gubbins 1977;Thurber and Ellsworth 1980;Aki and Richards 1980). Based on the bending method, Um and Thurber (1987) developed the pseudo-bending method, which was used to solve the twopoint ray tracing in continuous media (Um and Thurber 1987;Pereyra 1992;Xu et al. 2006), while is not suitable for the presence of strong velocity discontinuities. To solve this problem, Zhao et al. (1992Zhao et al. ( , 1994 and Zhao and Lei (2004) developed a method using the Snell's law to correct the path points on discontinuous interfaces. Xu et al. (2010Xu et al. ( , 2014) developed a ray tracing perturbation scheme of combination of pseudo-bending methods and segmentally iterative methods. In recent decades, the ray tracing methods developed include wavefront reconstruction methods (Vinje et al. 1993(Vinje et al. , 1996a(Vinje et al. , b, 1999, slowness matching methods (Symes 1996;Symes and Qian 2003), Huygens wavefront tracing methods (Sava and Fomel 2001), shortest path methods (Moser 1991;Fischer and Lees 1993;Zhou and Greenhalgh 2005) and simulated annealing methods (Bona et al. 2009).
Rough topography is very common, and we have to deal with it during the acquisition, processing and interpretation of seismic data (Neuberg and Pointer 2000;Bean et al. 2008;Lan and Zhang 2011;Lan et al. 2012;Bevc 2012). How to accurately describe the geological model with undulating interfaces is very important for solving the seismic wave propagation. The methods always used include the approximation method with ladderlike grids , the model expansion method with regular grids (Hole 1992;Ma and Zhang 2014), the nonuniform grid spacing method (Sun et al. 2012a, b), the hybrid grid spacing method (Sun et al. 2009(Sun et al. , 2012aBai et al. 2010Bai et al. , 2013Li et al. 2013), the curved grid spacing method (Thompson et al. 1985;Ruud 1994, 1998;Dong 2005;Wang and Liu 2006;Hestholm et al. 2006;Lan et al. 2012;Zhang 2011, 2013a, b;Ma and Zhang 2014) and the triangulated meshing method (Fomel 1997;Sethian 1999;Xu et al. 2006Xu et al. , 2010Kao et al. 2008;Yu et al. 2010;Bai et al. 2012). The triangular grids are simple and flexible in modeling and can describe any undulating terrain with high degree of approximation.
The elliptically anisotropic media is the anisotropic media with only four independent elastic parameters. The propagation of seismic waves can be described only by the vertical phase velocity and the horizontal phase velocity. In this paper, we construct the medium models with irregular surface parameterized by the flexible triangular grids and then develop the ray tracing method of turning waves in the elliptically anisotropic media with irregular surface, which can be regarded as the forward modeling of the further first-arrival tomography in the media.

Elliptically anisotropic media
We often use elastic parameters (or elastic tensor) to describe seismic anisotropic media (Cerveny 2001). The elliptically anisotropic media have only four independent elastic parameters (Schleicher and Aleixo 2007), and its elastic tensor has one more condition than that of the transversely isotropic media with vertical symmetry axis (VTI). If C ik denotes an element of the elastic tensor and A ik denotes an elastic tensor element with density normalization, that is, A ik = C ik /q, the density normalized elastic tensor matrix A in the elliptically anisotropic media is expressed as where the elastic tensor satisfies the following conditions (Schleicher and Aleixo 2007): Condition (2) is represented by the Thomsen parameters (Thomsen 1986) as e = d. Thomsen parameters were proposed and deduced to characterize the elastic properties of transversely isotropic media and to analyze the propagation characteristics of waves in weakly anisotropic media. e is used to measure the anisotropic intensity of P wave, which reflects the difference between the horizontal phase velocity and vertical phase velocity of P wave, whereas d is a transitional parameter associated with vertical phase velocity and horizontal phase velocity of P wave, reflecting the magnitude of anisotropy of phase velocity near the vertical direction of P wave (Niu et al. 2002). Rasolofosaon (1998) suggested that the anisotropy induced by stress may be elliptical anisotropy. Therefore, the elliptically anisotropic media is common in the ground which is the anisotropic media with the least parameters observed so far (Gurvich 1940;Kleyn 1956;Levin 1978;Daley and Hron 1979;Rogister and Slawinski 2005;Grechka 2009). The expression of the wavefront velocity (phase velocity) in the elliptically anisotropic media (Rogister and Slawinski 2005) where v H and v V denote the horizontal and vertical component of the wavefront velocity, respectively. h is the wavefront angle (phase angle) which is the angle between the normal direction of the seismic wavefront and the symmetry axis of the medium. It can be seen from the expression of the phase velocity (3) that the propagation of seismic waves in the elliptically anisotropic media is described only by the vertical and horizontal phase velocities. The ratio of these two velocities is called elliptical coefficient (Schleicher and Aleixo 2007;Grechka 2009). Compared with the isotropic media, the elliptically anisotropic media has only one more parameter to describe the propagation of seismic waves.

Model parameterization
Model parameterization is the first step of ray tracing. We construct models with irregular surface and then carry out the parameterization in two steps. Firstly, we mesh the model using the triangular grids ( Fig. 1) and define the vertical phase velocity v V and the horizontal isotropic plane phase velocity v H on the grid nodes, respectively. Here, we use the finite element pre-and post-processing software GiD 9.0.2 (Otin et al. 2005) to mesh the medium models with the unstructured triangles.
Then, we use the linear interpolation method to obtain the vertical and horizontal velocities everywhere inside all triangulated grids, that is, where v Vi and v Hi are the vertical and horizontal phase velocities of the three vertices in a triangular grid, respectively; u i is the corresponding area coordinate in a triangle (Xu et al. 2004(Xu et al. , 2005(Xu et al. , 2006; and v V and v H are vertical and horizontal phase velocities of an given point, respectively. After these two steps, we can obtain the velocity distribution within the whole model space (Fig. 2).
4 Ray tracing method

Ray equation
The eikonal equation is a nonlinear partial differential equation describing the propagation characteristics of wavefront in nonuniform continuous anisotropic media (Kravtsov and Orlov 1990;Cerveny 2001;Slawinski 2003). We introduce the ray equation in the elliptically anisotropic media from the eikonal equation.
The general form of the eikonal equation is expressed as (Slawinski 2003): where x ¼ ðx; zÞ is the position vector and p is the wavefront slowness which is normal to the wavefront. Phase velocity expression (3) is substituted into Eq. (5): where p x and p z are the horizontal and vertical component of the wavefront slowness vector p, respectively. h is the phase angle which is positive in the counterclockwise direction from the axis of symmetry. Equation (7) is the eikonal equation of elliptical anisotropy. The eikonal equation can be considered as a nonlinear partial differential equation, and then, ray equation is obtained by solving the eikonal equation using the method of characteristics (Daley and Hron 1979;Cerveny 2001). Considering the general form of the Hamiltonian function: where H = H(x i , p i ), i = 1, 2, and x 1 = x, x 2 = z, p 1 = p x , p 2 = p z . The Hamiltonian function containing the eikonal equation information in the elliptically anisotropic media can be obtained from Eq. (7) (Slawinski 2003).
The ray equations are obtained by substituting formula (9) into formula (8): where t is the traveltime along ray trajectories. The initial conditions of the equations above are: Based on the ray Eqs. (10), the Runge-Kutta method (Press et al. 2007) is used to calculate the numerical solutions of the shooting ray paths.

Linearly modified shooting method
The linearly modified shooting method (Xu et al. 2007) is used to trace the ray paths. The specific process is introduced as follows.
Firstly, given several departure group angles (also called ray angles, that is, the angles between the vectors of the group velocities and the symmetries of the media which are positive in the counterclockwise direction from the axis of symmetry) and the known velocity distribution, we can obtain the departure phase angles and the shooting ray trajectories using the Runge-Kutta method. The relationship of phase angle and group angle is expressed by Eq. (12) (Thomsen 1986). To substitute formula (6) into formula (12), we can get the departure phase angles expressed by group angles in Eq. (13): Secondly, we search for two ray paths adjacent to each side of the known receiver, and then, we calculate the departure group angle of a new ray between the two adjacent rays by linear interpretation (Xu et al. 2007). The equation of the linear interpretation is: where x a and x b are the horizontal coordinates of the detection points of the two rays adjacent to each side of the known receiver, respectively, and u a and u b are the corresponding departure group angles of the two rays, respectively. Finally, we calculate the departure phase angle from the new departure group angle, and a new ray path is obtained by shooting using the Runge-Kutta method. If the difference between the emergent position of the new path and the position of the given receiver satisfies the accuracy requirement, tracing is over, otherwise, the previous procedure is repeated to find the final ray path.

Ray tracing in transversely homogeneous elliptically anisotropic media
We construct a model with a range of 10 km in horizontal direction and 5 km in depth direction and with a gently undulating surface (Figs. 3, 5). The transversely homogeneous elliptically anisotropic velocity model is defined by three parameters: a, b and v, where it assumes that the wavefront velocity is only related to the depth and the propagation direction of the vertical plane (Slawinski et al. 2004): where a is the initial vertical velocity, b is the vertical velocity gradient, z is the depth and v is a dimensionless constant which indicates the relationship between the vertical phase velocity and horizontal phase velocity.
In the transversely homogeneous elliptically anisotropic velocity model defined by Eq. (15), ray trajectories have analytical solutions. Thus, it is useful to examine our ray tracing method.
Given that the source is located at the point (x 0 , z 0 ), and the departure group angle is / 0 , the initial phase angle h 0 can be calculated by Eq. (13). The analytical expression of the ray path can be derived as (Rogister and Slawinski 2005): The traveltime expressed by the offset x on the surface is: We set the values of each parameter as: a = 2.5 km/s, b = 0.7 s -1 , v = 0.06 and the source is set at the point with position of (0, 0). Then, three ray paths (blue curves in Fig. 3) are shooting by the Runge-Kutta numerical method with the time step Dt being 0.01 s. The corresponding theoretical ray paths are obtained by Eq. (16) (red curves in Fig. 3). Meanwhile, we also compare the traveltime of the numerical results with the traveltime of the corresponding theoretical ray paths calculated by Eq. (17) (shown in Fig. 4). We carry out another test with the time step being  0.001 s and then get the numerical (blue curves) and theoretical (red curves) ray paths shown in Fig. 5, respectively. The comparison results of the traveltime between the numerical and theoretical results are shown in Fig. 6. It can be seen that the precision of the numerical results is higher with the time step reducing. When the time step is small enough, the numerical results gradually approach to the theoretical ones, but the shooting process will become more time-consuming. As a result, an appropriate time step should be chosen for the balance between the accuracy and efficiency.

Ray tracing in transversely inhomogeneous elliptically anisotropic media
We construct a medium model with a range of 10 km in horizontal direction and 5 km in depth direction. The surface undulation is controlled by formula (18).
We place the source on the surface with a horizontal position of 0.5 km. The receivers are also placed on the surface. We set 18 receivers whose horizontal positions start from 1.5 to 10 km with the spacing of 0.5 km.
Firstly, we mesh the model by the triangular grids and then define the vertical phase velocity v V and the horizontal phase velocity v H in all grid nodes. The velocity model we use is: where the values of all parameters are a = 2.5 km/s, b = 0.6 s -1 , v = 0.1, l = 10 km, k = 0.12, respectively. It can be seen from Eq. (19) that the vertical phase velocity field is transversely homogeneous, but the horizontal phase velocity field is transversely inhomogeneous. The velocity Fig. 5 Ray paths in the horizontal homogeneous elliptically anisotropic medium with the time step Dt being 0.001 s. a Vertical phase velocity field; b horizontal phase velocity field; blue curves denote the numerical solutions, while red curves denote the analytical solutions. Note that red curves almost coincide with blue ones Fig. 6 Comparison of the traveltime between the numerical and theoretical results of the test shown as Fig. 5. Blue crosses denote the traveltimes of the numerical results in Fig. 5, while red crosses denote the corresponding theoretical traveltimes. Note that red crosses almost coincide with blue ones due to small traveltime difference distributions in both directions everywhere inside the triangulated grids are obtained by linear interpolation. We try-shoot 16 rays (Fig. 7) by the Runge-Kutta method with the departure group angle of the first ray being p/4 and the departure angle of each subsequent ray being increased by p/100. Then, the final ray tracing results (Fig. 8) are obtained by the shooting method as described above. Here, the precision of shooting ray tracing is defined as the distance between the position of emergent points of ray path and the receiver locations, which is typically equal to 15 m in this paper. Figure 7 shows that triangular grids are flexible with high degree of approximation to describe a strong undulating terrain, and ray tracing of turning waves in the transversely inhomogeneous elliptically anisotropic media is achieved. Note that some receivers are located in the shadow zone and cannot be traced due to sharp undulations of the terrain.

Conclusions
We employ vertical and horizontal phase velocity for seismic wave propagation in elliptically anisotropic media. Models with irregular surface are constructed, and model parameterization is realized by applying flexible triangular grids with high degree of approximation to the irregular surface. The vertical and horizontal phase velocities are defined in the triangular grids, respectively, and then, the phase velocity distribution everywhere in the model is obtained by linear interpolation. A shooting method for turning wave in the elliptically anisotropic media with irregular surface has developed. The partial differential equation of seismic ray in elliptically anisotropic media is solved by Runge-Kutta method, and the linearly modified method is used for adjusting departure phase angles in the shooting scheme. Numerical tests show that seismic tracing rays coincide well with analytical trajectories in transversely homogeneous elliptically anisotropic media. Seismic ray tracing results in transversely inhomogeneous elliptically anisotropic media demonstrate that our method is effective for further first-arrival tomography in elliptically anisotropic media with an irregular surface.