**1. Introduction**

104 Earthquake Engineering

Iran, 9(3).

Buffalo, NY.

and Stokoe, K.H. 2001. Liquefaction resistance of soils. Summary report from the 1996 NCEER and 1998 NCEER/NSF workshops on evaluation of liquefaction resistance of

[76] Boulanger, R. and Idriss, I.M. 2004. State normalization of penetration resistance and the effect of overburden stress on liquefaction resistance. Proc. 11th International Conf. on Soil Dynamics and Earthquake Engineering and 3rd International Conference on

[77] Fatemi-Agdha, S.M., Teshnehlab, M., Suzuki, A., Akiyoshi, T., and Kitazono, Y. 1998. Liquefaction potential assesment using multilayer artificial neural network. J. Sci. I.R.

[78] Juang, C. H., Chen, C. J., and Tien, Y. M. 1999. Appraising cone penetration test based liquefaction resistance evaluation methods: Artificial neural networks approach.

[79] Juang, C. H., Yuan, H. M., Lee, D. H., and Lin 2003, P. S., "Simplified cone penetration test-based method for evaluating liquefaction resistance of soils," Journal of

[83] Bartlett, S. F. and Youd, T. L., 1992, "Empirical Analysis of Horizontal Ground Displacement Generated by Liquefaction Induced Lateral Spreads", Tech. Rept. NCEER 92 - 0021, National Center for Earthquake Engineering Research, SUNY - Buffalo,

[84] Youd T.L., Hansen C.M. and Bartlett S.F. (2002) 'Revised Multilinear Regression Equations for Prediction of Lateral Spread Displacement' Journal of Geotechnical and

Geotechnical and Geoenvironmental Engineering, Vol. 129, No. 1, pp. 66-80. [80] Baziar, M. H. and Nilipour, N. 2003. Evaluation of liquefaction potential using neuralnetworks and CPT results. Soil Dynamics and Earthquake Engineering, 23(7) 631-636. [81] Chern, S. and Lee, C. 2009. CPT-based simplified liquefaction assessment by using fuzzy-neural network. Journal of Marine Science and Technology, 17(4) 326-331 [82] García S.R. and Romo M.P. 2007. GENES: Genetic Generator of Signals, a Synthetic

Earthquake Geotechnical Engineering, Univ. of California, Berkeley, CA.

soils. J. Geotech. Geoenviron. Eng., 127(10), 817–833.

Canadian Geotechnical Journal, 36(3) 443-454.

Accelerograms Application. , Proc. of the SEE5, SM-80.

Geoenvironmental Engineering, Vol. 128, No. 12, pp. 1007-1017.

Rough topography is very common and we have to deal with it during the acquisition, processing and interpretation of seismic data. For example, in the context of the deep seismic soundings to explore the crustal structure, seismic experiments are usually carried out across: (a) orogenic belts for understanding the mechanisms; (b) basins to understand the formation mechanisms; (c) transition zones for the study of its interaction (Al-Shukri et al., 1995; Ashford et al., 1997; Boore, 1972; Jih et al., 1988; Levander, 1990; Robertsson, 1996; Zhang et al., 2010). In oil/gas seismic exploration, seismologists also have a similar problem with the undulating topography along the survey line.

In the last two decades, several approaches have been proposed to simulate wave propagation in heterogeneous medium with irregular topography. These schemes include finite element method (Rial et al., 1992; Toshinawa and Ohmachi, 1992), spectral element method (Komatitsch and Tromp, 1999, 2002), pseudo-spectral method (Nielsen et al., 1994; Tessmer et al., 1992; Tessmer and Kosloff, 1994), boundary element method (Bouchon et al., 1989; Campillo and Bouchon, 1985; Sánchez-Sesma and Campillo, 1993; Sánchez-Sesma et al., 2006), finite difference method (Frankel and Vidale, 1992; Gao and Zhang, 2006; Hestholm and Ruud, 1994, 1998; Jih et al., 1988; Lombard et al., 2008; Robertsson, 1996; Zhang and Chen, 2006), and also a hybrid approach which combines the staggered-grid finite difference scheme with the finite element method (Galis et al., 2008; Moczo et al., 1997). Both the spectral element and the finite element methods satisfy boundary conditions on the free surface naturally. 3D surface and interface topographies can be modeled using curved piecewise elements. However, the classical finite element method suffers from a high computational cost, and, on the other hand, a smaller spectral element than the one required by numerical dispersion is required to describe a highly curved topography, as

