Three-Dimensional Wavefield Simulation in Heterogeneous Transversely Isotropic Medium with Irregular Free Surface

1. Introduction

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 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.

Nevertheless, the earth is often seismically anisotropic resulting from fractured rocks, fluid-filled 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:

$$x=x(q,r,s)$$$$y=y(q,r,s)$$

$$z=z(q,r,s)$$E1

(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 }_x=q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s$$$${\partial }_y=q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s$$

$${\partial }_z=q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s$$E2

(2)

and similarly in other cases

$${\partial }_q=x_q{\partial }_x+y_q{\partial }_y+z_q{\partial }_z$$$${\partial }_r=x_r{\partial }_x+y_r{\partial }_y+z_r{\partial }_z$$

$${\partial }_s=x_s{\partial }_x+y_s{\partial }_y+z_s{\partial }_z$$E3

(3)

where $$q_x$$ denotes $$\partial q(x,y,z)/\partial x$$ and the similar in other cases. These derivatives are called metric derivatives or simply the metric. We can also find the metric derivatives

$$\begin{array}{l}{q}_x=\frac{1}{J}(y_r{z}_s-z_r{y}_s)\\ r_x=\frac{1}{J}(y_s{z}_q-z_s{y}_q)\\ s_x=\frac{1}{J}(y_q{z}_r-z_q{y}_r)\end{array}$$$$\begin{array}{l}{q}_y=\frac{1}{J}(z_r{x}_s-x_r{z}_s)\\ r_y=\frac{1}{J}(z_s{x}_q-x_s{z}_q)\\ s_y=\frac{1}{J}(z_q{x}_r-x_q{z}_r)\end{array}$$

$$\begin{array}{l}{q}_z=\frac{1}{J}(x_r{y}_s-y_r{x}_s)\\ r_z=\frac{1}{J}(x_s{y}_q-y_s{x}_q)\\ s_z=\frac{1}{J}(x_q{y}_r-y_q{x}_r)\end{array}$$E4

(4)

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

$$J=x_q{y}_r{z}_s-x_q{y}_s{z}_r-x_r{y}_q{z}_s+x_r{y}_s{z}_q+x_s{y}_q{z}_r-x_s{y}_r{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:

$$\begin{array}{l}\rho \frac{{\partial }^{2}u}{\partial t^2}=\frac{\partial }{\partial x}(c_{11}\frac{\partial u}{\partial x}+c_{12}\frac{\partial v}{\partial y}+c_{13}\frac{\partial w}{\partial z})+\frac{\partial }{\partial y}(c_{66}\frac{\partial u}{\partial y}+c_{66}\frac{\partial v}{\partial x})+\frac{\partial }{\partial z}(c_{44}\frac{\partial u}{\partial z}+c_{44}\frac{\partial w}{\partial x})(a)\\ \rho \frac{{\partial }^{2}v}{\partial t^2}=\frac{\partial }{\partial x}(c_{66}\frac{\partial v}{\partial x}+c_{66}\frac{\partial u}{\partial y})+\frac{\partial }{\partial y}(c_{11}\frac{\partial v}{\partial y}+c_{12}\frac{\partial u}{\partial x}+c_{13}\frac{\partial w}{\partial z})+\frac{\partial }{\partial z}(c_{44}\frac{\partial v}{\partial z}+c_{44}\frac{\partial w}{\partial y})(b)\\ \rho \frac{{\partial }^{2}w}{\partial t^2}=\frac{\partial }{\partial x}(c_{44}\frac{\partial w}{\partial x}+c_{44}\frac{\partial u}{\partial z})+\frac{\partial }{\partial y}(c_{44}\frac{\partial w}{\partial y}+c_{44}\frac{\partial v}{\partial z})+\frac{\partial }{\partial z}(c_{33}\frac{\partial w}{\partial z}+c_{13}\frac{\partial u}{\partial x}+c_{13}\frac{\partial v}{\partial y})(c)\end{array}$$E5

(5)

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

$$\begin{array}{l}u(x,y,z,0)=u_0(x,y,z)\\ \frac{\partial u(x,y,z,0)}{\partial t}=u_1(x,y,z)\end{array}$$$$\begin{array}{l}v(x,y,z,0)=v_0(x,y,z)\\ \frac{\partial v(x,y,z,0)}{\partial t}=v_1(x,y,z)\end{array}$$

$$\begin{array}{l}w(x,y,z,0)=w_0(x,y,z)\\ \frac{\partial w(x,y,z,0)}{\partial t}=w_1(x,y,z)\end{array}$$E6

(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):

$$\begin{array}{c}J\rho \frac{{\partial }^{2}u}{\partial t^2}=\frac{\partial }{\partial q}\{J{q}_x[c_{11}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)u+c_{12}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)v+c_{13}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)w]\\ +J{q}_y[c_{66}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)v+c_{66}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)u]\\ +J{q}_z[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)u+c_{44}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)w]\}\\ +\frac{\partial }{\partial r}\{J{r}_x[c_{11}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)u+c_{12}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)v+c_{13}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)w]\\ +J{r}_y[c_{66}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)v+c_{66}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)u]\\ +J{r}_z[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)u+c_{44}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)w]\}\\ +\frac{\partial }{\partial s}\{J{s}_x[c_{11}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)u+c_{12}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)v+c_{13}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)w]\\ +J{s}_y[c_{66}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)v+c_{66}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)u]\\ +J{s}_z[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)u+c_{44}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)w]\}\end{array}$$E7

(7)

