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

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

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:

and similarly in other cases

where

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

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:

(5) |

where

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

(7) |

(8) |

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

(10) |

Here

(11) |

(12) |

(13) |

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

## 5. A discretization scheme on curvilinear grid

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

where l, w, h are the length of the rectangular solid in q-, r- and s-directions, respectively;

and the derivation operators as

(15) |

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

Here

(17) |

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

(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

(19) |

(20) |

(21) |

in the grid points

(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

(23) |

where

### 5.2. A discretization on curvilinear grid: Free boundary conditions

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

(24) |

(25) |

(26) |

The key step in obtaining a stable explicit discretization is to use the operator

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

(27)with

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.

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

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

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

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.

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

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.

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.

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.

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