© 2012 Lan and Zhongjie Zhang, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Lan and Zhongjie Zhang, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

demonstrated in seismic modeling of a hemispherical crater (Komatitsc and Tromp, 1999). The pseudo-spectral method is limited to a free surface with smoothly varying topography and leads to inaccuracies for models with strong heterogeneity or sharp boundaries (Tessmer et al., 1992). The boundary integral equation and boundary element methods are not suitable for near-surface regions with large velocity contrasts (Bouchon et al., 1995). The finite difference method is one of the most popular numerical methods used in computational seismology. In comparison to other methods, the finite difference method is simpler and more flexible, although it has some difficulty in dealing with surface topography. The situation has improved recently. For rectangular domains, a stable and explicit discretization of the free surface boundary conditions has been presented by Nilsson et al. (2007). By using boundary-modified difference operators, Nilsson et al. (2007) introduce a discretization of the mixed derivatives in the governing equations; they also show that the method is second order accurate for problems with smoothly varying material properties and stable under standard Courant-Friedrichs-Lewy (CFL) constraints, for arbitrarily varying material properties. We have investigated 6 free-surface boundary condition approximate schemes in seismic wavefield modelling and evaluated their stability and applicability by comparing with corresponding analytical solutions, the results reveal that Nilsson et al.'s method is more effective than others (Lan & Zhang, 2011a). Recently, Appelo and Petersson (2009) have generalized the results of Nilsson et al. (2007) to curvilinear coordinate systems, allowing for simulations on non-rectangular domains. They construct a stable discretization of the free surface boundary conditions on curvilinear grids, and they prove that the strengths of the proposed method are its ease of implementation, efficiency (relative to low-order unstructured grid methods), geometric flexibility, and, most importantly, the "bullet-proof" stability (Appelo and Petersson, 2009), even though they deal with 2D isotropic medium.

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 107

Thompson et al., 1985), and it was early used by Fornberg (1988) in seismic wave simulation with the pseudo-spectral method. A grid of this type is achieved by carrying out a transformation between the (curvilinear) computational space and the (Cartesian) physical space as illustrated in Figure 1. By means of this transformation, the curvilinear coordinates q, r and s are mapped into Cartesian coordinates within the physical space, where both systems have positive direction downward for the vertical coordinate. A boundary in the physical space presents a constant value of one of the curvilinear coordinates-be it a curve in

Boundary conforming grids may be of two fundamentally different types: structured and unstructured (or irregular) grids. A structured type grid (Figure 1) is characterized by having a fixed number of elements along each of the coordinate directions. The general element is a hexahedron in 3D, just as in the left panel of Figure 1. Neighboring elements in the physical space are also adjacent in the computational space, which is one of the great advantages of this type of grid. This property makes it relatively simple to implement in a computer. Structured grids are mainly used in finite difference and finite volume solvers. Here, we focus on structured boundary conforming grids. Several methods may be used to generate these grids, namely: Partial Differential Equation (PDE) methods, algebraic methods, co-normal mapping methods and variational methods. Here we use PDE methods

**Figure 1.** Mapping between computational and physical space in three dimensions (after Hvid, S.,

point can be determined from the curvilinear coordinates through the equations:

the curvilinear coordinate system (q, r, s) following the chain rule:

After generating the boundary-conforming grid, the Cartesian coordinates of every grid

then, we can express the spatial derivatives in the Cartesian coordinate system (x, y, z) from

x xq xr xs =q +r +s y y <sup>q</sup> <sup>y</sup> <sup>r</sup> <sup>y</sup> <sup>s</sup> =q +r +s z zq zr zs =q +r +s (2)

q qx q <sup>y</sup> q z =x +y +z r rx r <sup>y</sup> r z =x +y +z s sx s <sup>y</sup> s z =x + y + z (3)

x = x(q,r,s) y = y(q,r,s) z = z(q,r,s) (1)

two dimensions or a surface in three dimensions.

(see Hvid, 1994; and Thompson et al., 1985 for details).

1994).