$$\begin{array}{c}J\rho \frac{{\partial }^{2}v}{\partial t^2}=\frac{\partial }{\partial q}\{J{q}_x[c_{66}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)v+c_{66}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)u]\\ +J{q}_y[c_{11}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)v+c_{12}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)u+c_{13}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)w]\\ +J{q}_z[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)v+c_{44}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)w]\}\\ +\frac{\partial }{\partial r}\{J{r}_x[c_{66}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)v+c_{66}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)u]\\ +J{r}_y[c_{11}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)v+c_{12}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)u+c_{13}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)w]\\ +J{r}_z[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)v+c_{44}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)w]\}\\ +\frac{\partial }{\partial s}\{J{s}_x[c_{66}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)v+c_{66}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)u]\\ +J{s}_y[c_{11}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)v+c_{12}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)u+c_{13}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)w]\\ +J{s}_z[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)v+c_{44}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)w]\}\end{array}$$E8

(8)

$$\begin{array}{c}J\rho \frac{{\partial }^{2}w}{\partial t^2}=\frac{\partial }{\partial q}\{J{q}_x[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)u+c_{44}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)w]\\ +J{q}_y[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)v+c_{44}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)w]\\ +J{q}_z[c_{33}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)w+c_{13}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)u+c_{13}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)v]\}\\ +\frac{\partial }{\partial r}\{J{r}_x[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)u+c_{44}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)w]\\ +J{r}_y[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)v+c_{44}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)w]\\ +J{r}_z[c_{33}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)w+c_{13}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)u+c_{13}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)v]\}\\ +\frac{\partial }{\partial s}\{J{s}_x[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)u+c_{44}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)w]\\ +J{s}_y[c_{44}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)v+c_{44}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)w]\\ +J{s}_z[c_{33}(q_z{\partial }_q+r_z{\partial }_r+s_z{\partial }_s)w+c_{13}(q_x{\partial }_q+r_x{\partial }_r+s_x{\partial }_s)u+c_{13}(q_y{\partial }_q+r_y{\partial }_r+s_y{\partial }_s)v]\}\end{array}$$E9

(9)

4. Free boundary conditions in the Cartesian and curvilinear coordinate systems

At the free surface, the boundary conditions in the Cartesian coordinates are given by:

$$\left[\begin{array}{ccc}{c}_{11}\frac{\partial u}{\partial x}+c_{12}\frac{\partial v}{\partial y}+c_{13}\frac{\partial w}{\partial z}& c_{66}\frac{\partial u}{\partial y}+c_{66}\frac{\partial v}{\partial x}& c_{44}\frac{\partial u}{\partial z}+c_{44}\frac{\partial w}{\partial x}\\ c_{66}\frac{\partial v}{\partial x}+c_{66}\frac{\partial u}{\partial y}& c_{11}\frac{\partial v}{\partial y}+c_{12}\frac{\partial u}{\partial x}+c_{13}\frac{\partial w}{\partial z}& c_{44}\frac{\partial v}{\partial z}+c_{44}\frac{\partial w}{\partial y}\\ c_{44}\frac{\partial w}{\partial x}+c_{44}\frac{\partial u}{\partial z}& c_{44}\frac{\partial w}{\partial y}+c_{44}\frac{\partial v}{\partial z}& c_{33}\frac{\partial w}{\partial z}+c_{13}\frac{\partial u}{\partial x}+c_{13}\frac{\partial v}{\partial y}\end{array}\right]\left[\begin{array}{c}{n}_x\\ n_y\\ n_z\end{array}\right]=0$$E10

(10)

Here $${\left[n_x,n_y,n_z\right]}^T$$ is the inward normal of the free surface. Using relationships (2), the above boundary conditions in the curvilinear coordinates can be re-written as

$$\begin{array}{l}{\stackrel{-}{s}}_x\left[c_{11}(q_x{u}_q+r_x{u}_r+s_x{u}_s)+c_{12}(q_y{v}_q+r_y{v}_r+s_y{v}_s)+c_{13}(q_z{w}_q+r_z{w}_r+s_z{w}_s)\right]\\ +{\stackrel{-}{s}}_y\left[c_{66}(q_y{u}_q+r_y{u}_r+s_y{u}_s)+c_{66}(q_x{v}_q+r_x{v}_r+s_x{v}_s)\right]\\ +{\stackrel{-}{s}}_z\left[c_{44}(q_z{u}_q+r_z{u}_r+s_z{u}_s)+c_{44}(q_x{w}_q+r_x{w}_r+s_x{w}_s)\right]=0\end{array}$$E11

(11)

$$\begin{array}{l}{\stackrel{-}{s}}_x\left[c_{66}(q_x{v}_q+r_x{v}_r+s_x{v}_s)+c_{66}(q_y{u}_q+r_y{u}_r+s_y{u}_s)\right]\\ +{\stackrel{-}{s}}_y\left[c_{11}(q_y{v}_q+r_y{v}_r+s_y{v}_s)+c_{12}(q_x{u}_q+r_x{u}_r+s_x{u}_s)+c_{13}(q_z{w}_q+r_z{w}_r+s_z{w}_s)\right]\\ +{\stackrel{-}{s}}_z\left[c_{44}(q_z{v}_q+r_z{v}_r+s_z{v}_s)+c_{44}(q_y{w}_q+r_y{w}_r+s_y{w}_s)\right]=0,\end{array}$$E12

(12)

