# Discrete wavenumber solutions to numerical wave propagation in piecewise heterogeneous media - I. Theory of two-dimensional SH case

Li-Yun Fu, Michel Bouchon

► To cite this version:

Li-Yun Fu, Michel Bouchon. Discrete wavenumber solutions to numerical wave propagation in piecewise heterogeneous media - I. Theory of two-dimensional SH case. Geophysical Journal International, 2004, 157, pp.481-498. 10.1111/j.1365-246X.2004.02135.x. insu-03607131

HAL Id: insu-03607131

https://insu.hal.science/insu-03607131

Submitted on 13 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
relatively consistent lithologic constitution, spreading over a relatively large tectonic region. Heterogeneities inside the geological sequence can be described by a reference medium property accompanied with random medium perturbations of different scales. Wave scattering by such small-scale volume heterogeneities, leading to amplitude/phase fluctuations and wave attenuation, has been extensively studied (e.g. Chernov 1960; Aki 1980; Wu 1989a). Numerical modelling schemes have to be highly accurate to yield reliable results about subtle scattering effects by volume heterogeneities. On the other hand, the irregular interfaces that separate different geological formations, constitute a large-scale geometric system with strong-contrast impedance boundaries to control the principal characteristics of wave propagation (e.g. Kennett 1984, 1986; Campillo 1987a, 1990; Maupin 1989). It poses a special requirement for seismic modelling schemes to simulate subtle reflection/transmission effects across interfaces. Clearly, such inherent piecewise heterogeneous media of sedimentary basins challenge numerical modelling techniques being able to separate large-scale boundary waves and small-scale volume-scattering waves, and subsequently to handle them separately with scales of approximate accuracies in terms of incident wavelengths. The purpose of this paper is to develop such media-oriented numerical techniques for wave propagation in piecewise heterogeneous media.

The finite-difference and finite-element methods are universal numerical techniques to simulate wave propagation in complex media. However, these methods are characteristic of an implicit use of boundary continuity conditions across interfaces and hence cannot separate wavefields that come from different types of media. The explicit use of boundary continuity conditions will lead to a category of semi-analytical, semi-numerical methods that are more accurate in the simulation of reflection/transmission across irregular interfaces. These semi-analytical, semi-numerical methods either directly solve the boundary integral equation using boundary conditions or expand the solution into a set of plane waves with coefficients to be determined using boundary conditions. The boundary integral equation discrete wavenumber method (Bouchon 1982; Campillo & Bouchon 1985; Bouchon et al. 1989) incorporates the discrete wavenumber Green's function representation into the boundary integral equation techniques to model multilayered media having irregular interfaces. The discrete wavenumber decomposition allows the consideration of more complex structures (Campillo 1987b) and has been used for regional seismograms (Campillo & Paul 1992). The method has been extended to P-SV waves for 2-D (Gaffet & Bouchon 1989; Kawase & Aki 1989) and 3-D problems (Bouchon et al. 1996) for studying the effect of topography on surface motion. Chen (1990) combined the wavefield plane wave expansion, the discrete wavenumber Green's function representation and the boundary integral equation to develop a global generalized reflection/transmission matrices method for irregular layered media. The application to the 1987 Whittier–Narrows earthquake was demonstrated (Chen 1995) for interpretation of the observed records. This method can be viewed as an extension of the generalized R/T coefficients method for horizontally layered media (Kennett 1983; Luco & Apsel 1983) to irregular layered cases. The coupled mode method (Kennett 1984; Maupin & Kennett 1987; Kennett 1998) represents the wave train as a superposition of modal contributions for the flat-layered reference structure and incorporates the effect of lateral variations at irregular interfaces into the matrices of reflection and transmission coefficients of the different modes.

An efficient alternative with an explicit use of boundary continuity conditions is to directly solve the boundary integral equation using the boundary-element (BE) method. The method has been extensively used to study the effects of topography or sedimentary basin structures on ground motions at the surface during the last several decades (e.g. Wong & Jennings 1975; Sánchez-Sesma & Esquivel 1979; Dravinksi 1983; Lee & Langston 1983; Sánchez-Sesma & Campillo 1991; Gibson & Campillo 1994) and has gained popularity among theoretical seismologists. The BE method can be directly extended to complex geological structures (Fu 1996) for exploration-oriented seismic modelling. Because of the computational intensity, the BE method is often limited to cases of low frequencies for regional wave propagation at regional distances. At higher frequencies, the computation becomes extremely prohibitive because of the very large size of the resulting matrices to be inverted. Using a boundary connection technique, an efficient modification to the BE method has been made to investigate regional wave propagation over more than 1000 km in the Tibet area (Fu & Wu 2001) and for theoretical study of the scattering effects of random topography on regional wave attenuation and amplification (Fu et al. 2002).

The major advantages of the plane wave decomposition and BE methods mentioned above are:

(i) the geometrically accurate description of irregular interfaces;

(ii) the explicit use of the boundary conditions of continuities for displacement and traction across interfaces; and

(iii) the ability to separate wavefields generated by different parts of the medium.

However, these boundary methods are limited by their abilities to handle volume heterogeneities. The presence of volume heterogeneities leads to a coupled boundary–volume interaction that requires an extension of the traditional boundary methods to accounting for volume heterogeneities. To pursue an explicit use of boundary conditions when modelling piecewise heterogeneous media, we must formulate wave propagation involving volume heterogeneities in terms of integral equations. Some progress has been made to model piecewise heterogeneous media using integral equation numerical techniques. Fu (2003) incorporates the boundary integral representation into the Lippmann–Schwinger integral equation to model rough topography in complex near-surface areas. The former is very convenient to deal with the traction-free condition and can give a full description to arbitrarily rugged topography with irregular acquisition geometry. The latter can be used to handle generally heterogeneous media in the near-surface region and has been accurately solved using a velocity-weighted wavefield (Fu et al. 1997). In the generalized Lippmann–Schwinger integral (GLSI) equation, all the integral kernels are related to the Green’s function in the background medium, which avoids the necessity of the Green’s function of a heterogeneous medium. The numerical method is used to study the influence of widely observed strong scattered noises caused by rugged topographies and strong volume heterogeneities in complex near-surface areas in exploration seismology. Fu (2002a) extends the GLSI equation numerical method to heterogeneous layered media having irregular interfaces.

In the present work, the discrete wavenumber representation is incorporated into the boundary–volume integral equation to simulate 2-D SH-wave propagation in multilayered media with volume heterogeneities. The incident, boundary-scattering and volume-scattering waves
<figure><img src="image_2.png" /><figcaption>Figure 1. Configuration of the problem considered.</figcaption></figure>

are respectively formulated in the discrete wavenumber domain and subsequently handled flexibly at various accuracies using approximation methods. In this context, the well-known singularity problem of the Green's function is avoided by choosing the frequency to be complex and using truncated series in place of actual Green's functions. The full-waveform method presented in the paper is able to deal with fast spatial fluctuations of media properties. It is highly accurate to be sensitive to subtle effects of arbitrary volume heterogeneities on wave propagation. The major disadvantage of the numerical method is the considerable computer time and memory requirements because of the prohibitively expensive boundary-volume integral equation system. In particular, inverting the matrix generated by the volume integral equation takes a very large part of computational time. An efficient approach to this problem is to transform the implicit volume integral equation into an explicit integration expression. The approximation is related to our primary concern of how much the volume heterogeneities contribute to the total wavefield in comparison with the boundary-scattering wave. In practical applications, the volume heterogeneities inside each geological formation might be relatively smooth spatially at seismic wavelengths. Wave propagation in piecewise heterogeneous media is actually influenced profoundly by the large-scale geometric system. Therefore, we develop approximate solutions to the boundary-volume integral equation. The solutions are approximated from the semi-explicit volume integral equation by the Fresnel-zone aperture theory to the full explicit volume integral equation by the Born series approximation. These approximations are scaled to adapt the numerical method to the smoothness of the volume media in terms of incident wavelength. It is worth mentioning that at crustal seismic frequencies the heterogeneous crustal waveguide might be so smooth that one-way approximation to regional wave propagation can be applicable (Wu 1996).

PROBLEM DEFINITION AND SCATTERING INTEGRAL REPRESENTATION

Piecewise heterogeneous media that the earth presents consist of a series of geological formations with different scales of internal volume heterogeneities, separated by irregular interfaces over a half-space. Wave propagation in such media can be formulated as the superposition of incident, boundary-scattering and volume-scattering waves. The problem configuration for the jth formation $\Omega^{(j)}$ is depicted in Fig. 1. It is bounded by the top interface $\partial\Omega^{(j-1)}$ and the bottom interface $\partial\Omega^{(j)}$, expressed respectively by $z = \xi^{(j-1)}(x)$ and $z = \xi^{(j)}(x)$. The side boundaries at two edges of $\Omega^{(j)}$ are assumed to extend to infinity. The uppermost interface $\partial\Omega^{(0)}$ is a free surface and an arbitrary seismic source is embedded in $\Omega^{(j)}$. For simplicity, the present study is restricted to the 2-D SH problem (or acoustic problem). The elastic properties in $\Omega^{(j)}$ are described by the shear modulus $\mu^{(j)}(\mathbf{r})$ and density $\rho^{(j)}(\mathbf{r})$ with the corresponding reference values $\mu_0^{(j)}$ and $\rho_0^{(j)}$. $\mathbf{r} = (x, z)$ is the position vector. The solution domain of the problem is defined as $\tilde{\Omega}^{(j)} \in \Omega^{(j)} + \partial\Omega^{(j-1)} + \partial\Omega^{(j)}$. Seismic response $u(\mathbf{r})$ for steady-state scalar wave propagation in $\tilde{\Omega}^{(j)}$ satisfies the following scalar equation

$$ \nabla^2 u(\mathbf{r}) + [K^{(j)}(\mathbf{r})]^2 u(\mathbf{r}) = -s^{(j)}(\mathbf{r}, \omega), \quad \mathbf{r} \in \tilde{\Omega}^{(j)}, \qquad (1) $$

where the wavenumber $[K^{(j)}(\mathbf{r})]^2 = \omega^2 \rho^{(j)}(\mathbf{r}) / \mu^{(j)}(\mathbf{r})$ with the corresponding reference wavenumber $[K_0^{(j)}]^2 = \omega^2 \rho_0^{(j)} / \mu_0^{(j)}$ and $s(\mathbf{r}, \omega)$ is the body force occupying a region $\Omega_s^{(j)}$. The seismic response $u(\mathbf{r})$ also satisfies the following boundary conditions:

(i) The traction-free condition on the free surface, $\partial u(\mathbf{r}) / \partial n = 0$ at $\mathbf{r} \in \partial \Omega^{(0)}$.

(ii) The continuities of displacement and traction at the interface $\partial \Omega^{(j)}$,