and similarly in other cases

Nevertheless, the earth is often seismically anisotropic resulting from fractured rocks, fluidfilled cracks (Crampin, 1981; Hudson, 1981; Liu et al., 1993; Schoenberg and Muir, 1989; Zhang et al., 1999), thin isotropic layering (Backus, 1962; Helbig, 1984), lack of homogeneity (Grechka and McMechan, 2005), or even preferential orientation of olivine (Dziewonski and Anderson, 1981; Forsyth, 1975). Here, we give an introduction of the method in 3D case with the purpose of simulating seismic wave propagation in 3D heterogeneous anisotropic medium with non-flat surface topography. The chapter is organized as follows: firstly, we give a brief introduction on the boundary-conforming grid and the transformation between curvilinear coordinates and Cartesian coordinates; then we write the wave equations and free boundary conditions in these two coordinate systems; after that we introduce a numerical method to discretize both the wave equations and the free surface boundary conditions. Finally, several numerical examples are presented to demonstrate the accuracy and efficiency of the method.

#### **2. Transformation between curvilinear and Cartesian coordinates**

As to the topographic surface, the discrete grid must conform to the free surface to suppress artificial scattered waves. Such grid is named boundary-conforming grid (Hvid, 1994; Thompson et al., 1985), and it was early used by Fornberg (1988) in seismic wave simulation with the pseudo-spectral method. A grid of this type is achieved by carrying out a transformation between the (curvilinear) computational space and the (Cartesian) physical space as illustrated in Figure 1. By means of this transformation, the curvilinear coordinates q, r and s are mapped into Cartesian coordinates within the physical space, where both systems have positive direction downward for the vertical coordinate. A boundary in the physical space presents a constant value of one of the curvilinear coordinates-be it a curve in two dimensions or a surface in three dimensions.

Boundary conforming grids may be of two fundamentally different types: structured and unstructured (or irregular) grids. A structured type grid (Figure 1) is characterized by having a fixed number of elements along each of the coordinate directions. The general element is a hexahedron in 3D, just as in the left panel of Figure 1. Neighboring elements in the physical space are also adjacent in the computational space, which is one of the great advantages of this type of grid. This property makes it relatively simple to implement in a computer. Structured grids are mainly used in finite difference and finite volume solvers. Here, we focus on structured boundary conforming grids. Several methods may be used to generate these grids, namely: Partial Differential Equation (PDE) methods, algebraic methods, co-normal mapping methods and variational methods. Here we use PDE methods (see Hvid, 1994; and Thompson et al., 1985 for details).

**Figure 1.** Mapping between computational and physical space in three dimensions (after Hvid, S., 1994).

After generating the boundary-conforming grid, the Cartesian coordinates of every grid point can be determined from the curvilinear coordinates through the equations:

$$\mathbf{x} = \mathbf{x}(\mathbf{q}\_l \mathbf{r}, \mathbf{s}) \qquad \mathbf{y} = \mathbf{y}(\mathbf{q}\_l \mathbf{r}, \mathbf{s}) \qquad \mathbf{z} = \mathbf{z}(\mathbf{q}\_l \mathbf{r}, \mathbf{s}) \tag{1}$$

then, we can express the spatial derivatives in the Cartesian coordinate system (x, y, z) from the curvilinear coordinate system (q, r, s) following the chain rule:

$$
\partial\_{\mathbf{x}} \cdot \mathbf{=} \mathbf{q}\_{\mathbf{x}} \partial\_{\mathbf{q}} + \mathbf{r}\_{\mathbf{x}} \partial\_{\mathbf{r}} + \mathbf{s}\_{\mathbf{x}} \partial\_{\mathbf{s}} \qquad \partial\_{\mathbf{y}} = \mathbf{q}\_{\mathbf{y}} \partial\_{\mathbf{q}} + \mathbf{r}\_{\mathbf{y}} \partial\_{\mathbf{r}} + \mathbf{s}\_{\mathbf{y}} \partial\_{\mathbf{s}} \qquad \partial\_{\mathbf{z}} = \mathbf{q}\_{\mathbf{z}} \partial\_{\mathbf{q}} + \mathbf{r}\_{\mathbf{z}} \partial\_{\mathbf{r}} + \mathbf{s}\_{\mathbf{z}} \partial\_{\mathbf{s}} \tag{2}
$$