$$\begin{array}{l}{\stackrel{-}{s}}_x\left[c_{44}(q_x{w}_q+r_x{w}_r+s_x{w}_s)+c_{44}(q_z{u}_q+r_z{u}_r+s_z{u}_s)\right]\\ +{\stackrel{-}{s}}_y\left[c_{44}(q_y{w}_q+r_y{w}_r+s_y{w}_s)+c_{44}(q_z{v}_q+r_z{v}_r+s_z{v}_s)\right]\\ +{\stackrel{-}{s}}_z\left[c_{33}(q_z{w}_q+r_z{w}_r+s_z{w}_s)+c_{13}(q_x{u}_q+r_x{u}_r+s_x{u}_s)+c_{13}(q_y{v}_q+r_y{v}_r+s_y{v}_s)\right]=0\end{array}$$E13

(13)

Note that here the normal is represented by the normalized metric (evaluated along the free surface)

$${\stackrel{-}{s}}_x=\frac{s_x}{\sqrt{s_x^2+s_y^2+s_z^2}},$$$${\stackrel{-}{s}}_y=\frac{s_y}{\sqrt{s_x^2+s_y^2+s_z^2}},$$$${\stackrel{-}{s}}_z=\frac{s_z}{\sqrt{s_x^2+s_y^2+s_z^2}}.$$

5. A discretization scheme on curvilinear grid

To approximate (7)-(9) we discretize the rectangular solid (Figure 2)

Figure 2.

Grids distributions in curvilinear coordinate.

$$\begin{array}{l}{q}_i=(i-1)h_q\\ r_j=(j-1)h_r\\ s_k=(k-1)h_s\end{array}$$$$\begin{array}{l}i=1,\mathrm{...},N_q\\ j=1,\mathrm{...},N_r\\ k=1,\mathrm{...},N_s\end{array}$$

$$\begin{array}{l}{h}_q=l/(N_q-1);\\ h_r=w/(N_r-1);\\ h_s=h/(N_s-1)\end{array}$$E14

(14)

where l, w, h are the length of the rectangular solid in q-, r- and s-directions, respectively; $$h_q,h_r,h_s>0$$define the grid size in q-, r- and s-directions, respectively. The three components of the wave field are given by

$$[u_{i,j,k}(t),v_{i,j,k}(t),w_{i,j,k}(t)]=[u(q_i,r_j,s_k,t),v(q_i,r_j,s_k,t),w(q_i,r_j,s_k,t)]$$

and the derivation operators as

$$\begin{array}{l}{D}_{+}^q{u}_{i,j,k}=\frac{u_{i+1,j,k}-u_{i,j,k}}{h_q}\\ D_{+}^r{u}_{i,j,k}=\frac{u_{i,j+1,k}-u_{i,j,k}}{h_r}\\ D_{+}^s{u}_{i,j,k}=\frac{u_{i,j,k+1}-u_{i,j,k}}{h_s}\end{array}$$$$\begin{array}{l}{D}_{-}^q{u}_{i,j,k}=D_{+}^q{u}_{i-1,j,k}\\ D_{-}^r{u}_{i,j,k}=D_{+}^r{u}_{i,j-1,k}\\ D_{-}^s{u}_{i,j,k}=D_{+}^s{u}_{i,j,k-1}\end{array}$$

$$\begin{array}{l}{D}_0^q{u}_{i,j,k}=\frac{1}{2}(D_{+}^q{u}_{i,j,k}+D_{-}^q{u}_{i,j,k})\\ D_0^r{u}_{i,j,k}=\frac{1}{2}(D_{+}^r{u}_{i,j,k}+D_{-}^r{u}_{i,j,k})\\ D_0^s{u}_{i,j,k}=\frac{1}{2}(D_{+}^s{u}_{i,j,k}+D_{-}^s{u}_{i,j,k})\end{array}$$E15

(15)

The right hand sides of eqs. (7)-(9) contain spatial derivatives of nine basic types, which are discretized according to the following equations

$$\begin{array}{l}\frac{\partial }{\partial q}(a{\omega }_q)\approx D_{-}^q(E_{1/2}^q(a)D_{+}^q\omega )\\ \frac{\partial }{\partial r}(d{\omega }_q)\approx D_0^r(d{D}_0^q\omega )\\ \frac{\partial }{\partial s}(g{\omega }_q)\approx \widetilde{D_0^s}(g{D}_0^q\omega )\end{array}$$$$\begin{array}{l}\frac{\partial }{\partial q}(b{\omega }_r)\approx D_0^q(b{D}_0^r\omega )\\ \frac{\partial }{\partial r}(e{\omega }_r)\approx D_{-}^r(E_{1/2}^r(e)D_{+}^r\omega )\\ \frac{\partial }{\partial s}(m{\omega }_r)\approx \widetilde{D_0^s}(m{D}_0^r\omega )\end{array}$$

$$\begin{array}{l}\frac{\partial }{\partial q}(c{\omega }_s)\approx D_0^q(c\widetilde{D_0^s}\omega )\\ \frac{\partial }{\partial r}(f{\omega }_s)\approx D_0^r(f\widetilde{D_0^s}\omega )\\ \frac{\partial }{\partial s}(p{\omega }_s)\approx D_{-}^s(E_{1/2}^s(p)D_{+}^s\omega )\end{array}$$E16

(16)

Here $$\omega$$ represents u, v or w; a, b, c, d, e, f, g, m and p are combinations of metric and material coefficients. We introduce the following averaging operators:

$$\begin{array}{l}{E}_{1/2}^q({\gamma }_{i,j,k})={\gamma }_{i+1/2,j,k}=\frac{{\gamma }_{i,j,k}+{\gamma }_{i+1,j,k}}{2}\\ E_{1/2}^r({\gamma }_{i,j,k})={\gamma }_{i,j+1/2,k}=\frac{{\gamma }_{i,j,k}+{\gamma }_{i,j+1,k}}{2}\\ E_{1/2}^s({\gamma }_{i,j,k})={\gamma }_{i,j,k+1/2}=\frac{{\gamma }_{i,j,k}+{\gamma }_{i,j,k+1}}{2}\end{array}$$E17