$$ \left\{ \begin{array}{l} u_{-}^{(j)}(\mathbf{r}) = u_{+}^{(j)}(\mathbf{r}) \\ \mu^{(j)} \frac{\partial u_{-}^{(j)}(\mathbf{r})}{\partial n} = \mu^{(j+1)} \frac{\partial u_{+}^{(j)}(\mathbf{r})}{\partial n}, \quad \mathbf{r} \in \partial \Omega^{(j)}, \end{array} \right. \qquad (2) $$

where ‘-’ denotes the top side of $\partial \Omega^{(j)}$ toward $\Omega^{(j)}$ and ‘+’ denotes the underside of $\partial \Omega^{(j)}$ toward $\Omega^{(j+1)}$.
(iii) The radiation boundary condition imposed on the far-field behavior at infinity,

$$ \left\{ \begin{array}{l} \lim_{|\mathbf{r}| \to \infty} u(\mathbf{r}) = 0 \\ \lim_{|\mathbf{r}| \to \infty} \frac{\partial u(\mathbf{r})}{\partial r} = i K_0 u(\mathbf{r}) \end{array} \right. \qquad (3) $$

To explicitly use the boundary condition defined by eq. (2) in the solution of the problem, we need to transform eq. (1) into an integral equation in terms of boundary and volume integrals of the solution domain. This integral representation for seismic response $u(\mathbf{r})$ enables us to separate the wavefields contributed from different scattering sources and subsequently handle them respectively for different purposes. The scattering sources in the solution domain $\tilde{\Omega}^{(j)}$ include the top boundary $\partial\Omega^{(j-1)}$, the bottom boundary $\partial\Omega^{(j)}$, the seismic source $s(\mathbf{r}, \omega)$ and the volume heterogeneities in $\Omega^{(j)}$. That is, the total seismic response $u(\mathbf{r})$ at a location $\mathbf{r} \in \tilde{\Omega}^{(j)}$ is composed of

$$ u(\mathbf{r}) = u_0(\mathbf{r}) + u_1(\mathbf{r}) + u_2(\mathbf{r}) + u_3(\mathbf{r}). \qquad (4) $$

$u_0(\mathbf{r})$ is the incident field in the background medium and can be represented as

$$ u_0(\mathbf{r}) = \int_{\tilde{\Omega}^{(j)}} s^{(j)}(\mathbf{r}', \omega) G^{(j)}(\mathbf{r}, \mathbf{r}') d\mathbf{r}'. \qquad (5) $$

$u_1(\mathbf{r})$ is the boundary field scattered by the top boundary $\partial\Omega^{(j-1)}$ and satisfies the following boundary integral equation:

$$ u_1(\mathbf{r}) = \int_{\partial\Omega^{(j-1)}} \left[ G^{(j)}(\mathbf{r}, \mathbf{r}') t^{(j-1)}(\mathbf{r}') - u^{(j-1)}(\mathbf{r}') \frac{\partial G^{(j)}(\mathbf{r}, \mathbf{r}')}{\partial n} \right] d\mathbf{r}', \qquad (6) $$

where $\partial/\partial n$ denotes differentiation with respect to the outward normal of the boundary $\partial\Omega^{(j-1)}$ and $t^{(j-1)}(\mathbf{r}) = \partial u^{(j-1)}(\mathbf{r})/\partial n$ the normal gradient of the displacement on $\partial\Omega^{(j-1)}$. Similarly, $u_2(\mathbf{r})$ is the boundary field scattered by the bottom boundary $\partial\Omega^{(j)}$ and satisfies the following boundary integral equation:

$$ u_2(\mathbf{r}) = \int_{\partial\Omega^{(j)}} \left[ G^{(j)}(\mathbf{r}, \mathbf{r}') t^{(j)}(\mathbf{r}') - u^{(j)}(\mathbf{r}') \frac{\partial G^{(j)}(\mathbf{r}, \mathbf{r}')}{\partial n} \right] d\mathbf{r}'. \qquad (7) $$

$u_3(\mathbf{r})$ is the volume field scattered by the volume heterogeneities in $\Omega^{(j)}$ and satisfies the following Lippmann-Schwinger integral equation:

$$ u_3(\mathbf{r}) = [K_0^{(j)}]^2 \int_{\Omega^{(j)}} O^{(j)}(\mathbf{r}') w^{(j)}(\mathbf{r}') G^{(j)}(\mathbf{r}, \mathbf{r}') d\mathbf{r}', \qquad (8) $$

where the relative slowness perturbation $O^{(j)}(\mathbf{r}) = \rho^{(j)}(\mathbf{r})\mu_0^{(j)}/\rho_0^{(j)}\mu^{(j)}(\mathbf{r}) - 1$. $w^{(j)}(\mathbf{r})$ is used to denote the displacement at internal points inside $\Omega^{(j)}$, while $u^{(j-1)}(\mathbf{r})$ and $u^{(j)}(\mathbf{r})$ in eqs (6) and (7) denotes the displacement at boundary points of $\Omega^{(j)}$. The causal Green's function is defined everywhere in the free space, relating an observation point $\mathbf{r}$ to a scattering point $\mathbf{r}'$. It satisfies the homogeneous Helmholtz equation in the reference medium:

$$ \nabla^2 G^{(j)}(\mathbf{r}, \mathbf{r}') + [K_0^{(j)}]^2 G^{(j)}(\mathbf{r}, \mathbf{r}') = -\delta(\mathbf{r} - \mathbf{r}'), \qquad (9) $$

for all $\mathbf{r}, \mathbf{r}' \in \tilde{\Omega}^{(j)}$. For 2-D problems, the Green's function is given by

$$ G^{(j)}(\mathbf{r}, \mathbf{r}') = \frac{i H_0^{(1)}[K_0^{(j)} |\mathbf{r}' - \mathbf{r}|]}{4}, \qquad (10) $$

where $i = \sqrt{-1}$ and $H_0^{(1)}$ is the Hankel function of the first kind and of zeroth order. These integral equations naturally satisfy Sommerfeld non-reflecting and decay boundary conditions defined by eq. (3).

For piecewise homogeneous media [i.e. $O^{(j)}(\mathbf{r}) = 0$], $u_3(\mathbf{r}) = 0$ and eq. (4) reduces to a standard boundary integral equation. In the next section, these integrals in eqs (5) to (8) for piecewise heterogeneous media are formulated by the discrete wavenumber Green's function representation and then assembled into a set of simultaneous matrix equations using the boundary conditions of continuity for displacement and traction across all interfaces. This global matrix presents a recursive relation for wave propagation in piecewise heterogeneous media.

## DISCRETE WAVENUMBER REPRESENTATION

A pioneering work on the discrete wavenumber representation of wavefield can be found in Larner (1970). A discrete wavenumber numerical method was developed to study the effect of topography (Bouchon 1973) and sediment-filled valleys (Bard & Bouchon 1980) on surface motion. To use the discrete wavenumber method to discretize the boundary/volume integrals, we first need to express the Green's function $G^{(j)}(\mathbf{r}, \mathbf{r}')$ in the form of a horizontal wavenumber integral (Morse & Feshbach 1953):

$$ G^{(j)}(\mathbf{r}, \mathbf{r}') = \frac{i}{4\pi} \int_{-\infty}^{+\infty} [\gamma^{(j)}]^{-1} \exp[ik(x - x') + i\gamma^{(j)}|z - z'|] dk, \qquad (11) $$

where $k$ is the horizontal wavenumber and the vertical wavenumber $\gamma^{(j)} = \sqrt{[K_0^{(j)}]^2 - k^2}$. The discretization of the wavenumber integral requires a spatial periodicity in the $x$ direction. We assume that the interface irregularity and volume heterogeneity are localized, resulting in
periodic interfaces and media. We also assume the distribution of sources is periodic. That is,

$$
\left\{
\begin{aligned}
\xi^{(j)}(x) &= \xi^{(j)}(x + nL) \\
O^{(j)}(x, z) &= O^{(j)}(x + nL, z) \\
s^{(j)}(x, z, \omega) &= s^{(j)}(x + nL, z, \omega),
\end{aligned}
\right.
\qquad (12)
$$

where $n = 0, \pm 1, \pm 2, \dots$, and $L$ is the spatial periodicity length in the $x$ direction. The periodicities of both the medium and source will produce the following periodic wavefields:

$$
\left\{
\begin{aligned}
u(x, z) &= u(x + nL, z) \\
t(x, z) &= t(x + nL, z) \\
w(x, z) &= w(x + nL, z),
\end{aligned}
\right.
\qquad (13)
$$

which allows a plane wave expansion of these wavefields (Waterman 1975). Eq. (11) is then replaced by

$$
G^{(j)}(\mathbf{r}, \mathbf{r}') = \frac{i}{2L} \sum_{n=-\infty}^{+\infty} [\gamma_n^{(j)}]^{-1} \exp[-ik_n x' + i\gamma_n^{(j)}|z-z'|] \exp(ik_n x), \quad (14)
$$

where $k_n = 2\pi n/L$ and $\gamma_n^{(j)} = \sqrt{[K_0^{(j)}]^2 - k_n^2}$ with $\text{Im}[\gamma_n^{(j)}] > 0$.

The integrations over the periodic interfaces $[\partial\Omega^{(j-1)}$ and $\partial\Omega^{(j)}]$ and the periodic media ($\Omega^{(j)}$ and $\Omega_s$) can be broken up into a sum of integrals over individual periods. Considering $d\mathbf{r}' = dx' dz'$ for volume integrals, we substitute eqs (11) in (5) to obtain

$$
u_0(x, z) = \frac{i}{4\pi} \int_{-\infty}^{+\infty} \frac{\exp(ikx)dk}{\gamma^{(j)}} \int_{-\infty}^{+\infty} dx' \int_{\xi_{-}^{s}(x')}^{\xi_{+}^{s}(x')} s^{(j)}(x', z') \exp[-ikx' + i\gamma^{(j)}|z-z'|] dz', \quad (15)
$$

where $\xi_+^s$ and $\xi_-^s$ are the top and bottom boundaries of the source medium $\Omega_s^{(j)}$. Using the periodic condition in eq. (12), we have

$$
u_0(x, z) = \sum_{n=-\infty}^{+\infty} \frac{i \exp(ik_n x)}{2L \gamma_n^{(j)}} \int_0^L dx' \int_{\xi_n^{(j)}(x')}^{\xi_n^{(j)}(x')} s^{(j)}(x', z') \exp[-ik_n x' + i\gamma_n^{(j)}|z-z'|] dz'. \quad (16)
$$

For the boundary integrals, we express $\partial G^{(j)}(\mathbf{r}, \mathbf{r}')/\partial n = \mathbf{n} \cdot \nabla G^{(j)}(\mathbf{r}, \mathbf{r}')$ with $\mathbf{n} = (n_x, n_z)$ and $\nabla G^{(j)}(\mathbf{r}, \mathbf{r}') = [\partial G^{(j)}/\partial(x - x'), \partial G^{(j)}/\partial(z - z')]$. For the top boundary $\partial\Omega^{(j-1)}$, we have

$$
\left\{
\begin{aligned}
n_x[x', \xi^{(j-1)}(x')] &= -\xi_x^{(j-1)}(x')/J^{(j-1)}(x') \\
n_z[x', \xi^{(j-1)}(x')] &= 1/J^{(j-1)}(x') \\
d\mathbf{r}' &= J^{(j-1)}(x')dx' \\
J^{(j-1)}(x') &= \sqrt{1 + [\xi_x^{(j-1)}(x')]^2},
\end{aligned}
\right.
\qquad (17)
$$

where $\xi_x^{(j-1)}(x') = \partial\xi^{(j-1)}(x')/\partial x'$. Similarly for the bottom boundary $\partial\Omega^{(j)}$, we have

$$
\left\{
\begin{aligned}
n_x[x', \xi^{(j)}(x')] &= -\xi_x^{(j)}(x')/J^{(j)}(x') \\
n_z[x', \xi^{(j)}(x')] &= 1/J^{(j)}(x') \\
d\mathbf{r}' &= J^{(j)}(x') dx' \\
J^{(j)}(x') &= \sqrt{1 + [\xi_x^{(j)}(x')]^2}.
\end{aligned}
\right.
\qquad (18)
$$

Using the above expressions, eq. (6) can be transformed into