and similarly in other cases

106 Earthquake Engineering

deal with 2D isotropic medium.

and efficiency of the method.

demonstrated in seismic modeling of a hemispherical crater (Komatitsc and Tromp, 1999). The pseudo-spectral method is limited to a free surface with smoothly varying topography and leads to inaccuracies for models with strong heterogeneity or sharp boundaries (Tessmer et al., 1992). The boundary integral equation and boundary element methods are not suitable for near-surface regions with large velocity contrasts (Bouchon et al., 1995). The finite difference method is one of the most popular numerical methods used in computational seismology. In comparison to other methods, the finite difference method is simpler and more flexible, although it has some difficulty in dealing with surface topography. The situation has improved recently. For rectangular domains, a stable and explicit discretization of the free surface boundary conditions has been presented by Nilsson et al. (2007). By using boundary-modified difference operators, Nilsson et al. (2007) introduce a discretization of the mixed derivatives in the governing equations; they also show that the method is second order accurate for problems with smoothly varying material properties and stable under standard Courant-Friedrichs-Lewy (CFL) constraints, for arbitrarily varying material properties. We have investigated 6 free-surface boundary condition approximate schemes in seismic wavefield modelling and evaluated their stability and applicability by comparing with corresponding analytical solutions, the results reveal that Nilsson et al.'s method is more effective than others (Lan & Zhang, 2011a). Recently, Appelo and Petersson (2009) have generalized the results of Nilsson et al. (2007) to curvilinear coordinate systems, allowing for simulations on non-rectangular domains. They construct a stable discretization of the free surface boundary conditions on curvilinear grids, and they prove that the strengths of the proposed method are its ease of implementation, efficiency (relative to low-order unstructured grid methods), geometric flexibility, and, most importantly, the "bullet-proof" stability (Appelo and Petersson, 2009), even though they

Nevertheless, the earth is often seismically anisotropic resulting from fractured rocks, fluidfilled cracks (Crampin, 1981; Hudson, 1981; Liu et al., 1993; Schoenberg and Muir, 1989; Zhang et al., 1999), thin isotropic layering (Backus, 1962; Helbig, 1984), lack of homogeneity (Grechka and McMechan, 2005), or even preferential orientation of olivine (Dziewonski and Anderson, 1981; Forsyth, 1975). Here, we give an introduction of the method in 3D case with the purpose of simulating seismic wave propagation in 3D heterogeneous anisotropic medium with non-flat surface topography. The chapter is organized as follows: firstly, we give a brief introduction on the boundary-conforming grid and the transformation between curvilinear coordinates and Cartesian coordinates; then we write the wave equations and free boundary conditions in these two coordinate systems; after that we introduce a numerical method to discretize both the wave equations and the free surface boundary conditions. Finally, several numerical examples are presented to demonstrate the accuracy

**2. Transformation between curvilinear and Cartesian coordinates** 