(17)

The cross terms which contain a normal derivative on the boundary are discretized one-sided in the direction normal to the boundary:

$$\widetilde{D_0^s}{u}_{i,j,k}=\{\begin{array}{cc}{D}_{+}^s{u}_{i,j,k},& k=1,\\ D_0^s{u}_{i,j,k},& k\ge 2.\end{array}$$E18

(18)

5.1. A discretization on curvilinear grid: Elastic wave equations

We approximate the spatial operators in eqs. (7)-(9) by (16). After suppressing grid indexes, this leads to

$$\begin{array}{c}J\rho \frac{{\partial }^{2}u}{\partial t^2}=D_{-}^q[E_{1/2}^q(M_1^{qq})D_{+}^{q}u+E_{1/2}^q(M_2^{qq})D_{+}^{q}v+E_{1/2}^q(M_3^{qq})D_{+}^{q}w]\\ +D_0^q[M_1^{qs}\widetilde{D_0^s}u+M_2^{qs}\widetilde{D_0^s}v+M_3^{qs}\widetilde{D_0^s}w]\\ +D_0^r[M_1^{rs}\widetilde{D_0^s}u+M_2^{rs}\widetilde{D_0^s}v+M_3^{rs}\widetilde{D_0^s}w]\\ +\widetilde{D_0^s}[M_1^{sq}{D}_0^{q}u+M_2^{sq}{D}_0^{q}v+M_3^{sq}{D}_0^{q}w]\\ +\widetilde{D_0^s}[M_1^{sr}{D}_0^{r}u+M_2^{sr}{D}_0^{r}v+M_3^{sr}{D}_0^{r}w]\\ +D_0^q[M_1^{qr}{D}_0^{r}u+M_2^{qr}{D}_0^{r}v+M_3^{qr}{D}_0^{r}w]\\ +D_0^r[M_1^{rq}{D}_0^{q}u+M_2^{rq}{D}_0^{q}v+M_3^{rq}{D}_0^{q}w]\\ +D_{-}^r[E_{1/2}^r(M_1^{rr})D_{+}^{r}u+E_{1/2}^r(M_2^{rr})D_{+}^{r}v+E_{1/2}^r(M_3^{rr})D_{+}^{r}w]\\ +D_{-}^s[E_{1/2}^s(M_1^{ss})D_{+}^{s}u+E_{1/2}^s(M_2^{ss})D_{+}^{s}v+E_{1/2}^s(M_3^{ss})D_{+}^{s}w]\\ \equiv L^{(u)}(u,v,w)\end{array}$$E19

(19)

$$\begin{array}{c}J\rho \frac{{\partial }^{2}v}{\partial t^2}=D_{-}^q[E_{1/2}^q(M_5^{qq})D_{+}^{q}v+E_{1/2}^q(M_2^{qq})D_{+}^{q}u+E_{1/2}^q(M_4^{qq})D_{+}^{q}w]\\ +D_0^q[M_5^{qs}\widetilde{D_0^s}v+M_2^{sq}\widetilde{D_0^s}u+M_4^{qs}\widetilde{D_0^s}w]\\ +D_0^r[M_5^{rs}\widetilde{D_0^s}v+M_2^{sr}\widetilde{D_0^s}u+M_4^{rs}\widetilde{D_0^s}w]\\ +\widetilde{D_0^s}[M_5^{sq}{D}_0^{q}v+M_2^{qs}{D}_0^{q}u+M_4^{sq}{D}_0^{q}w]\\ +\widetilde{D_0^s}[M_5^{sr}{D}_0^{r}v+M_2^{rs}{D}_0^{r}u+M_4^{sr}{D}_0^{r}w]\\ +D_0^q[M_5^{qr}{D}_0^{r}v+M_2^{rq}{D}_0^{r}u+M_4^{qr}{D}_0^{r}w]\\ +D_0^r[M_5^{rq}{D}_0^{q}v+M_2^{qr}{D}_0^{q}u+M_4^{rq}{D}_0^{q}w]\\ +D_{-}^r[E_{1/2}^r(M_5^{rr})D_{+}^{r}v+E_{1/2}^r(M_2^{rr})D_{+}^{r}u+E_{1/2}^r(M_4^{rr})D_{+}^{r}w]\\ +D_{-}^s[E_{1/2}^s(M_5^{ss})D_{+}^{s}v+E_{1/2}^s(M_2^{ss})D_{+}^{s}u+E_{1/2}^s(M_4^{ss})D_{+}^{s}w]\\ \equiv L^{(v)}(u,v,w)\end{array}$$E20

(20)

$$\begin{array}{c}J\rho \frac{{\partial }^{2}w}{\partial t^2}=D_{-}^q[E_{1/2}^q(M_3^{qq})D_{+}^{q}u+E_{1/2}^q(M_4^{qq})D_{+}^{q}v+E_{1/2}^q(M_6^{qq})D_{+}^{q}w]\\ +D_0^q[M_3^{sq}\widetilde{D_0^s}u+M_4^{sq}\widetilde{D_0^s}v+M_6^{qs}\widetilde{D_0^s}w]\\ +D_0^r[M_3^{sr}\widetilde{D_0^s}u+M_4^{sr}\widetilde{D_0^s}v+M_6^{rs}\widetilde{D_0^s}w]\\ +\widetilde{D_0^s}[M_3^{qs}{D}_0^{q}u+M_4^{qs}{D}_0^{q}v+M_6^{sq}{D}_0^{q}w]\\ +\widetilde{D_0^s}[M_3^{rs}{D}_0^{r}u+M_4^{rs}{D}_0^{r}v+M_6^{sr}{D}_0^{r}w]\\ +D_0^q[M_3^{rq}{D}_0^{r}u+M_4^{rq}{D}_0^{r}v+M_6^{qr}{D}_0^{r}w]\\ +D_0^r[M_3^{qr}{D}_0^{q}u+M_4^{qr}{D}_0^{q}v+M_6^{rq}{D}_0^{q}w]\\ +D_{-}^r[E_{1/2}^r(M_3^{rr})D_{+}^{r}u+E_{1/2}^r(M_4^{rr})D_{+}^{r}v+E_{1/2}^r(M_6^{rr})D_{+}^{r}w]\\ +D_{-}^s[E_{1/2}^s(M_3^{ss})D_{+}^{s}u+E_{1/2}^s(M_4^{ss})D_{+}^{s}v+E_{1/2}^s(M_6^{ss})D_{+}^{s}w]\\ \equiv L^{(w)}(u,v,w)\end{array}$$E21

(21)

in the grid points$$(q_i,r_j,s_k)$$,$$(i,j,k)\in [1,N_q]\times [1,N_r]\times [1,N_s]$$. We have introduced the following notations for the material and metric terms in order to express the discretized equations in a more compact form:

$$\begin{array}{l}{M}_1^{kl}=J{k}_x{l}_x{c}_{11}+J{k}_y{l}_y{c}_{66}+J{k}_z{l}_z{c}_{44}{M}_2^{kl}=J{k}_x{l}_y{c}_{12}+J{k}_y{l}_x{c}_{66}\\ M_3^{kl}=J{k}_x{l}_z{c}_{13}+J{k}_z{l}_x{c}_{44}{M}_4^{kl}=J{k}_y{l}_z{c}_{13}+J{k}_z{l}_y{c}_{44}\\ M_5^{kl}=J{k}_x{l}_x{c}_{66}+J{k}_y{l}_y{c}_{11}+J{k}_z{l}_z{c}_{44}{M}_6^{kl}=J{k}_x{l}_x{c}_{44}+J{k}_y{l}_y{c}_{44}+J{k}_z{l}_z{c}_{33}\end{array}$$E22

(22)

where k and l represent the metric coefficients q, r or s.

We discretize in time using second-order accurate centered differences. The full set of discretized equations is

$$\begin{array}{l}\rho (\frac{u^{n+1}-2u^n+u^{n-1}}{{\delta }_t^2})=L^{(u)}(u^n,v^n,w^n)\\ \rho (\frac{v^{n+1}-2v^n+v^{n-1}}{{\delta }_t^2})=L^{(v)}(u^n,v^n,w^n)\\ \rho (\frac{w^{n+1}-2w^n+w^{n-1}}{{\delta }_t^2})=L^{(w)}(u^n,v^n,w^n)\end{array}$$E23

(23)

where $${\delta }_t$$ represents the time step.

5.2. A discretization on curvilinear grid: Free boundary conditions

The boundary conditions (11)-(13) are discretized by

$$\begin{array}{l}\frac{1}{2}[{(}_{M_1^{ss}}{D}_{+}^s{u}_{i,j,1}+{(}_{M_1^{ss}}{D}_{+}^s{u}_{i,j,0}]+{(}_{M_1^{sq}}{D}_0^q{u}_{i,j,1}+{(}_{M_2^{sq}}{D}_0^q{v}_{i,j,1}\\ +{(}_{M_3^{sq}}{D}_0^q{w}_{i,j,1}+\frac{1}{2}[{(}_{M_2^{ss}}{D}_{+}^s{v}_{i,j,1}+{(}_{M_2^{ss}}{D}_{+}^s{v}_{i,j,0}]+{(}_{M_1^{sr}}{D}_0^r{u}_{i,j,1}\\ +{(}_{M_2^{sr}}{D}_0^r{v}_{i,j,1}+{(}_{M_3^{sr}}{D}_0^r{w}_{i,j,1}+\frac{1}{2}[{(}_{M_3^{ss}}{D}_{+}^s{w}_{i,j,1}+{(}_{M_3^{ss}}{D}_{+}^s{w}_{i,j,0}]=0\\ \end{array}$$E24

(24)

$$\begin{array}{l}\frac{1}{2}[{(}_{M_5^{ss}}{D}_{+}^s{v}_{i,j,1}+{(}_{M_5^{ss}}{D}_{+}^s{v}_{i,j,0}]+{(}_{M_5^{sq}}{D}_0^q{v}_{i,j,1}+{(}_{M_2^{qs}}{D}_0^q{u}_{i,j,1}\\ +{(}_{M_4^{sq}}{D}_0^q{w}_{i,j,1}+\frac{1}{2}[{(}_{M_2^{ss}}{D}_{+}^s{u}_{i,j,1}+{(}_{M_2^{ss}}{D}_{+}^s{u}_{i,j,0}]+{(}_{M_5^{sr}}{D}_0^r{v}_{i,j,1}\\ +{(}_{M_2^{rs}}{D}_0^r{u}_{i,j,1}+{(}_{M_4^{sr}}{D}_0^r{w}_{i,j,1}+\frac{1}{2}[{(}_{M_4^{ss}}{D}_{+}^s{w}_{i,j,1}+{(}_{M_4^{ss}}{D}_{+}^s{w}_{i,j,0}]=0\\ \end{array}$$E25

(25)

$$\begin{array}{l}\frac{1}{2}[{(}_{M_3^{ss}}{D}_{+}^s{u}_{i,j,1}+{(}_{M_3^{ss}}{D}_{+}^s{u}_{i,j,0}]+{(}_{M_3^{qs}}{D}_0^q{u}_{i,j,1}+{(}_{M_4^{qs}}{D}_0^q{v}_{i,j,1}\\ +{(}_{M_6^{sq}}{D}_0^q{w}_{i,j,1}+\frac{1}{2}[{(}_{M_4^{ss}}{D}_{+}^s{v}_{i,j,1}+{(}_{M_4^{ss}}{D}_{+}^s{v}_{i,j,0}]+{(}_{M_3^{rs}}{D}_0^r{u}_{i,j,1}\\ +{(}_{M_4^{rs}}{D}_0^r{v}_{i,j,1}+{(}_{M_6^{sr}}{D}_0^r{w}_{i,j,1}+\frac{1}{2}[{(}_{M_6^{ss}}{D}_{+}^s{w}_{i,j,1}+{(}_{M_6^{ss}}{D}_{+}^s{w}_{i,j,0}]=0\end{array}$$E26

(26)$$i=1,\mathrm{...},N_q;j=1,\mathrm{...},N_r.$$

The key step in obtaining a stable explicit discretization is to use the operator$$\widetilde{D_0^s}$$ ( which is one-sided on the boundary) for the approximation of the normal derivative in $${\partial }_q{\partial }_s,{\partial }_r{\partial }_s,{\partial }_s{\partial }_q$$ and $${\partial }_s{\partial }_r$$ cross derivatives. At first glance, it may appear that using a one-sided operator the accuracy of the method would be reduced to the first-order. However, as it was theoretically shown by Nilsson et al. (2007) (for a Cartesian discretization), a first-order error on the boundary in the differential equations (19)-(21) can be absorbed as a second-order perturbation of the boundary conditions (24)-(26).

Figure 3.

Model of a half-space with a planar free surface.

6. Accuracy and efficiency tests

6.1. Accuracy

The accuracy of the proposed method is examined by comparing numerical results with the analytical solution of the Lamb’s problem, for a transversely isotropic medium with a vertical symmetry axis (VTI medium). The elastic parameters describing the VTI medium are given in Table 1. The analytical solution is obtained by convolving the free-surface Green-function with the source function (Payton, 1983). A vertical point source of the type

$$c_{11}$$; (GPa)	$$c_{12}$$; (GPa)	$$c_{13}$$; (GPa)	$$c_{33}$$; (GPa)	$$c_{44}$$; (GPa)	$$\rho$$; $$(g/c{m}^3)$$
25.5	2.0	14.0	18.4	5.6	2.4

Table 1.

Medium parameters in the homogeneous half-space

$$f(t)=e^{-0.5f_0^2{(t-t_0)}^2}cos\pi f_0(t-t_0)$$E27

(27)

with $$t_0$$= 0.5 s and a high cut-off frequency $$f_0$$= 10 Hz, is assumed to be located at (300 m, 2000 m) at the surface, which is marked as an asterisk in Figure 3. It should be mentioned that Carcione et al. (1992) and Carcione (2000) presented an analytical comparison of the point-source response in a 3-D VTI medium in the absence of the free surface. The comparisons are performed by first transforming the 3-D numerical results into a line-source response by carrying out an integration along the receiver line (Wapenaar et al., 1992) and then comparing the emerging results with the 2-D Lamb’s analytical solutions. The numerical model contains 401 x 401 x 191 grid nodes in the x-, y- and z-direction, respectively. The grid spacings are 10 m in all directions. The solution is advanced using a time step of 1.25 ms for 3.5 s.

Figure 4.

Comparison between numerical and analytical vertical components of the displacement for the VTI medium.

Three receiver lines are positioned on the free surface, two of which are parallel to the y-direction with respective normal distances of 130 (Line 1) and 1000 m (Line 2) away from the point source, the other crosses the source location and parallels to the x-direction (Line 3). The integrations are performed along the first two receiver lines, these represent 2-D results of 130 and 1000 m away from the source, respectively. Figure 4 shows the comparisons between the resulting numerical and 2-D analytical z-components of the displacement for the VTI medium. In spite of the errors resulting from the transformation of the point-source response into the line-source one, numerical and analytical results agree well in Figure 4. These comparisons demonstrate the accuracy of our corresponding algorithm.

Figure 5.

Seismogram sections at Line 3 for the planar surface model. Symbol qP indicates the qP wave and R indicates Rayleigh wave.

Figure 6.

Snapshots of the vertical component of the wavefield at the surface (xy-plane) of the planar surface model.

Synthetic seismograms are computed at Line 3. The seismograms in Figure 5 show the direct quasi-P wave (qP) and a high-amplitude Rayleigh wave (R). Snapshots of the vertical component of the wavefield in horizontal (xy-) plane at the propagation time of 1.4 s are displayed in Figure 6. We define the incidence plane by the propagation direction and the z-axis, quasi-P wave and quasi-SV wave (qSV) motions lie in this plane, while SH motion is normal to the plane. Hence, the z-component does not contain SH motion. The xy-plane of a transversely isotropic medium is a plane of isotropy, where the velocity of the qP wave is about 3260 $$m{s}^{-1}$$ and the velocity of the qSV wave is about 1528$$m{s}^{-1}$$. The amplitude of the qP wave is so weak compared with that of the Rayleigh and qSV wave that one can hardly identify it in the snapshot (Figure 6a). In order to observe the qP wave, a gain has been given to the amplitudes of the wavefield. Owing to this, side reflections also appear in the photo, as shown in Figure 6 b. As the velocity of the Rayleigh wave is very close to that of the qSV wave, the two waves are almost superimposed and it is difficult to distinguish between the two in synthetic seismograms and snapshots.

Figure 7.

Snapshots of the x-component of the wavefield in the vertical (xz-) plane which contains the receiver line and the source at 1.4 (a) and 2.3 s (b) propagation times for the planar surface model.

Figure 7 shows the x-component of the wavefield in the vertical (xz-) plane at 1.4 and 2.3 s propagation times. The xz-plane contains the receiver line (Line 3) and the source location. Both snapshots show the wave front of the qP-wave and the qSV-wave. The former snapshot (1.4 s) shows the qSV-wave with the cusps. A headwave (H) can also be found in the photos, the headwave is a quasi-shear wave and is guided along the surface by the qP-wave.

6.2. Numerical simulations on an irregular (non-flat) free surface

Three numerical experiments with irregular free surfaces are now investigated. The first example is a test on smooth boundaries, while the second example consists of a hemispherical depression to test the ability of the method for non-smooth topography. For sake of simplicity both these examples are based on homogeneous half-spaces, i.e., the medium parameters are the same as in the case of flat surface (Table 1). The same source is located at the same place as in the planar surface model, the time step is 0.8 ms. The total propagation time is 3.5 s for the two models. Finally, we consider a two-layered model with a realistic topography.

6.2.1. Topography simulating a shaped Gaussian hill

The first model considered here is a half-space whose free surface is a hill-like feature (Figure 8). The shape of the hill resembles a Gaussian curved surface given by the function

$$z(x,y)=-150exp(-{(\frac{x-2000}{150})}^2-{(\frac{y-2000}{150})}^2)m$$

$$(x,y)\in {[0m,4000m]}^2$$E28

(28)

Figure 8.

Model of a half-space with Gaussian shape hill topography.

The computational domain extends to depth z(x, y)=2000 m. The volume is discretized with equal grid nodes in each direction as in the planar surface one. The grid spacings are 10 m in the x, y-directions and about 10.5 m in the z-direction for average. The vertical spacing varies with depth, it is smaller toward the free surface and larger toward the bottom of the model. The minimum and maximum of the vertical spacings are 6 and 12 m, respectively.

The gridding scheme which shows the detailed cross-section of the grids along Line3 is shown in Figure 9. Synthetic seismograms are also computed at Line3 (Figure 10). As a result of the hill-shaped free surface (and compare with the synthetic seismograms in Figure 5), the amplitudes of the quasi-P wave and Rayleigh wave are reduced in the right-hand part of the sections. In addition, after the ordinary quasi-P wave a secondary quasi-P wave (RqPf) induced by the scattering of the direct Rayleigh wave can be observed. Similarly, a secondary Rayleigh wave (qPRf) which travels in front of the ordinary Rayleigh wave induced by the scattering of the direct quasi-P wave can also be distinguished. Some energy is scattered back to the left-hand side as a Rayleigh wave (qPRb, RR) and a quasi-P wave (RqPb).

Figure 9.

The gridding scheme which shows the detailed cross-section of the grids along Line3 in the Gaussian shape hill topography model. For clarity, the grids are displayed with a reducing density factor of 3.

Figure 10.

Seismograms along the receiver line for the Gaussian shape hill topography model: (a) x-component (horizontal) of the displacement; (b) z-component (vertical). Symbols mean the following: (qPd) qP wave diffracts to qP wave; (Rd) Rayleigh wave diffracts to Rayleigh wave; (qPRf) qP wave scatters to Rayleigh wave and propagates forward; (qPRb) qP wave scatters to Rayleigh wave and propagates backward; (RqPf) Rayleigh wave scatters to qP wave and propagates forward; (RqPb) Rayleigh wave scatters to qP wave and propagates backward; (RR) Rayleigh wave reflectes to Rayleigh wave.

Figure 11.

Snapshots of the vertical component of the wavefield at the surface (xy-plane) at different propagation times for the Gaussian shape hill topography model.

Snapshots of the wavefield in the horizontal (xy-) plane at different propagation times are displayed in Figures 11. The amplitudes are also gained. In the beginning the wavefield propagates undisturbed along the free surface. At 1.1 s the direct quasi-P wave hits the hill and generates a circular diffracted wave. This wave is a Rayleigh wave, which is marked as two parts, one travels forward (qPRf), and the other travels backward (qPRb). These can be seen clearly in the later snapshots (1.4 - 2.3 s). In addition, a reflected Rayleigh wave (RR) can be observed. The direct quasi-P wave (qP) and Rayleigh wave (R) are also marked in the figure. By the way, side reflections from the boundaries can also be noted in the plane. Figure 12 shows the x-component of the wavefield in the vertical (xz-) plane. The xz-plane contains the receiver line and source location. The snapshots show the diffracted quasi-P and quasi-SV waves clearly in the vertical plane.

Figure 12.

Snapshots of the x-component of the wavefield in the vertical (xz-) plane which contains the receiver line and the source at different propagation times for the Gaussian shape hill topography model.

6.2.2. Topography simulating a shaped hemispherical depression

In the second model, we consider a hemispherical depression model as illustrated in Figure 13. The first model that we have considered is of smooth topography, that is, with continuous and finite slopes everywhere. However, the shaped hemispherical depression here taken as reference is a case of extreme topography, such that the vertical-to-horizontal ratio of the depression is very large (1:2) and the slopes of the edges tend to infinity. The hemispherical depression is at the middle of the free surface and the radius is 150 m.

Figure 13.

Model of a half-space with hemispherical shape depression topography.

Figure 14.

The gridding scheme which shows the detailed cross-section of the grids along Line3 in the hemispherical shape depression topography model. For clarity, the grids are displayed with a reducing density factor of 3.

The numerical model is discretized in the same way as in the hill topography model. The gridding scheme which shows the detailed cross-section of the grids along Line3 is shown in Figure 14. Owing to the existence of model edges with strong slopes at x=1850 and x=2150 m along the receiver line, both body and Rayleigh waves scattered by sharp changes in the topography can be clearly observed on the synthetic seismograms shown in Figure 15. Owing to its shorter wavelength, the scattering of Rayleigh wave is much stronger than that of the body wave when propagating through the hemispherical depression, this indicating that such sharp depression can affect the propagation of Rayleigh wave significantly.

Figure 15.

Seismograms at the receiver line for the hemispherical shape depression topography model: (a) x-component (horizontal) of the displacement; (b) z-component (vertical). Symbols mean the same as in Fig. 10.

The photos in Figure 16 show the vertical component of the wavefield in the horizontal (xy-) plane. Compared with the photos of the hill topography model, we can see the Rayleigh wave scattering at the edges of the hemispherical depression; it seems as if the reflected Rayleigh wave propagating faster in the hemispherical depression model than that in the hill topography model. What’s more, the back scattered waves of Rayleigh wave in the hemispherical depression model are much stronger, this may also indicating that such sharp depression blocks the propagation of Rayleigh wave more significantly.

6.2.3. Real topography simulating

It is also interesting to study a realistic example. We consider a model in Tibet (Figure 17). The length and width of the model are 21.6 km, and the “average” height of the topography is roughly -3560 m (3560 m in the geodetic coordinate system). The computational domain is extended to depth z(x, y)= 7200 m. For simplicity we use a two-layered model with parameters given in the model sketch (Figure 17) instead of the “real” velocity structure under the realistic topography. It consists of 241×241×121 grid nodes in the x-, y-, and z-direction, respectively, with equal vertical grid nodes in each layer. A vertical point source like the used in previous models is loaded in the middle of the free surface (indicated by the asterisk in Fig. 17), where the high cut-off frequency has been changed to 2.7 Hz and the time-shift is 1.5 s. Two lines of receivers crossing the source location and paralleling to the x- and y-direction respectively are placed at the free surface. The time step is 5 ms, and the total propagation time is 8 s.

Figure 16.

Snapshots of the vertical component of the wavefield at the surface (xy-plane) at different propagation times for the hemispherical shape depression topography model.

Figure 17.

A two-layered model with a realistic topography. The medium parameters of each layer are also given in the figure. The units for the elasticity and density are GPa and$$g/c{m}^3$$, respectively.

Snapshots of the z-component of the wavefield in the vertical plane which contains receiver line Line1 and the source location are presented on Figure 18, and the seismograms of the z-component are also computed at the two receiver lines (Figure 19). We can see that the effect of the topography is very important, with strong scattered phases that are superimposed to the direct and reflected waves, which make it very difficult to identify effective reflections from subsurface interface. The scattering in the seismograms also reflect different features of the surface. The scattering in the seismograms at Line 1 (Figure 19a) is much stronger than that in the seismogram at Line 2 (Figure 19b), indicating that the surface along Line 1 is much rougher than that along Line 2, which also can be observed in Figure 17. What’s more, the scattering in Figure 19a is approximately uniformly distributed while in Figure 19b it is mostly distributed in the vicinity of the shot. These may due to different distributions of the surface topography along these two lines.

Figure 18.

Snapshots of the vertical component of the wavefield in the vertical (xz-) plane along Line 1 at different propagation times for the two-layered model with a realistic topography.

Figure 19.

Vertical-component synthetic seismograms coming from the two-layered model with real topography represented in Figure 19: (a) at the receiver line (Line 1) that crosses the source location and is parallel to x-direction (Fig. 19); (b) at the receiver line (Line 2) that crosses the source location and is parallel to y-direction.

7. Conclusion

We propose a stable and explicit finite difference method to simulate with second-order accuracy the propagation of seismic waves in a 3D heterogeneous transversely isotropic medium with non-flat free surface. The method is an extension of the 2D method proposed by Appelo and Petersson (2009) to the 3D anisotropic case. The surface topography is introduced via mapping rectangular grids to curved grids. The accurate application of the free surface boundary conditions is done by using boundary-modified difference operators to discretize the mixed derivatives in the governing equations of the problem. Several numerical examples under different assumptions of free surface are given to highlight the complications of realistic seismic wave propagation in the vicinity of the earth surface. Synthetic seismograms and snapshots explain diffractions, scattering, multiple reflections, and converted waves provoked by the features of the free surface topography. The typical cuspidal triangles of the quasi-transverse (qS) mode also appear in the snapshots of the anisotropic medium.

The future directions of our research will include an extension of the schemes to the viscoelastic case. This will allow a realistic attenuation of the seismic waves due to the presence of a weathered layer to be included (2000).

Acknowledgement

We are grateful to Jinlai Hao and Jinhai Zhang for their assistance and the facilities given in the course of this work. Fruitful discussions with Tao Xu, Kai Liang, Wei Zhang are also greatly appreciated. The Ministry of Science and Technology of China (SINOPROBE-02-02) supported this research.