$$
\begin{align}
u_1(x, z) = \int_{-\infty}^{+\infty} \frac{dk}{4\pi \gamma^{(j)}} \int_{-\infty}^{+\infty} & \left[ i J^{(j-1)}(x') t^{(j-1)}(x', \xi^{(j-1)}(x')) - (\xi_x^{(j-1)}(x') k - \gamma^{(j)} \operatorname{sgn}(z - \xi^{(j-1)}(x')) \right] \nonumber \\
& \cdot u^{(j-1)}(x', \xi^{(j-1)}(x')) \exp \left[ ik(x - x') + i\gamma^{(j)} |z - \xi^{(j-1)}(x')| \right] dx'. \tag{19}
\end{align}
$$

Considering the interface periodicity, the above expression reduces to

$$
\begin{equation}
\begin{split}
u_1(x, z) = \sum_{n=-\infty}^{+\infty} \frac{\exp(ik_n x)}{2L\gamma_n^{(j)}} & \int_0^L \left[ iJ^{(j-1)}(x')t^{(j-1)}(x', \xi^{(j-1)}(x')) - (\xi_x^{(j-1)}(x')k_n - \gamma_n^{(j)}\operatorname{sgn}(z-\xi^{(j-1)}(x')) \right] \\
& \cdot u^{(j-1)}(x', \xi^{(j-1)}(x')) \exp[-ik_n x' + i\gamma_n^{(j)}|z-\xi^{(j-1)}(x')|] dx'.
\end{split}
\tag{20}
\end{equation}
$$

Similarly for the bottom boundary $\partial\Omega^{(j)}$, eq. (7) can be rewritten in the form

$$
\begin{equation}
\begin{split}
u_2(x, z) = \sum_{n=-\infty}^{+\infty} \frac{\exp(ik_n x)}{2L\gamma_n^{(j)}} & \int_0^L \left[ i J^{(j)}(x')t^{(j)}(x', \xi^{(j)}(x')) - (\xi_x^{(j)}(x')k_n - \gamma_n^{(j)}\operatorname{sgn}(z-\xi^{(j)}(x')) \right] \\
& \cdot u^{(j)}(x', \xi^{(j)}(x')) \exp[-ik_n x' + i\gamma_n^{(j)}|z-\xi^{(j)}(x')|] dx'.
\end{split}
\tag{21}
\end{equation}
$$
For the volume integral over $\Omega^{(j)}$, eq. (8) becomes

$$u_3(x, z) = \frac{i [K_0^{(j)}]^2}{4\pi} \int_{-\infty}^{+\infty} \frac{\exp(ikx)dk}{\gamma^{(j)}} \int_{-\infty}^{+\infty} dx' \int_{\xi^{(j-1)}(x')}^{\xi^{(j)}(x')} O^{(j)}(x', z') w^{(j)}(x', z') \exp[-ikx' + i\gamma^{(j)}|z - z'|] dz'. \quad (22)$$

Considering the medium periodicity over $\Omega^{(j)}$, the above equation reduces to

$$u_3(x, z) = \sum_{n=-\infty}^{+\infty} \frac{i (K_0^{(j)})^2 \exp(ik_n x)}{2L \gamma_n^{(j)}} \int_0^L dx' \int_{\xi^{(j-1)}(x')}^{\xi^{(j)}(x')} O^{(j)}(x', z') w^{(j)}(x', z') \exp[-ik_n x' + i\gamma_n^{(j)}|z - z'|] dz'. \quad (23)$$

To obtain the total wavefield $u(\mathbf{r})$ using eq. (4), we need to compute the individual scattering components. Using eq. (16), we can directly calculate the source term $u_0(x, z)$ from the known source distribution. However, the terms $u_1(x, z)$, $u_2(x, z)$ and $u_3(x, z)$ cannot be determined as a result of the unknowns in the integrals of eqs (20), (21) and (23). It is necessary to use the boundary conditions for calculation of these unknowns.

## SIMULTANEOUS MATRIX EQUATIONS AND FULL-WAVEFORM SOLUTIONS

We assume that the infinite sum over the horizontal wavenumber could be properly truncated:

$$G^{(j)}(\mathbf{r}, \mathbf{r}') = \frac{i}{2L} \sum_{n=-M}^{+M} (\gamma_n^{(j)})^{-1} \exp[-ik_n x' + i\gamma_n^{(j)}|z-z'|] \exp(ik_n x), \quad (24)$$

where $M$ is an integer large enough to ensure the convergence of the series and accuracy requirement. To avoid numerical singularity in eq. (24) when $\gamma_n^{(j)} \to 0$, the frequency is chosen to be complex. The choice of the imaginary part of the frequency has been discussed in Bouchon & Aki (1977). The entire periodic media, including interfaces, volume medium and source medium, are then discretized in $x$ at equal $\Delta x$ intervals, resulting in an odd number of points $N_x$ along the $x$ direction. We set $L = N_x \Delta x$ and $N_x = 2M + 1$ that defines the sampling of the wavenumber space. That is, the wavenumber summation in eq. (24) is restricted to the intervals (Bouchon et al. 1989):

$$[-\frac{\pi}{\Delta x}, \frac{\pi}{\Delta x}] \text{ or } [-\frac{\pi}{L}N_x, \frac{\pi}{L}N_x].$$

In addition, the discretization interval $\Delta x$ should be selected to ensure an enough number of points per seismic wavelength. Using $\Delta x$ determined above, the volume $\Omega^{(j)}$ is discretized into $N_p^{(j)}$ points and $\Omega_s^{(j)}$ into $N_s^{(j)}$.

We first evaluate the integral over the source domain $\Omega_s^{(j)}$ by discretizing eq. (16):

$$\left\{ \begin{aligned} u_0(x_i, z_i) &= \sum_{l=1}^{N_x^{(j)}} b_{il}^{(j)} s_l^{(j)}, & i = 1, 2, \dots, (2N_x + N_p^{(j)}), \\ b_{il}^{(j)} &= \frac{i \Delta z}{2N_x} \sum_{n=-M}^{M} \frac{\exp(ik_n x_i)}{\gamma_n^{(j)}} \exp[-ik_n x'_i + i\gamma_n^{(j)}|z_i - z'_i|]. \end{aligned} \right. \quad (25)$$

We then evaluate the integral over the top boundary $\partial\Omega^{(j-1)}$ by discretizing eq. (20):

$$\left\{ \begin{aligned} u_1(x_i, z_i) &= \sum_{l=1}^{N_x} [g_{il}^{(j,1)} t_l^{(j-1)} + \bar{h}_{il}^{(j,1)} u_l^{(j-1)}], & i = 1, 2, \dots, (2N_x + N_p^{(j)}), \\ g_{il}^{(j,1)} &= \frac{i J^{(j-1)}(x'_i)}{2N_x} \sum_{n=-M}^{M} \frac{\exp(ik_n x_i)}{\gamma_n^{(j)}} \exp[-ik_n x'_i + i\gamma_n^{(j)}|z_i - \xi^{(j-1)}(x'_i)|], \\ \bar{h}_{il}^{(j,1)} &= \frac{-1}{2N_x} \sum_{n=-M}^{M} \frac{\exp(ik_n x_i)}{\gamma_n^{(j)}} [\xi_n^{(j-1)}(x'_i)k_n - \gamma_n^{(j)}\operatorname{sgn}(z_i - \xi^{(j-1)}(x'_i))], \\ & \cdot \exp[-ik_n x'_i + i\gamma_n^{(j)}|z_i - \xi^{(j-1)}(x'_i)|]. \end{aligned} \right. \quad (26)$$

Similarly for the integral over the bottom boundary $\partial\Omega^{(j)}$, we have

$$\left\{ \begin{aligned} u_2(x_i, z_i) &= \sum_{l=1}^{N_x} [g_{il}^{(j,2)} t_l^{(j)} + \bar{h}_{il}^{(j,2)} u_l^{(j)}], & i = 1, 2, \dots, (2N_x + N_p^{(j)}), \\ g_{il}^{(j,2)} &= \frac{i J^{(j)}(x'_i)}{2N_x} \sum_{n=-M}^{M} \frac{\exp(ik_n x_i)}{\gamma_n^{(j)}} \exp[-ik_n x'_i + i\gamma_n^{(j)}|z_i - \xi^{(j)}(x'_i)|], \\ \bar{h}_{il}^{(j,2)} &= \frac{-1}{2N_x} \sum_{n=-M}^{M} \frac{\exp(ik_n x_i)}{\gamma_n^{(j)}} [\xi_n^{(j)}(x'_i)k_n - \gamma_n^{(j)}\operatorname{sgn}(z_i - \xi^{(j)}(x'_i))], \\ & \cdot \exp[-ik_n x'_i + i\gamma_n^{(j)}|z_i - \xi^{(j)}(x'_i)|]. \end{aligned} \right. \quad (27)$$
For the volume integral over $\Omega^{(j)}$, we can discretize eq. (23) as

$$ egin{cases} u_3(x_i, z_i) = \sum_{l=1}^{N_p^{(j)}} \bar{d}_{il}^{(j)} w_l^{(j)}, & i = 1, 2, \dots, [2N_x + N_p^{(j)}], \\ \bar{d}_{il}^{(j)} = \frac{i \Delta z O_l^{(j)} (K_0^{(j)})^2}{2N_x} \sum_{n=-M}^{M} \frac{\exp(ik_n x_i)}{\gamma_n^{(j)}} \exp[-i k_n x'_i + i \gamma_n^{(j)} |z_i - z'_i|]. \end{cases} \quad (28) $$

Substituting eqs (25)–(28) into eq. (4), we obtain

$$ u(x_i, z_i) = \sum_{l=1}^{N_p^{(j)}} b_{il}^{(j)} s_l^{(j)} + \sum_{l=1}^{N_x} \left[ g_{il}^{(j,1)} t_l^{(j-1)} + \bar{h}_{il}^{(j,1)} u_l^{(j-1)} + g_{il}^{(j,2)} t_l^{(j)} + \bar{h}_{il}^{(j,2)} u_l^{(j)} \right] + \sum_{l=1}^{N_p^{(j)}} \bar{d}_{il}^{(j)} w_l^{(j)}. \quad (29) $$

Letting

$$ egin{cases} h_{il}^{(j,1)} = \bar{h}_{il}^{(j,1)} - \delta_{il}, \\ h_{il}^{(j,2)} = \bar{h}_{il}^{(j,2)} - \delta_{il}, \\ d_{il}^{(j)} = \bar{d}_{il}^{(j)} - \delta_{il}, \end{cases} \quad (30) $$

where $\delta_{il}$ is the Kronecker delta function, eq. (29) for $i = 1, 2, \dots, [2N_x + N_p^{(j)}]$ can be compacted as a matrix equation:

$$ \mathbf{g}^{(j,1)}\mathbf{t}^{(j-1)} + \mathbf{h}^{(j,1)}\mathbf{u}^{(j-1)} + \mathbf{g}^{(j,2)}\mathbf{t}^{(j)} + \mathbf{h}^{(j,2)}\mathbf{u}^{(j)} + \mathbf{d}^{(j)}\mathbf{w}^{(j)} = -\mathbf{s}^{(j)}, \quad (31) $$

where $\mathbf{s}^{(j)}$ is the $[2N_x + N_p^{(j)}] \times [N_s^{(j)}]$ source matrix with its elements calculated from eq. (25), $\mathbf{g}^{(j,1)}$ and $\mathbf{h}^{(j,1)}$ are the $[2N_x + N_p^{(j)}] \times (N_x)$ top-boundary coefficients matrices, $\mathbf{g}^{(j,2)}$ and $\mathbf{h}^{(j,2)}$ are the $[2N_x + N_p^{(j)}] \times (N_x)$ bottom-boundary coefficients matrices and $\mathbf{d}^{(j)}$ is the $[2N_x + N_p^{(j)}] \times [N_p^{(j)}]$ volume coefficients matrix. Assuming that the source is located in the shallowest volume $\Omega^{(1)}$ as is the case in seismic exploration and using the traction-free condition $\mathbf{t}^{(0)} = 0$ on the free surface $\partial\Omega^{(0)}$, eq. (31) for $\Omega^{(1)}$ becomes

$$ \mathbf{h}^{(1,1)}\mathbf{u}^{(0)} + \mathbf{g}^{(1,2)}\mathbf{t}^{(1)} + \mathbf{h}^{(1,2)}\mathbf{u}^{(1)} + \mathbf{d}^{(1)}\mathbf{w}^{(1)} = -\mathbf{s}^{(1)}. \quad (32) $$

We assume that the medium in the lowermost formation $\Omega^{(N+1)}$ is homogeneous bounded by the lowermost interface $\partial\Omega^{(N)}$ and an arc at infinity. In this case, eq. (31) for $\Omega^{(N+1)}$ reduces to

$$ \mathbf{g}^{(N+1,1)}\mathbf{t}^{(N)} + \mathbf{h}^{(N+1,1)}\mathbf{u}^{(N)} = 0. \quad (33) $$

Defining the boundary coefficient matrices and the unknown boundary displacement-traction vector as

$$ egin{cases} A_1^{(j)} = [\mathbf{g}^{(j,1)}, \mathbf{h}^{(j,1)}], \\ A_2^{(j)} = [\mathbf{g}^{(j,2)}, \mathbf{h}^{(j,2)}], \\ \mathbf{Q}^{(j)} = [\mathbf{t}^{(j)}, \mathbf{u}^{(j)}], \end{cases} \quad (34) $$

eq. (28) can be further compacted as

$$ A_1^{(j)} \mathbf{Q}^{(j-1)} + A_2^{(j)} \mathbf{Q}^{(j)} + \mathbf{d}^{(j)} \mathbf{w}^{(j)} = -\mathbf{s}^{(j)}. \quad (35) $$

For each subregion $\Omega^{(j)}$, the matrices $A_1^{(j)}$, $A_2^{(j)}$ and $\mathbf{d}^{(j)}$ are full with the complex coefficients being functions of frequency, material property and geometry. Using the continuity of the displacement-traction vector across interfaces, we can build the following global matrix equation that describes wave propagation in the whole model.

$$ egin{cases} A_1^{(N+1)} \mathbf{Q}^{(N)} = 0, \\ A_1^{(N)} \mathbf{Q}^{(N-1)} + A_2^{(N)} \mathbf{Q}^{(N)} + \mathbf{d}^{(N)} \mathbf{w}^{(N)} = 0, \\ \vdots \\ A_1^{(j)} \mathbf{Q}^{(j-1)} + A_2^{(j)} \mathbf{Q}^{(j)} + \mathbf{d}^{(j)} \mathbf{w}^{(j)} = 0, \\ A_1^{(j-1)} \mathbf{Q}^{(j-2)} + A_2^{(j-1)} \mathbf{Q}^{(j-1)} + \mathbf{d}^{(j-1)} \mathbf{w}^{(j-1)} = 0, \\ \vdots \\ A_1^{(2)} \mathbf{Q}^{(1)} + A_2^{(2)} \mathbf{Q}^{(2)} + \mathbf{d}^{(2)} \mathbf{w}^{(2)} = 0, \\ A_1^{(1)} \mathbf{u}^{(0)} + A_2^{(1)} \mathbf{Q}^{(1)} + \mathbf{d}^{(1)} \mathbf{w}^{(1)} = -\mathbf{s}^{(1)}. \end{cases} \quad (36) $$

We see that these matrix equations are coupled in the manner of a Markovian chain. Solving the linear equation system of eq. (36) results in seismic responses $\mathbf{u}(\mathbf{r})$ for all nodes in the medium. The Gaussian elimination algorithms can be used for small-scale problems. For large-scale problems or models with complex geometry, the resultant total coefficient matrix is sparse and the corresponding equation is better to be solved by (i) an improved block Gaussian elimination algorithm if seismic survey is set at the surface (Fu 2002a), or (ii) conjugate gradient algorithms. The condition number of the system guarantees fast algorithm convergence because the boundary–volume integral formulations often give rise to coefficient matrices such that their magnitudes are always larger for the main diagonal terms. In addition, our numerical calculations are performed in the frequency domain, which can be easily vectorized and parallelized. For some problems of practical
interest, it contains many more non-zero elements and can be too large to solve. We can reduce the numerical burden by Bouchon's sparsity approximation using a threshold criterion to remove very small entries in the coefficient matrices. The approach is based on the significant spatial decay of Green's functions. Using various threshold criteria, Ortiz-Alemán et al. (1998) demonstrated encouraging results with great savings on both computing time and memory requirements.

## APPROXIMATE SOLUTIONS

The efficiency of the seismic modelling presented above relies on:

(i) the maximum frequency to be calculated,

(ii) the total number of nodes that depends on points per wavelength, and

(iii) the sparsity of the coefficient matrices depending on the geometrical complexity of the used model.

The major drawback of the numerical method is the considerable computer time and memory requirements as indicated clearly by eq. (36). For each $\Omega^{(j)}$, we have total $[2N_x + N_p^{(j)}]$ nodes and the same number of independent equations. Solving these equations involves inverting a matrix of the order $[2N_x + N_p^{(j)}] \times [2N_x + N_p^{(j)}]$. Up to several orders of magnitude less computations can be made to the pure full-wave numerical method by a variety of approximations, such as Born, Rytov, Fresnel and others, to the prohibitively expensive boundary and volume integrals in eq. (36).

Some surveys on the Born approximation to topography and irregular interfaces for the boundary-scattering wave can be found in Levander (1990). Kennett (1972) used the Born approximation in the propagator matrix formalism to model smooth irregular surfaces with small height and small slope. Hudson et al. (1973) studied the accuracy of perturbation theory (Gilbert & Knopoff 1960) applied to Rayleigh wave scattering from triangular grooved surfaces. Perturbation theory to simulate the boundary-scattering wave is usually appropriate for slightly rough surfaces where the scattered field is regarded as only slightly altered by the presence of roughness. Fu et al. (2002) carried out numerical experiments of wave scattering from gentle to rugged surfaces for regional phases. They demonstrate that for gentle rough surfaces the Born approximation can model, to a large degree, both the large-scale and small-scale energy components. However, the error in the small-scale component increases rapidly with propagation distance where the energy fluctuations are obviously flattened because of the Born approximation. For rugged surfaces, the Born approximation significantly destroys the large-scale energy components that are measured by a scattering $Q$. However, it can preserve a relative variation in the small-scale energy components for a certain propagation distance. The above review suggests that the large-scale boundary structures control the principal characteristics of wave propagation. Perturbation theory solutions are restricted to smooth roughness and weak-contrast situations. The full-waveform solutions are required to accurately simulate the reflection/transmission across strong-contrast irregular boundaries.

In this paper, the boundary-scattering wave is simulated in a full implicit way, but approximate solutions are pursued for the volume-scattering wave. The approximation is based on the fact that the volume heterogeneities within each geological formation may be relatively smooth spatially at seismic wavelengths. As indicated from theoretical studies (Aki 1973) on scattering in Gaussian random media, traveltime anomalies across seismic arrays require the presence of large-scale velocity fluctuations in the lithosphere. This result has been quantified by numerical experiments (Frankel & Clayton 1986). The small-scale fluctuations are too smooth to produce traveltime anomalies comparable in size to those reported for actual arrays. However, the small-scale heterogeneities in the lithosphere generate the high-frequency seismic coda. Wu & Aki (1985) characterized the lithosphere by two scale lengths (ca 10 and 1 km) based on theoretical analysis of traveltime variations and coda amplitude. Levander & Holliger (1992) used the finite-difference synthetics to characterize small-scale heterogeneity and large-scale velocity structure of the continental crust. These previous studies imply that the volume heterogeneities within a geological formation are smooth and could be accurately modelled using approximate methods. Literature is abound with studies aimed at describing the effect of volume scattering in the crust and upper mantle on traveltime, amplitude and waveform variations. Much of the theoretical work is based on weak-scattering theory (e.g. Chernov 1960; Sato 1982; Wu 1982) considering single scattering. More accurate calculations can be obtained using energy transfer theory considering the effect of multiple scattering (Wu 1985; Sato & Fehler 2000). An efficient alternative that includes multiple scattering and is more flexible in both the time and frequency domains is the Neumann series (or the Born series) originated in the context of classical integral equation theory. Schuster (1985) incorporated the Born series into the boundary integral equation in solving multibody scattering problems. The method is based on perturbing the boundary integral equation matrix into two parts, with one being inverted and the other expressed by the Born series. In this paper, the Born series is used to approximate the volume integral equation matrix, but the boundary integral equation matrix is kept to be inverted. The accuracies of approximate solutions are quantified by applying them to dimensionless frequency responses of heterogeneous alluvial valleys with a vertical incident SH wave.

Eq. (8) is approximated by the Born series. Rewriting eq. (4) for $\mathbf{r} \in \Omega^{(j)}$:

$$w^{(j)}(\mathbf{r}) = u_0^{(j)}(\mathbf{r}) + u_1^{(j)}(\mathbf{r}) + u_2^{(j)}(\mathbf{r}) + [K_0^{(j)}]^2 \int_{\Omega^{(j)}} O^{(j)}(\mathbf{r}')w^{(j)}(\mathbf{r}')G^{(j)}(\mathbf{r}, \mathbf{r}')d\mathbf{r}', \quad (37)$$

where the interaction between individual volume scatterers in $\Omega^{(j)}$ is handled with a fully implicit manner that accounts for arbitrary multi-scatterings between the volume scatterers. The above equation is a non-linear Fredholm integral equation in the sense that the field-dependent scattering sources are involved, i.e. adding scatterers does not linearly add their scattered fields but changes the distribution of the whole field-dependent scattering sources. For layered geological formations, the contribution to $w^{(j)}(\mathbf{r})$, however, comes mainly from the seismic
source and the strong-contrast boundaries [$\partial\Omega^{(j-1)}$ and $\partial\Omega^{(j)}$], which can be expressed as

$$f^{(j)}(\mathbf{r}) = u_0^{(j)}(\mathbf{r}) + u_1^{(j)}(\mathbf{r}) + u_2^{(j)}(\mathbf{r}). \qquad (38)$$

It is the background field calculated from piecewise homogeneous media. The contribution to $w^{(j)}(\mathbf{r})$ from the volume scatterers in $\Omega^{(j)}$ can be regarded as fluctuations in amplitude and phase over the background field. The strength of the volume-scattering field depends on both the interaction between volume heterogeneities determined by the two-point kernel function $G^{(j)}(\mathbf{r}, \mathbf{r}')$ and the scale of the volume heterogeneity expressed by $O^{(j)}(\mathbf{r}')$. Letting $\lambda^{(j)} = [K_0^{(j)}]^2$ and $\kappa^{(j)}(\mathbf{r}, \mathbf{r}') = G^{(j)}(\mathbf{r}, \mathbf{r}')O^{(j)}(\mathbf{r}')$, eq. (37) can be rewritten in the standard form

$$w^{(j)}(\mathbf{r}) = f^{(j)}(\mathbf{r}) + \lambda^{(j)} \int_{\Omega^{(j)}} \kappa^{(j)}(\mathbf{r}, \mathbf{r}') w^{(j)}(\mathbf{r}') d\mathbf{r}'. \qquad (39)$$

As might be deduced from the $1/r$ decay in $G^{(j)}(\mathbf{r}, \mathbf{r}')$, the interaction between volume heterogeneities is only important if the scatterers are fairly close (Domany et al. 1984). Therefore, the minimum aperture theory using the average Fresnel radius can be employed to design an optimal integration region $\Omega_F^{(j)}$, centered at $\mathbf{r}$. Assuming the volume $\Omega_F^{(j)}$ is discretized by $\Delta x$ into $N_F^{(j)}$ points, we generally have $\Omega_F^{(j)} \ll \Omega^{(j)}$. Therefore, this Fresnel-zone approximation to volume integrations can reduce numerical burden by making the coefficient matrix sparser:

$$w^{(j)}(\mathbf{r}) \approx f^{(j)}(\mathbf{r}) + \lambda^{(j)} \int_{\Omega_F^{(j)}} \kappa^{(j)}(\mathbf{r}, \mathbf{r}') w^{(j)}(\mathbf{r}') d\mathbf{r}', \qquad (40)$$

which functions somewhat analogous to the Bouchon's sparsity approximation to boundary integrations as mentioned previously.

We define an integral operator $\mathbf{K}^{(j)}$ by

$$\mathbf{K}^{(j)} w^{(j)}(\mathbf{r}) = \int_{\Omega^{(j)}} \kappa^{(j)}(\mathbf{r}, \mathbf{r}') w^{(j)}(\mathbf{r}') d\mathbf{r}', \quad \mathbf{r} \in \Omega^{(j)}. \qquad (41)$$

It is easy to check $\mathbf{K}^{(j)}(\phi_1 + \phi_2) = \mathbf{K}^{(j)}\phi_1 + \mathbf{K}^{(j)}\phi_2$ and $\mathbf{K}^{(j)}(\lambda, \phi) = \lambda \mathbf{K}^{(j)}\phi$ if $\lambda$ is a constant. That is, $\mathbf{K}^{(j)}$ is a linear mapping transformation between suitable vector spaces of functions. In operator form, eq. (39) can be rewritten

$$w^{(j)} = f^{(j)} + \lambda^{(j)} \mathbf{K}^{(j)} w^{(j)}, \qquad (42)$$

or more standard

$$[\mathbf{I} - \lambda^{(j)} \mathbf{K}^{(j)}] w^{(j)} = f^{(j)}. \qquad (43)$$

$\mathbf{K}^{(j)}w^{(j)}$ is continuous on $\Omega^{(j)}$ as a result of the continuity of the Green's function, so the linear map $\mathbf{I} - \lambda^{(j)}\mathbf{K}^{(j)}$ has an inverse and eq. (43) has the unique solution $w^{(j)} = [\mathbf{I} - \lambda^{(j)}\mathbf{K}^{(j)}]^{-1}f^{(j)}$. If $\|\lambda^{(j)}\mathbf{K}^{(j)}\| < 1$, the iterative solution is possible with eq. (42). We start with eq. (42) and substitute $f^{(j)} + \lambda^{(j)}\mathbf{K}^{(j)}w^{(j)}$ for $w^{(j)}$ on the right side of this equation. We obtain $w^{(j)} = f^{(j)} + \lambda^{(j)}\mathbf{K}^{(j)}f^{(j)} + [\lambda^{(j)}]^2[\mathbf{K}^{(j)}]^2w^{(j)}$ and repeated application yields

$$w^{(j)} = \sum_{n=0}^{N-1} [\lambda^{(j)}]^n [\mathbf{K}^{(j)}]^n f^{(j)} + [\lambda^{(j)}]^N [\mathbf{K}^{(j)}]^N w^{(j)}. \qquad (44)$$

Here we define $[\mathbf{K}^{(j)}]^{n+1} f^{(j)} = \mathbf{K}^{(j)} [\mathbf{K}^{(j)}]^n f^{(j)}$. Eq. (44) will tend to $[\mathbf{I} - \lambda^{(j)}\mathbf{K}^{(j)}]^{-1}f^{(j)}$ provided the series $\sum [\lambda^{(j)}]^n [\mathbf{K}^{(j)}]^n f^{(j)}$ converges and $\|[\lambda^{(j)}]^N[\mathbf{K}^{(j)}]^N w^{(j)}\| \to 0$, both being true under the assumption $\|\lambda^{(j)}\mathbf{K}^{(j)}\| < 1$. The series expansion $\sum [\lambda^{(j)}]^n[\mathbf{K}^{(j)}]^n$ for $[\mathbf{I} - \lambda^{(j)}\mathbf{K}^{(j)}]^{-1}$ is called the Neumann series (or the Born series). Using the Born series approximation, the implicit eq. (43) becomes an explicit summation:

$$w^{(j)} = \sum_{n=0}^{\infty} [\lambda^{(j)}]^n [\mathbf{K}^{(j)}]^n f^{(j)}, \qquad (45)$$

in which the first-order Born approximation ($n=1$) has been used extensively and successfully in both acoustic and elastic wave scattering if the scattered field is comparatively weak (Wu 1982, 1989b). For strong-contrast volume heterogeneities, adding more terms of the Born series would provide a satisfactory approximation solution.

Using eq. (45), the boundary-volume integral equation numerical method reduces to a relatively inexpensive boundary integral equation method. The matrix eq. (35) will be reduced from the order $[2N_x + N_p^{(j)}] \times [2N_x + N_p^{(j)}]$ to $(2N_x) \times (2N_x)$, which is expected with a great saving of computing time and memory by several orders. The mathematical formulations are demonstrated as follows. Using the first-order Born approximation, eq. (4) for $\mathbf{r} \in \partial\Omega^{(j-1)} + \partial\Omega^{(j)}$ can be expressed as

$$u^{(j)}(\mathbf{r}) = u_0^{(j)}(\mathbf{r}) + u_1^{(j)}(\mathbf{r}) + u_2^{(j)}(\mathbf{r}) + [K_0^{(j)}]^2 \int_{\Omega^{(j)}} O^{(j)}(\mathbf{r}')G^{(j)}(\mathbf{r}, \mathbf{r}') [u_0^{(j)}(\mathbf{r}') + u_1^{(j)}(\mathbf{r}') + u_2^{(j)}(\mathbf{r}')] d\mathbf{r}'. \qquad (46)$$

Considering eqs (25) to (27), the integral in the above expression can be discretized as

$$I(\mathbf{r}) = [K_0^{(j)}]^2 \Delta x \Delta z \sum_{k=1}^{N_p^{(j)}} O^{(j)}(\mathbf{r}_k') G^{(j)}(\mathbf{r}, \mathbf{r}_k') \left\{ \sum_{l=1}^{N_s^{(j)}} b_{kl}^{(j)} s_l^{(j)} + \sum_{l=1}^{N_x} \left[ g_{kl}^{(j,1)} t_l^{(j-1)} + \bar{h}_{kl}^{(j,1)} u_l^{(j-1)} + g_{kl}^{(j,2)} t_l^{(j)} + \bar{h}_{kl}^{(j,2)} u_l^{(j)} \right] \right\}. \qquad (47)$$

Using eqs (28) and (31), the above equation can be rewritten in the matrix form:

$$I(\mathbf{r}_i) = \sum_{k=1}^{N_p^{(j)}} \bar{d}_{ik}^{(j)} \left[ \mathbf{s}_k^{(j)} + \mathbf{g}_k^{(j,1)} \mathbf{t}^{(j-1)} + \bar{h}_k^{(j,1)} \mathbf{u}^{(j-1)} + \mathbf{g}_k^{(j,2)} \mathbf{t}^{(j)} + \bar{h}_k^{(j,2)} \mathbf{u}^{(j)} \right]. \qquad (48)$$
Letting

$$
\left\{
\begin{aligned}
\tilde{s}_i^{(j)} &= \sum_{k=1}^{N_p^{(j)}} \tilde{d}_{ik}^{(j)} s_k^{(j)}, \\
\tilde{g}_i^{(j,1)} &= \sum_{k=1}^{N_p^{(j)}} \tilde{d}_{ik}^{(j)} g_k^{(j,1)}, \quad \tilde{g}_i^{(j,2)} = \sum_{k=1}^{N_p^{(j)}} \tilde{d}_{ik}^{(j)} g_k^{(j,2)}, \\
\tilde{h}_i^{(j,1)} &= \sum_{k=1}^{N_p^{(j)}} \tilde{d}_{ik}^{(j)} \tilde{h}_k^{(j,1)}, \quad \tilde{h}_i^{(j,2)} = \sum_{k=1}^{N_p^{(j)}} \tilde{d}_{ik}^{(j)} \tilde{h}_k^{(j,2)},
\end{aligned}
\right.
\qquad (49)
$$

eq. (48) for $i = 1, 2, \dots, 2N_x$ can be compacted as

$$
I = \tilde{\mathbf{s}}^{(j)} + \tilde{\mathbf{g}}^{(j,1)} \mathbf{t}^{(j-1)} + \tilde{\mathbf{h}}^{(j,1)} \mathbf{u}^{(j-1)} + \tilde{\mathbf{g}}^{(j,2)} \mathbf{t}^{(j)} + \tilde{\mathbf{h}}^{(j,2)} \mathbf{u}^{(j)}. \quad (50)
$$

Substituting in eq. (46) and considering eq. (31) for the background field, we have

$$
[\mathbf{g}^{(j,1)} + \tilde{\mathbf{g}}^{(j,1)}] \mathbf{t}^{(j-1)} + [\mathbf{h}^{(j,1)} + \tilde{\mathbf{h}}^{(j,1)}] \mathbf{u}^{(j-1)} + [\mathbf{g}^{(j,2)} + \tilde{\mathbf{g}}^{(j,2)}] \mathbf{t}^{(j)} + [\mathbf{h}^{(j,2)} + \tilde{\mathbf{h}}^{(j,2)}] \mathbf{u}^{(j)} = -\mathbf{s}^{(j)} - \tilde{\mathbf{s}}^{(j)}. \quad (51)
$$

Setting the source vector and the boundary coefficient matrices as

$$
\left\{
\begin{aligned}
\tilde{\mathbf{S}}^{(j)} &= \mathbf{s}^{(j)} + \tilde{\mathbf{s}}^{(j)}, \\
\tilde{\mathbf{A}}_1^{(j)} &= [\mathbf{g}^{(j,1)} + \tilde{\mathbf{g}}^{(j,1)}, \mathbf{h}^{(j,1)} + \tilde{\mathbf{h}}^{(j,1)}], \\
\tilde{\mathbf{A}}_2^{(j)} &= [\mathbf{g}^{(j,2)} + \tilde{\mathbf{g}}^{(j,2)}, \mathbf{h}^{(j,2)} + \tilde{\mathbf{h}}^{(j,2)}],
\end{aligned}
\right.
\qquad (52)
$$

eq. (51) can be further compacted as

$$
\tilde{\mathbf{A}}_1^{(j)} \mathbf{Q}^{(j-1)} + \tilde{\mathbf{A}}_2^{(j)} \mathbf{Q}^{(j)} = -\tilde{\mathbf{S}}^{(j)}. \qquad (53)
$$

Therefore, the boundary–volume matrix equation system of eq. (36) reduces to a boundary matrix system:

$$
\left\{
\begin{array}{l}
\tilde{\mathbf{A}}_1^{(N+1)} \mathbf{Q}^{(N)} = 0, \\
\tilde{\mathbf{A}}_1^{(N)} \mathbf{Q}^{(N-1)} + \tilde{\mathbf{A}}_2^{(N)} \mathbf{Q}^{(N)} = 0, \\
\vdots \\
\tilde{\mathbf{A}}_1^{(i)} \mathbf{Q}^{(i-1)} + \tilde{\mathbf{A}}_2^{(i)} \mathbf{Q}^{(i)} = 0, \\
\tilde{\mathbf{A}}_1^{(i-1)} \mathbf{Q}^{(i-2)} + \tilde{\mathbf{A}}_2^{(i-1)} \mathbf{Q}^{(i-1)} = 0, \\
\vdots \\
\tilde{\mathbf{A}}_1^{(2)} \mathbf{Q}^{(1)} + \tilde{\mathbf{A}}_2^{(2)} \mathbf{Q}^{(2)} = 0, \\
\tilde{\mathbf{A}}_1^{(1)} \mathbf{u}^{(0)} + \tilde{\mathbf{A}}_2^{(1)} \mathbf{Q}^{(1)} = -\tilde{\mathbf{S}}^{(1)}.
\end{array}
\right.
\qquad (54)
$$

The first-order Born approximation considers only single scattering between a boundary point $\mathbf{r} \in \partial\Omega^{(j-1)} + \partial\Omega^{(j)}$ and an internal point $\mathbf{r}' \in \Omega^{(j)}$. It neglects multiple-scattering energy. Similarly for the second-order Born approximation, using eq. (38) we can express eq. (4) for $\mathbf{r} \in \partial\Omega^{(j-1)} + \partial\Omega^{(j)}$ as

$$
u^{(j)}(\mathbf{r}) = f^{(j)}(\mathbf{r}) + [K_0^{(j)}]^2 \int_{\Omega^{(j)}} O^{(j)}(\mathbf{r}') G^{(j)}(\mathbf{r}, \mathbf{r}') f^{(j)}(\mathbf{r}') d\mathbf{r}' \\
+ [K_0^{(j)}]^4 \int_{\Omega^{(j)}} d\mathbf{r}' \int_{\Omega^{(j)}} O^{(j)}(\mathbf{r}') O^{(j)}(\mathbf{r}'') G^{(j)}(\mathbf{r}, \mathbf{r}') G^{(j)}(\mathbf{r}, \mathbf{r}'') f^{(j)}(\mathbf{r}'') d\mathbf{r}''. \quad (55)
$$

The first two terms on the right side of this equation are the first-order Born approximation. The third term can be treated in a similar way:

$$
I(\mathbf{r}) = [K_0^{(j)}]^2 \Delta x \Delta z \sum_{k=1}^{N_p^{(j)}} O^{(j)}(\mathbf{r}_k') G^{(j)}(\mathbf{r}, \mathbf{r}_k') [\tilde{\mathbf{s}}_k^{(j)} + \tilde{\mathbf{g}}_k^{(j,1)} \mathbf{t}^{(j-1)} + \tilde{\mathbf{h}}_k^{(j,1)} \mathbf{u}^{(j-1)} + \tilde{\mathbf{g}}_k^{(j,2)} \mathbf{t}^{(j)} + \tilde{\mathbf{h}}_k^{(j,2)} \mathbf{u}^{(j)}]. \quad (56)
$$

Using eq. (28), the above equation can be further discretized in the form

$$
I(\mathbf{r}_i) = \sum_{k=1}^{N_p^{(j)}} \tilde{d}_{ik}^{(j)} [\tilde{\mathbf{s}}_k^{(j)} + \tilde{\mathbf{g}}_k^{(j,1)} \mathbf{t}^{(j-1)} + \tilde{\mathbf{h}}_k^{(j,1)} \mathbf{u}^{(j-1)} + \tilde{\mathbf{g}}_k^{(j,2)} \mathbf{t}^{(j)} + \tilde{\mathbf{h}}_k^{(j,2)} \mathbf{u}^{(j)}]. \quad (57)
$$
Letting

$$
\left\{
\begin{aligned}
\tilde{\mathbf{s}}_i^{(j)} &= \sum_{k=1}^{N_p^{(j)}} \bar{d}_{ik}^{(j)} \tilde{\mathbf{s}}_k^{(j)}, \\
\tilde{\mathbf{g}}_i^{(j,1)} &= \sum_{k=1}^{N_p^{(j)}} \bar{d}_{ik}^{(j)} \tilde{\mathbf{g}}_k^{(j,1)}, \quad \tilde{\mathbf{g}}_i^{(j,2)} = \sum_{k=1}^{N_p^{(j)}} \bar{d}_{ik}^{(j)} \tilde{\mathbf{g}}_k^{(j,2)}, \\
\tilde{\mathbf{h}}_i^{(j,1)} &= \sum_{k=1}^{N_p^{(j)}} \bar{d}_{ik}^{(j)} \tilde{\mathbf{h}}_k^{(j,1)}, \quad \tilde{\mathbf{h}}_i^{(j,2)} = \sum_{k=1}^{N_p^{(j)}} \bar{d}_{ik}^{(j)} \tilde{\mathbf{h}}_k^{(j,2)},
\end{aligned}
\right.
\qquad (58)
$$

eq. (57) for $i = 1, 2, \dots, 2N_x$ can be compacted as

$$
I = \tilde{\mathbf{s}}^{(j)} + \tilde{\mathbf{g}}^{(j,1)} t^{(j-1)} + \tilde{\mathbf{h}}^{(j,1)} \mathbf{u}^{(j-1)} + \tilde{\mathbf{g}}^{(j,2)} t^{(j)} + \tilde{\mathbf{h}}^{(j,2)} \mathbf{u}^{(j)}. \quad (59)
$$

Substituting in eq. (55) and letting the boundary coefficient matrices and the source matrix be

$$
\left\{
\begin{aligned}
\tilde{\mathbf{S}}^{(j)} &= \mathbf{s}^{(j)} + \tilde{\mathbf{s}}^{(j)} + \tilde{\mathbf{s}}^{(j)}, \\
\tilde{\mathbf{A}}_1^{(j)} &= \left[ \mathbf{g}^{(j,1)} + \tilde{\mathbf{g}}^{(j,1)} + \tilde{\mathbf{g}}^{(j,1)}, \mathbf{h}^{(j,1)} + \tilde{\mathbf{h}}^{(j,1)} + \tilde{\mathbf{h}}^{(j,1)} \right], \\
\tilde{\mathbf{A}}_2^{(j)} &= \left[ \mathbf{g}^{(j,2)} + \tilde{\mathbf{g}}^{(j,2)} + \tilde{\mathbf{g}}^{(j,2)}, \mathbf{h}^{(j,2)} + \tilde{\mathbf{h}}^{(j,2)} + \tilde{\mathbf{h}}^{(j,2)} \right],
\end{aligned}
\right.
\qquad (60)
$$

we obtain a concise matrix equation for $\Omega^{(j)}$:

$$
\tilde{\mathbf{A}}_1^{(j)} \mathbf{Q}^{(j-1)} + \tilde{\mathbf{A}}_2^{(j)} \mathbf{Q}^{(j)} = -\tilde{\mathbf{S}}^{(j)}. \qquad (61)
$$

Therefore, the global matrix equation system for the second-order Born approximation is

$$
\left\{
\begin{array}{l}
\tilde{\tilde{\mathbf{A}}}_1^{(N+1)} \mathbf{Q}^{(N)} = 0, \\
\tilde{\tilde{\mathbf{A}}}_1^{(N)} \mathbf{Q}^{(N-1)} + \tilde{\tilde{\mathbf{A}}}_2^{(N)} \mathbf{Q}^{(N)} = 0, \\
\vdots \\
\tilde{\tilde{\mathbf{A}}}_1^{(i)} \mathbf{Q}^{(i-1)} + \tilde{\tilde{\mathbf{A}}}_2^{(i)} \mathbf{Q}^{(i)} = 0, \\
\tilde{\tilde{\mathbf{A}}}_1^{(i-1)} \mathbf{Q}^{(i-2)} + \tilde{\tilde{\mathbf{A}}}_2^{(i-1)} \mathbf{Q}^{(i-1)} = 0, \\
\vdots \\
\tilde{\tilde{\mathbf{A}}}_1^{(2)} \mathbf{Q}^{(1)} + \tilde{\tilde{\mathbf{A}}}_2^{(2)} \mathbf{Q}^{(2)} = 0, \\
\tilde{\tilde{\mathbf{A}}}_1^{(1)} \mathbf{u}^{(0)} + \tilde{\tilde{\mathbf{A}}}_2^{(1)} \mathbf{Q}^{(1)} = -\tilde{\mathbf{S}}^{(1)}.
\end{array}
\right.
\qquad (62)
$$

The second-order Born approximation considers two scatterings, one between a boundary point $\mathbf{r} \in \partial\Omega^{(j-1)} + \partial\Omega^{(j)}$ and an internal point $\mathbf{r}' \in \Omega^{(j)}$, the other between two internal points $\mathbf{r}' \in \Omega^{(j)}$ and $\mathbf{r}'' \in \Omega^{(j)}$. We see that the generalization of the scheme to multi-order Born approximation is straightforward for multiple scatterings.

NUMERICAL TESTS USING SEMICIRCULAR CANYON AND ALLUVIAL VALLEYS

The proposed numerical method and the implementing computation program are tested by modelling a semicircular canyon topography, a semicircular homogeneous valley and a semicircular heterogeneous valley. Previously published results for these typical topography structures are used for comparison. The semicircular canyon topography of radius *a* is shown in Fig. 2(a) in an elastic homogeneous half-space. The frequency response distribution along the horizontal surface and the canyon topography to vertical incident *SH* wave is shown in Fig. 2(b) for the dimensionless frequency *η* = 1.0 and compared with those by Sánchez-Sesma & Campillo (1991) and Fu & Wu (2001). The latter has been verified by Trifunac's analytical solution (Trifunac 1973). The dimensionless frequency is defined as *η* = 2*a*/λ = *a*ω/π*β*, where *ω* is the circular frequency and *β* is the velocity of shear-wave propagation. We see an excellent agreement between these three results, particularly around the two sharp edges at *x* = ±*a*. Some minor departures between the solid line and the dotted line are observed along the surface outside the canyon possibly because of the imaginary part of the frequency introduced to the discrete wavenumber Green's function computation. This artificial damping to suppress the fictitious waves affects the dimensionless frequency responses considerably. It seems to indicate the difficulty of the same amount of imaginary part of the frequency to ensure accurate amplitude responses for both parts of the horizontal surface and the semicircular canyon. The time-domain responses of the semicircular canyon topography to vertical incident *SH* wave are shown in Fig. 2(c) in comparison to Fig. 2(d) computed by Fu (1996). The remarkable agreement can be seen between the two solutions. Synthetic seismograms with the time axis shifted properly are calculated using a Ricker wavelet in the frequency range of *ca* 0.0–4.0 Hz.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison with other methods for a semicircular canyon topography with vertical incident *SH* wave. (a) Semicircular canyon topography of radius *a* and the surrounding homogeneous half-space with the density *ρ* and shear wave velocity *β*. (b) Frequency responses by our method (dotted line), Fu & Wu (2001) (solid line) and Sánchez-Sesma & Campillo (1991) (open circles) for the dimensionless frequency η = 1.0. (c) Time responses along the topography in the frequency range of 0–4 Hz. (d) Time responses by Fu (1996).</figcaption></figure>

Observation points are set along the surface both inside and outside the canyon and the shear wave velocity is 2000 m s⁻¹. Although the semicircular canyon is too simple for realistic topography, the reflection and diffraction along the canyon irregularity are quite complex in the time domain. We see that the amplitude distribution of the seismograms in the time domain is consistent with that in the frequency domain (Fig. 2b). However, the destructive/constructive interferences between the direct, reflected and diffracted waves are much more clearly seen in the time domain, particularly near the edges and inside the canyon. In summary, the comparison in Fig. 2 shows the validity of the present method and computation program.

The seismic response of alluvial valleys has been extensively studied in the seismological and engineering perspectives to gain physical insight into ground motion behaviour of alluvial valleys. For instance, Bard & Bouchon (1985) demonstrated the existence of resonance patterns from the numerical results on seismic responses of symmetrical and homogeneous sediment-filled valleys. Sánchez-Sesma (1983) studied the elastic response of 3-D axisymmetric canyons and valleys. Pedersen et al. (1995) modelled azimuth dependence of wave amplification in alluvial valleys. However, numerical modelling of the site amplification has often been found to systematically underestimate the actual amplification observed in the field, implying that modelling should consider the more complex nature of the earth. Our primary concern here is how much the localized volume heterogeneities inside an alluvial valley affect wave amplification. We first calculate the frequency and time responses of a semicircular homogeneous valley to verify our computation program and also to provide a comparison for various heterogeneous valleys. Fig. 3(a) shows a semicircular homogeneous valley of radius *a*, with the density $\bar{\rho}$ and the shear wave velocity $\bar{\beta}$, in the surrounding homogeneous half-space with the density *ρ* and shear wave velocity *β*. The frequency response distribution along the horizontal topography to vertical incident *SH* wave is shown in Fig. 3(b) for the dimensionless frequency η = 1.0 and compared well to that by Fu (2002a). The latter has been verified by Trifunac's analytical solution (Trifunac 1971). Small discrepancies, however, are observed along the surface outside the valley. Again, this may indicate the difficulty of the choice of the imaginary part of the frequency for both inside and outside the valley. The time-domain response of the semicircular homogeneous valley to vertical incident *SH* wave is shown in Fig. 3(c) in good comparison to Fig. 3(d) computed by Fu (2002a). Synthetic seismograms are calculated in the frequency range of *ca* 0.0–4.0 Hz, with observation points setting along the surface. The shear wave velocity inside the valley is 1500 m s⁻¹ in the background medium of 3000 m s⁻¹.

To improve the quantitative account of site responses, the modelling is extended to heterogeneous valleys by perturbing the velocity inside the valley randomly by *ca* 5–20 per cent. For example, the 20 per cent perturbation of the 3000 m s⁻¹ reference velocity indicates that the valley velocities will change randomly from 2400 to 3600 m s⁻¹. The top panel in Fig.4 shows the semicircular heterogeneous valley of radius *a*, with the reference density $\bar{\rho}_0$, the reference shear wave velocity $\bar{\beta}_0$ and the velocity perturbation $\delta$, in the surrounding homogeneous half-space. First the valley velocity is perturbed by 5 and 10 per cent and frequency responses are computed along the valley surface resulting from a vertical incident *SH* wave for the dimensionless frequency η = 1.0. The results are presented in the middle panel of Fig. 4. We see
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison with Fu (2002a) for a semicircular homogeneous valley with vertical incident *SH* wave. (a) Semicircular valley with the radius *a*, the density *ρ̄* and the shear wave velocity *β̄* in the surrounding homogeneous half-space with the density *ρ* and shear wave velocity *β*. (b) Frequency responses by our method (dotted line) and Fu (2002a) (solid line) for the dimensionless frequency *η* = 1.0. (c) Time responses along the topography in the frequency range of 0–4 Hz. (d) Time responses by Fu (2002a).</figcaption></figure>

relatively moderate amplitude fluctuations along the reference response curve, mainly occurring at the two edges and in the middle of the valley. Remarkable amplitude fluctuations can be expected when the velocity perturbation increases to 15 and 20 per cent, shown in the bottom panel. We see several times of amplification as a result of heterogeneities introduced in the homogeneous valley. The *ca* 15–20 per cent velocity perturbations are not rare in most complex near-surface areas.

It is well known that the first-order Born approximation is theoretically restricted to weak heterogeneities. The small-perturbation assumption and its associated effects on wave phenomena have been extensively studied, leading to some qualitative and semi-quantitative results. For example, velocity perturbation should be less than 10 per cent. This result can be verified exactly by our method through dimensionless frequency responses of heterogeneous valleys to incident *SH* plane wave. The valley velocity is perturbed by 10, 15 and 20 per cent and then the first-order Born approximation (eq. 54) is used to compute frequency responses along the valley surface resulting from vertical incident *SH* wave for the dimensionless frequency *η* = 1.0. The results are presented respectively in the top, middle and bottom panels of Fig. 5. The remarkable agreement between the accurate solution (solid line, from eq. 36) and the first-order Born approximation solution (dotted line, from eq. 54) at *δ* = 10 per cent confirms the applicable condition of the first-order Born approximation. We also see that the first-order Born approximation solution might be applicable up to *δ* = 15 per cent for some applications, but significantly underestimates amplitude responses at *δ* = 20 per cent. It must be stressed that in this computation the full-waveform implicit solution is used for the boundary-scattering waves inside and around the valley. The approximation is made only at points inside the valley for the volume-scattering wave and compared with the full-waveform implicit solutions for both the boundary-scattering and volume-scattering waves. The first-order Born approximation accounts for only single scattering between a boundary point and an internal point. It neglects multiple scattering between boundary and internal points and also neglects any scattering between internal points. It is worth a mention that Fu (2002b) formulated the first-order Born approximation into the Born dispersion equation in order to analytically evaluate the small-perturbation assumption. The Born dispersion equation can be modified more accurately for one-way wave propagation.

The second-order Born approximation considers single scattering between boundary points and internal points and single scattering between two internal points. Its accuracy can be studied through dimensionless frequency responses of heterogeneous valleys to incident *SH* plane wave. Fig. 6 presents a comparison between the full-waveform implicit solutions (solid lines, from eq. 36) for both the boundary-scattering and volume-scattering waves and the second-order Born approximation solution (dotted lines, from eq. 62) for the volume-scattering wave, respectively for velocity perturbations of 15 (top panel), 20 (middle panel) and 25 per cent (bottom panel). We see that the second-order Born approximation solution is valid strictly for velocity perturbation of less than 15 per cent and approximately applicable up to 20 per cent for some applications. It significantly overestimates amplitude responses at *δ* = 25 per cent possibly as a result of single scattering approximation between two internal points.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Frequency responses resulting from vertical incident SH wave with $\eta = 1.0$ to a semicircular heterogeneous valley with various velocity perturbations. The top panel shows the semicircular heterogeneous valley with the radius $a$, the reference density $\bar{\rho}_0$, the reference shear-wave velocity $\bar{\beta}_0$ and the velocity perturbation $\delta$ in the surrounding homogeneous half-space with the density $\rho$ and the shear wave velocity $\beta$. The middle panel shows frequency responses for random velocity perturbations of 0 (solid line), 5 (dash line) and 10 per cent (dotted line), respectively. The bottom panel shows frequency responses for random velocity perturbations of 0 (solid line), 15 (dash line) and 20 per cent (dotted line), respectively.</figcaption></figure>

DISCUSSIONS AND CONCLUSIONS

A new method for solving the problem of wave excitation and propagation in the multilayered elastic media with volume heterogeneities separated by irregular interfaces has been presented. The method combines the boundary–volume integral equation with the discrete wavenumber Green’s function representation to avoid the well-known singular problem of the Green’s function in the integral equation numerical techniques. The incident, boundary-scattering and volume-scattering waves are separately formulated in order to model different parts of the media with different accuracies. These waves are accurately superposed by the generalized Lippmann–Schwinger integral equation. The earth is assumed to present a system of piecewise heterogeneous media where the large-scale geometric structure controls the main features of wave propagation, whereas the volume heterogeneity in each geological formation varies at different scales, depending on the incident wavelength. Some flexible approaches have been developed in the seismic modelling scheme used here, with a great saving of computing time and memory. The full-waveform boundary method is used to accurately simulate the reflection/transmission across strong-contrast boundaries. Meanwhile for volume heterogeneities, the modelling method is designed with the flexibility for numerical wave propagation at scales of approximate accuracies. The media-oriented solutions are classified in the following four cases:

(i) Solution implicitly for both the boundary-scattering and volume-scattering waves is obtained for each $\Omega^{(j)}$ by

$$A_1^{(j)}Q^{(j-1)} + A_2^{(j)}Q^{(j)} + d^{(j)}w^{(j)} = -s^{(j)}, \qquad (63)$$

where $A_1^{(j)}$ and $A_2^{(j)}$ are the $[2N_x + N_p^{(j)}] \times N_x$ boundary coefficient matrices, respectively, and $d^{(j)}$ is the $[2N_x + N_p^{(j)}] \times N_p^{(j)}$ volume coefficient
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison between the accurate solution (solid lines) of eq. (36) and the first-order Born approximate solution (dotted lines) of eq. (54) for random velocity perturbations of 10 (top panel), 15 (middle panel) and 20 per cent (bottom panel), respectively. Frequency responses are computed resulting from vertical incident SH wave with $\eta = 1.0$ to the semicircular heterogeneous valley.</figcaption></figure>

matrix. Using boundary conditions, we have a $[2N_x + N_p^{(j)}] \times [2N_x + N_p^{(j)}]$ matrix to be inverted. The solution is of high accuracy to model subtle effects of both boundary and volume scatterings, depending on the sampling to the model media that is relative to the scale lengths of the media and the wavelength of incident waves. Numerical modelling for dimensionless frequency responses of heterogeneous valleys to incident SH plane wave indicates that several times of amplification can be expected as a result of heterogeneities introduced in a homogeneous valley.

(ii) Solution implicitly for the boundary-scattering waves but semi-explicitly for the volume-scattering waves is obtained for each $\Omega^{(j)}$ by
$$A_1^{(j)}Q^{(j-1)} + A_2^{(j)}Q^{(j)} + d_F^{(j)}w^{(j)} = -\tilde{s}^{(j)}, \quad (64)$$

where $d_F^{(j)}$ is the $[2N_x + N_p^{(j)}] \times N_p^{(j)}$ volume coefficient matrix but with only $N_F^{(j)}$ non-zero elements in each row, calculated using eq. (40). That is, the total coefficient matrix for $\Omega^{(j)}$ is sparse rather than full in eq. (63).

(iii) Solution implicitly for the boundary-scattering waves but explicitly for the volume-scattering waves using the first-order Born approximation is obtained for each $\Omega^{(j)}$ by

$$\tilde{A}_1^{(j)}Q^{(j-1)} + \tilde{A}_2^{(j)}Q^{(j)} = -\tilde{S}^{(j)}, \quad (65)$$

where $\tilde{A}_1^{(j)}$ and $\tilde{A}_2^{(j)}$ are the $2N_x \times N_x$ boundary coefficient matrices calculated by eq. (52). That is, the total coefficient matrix for $\Omega^{(j)}$ reduces to the order $2N_x \times 2N_x$. Numerical tests with dimensionless frequency responses of heterogeneous valleys confirm that the first-order Born approximation to the volume-scattering waves is strictly valid for velocity perturbation of less than 10 per cent and approximately used up to 15 per cent for general applications.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparison between the accurate solution (solid lines) of eq. (36) and the second-order Born approximate solution (dotted lines) of eq. (62) for random velocity perturbations of 15 (top panel), 20 (middle panel) and 25 per cent (bottom panel), respectively. Frequency responses are computed resulting from vertical incident SH wave with η = 1.0 to the semicircular heterogeneous valley.</figcaption></figure>

(iv) Solution implicitly for the boundary-scattering waves but explicitly for the volume-scattering waves using the second-order Born approximation is obtained for each $\Omega^{(j)}$ by

$$ \tilde{\mathbf{A}}_1^{(j)} \mathbf{Q}^{(j-1)} + \tilde{\mathbf{A}}_2^{(j)} \mathbf{Q}^{(j)} = -\tilde{\mathbf{S}}^{(j)}, \qquad (66) $$

where $\tilde{\mathbf{A}}_1^{(j)}$ and $\tilde{\mathbf{A}}_2^{(j)}$ are the $2N_x \times N_x$ boundary coefficient matrices calculated by eq. (60). That is, the total coefficient matrix for $\Omega^{(j)}$ reduces to the order $2N_x \times 2N_x$. Numerical tests with dimensionless frequency responses of heterogeneous valleys confirm that the second-order Born approximation to the volume-scattering waves is strictly valid for velocity perturbation of less than 15 per cent and approximately used up to 20 per cent for general applications.

The above solutions with various scales of approximation accuracies are expected to cover various complex media with piecewise heterogeneities. Generally, no method is shown to be the best in all possible modelling problems. The method used here might be an optimal modelling scheme in the sense that the designed modelling accuracy can be always satisfied with a balance between computational cost, wavelengths and heterogeneities. Two points must be stressed here:

(i) seismic modelling methods should be flexible at various scales in order to achieve an optimization among computer time, memory requirements, computational wavelengths, medium heterogeneities and modelling accuracies;

(ii) seismic numerical modellings do not only produce final synthetic seismograms but more importantly provide an ability to investigate subtle physical effects of different parts (or different types) of the subsurface media on wave propagation like attenuation and dispersion.
This requires that numerical methods be able to separate wavefields and handle them individually. Therefore, developments of the flexible numerical modelling schemes have crucial importance in future seismological research. In summary, we can say that the seismic modelling scheme presented in this paper is a promising candidate for such targets.

ACKNOWLEDGMENTS

LYF is grateful for helpful discussions of Ru-Shan Wu and Francisco J. Sánchez-Sesma. Francisco J. Sánchez-Sesma provided the IBEM SH-wave solution for comparison in Fig. 2. Part of this material is based on the research carried out by LYF whilst working at the Institute of Tectonics, University of California in Santa Cruz. Valuable comments by Associate Editor, Malcolm Cambridge and two anonymous reviewers have significantly improved this paper. This research was supported by Funds of a Hundred Outstanding Talents Plan, Chinese Academy of Sciences.

REFERENCES

Aki, K., 1973. Scattering of *P* waves under the Montana LASA, *J. geophys. Res.*, **78**, 1334–1346.

Aki, K., 1980. Scattering and attenuation of shear waves in the lithosphere, *J. geophys. Res.*, **85**, 6496–6504.

Bard, P.Y. & Bouchon, M., 1980. Seismic response of sediment-filled valleys. Part 1. The case of incident SH waves, *Bull. seism. Soc. Am.*, **70**, 1263–1286.

Bard, P.Y. & Bouchon, M., 1985. The two-dimensional resonance of sediment-filled valleys, *Bull. seism. Soc. Am.*, **75**, 519–541.

Bouchon, M., 1973. Effect of topography on surface motion, *Bull. seism. Soc. Am.*, **63**, 615–632.

Bouchon, M., 1982. The complete synthesis of seismic crustal phases at regional distances, *J. geophys. Res.*, **82**, 1735–1741.

Bouchon, M. & Aki, K., 1977. Discrete wavenumber representation of seismic source wave fields, *Bull. seism. Soc. Am.*, **67**, 259–277.

Bouchon, M., Schultz, C.A. & Toksöz, M.N., 1996. Effect of three-dimensional topography on seismic motion, *J. geophys. Res.*, **101**, 5835–5846.

Bouchon, M., Campillo, M. & Gaffet, S., 1989. A boundary integral equation-discrete wavenumber representation method to study wave propagation in multilayered media having irregular interfaces, *Geophysics*, **54**, 1134–1140.

Campillo, M., 1987a. Lg wave propagation in a laterally varying crust and the distribution of the apparent quality factor in Central France, *J. geophys. Res.*, **92**, 12 604–12 614.

Campillo, M., 1987b. Modelling of SH-wave propagation in an irregularly layered medium—Application to seismic profiles near a dome, *Geophys. Prospect.*, **35**, 236–249.

Campillo, M., 1990. Propagation and attenuation characteristics of the crustal phase Lg, *Pure appl. Geophys.*, **132**, 1–19.

Campillo, M. & Bouchon, M., 1985. Synthetic SH-seismograms in a laterally varying medium by the discrete wavenumber method, *Geophys. J. R. astr. Soc.*, **83**, 307–317.

Campillo, M. & Paul, A., 1992. Influence of the lower crustal structure on the earth coda of regional seismograms, *J. geophys. Res.*, **97**, 3405–3416.

Chen, X.F., 1990. Seismogram synthesis for multilayered media with irregular interfaces by global generalized reflection/transmission matrices method. I. Theory of two-dimensional SH case, *Bull. seism. Soc. Am.*, **80**, 1696–1724.

Chen, X.F., 1995. Seismogram synthesis for multi-layered media with irregular interfaces by the global generalized reflection/transmission matrices method—Part II. Applications of 2-D SH case, *Bull. seism. Soc. Am.*, **85**, 1094–1106.

Chernov, L.A., 1960. *Wave propagation in a random medium*, McGraw-Hill, New York.

Domany, E., Entin-Wohlman, O. & Mizrachi, L., 1984. Multiple scattering formalism: Application to scattering by two spheres, *J. appl. Phys.*, **56**, 132–136

Dravinksi, M., 1983. Scattering of plane harmonic SH waves by dipping layers of arbitrary shape, *Bull. seism. Soc. Am.*, **73**, 1309–1319.

Frankel, A. & Clayton, R.W., 1986. Finite-difference simulations of seismic scattering: Implications for propagation of short-period seismic waves in the crust and models of crustal heterogeneity, *J. geophys. Res.*, **91**, 6465–6489.

Fu, L.Y., 1996. 3-D boundary element seismic modelling in complex geology. In *Expanded Abstracts of the Technical Program, SEG 66th Annual Meeting*, Denver, CO, USA, pp. 1239–1242.

Fu, L.Y., 2002a. Seismogram synthesis for piecewise heterogeneous media, *Geophys. J. Int.*, **150**, 800–808.

Fu, L.Y., 2002b. Born dispersion equation and Kirchhoff migration in laterally heterogeneous media. In *Expanded Abstracts of the Technical Program, SEG 72nd Annual Meeting*, Salt Lake City, UT, USA, pp. 895–898.

Fu, L.Y., 2003. Numerical study of generalized Lipmann–Schwinger integral equation including surface topography, *Geophysics*, **68**, 665–671.

Fu, L.Y. & Wu, R.S., 2001. A hybrid BE-GS method for modelling regional wave propagation, *Pure appl. Geophys.*, **158**, 1251–1277.

Fu, L.Y., Mu, Y.G. & Yang, H.J., 1997. Forward problem of nonlinear Fredholm integral equation in reference medium via velocity-weighted wavefield function, *Geophysics*, **62**, 650–656.

Fu, L.Y., Wu, R.S. & Campillo, M., 2002. Energy partition and attenuation of regional phases by random free surface, *Bull. seism. Soc. Am.*, **92**, 1992–2007.

Gaffet, S. & Bouchon, M., 1989. Effects of two-dimensional topographies using the discrete wavenumber-boundary integral equation method in *P-SV* cases, *J. acoust. Soc. Am.*, **85**, 2277–2283.

Gibson, R.L., J, Jr & Campillo, M., 1994. Numerical simulation of high- and low-frequency Lg-wave propagation, *Geophys. J. Int.*, **118**, 47–56.

Gilbert, F. & Knopoff, L., 1960. Seismic scattering from topographic irregularities, *J. geophys. Res.*, **65**, 3437–3444.

Hudson, J.A., Humphryes, R.F., Mason, I.M. & Kembhavi, V.K., 1973. The scattering of longitudinal elastic waves at a rough free surface, *J. Phys. D. Appl. Phys.*, **6**, 2174–2186.

Kawase, H. & Aki, K., 1989. A study on the response of a soft basin for incident S, P and Rayleigh waves, with special reference to the long duration observed in Mexico City, *Bull. seism. Soc. Am.*, **79**, 1361–1382.

Kennett, B.L.N., 1972. Seismic waves in laterally varying media, *Geophys. J. R. astr. Soc.*, **27**, 310–325.

Kennett, B.L.N., 1983. *Seismic wave propagation in stratified media*, Cambridge University Press, Cambridge.

Kennett, B.L.N., 1984. Guided wave propagation in laterally varying media-I: Theoretical development, *Geophys. J. R. astr. Soc.*, **79**, 235–255.

Kennett, B.L.N., 1986. Lg waves and structural boundaries, *Bull. seism. Soc. Am.*, **76**, 1133–1141.

Kennett, B.L.N., 1998. Guided waves in three-dimensional structures, *Geophys. J. Int.*, **133**, 159–174.

Larner, K.L., 1970. Near-receiver scattering of teleseismic body waves in layered crust-mantle models having irregular interfaces, *PhD thesis*, Mass. Inst. of Tech, Cambridge, MA, USA.

Lee, J.J. & Langston, C.A., 1983. Wave propagation in a three-dimensional circular basin, *Bull. seism. Soc. Am.*, **61**, 1637–1653.

Levander, A.R., 1990. Seismic scattering near the Earth's surface, *Pure appl. Geophys.*, **132**, 21–47.
Levander, A.R. & Holliger, K., 1992. Small-scale heterogeneity and large-scale velocity structure of the continental crust, *J. geophys. Res.*, **97**, 8797–8804.

Luco, J.E. & Apsel, R.J., 1983. On the Green's functions for a layered half-space. Part I, *Bull. seism. Soc. Am.*, **73**, 909–929.

Maupin, V., 1989. Numerical modelling of Lg wave propagation across the North Sea central graben, *Geophys. J. Int.*, **99**, 273–283.

Maupin, V. & Kennett, B.L.N., 1987. On the use of truncated model expansion in laterally varying media, *Geophys. J. R. astr. Soc.*, **91**, 837–851.

Morse, P.M. & Feshbach, H., 1953. *Methods of Theoretical Physics*, McGraw-Hill, New York.

Ortiz-Alemán, C., Sánchez-Sesma, F.J., Rodriguez-Zuniga, J.L. & Luzon, F., 1998. Computing topographical 3D site effects using a fast IBEM/conjugate gradient approach, *Bull. seism. Soc. Am.*, **88**, 393–399.

Pedersen, H.A., Campillo, M. & Sánchez-Sesma, F.J., 1995. Azimuth dependent wave amplification in alluvial valleys, *Soil Dyn. and Earthq. Eng.*, **14**, 289–300.

Sato, H., 1982. Attenuation of S waves in the lithosphere due to scattering by its random velocity structure, *J. geophys. Res.*, **87**, 7779–7785.

Sato, H. & Fehler, M.C., 2000. *Seismic wave propagation and scattering in the heterogeneous earth*, Springer, New York.

Schuster, G.T., 1985. A hybrid BIE+Born series modelling scheme: Generalized Born series, *J. acoust. Soc. Am.*, **77**, 865–879.

Sánchez-Sesma, F.J., 1983. Diffraction of elastic waves by three-dimensional surface irregularities, *Bull. seism. Soc. Am.*, **73**, 1621–1636.

Sánchez-Sesma, F.J. & Campillo, M., 1991. Diffraction of P, SV and Rayleigh waves by topographic features: a boundary integral formulation, *Bull. seism. Soc. Am.*, **81**, 2234–2253.

Sánchez-Sesma, F.J. & Esquivel, J., 1979. Ground motion on alluvial valleys under incident plane SH waves, *Bull. seism. Soc. Am.*, **69**, 1107–1120.

Trifunac, M.D., 1971. Surface motion of a semi-cylindrical alluvial valley for incident plane SH waves, *Bull. seism. Soc. Am.*, **61**, 1755–1770.

Trifunac, M.D., 1973. Scattering of plane SH waves by a semi-cylindrical canyon, *Earthquake Eng. and Struct. Dyn.*, **1**, 267–281.

Waterman, P.C., 1975. Scattering by periodic surfaces, *J. acoust. Soc. Am.*, **57**, 791–802.

Wong, H.L. & Jennings, P.C., 1975. Effect of canyon topographies on strong ground motion, *Bull. seism. Soc. Am.*, **65**, 1239–1257.

Wu, R.S., 1982. Mean field attenuation and amplitude attenuation due to wave scattering, *Wave Motion*, **4**, 305–316.

Wu, R.S., 1985. Multiple scattering and energy transfer of seismic waves—Separation of scattering effect from intrinsic attenuation—I. Theoretical modelling, *Geophys. J. R. astr. Soc.*, **82**, 57–80.

Wu, R.S., 1989a. Seismic wave scattering, in, *Encyclopedia of Geophysics*, pp. 1166–1187, ed. James, D.E., Van Nostrand Reinhold Co., New York.

Wu, R.S., 1989b. The perturbation method in elastic wave scattering, *Pageoph*, **131**, 605–637.

Wu, R.S., 1996. Synthetic seismogram in heterogeneous media by one-return approximation, *Pure appl. Geophys.*, **148**, 155–173.

Wu, R.S. & Aki, K., 1985. Elastic wave scattering by a random medium and the small-scale inhomogeneities in the Lithosphere, *J. geophys. Res.*, **90**, 10261–10273.