As to the topographic surface, the discrete grid must conform to the free surface to suppress artificial scattered waves. Such grid is named boundary-conforming grid (Hvid, 1994;

$$\boldsymbol{\partial}\_{\mathbf{q}} = \mathbf{x}\_{\mathbf{q}} \boldsymbol{\partial}\_{\mathbf{x}} + \mathbf{y}\_{\mathbf{q}} \boldsymbol{\partial}\_{\mathbf{y}} + \mathbf{z}\_{\mathbf{q}} \boldsymbol{\partial}\_{\mathbf{z}} \qquad \boldsymbol{\partial}\_{\mathbf{r}} = \mathbf{x}\_{\mathbf{r}} \boldsymbol{\partial}\_{\mathbf{x}} + \mathbf{y}\_{\mathbf{r}} \boldsymbol{\partial}\_{\mathbf{y}} + \mathbf{z}\_{\mathbf{r}} \boldsymbol{\partial}\_{\mathbf{z}} \qquad \boldsymbol{\partial}\_{\mathbf{s}} = \mathbf{x}\_{\mathbf{s}} \boldsymbol{\partial}\_{\mathbf{x}} + \mathbf{y}\_{\mathbf{s}} \boldsymbol{\partial}\_{\mathbf{y}} + \mathbf{z}\_{\mathbf{s}} \boldsymbol{\partial}\_{\mathbf{z}} \tag{3}$$

where *<sup>x</sup> q* denotes q(x,y,z) x and the similar in other cases. These derivatives are called metric derivatives or simply the metric. We can also find the metric derivatives

x rs r s x sq s q x qr q r <sup>1</sup> q= (y z -z <sup>y</sup> ) <sup>J</sup> <sup>1</sup> r= (y z -z <sup>y</sup> ) <sup>J</sup> <sup>1</sup> s= (y z -z <sup>y</sup> ) <sup>J</sup> y rs rs y sq sq y qr qr <sup>1</sup> q = (z x - x z ) <sup>J</sup> <sup>1</sup> r = (z x - x z ) <sup>J</sup> <sup>1</sup> s = (z x - x z ) <sup>J</sup> z rs rs z sq sq z qr qr <sup>1</sup> q = (x <sup>y</sup> - <sup>y</sup> x ) <sup>J</sup> <sup>1</sup> r = (x <sup>y</sup> - <sup>y</sup> x ) <sup>J</sup> <sup>1</sup> s = (x <sup>y</sup> - <sup>y</sup> x ) <sup>J</sup> (4)

Three-Dimensional Wavefield Simulation in

+ r + s )w]

(7)

(8)

(9)

q zr zs

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 109

s 12 y q y r y s 13 z q z r z s

)u + c (q + r + s )v + c (q + r + s )w]

<sup>x</sup> {

2 x 11 x q r x s 12 y q y r y s 13 z q z r z s

Jρ Jq [c (q + r + s )u + c (q + r + s )v + c (q + r + s )w]

 

 

 

x 66 x q x r x s 66 y q y r y s

x 44 z q z r z s 44 x q x r x s

(q + r + s )v + c (q + r + s )w]

+Js [c

+ {Js [c (q + r + s )u + c (q + r + s )w]

4 z q z r z s 44 y q y r y s

y 44 z q z r z s 44 y q y r y s

+Jr [c (q + r + s )v + c (q + r + s )w]

z 44 z q z r z s 44 y q y r y s

2 x 44 z q z r z s 44 x q x r x s

<sup>w</sup> Jρ = {Jq [c (q + r + s )u + c (q + r + s )w]

y 44 z q z r z s 44 y q y r y s

+Jq [c (q + r + s )v + c (q + r + s )w]

+Js [c (q + r + s )v + c (q + r + s )w]}

+ {Js [c (q + r + s )v + c (q + r + s )u]

z 44 z q z r z s 44 y q y r y s

+Jr [c (q + r + s )v + c (q + r + s )w]}

{Js [c (q + r + s )u + c (q + r + s )v + c (q

x 11 x q x r x s 12 y q y r y s 13 z

)v + c (q + r + s )u]

y 11 y q y r y s 12 x q x r x s 13 z q z r z s

y 11 y q y r y s 12 x q x r x s 13 z q z r z s

+Js [c

1 y q y r y s 12 x q x r x s 13 z q z r z s

(q + r + s )v + c (q + r + s )u + c (q + r + s )w]

)u + c (q + r + s )w]

z 33 z q z r z s 13 x q x r x s 13 y q y r y s

z 33 z q z r z s 13 x q x r x s 13 y q y r y s

z 33 z q z r z s 13 x q x r x s 13 y q y r y s

+Js [c (q + r + s )w + c (q + r + s )u + c (q + r + s )v]}

+Jr [c (q + r + s )w + c (q + r + s )u + c (q + r + s )v]}

s 44 x q x r x s

+Jq [c (q + r + s )w + c (q + r + s )u + c (q + r + s )v]}

+Jr [c (q + r + s )v + c (q + r + s )u + c (q + r + s )w]

s 66 y q y r y s

+Jq [c (q + r + s )v + c (q + r + s )u + c (q + r + s )w]

y 66 x q x r x s 66 y q y r y s z 44 z q z r z s 44 x q x r x s

+Jq [c (q + r + s )v + c (q + r + s )u] +Jq [c (q + r + s )u + c (q + r + s )w]}

y 66 q x r x s 66 y q y r y s z 44 z q z r z s 44 x q x r x s

y 66 x q x r x s 66 y q y r y s z 44 z q z r z s 44 x q x r x s

<sup>v</sup> Jρ = {Jq [c (q + r + s )v + c (q + r + s )u]

2 x 66 x q x r x s 66 y q y r y s

z 44 z q z r z s 44 y q y r y s

+Jq [c (q + r + s )v + c (q + r + s )w]}

+Js [c (q + r + s )v + c (q + r + s )u] +Js [c (q + r + s )u + c (q + r + s )w]}

+Jr [c (q + r + s )v + c (q + r + s )u] +Jr [c (q + r + s )u + c (q + r + s )w]}

2

*u*

t

t q

2

t q

2

s

y 4

r

s

y 1

r

*s*

*r*

*q*

*x*

+ {Jr [c (q + r + s *<sup>x</sup>*

+ {Jr [c (q + r + s

x 44 z q z r z

+ {Jr [c (q + r + s

x 66 x q x r x

x 11 q x r x

where J is the Jacobian of the transformation that is written as

$$\mathbf{J} = \mathbf{x\_q}\mathbf{y\_r}\mathbf{z\_s} - \mathbf{x\_q}\mathbf{y\_s}\mathbf{z\_r} - \mathbf{x\_r}\mathbf{y\_q}\mathbf{z\_s} + \mathbf{x\_r}\mathbf{y\_s}\mathbf{z\_q} + \mathbf{x\_s}\mathbf{y\_q}\mathbf{z\_r} - \mathbf{x\_s}\mathbf{y\_r}\mathbf{z\_q}$$

and whose detailed form can be found in Appendix A in Lan & Zhang (2011b).

It is worthful to note that even if the mapping equations (1) are given by an analytic function, the derivatives should still be calculated numerically to avoid spurious source terms due to the coefficients of the derivatives when the conservation form of the momentum equations are used (Thompson et al., 1985). In all examples presented in this paper the metric derivatives are computed numerically using second-order accurate finite difference approximations.

#### **3. Elastic wave equations in Cartesian and curvilinear coordinate systems**

In the following we consider a well-studied type of anisotropy in seismology, namely, a transversely isotropic medium. In the absence of external force, the elastic wave equations in the Cartesian coordinates are given by:

(b) 2 2 11 12 13 66 66 44 44 2 2 66 66 11 12 13 44 44 2 2 44 44 44 u uvw uv uw ρ = (c + c + c ) + (c + c ) + (c + c ) (a) t x x y z y y xz z x v vu vuw vw ρ = (c + c ) + (c + c + c ) + (c + c ) t x x yy y x z z z y w wu w ρ = (c + c ) + (c t x x zy (c) <sup>44</sup> 33 13 13 v wuv + c ) + (c + c + c ) y zz z x y (5)

where *ij c* (x, y, z) are elastic parameters and 66 11 12 c = 0.5(c - c ) ; u, v and w are the displacements in x-, y- and z-directions, respectively; *ρ* (x, y, z) is density. Equations (5a)- (5c) are complemented by the initial data:

$$\begin{aligned} \mathbf{u}(\mathbf{x}, \mathbf{y}, \mathbf{z}, 0) &= \mathbf{u}\_0(\mathbf{x}, \mathbf{y}, \mathbf{z}) & \mathbf{v}(\mathbf{x}, \mathbf{y}, \mathbf{z}, 0) &= \mathbf{v}\_0(\mathbf{x}, \mathbf{y}, \mathbf{z}) & \mathbf{w}(\mathbf{x}, \mathbf{y}, \mathbf{z}, 0) &= \mathbf{w}\_0(\mathbf{x}, \mathbf{y}, \mathbf{z})\\ \frac{\partial \mathbf{u}(\mathbf{x}, \mathbf{y}, \mathbf{z}, 0)}{\partial \mathbf{t}} &= \mathbf{u}\_1(\mathbf{x}, \mathbf{y}, \mathbf{z}) & \frac{\partial \mathbf{v}(\mathbf{x}, \mathbf{y}, \mathbf{z}, 0)}{\partial \mathbf{t}} &= \mathbf{v}\_1(\mathbf{x}, \mathbf{y}, \mathbf{z}) & \frac{\partial \mathbf{w}(\mathbf{x}, \mathbf{y}, \mathbf{z}, 0)}{\partial \mathbf{t}} = \mathbf{w}\_1(\mathbf{x}, \mathbf{y}, \mathbf{z}) \end{aligned} \tag{6}$$

Utilizing relationships (2), the wave equations (5a)-(5c) can be re-written in the curvilinear coordinate system in the following form (see Appendix B in Lan & Zhang, 2011b for details):

2 <sup>x</sup> { 2 x 11 x q r x s 12 y q y r y s 13 z q z r z s y 66 x q x r x s 66 y q y r y s z 44 z q z r z s 44 x q x r x s x 11 q x r x Jρ Jq [c (q + r + s )u + c (q + r + s )v + c (q + r + s )w] t +Jq [c (q + r + s )v + c (q + r + s )u] +Jq [c (q + r + s )u + c (q + r + s )w]} + {Jr [c (q + r + s *<sup>x</sup> u q r* s 12 y q y r y s 13 z q z r z s y 66 q x r x s 66 y q y r y s z 44 z q z r z s 44 x q x r x s x 11 x q x r x s 12 y q y r y s 13 z )u + c (q + r + s )v + c (q + r + s )w] +Jr [c (q + r + s )v + c (q + r + s )u] +Jr [c (q + r + s )u + c (q + r + s )w]} {Js [c (q + r + s )u + c (q + r + s )v + c (q *x s* q zr zs y 66 x q x r x s 66 y q y r y s z 44 z q z r z s 44 x q x r x s + r + s )w] +Js [c (q + r + s )v + c (q + r + s )u] +Js [c (q + r + s )u + c (q + r + s )w]} (7)

where *<sup>x</sup> q* denotes q(x,y,z) x and the similar in other cases. These derivatives are called

y rs rs

<sup>1</sup> q = (z x - x z ) <sup>J</sup>

y sq sq

<sup>1</sup> r = (z x - x z ) <sup>J</sup>

z rs rs

<sup>1</sup> q = (x <sup>y</sup> - <sup>y</sup> x ) <sup>J</sup>

z sq sq

<sup>1</sup> r = (x <sup>y</sup> - <sup>y</sup> x ) <sup>J</sup>

(4)

(b)

(5)

(6)

z qr qr

(c)

 <sup>44</sup> 33 13 13 v wuv + c ) + (c + c + c ) y zz z x y

0

w(x,y,z,0) = w (x,y,z)

w(x, y,z,0) = w (x,y,z) <sup>t</sup>

1

<sup>1</sup> s = (x <sup>y</sup> - <sup>y</sup> x ) <sup>J</sup>

y qr qr

q rs q sr r qs r sq s qr s rq J=x y z -x y z -x y z +x y z +x y z -x y z

It is worthful to note that even if the mapping equations (1) are given by an analytic function, the derivatives should still be calculated numerically to avoid spurious source terms due to the coefficients of the derivatives when the conservation form of the momentum equations are used (Thompson et al., 1985). In all examples presented in this paper the metric derivatives are computed numerically using second-order accurate finite

**3. Elastic wave equations in Cartesian and curvilinear coordinate systems** 

In the following we consider a well-studied type of anisotropy in seismology, namely, a transversely isotropic medium. In the absence of external force, the elastic wave equations in

 

u uvw uv uw ρ = (c + c + c ) + (c + c ) + (c + c ) (a) t x x y z y y xz z x v vu vuw vw ρ = (c + c ) + (c + c + c ) + (c + c ) t x x yy y x z z z y

2 11 12 13 66 66 44 44

2 66 66 11 12 13 44 44

where *ij c* (x, y, z) are elastic parameters and 66 11 12 c = 0.5(c - c ) ; u, v and w are the displacements in x-, y- and z-directions, respectively; *ρ* (x, y, z) is density. Equations (5a)-

v(x,y,z,0) = v (x,y,z)

v(x, y,z,0) = v (x,y,z) <sup>t</sup>

Utilizing relationships (2), the wave equations (5a)-(5c) can be re-written in the curvilinear coordinate system in the following form (see Appendix B in Lan & Zhang, 2011b for details):

0

1

<sup>1</sup> s = (z x - x z ) <sup>J</sup>

metric derivatives or simply the metric. We can also find the metric derivatives

and whose detailed form can be found in Appendix A in Lan & Zhang (2011b).

where J is the Jacobian of the transformation that is written as

x rs r s

<sup>1</sup> q= (y z -z <sup>y</sup> ) <sup>J</sup>

x sq s q

<sup>1</sup> r= (y z -z <sup>y</sup> ) <sup>J</sup>

x qr q r

difference approximations.

2

2

2

the Cartesian coordinates are given by:

 

ρ = (c + c ) + (c t x x zy

(5c) are complemented by the initial data:

u(x,y,z,0) = u (x,y,z)

u(x, y,z,0) = u (x,y,z) <sup>t</sup>

0

1

w wu w

2 44 44 44

<sup>1</sup> s= (y z -z <sup>y</sup> ) <sup>J</sup>

 2 2 x 66 x q x r x s 66 y q y r y s y 11 y q y r y s 12 x q x r x s 13 z q z r z s z 44 z q z r z s 44 y q y r y s x 66 x q x r x <sup>v</sup> Jρ = {Jq [c (q + r + s )v + c (q + r + s )u] t q +Jq [c (q + r + s )v + c (q + r + s )u + c (q + r + s )w] +Jq [c (q + r + s )v + c (q + r + s )w]} + {Jr [c (q + r + s r s 66 y q y r y s y 11 y q y r y s 12 x q x r x s 13 z q z r z s z 44 z q z r z s 44 y q y r y s x 66 x q x r x s 66 y q y r y s y 1 )v + c (q + r + s )u] +Jr [c (q + r + s )v + c (q + r + s )u + c (q + r + s )w] +Jr [c (q + r + s )v + c (q + r + s )w]} + {Js [c (q + r + s )v + c (q + r + s )u] s +Js [c 1 y q y r y s 12 x q x r x s 13 z q z r z s z 44 z q z r z s 44 y q y r y s (q + r + s )v + c (q + r + s )u + c (q + r + s )w] +Js [c (q + r + s )v + c (q + r + s )w]} (8)

 2 2 x 44 z q z r z s 44 x q x r x s y 44 z q z r z s 44 y q y r y s z 33 z q z r z s 13 x q x r x s 13 y q y r y s x 44 z q z r z <sup>w</sup> Jρ = {Jq [c (q + r + s )u + c (q + r + s )w] t q +Jq [c (q + r + s )v + c (q + r + s )w] +Jq [c (q + r + s )w + c (q + r + s )u + c (q + r + s )v]} + {Jr [c (q + r + s r s 44 x q x r x s y 44 z q z r z s 44 y q y r y s z 33 z q z r z s 13 x q x r x s 13 y q y r y s x 44 z q z r z s 44 x q x r x s y 4 )u + c (q + r + s )w] +Jr [c (q + r + s )v + c (q + r + s )w] +Jr [c (q + r + s )w + c (q + r + s )u + c (q + r + s )v]} + {Js [c (q + r + s )u + c (q + r + s )w] s +Js [c 4 z q z r z s 44 y q y r y s z 33 z q z r z s 13 x q x r x s 13 y q y r y s (q + r + s )v + c (q + r + s )w] +Js [c (q + r + s )w + c (q + r + s )u + c (q + r + s )v]} (9)
