**3. Seismoacoustic scattering modeling in mutilayered media**

## **3.1. Methodology statement**

The problem to be considered is illustrated in Figure 10. In this model, there are *L* homogeneous layers (*i*) (*i*=1, 2, …, *L*) over a half-space, among which the *i*th layer is bounded by two irregular interfaces (*<sup>i</sup>*-1) and (*<sup>i</sup>*) (Throughout the section, the superscripts with round brackets are used for layer index and the subscripts with round brackets for interface index). The uppermost layer is water which is assumed to be linear liquid without viscosity. An arbitrary seismic source is embedded in the *s*th layer. The properties of an arbitrary solid layer (*i*) are described by the elastic constants (*i*) , (*i*) and mass density (*i*) . This problem will be solved by a boundary element method with a global matrix propagator. In this Subsection, we will present the formulations of the method in the solid layers and the fluid layer respectively.

**Figure 10.** A multi-layered water/solid model with irregular interfaces

#### *3.1.1. Formulations in solid layers*

138 Wave Processes in Classical and New Solids

TIME (s)

**3.1. Methodology statement** 

homogeneous layers (*i*)

arbitrary solid layer (*i*)

layers and the fluid layer respectively.

**Figure 10.** A multi-layered water/solid model with irregular interfaces

**Figure 9.** Synthetic seismograms for the model in Figure 8. The characteristic frequency is 1.0 Hz.


OFFSET (km)

The problem to be considered is illustrated in Figure 10. In this model, there are *L*

bounded by two irregular interfaces (*<sup>i</sup>*-1) and (*<sup>i</sup>*) (Throughout the section, the superscripts with round brackets are used for layer index and the subscripts with round brackets for interface index). The uppermost layer is water which is assumed to be linear liquid without viscosity. An arbitrary seismic source is embedded in the *s*th layer. The properties of an

This problem will be solved by a boundary element method with a global matrix propagator. In this Subsection, we will present the formulations of the method in the solid

are described by the elastic constants

(*i*=1, 2, …, *L*) over a half-space, among which the *i*th layer is

(*i*) , (*i*)

and mass density

(*i*) .

**3. Seismoacoustic scattering modeling in mutilayered media** 

Consider the boundary element problem in a region, where the elastic wave equation in frequency domain is given by

$$
\mu \nabla^2 \mathbf{u} + (\lambda + \mu) \nabla \nabla \cdot \mathbf{u} + \rho \omega^2 \mathbf{u} = -\mathbf{f},
\tag{29}
$$

in which **u** is the seismic displacement response in frequency domain, and **f** is the seismic source. Suppose the seismic source distribution consists of a simple point source at a position vector **s**. With the aid of the free-space Green's function, Eq. (29) can thus be transformed into the following boundary integral equation for the displacement *uj*(**r**) on the boundary of the region [30]

$$
\Delta(\mathbf{r})u\_k(\mathbf{r}) + \int\_{\Gamma} u\_j(\mathbf{r'}) T\_{jk} \left(\mathbf{r}, \mathbf{r'}\right) d\Gamma(\mathbf{r'}) = \int\_{\Gamma} t\_j(\mathbf{r'}) U\_{jk} \left(\mathbf{r}, \mathbf{r'}\right) d\Gamma(\mathbf{r'}) + f\_j U\_{jk} \left(\mathbf{r}, \mathbf{s}\right), \qquad (30)
$$

where *j* and *k* can be 1, 2. **r** and **r** are the position vectors of "field point" and "source point" on the boundary , respectively. The coefficient (**r**) generally depends on the local geometry at **r**, *tj*(**r**) are the traction components, and *Ujk*(**r**, **r**) and *Tjk*(**r**, **r**) are the fundamental solutions for displacements and tractions, respectively. The source exciting direction can be controlled arbitrarily by changing one of two components *fj*. Application of Eq. (30) in domain (*i*) and discretization of interfaces (*<sup>i</sup>*) and (*<sup>i</sup>*+1) each into *N* elements (for simplicity, we use constant element) give rise to

$$\mathbf{H}^{(l)}\mathbf{u}^{(l)} = \mathbf{G}^{(l)}\mathbf{t}^{(l)} + \mathbf{F}\boldsymbol{\delta}\_{l\boldsymbol{\delta}\boldsymbol{\prime}}\tag{31}$$

where **H**(*i*) and **G**(*i*) are the coefficient matrices obtained by integrating the fundamental solutions *Tjk* and *Ujk* over elements at both interfaces (*<sup>i</sup>*) and (*<sup>i</sup>*+1), and **u**(*i*) and **t**(*i*) are vectors containing element displacements and tractions at both interfaces (*<sup>i</sup>*) and (*<sup>i</sup>*+1), respectively. Obviously, **H**(*i*) and **G**(*i*) are of the dimension 4*N*4*N*, and **u**(*i*) and **t**(*i*) are of the dimension 4*N*1. In the multilayered medium (as shown in Figure 10), considering the continuity conditions of displacements and tractions at inner interfaces, Eq. (31) needs to be rewritten into

$$
\begin{bmatrix}
\mathbf{H}\_{\langle l\rangle}^{\langle l\rangle} & \mathbf{H}\_{\langle l+1\rangle}^{\langle l\rangle}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}\_{\langle l\rangle}^{\langle l\rangle} \\
\mathbf{u}\_{\langle l+1\rangle}^{\langle l\rangle}
\end{bmatrix} = 
\begin{bmatrix}
\mathbf{G}\_{\langle l\rangle}^{\langle l\rangle} & \mathbf{G}\_{\langle l+1\rangle}^{\langle l\rangle}
\end{bmatrix}
\begin{bmatrix}
\mathbf{t}\_{\langle l\rangle}^{\langle l\rangle} \\
\mathbf{t}\_{\langle l+1\rangle}^{\langle l\rangle}
\end{bmatrix} + \mathbf{F}\delta\_{sl\nu} \tag{32}
$$

where **u**(*i*) and **t**(*i*) with subscripts (*i*) and (*i*+1) corresponds separately to the element displacements and tractions at (*<sup>i</sup>*) and (*<sup>i</sup>*+1). Similarly, the subscripts (*i*) and (*i*+1) in **H**(*i*) and **G**(*i*) indicate the integration over elements at (*<sup>i</sup>*) and (*<sup>i</sup>*+1), respectively. Obviously, the matrices **H**(*i*) and **G**(*i*) with subscripts (*i*) or (*i*+1) are of the dimension 4*N*2*N*.

For the same purpose as in Section 2, the upward direction is defined as the positive direction for tractions so that they have unique values at each interface [18]. Eq. (32) thus becomes

$$\mathbf{H}\_{\{l\}}^{\{l\}}\mathbf{u}\_{\{l\}} + \mathbf{H}\_{\{l+1\}}^{\{l\}}\mathbf{u}\_{\{l+1\}} = \mathbf{G}\_{\{l\}}^{\{l\}}\mathbf{t}\_{\{l\}} - \mathbf{G}\_{\{l+1\}}^{\{l\}}\mathbf{t}\_{\{l+1\}} + \mathbf{F}\delta\_{\text{sl}}\tag{33}$$

$$\begin{cases} \mathbf{u}\_{(l)} = \mathbf{D}\_{uu}^{(l)} \mathbf{u}\_{(l+1)}\\ \mathbf{t}\_{(l+1)} = \mathbf{D}\_{tu}^{(l)} \mathbf{u}\_{(l+1)} \end{cases}, \ t = 2, 3, \cdots, s - 1,\tag{34}$$

$$\mathbf{H}\_{\{l\}}^{\{l\}}\mathbf{D}\_{uu}^{\{l\}}\mathbf{u}\_{\{l+1\}} + \mathbf{H}\_{\{l+1\}}^{\{l\}}\mathbf{u}\_{\{l+1\}} = \mathbf{G}\_{\{l\}}^{\{l\}}\mathbf{D}\_{tu}^{\{l-1\}}\mathbf{D}\_{uu}^{\{l\}}\mathbf{u}\_{\{l+1\}} - \mathbf{G}\_{\{l+1\}}^{\{l\}}\mathbf{D}\_{tu}^{\{l\}}\mathbf{u}\_{\{l+1\}}.\tag{35}$$

$$
\begin{pmatrix} \mathbf{D}\_{uu}^{(l)} \\ \mathbf{D}\_{tu}^{(l)} \end{pmatrix} = - \begin{bmatrix} \mathbf{H}\_{(l)}^{(l)} - \mathbf{G}\_{(l)}^{(l)} \mathbf{D}\_{tu}^{(l-1)}, & \mathbf{G}\_{(l+1)}^{(l)} \end{bmatrix}^{-1} \mathbf{H}\_{(l+1)}^{(l)}, \quad i = 2, 3, \cdots, s - 1. \tag{36}
$$

$$\begin{cases} \mathbf{u}\_{(l+1)} = \mathbf{D}\_{uu}^{(l)} \mathbf{u}\_{(l)}\\ \mathbf{t}\_{(l)} = \mathbf{D}\_{tu}^{(l)} \mathbf{u}\_{(l)} \end{cases}, \ t = s+1, s+2, \cdots, L. \tag{37}$$

$$
\begin{pmatrix} \mathbf{D}\_{uu}^{(l)} \\ \mathbf{D}\_{tu}^{(l)} \end{pmatrix} = - \begin{bmatrix} \mathbf{H}\_{(l+1)}^{(l)} + \mathbf{G}\_{(l+1)}^{(l)} \mathbf{D}\_{tu}^{(l+1)} \end{bmatrix}, \quad -\mathbf{G}\_{(l)}^{(l)} \Big|^{-1} \mathbf{H}\_{(l)}^{(l)}, \ i = L - 1, L - 2, \cdots, s + 1,\tag{38}
$$

$$\begin{cases} \mathbf{D}\_{uu}^{(L)} = \mathbf{0} \\ \mathbf{D}\_{tu}^{(L)} = \left[ \mathbf{G}\_{(L)}^{(L)} \right]^{-1} \mathbf{H}\_{(L)}^{(L)} \end{cases} \tag{39}$$

$$
\nabla^2 p + k^2 p = 0,\tag{40}
$$

$$\mathbf{c}(\mathbf{r})p(\mathbf{r}) + \int\_{\Gamma} p(\mathbf{r'}) \frac{\partial \mathcal{G}(\mathbf{r}, \mathbf{r})}{\partial n\_{\mathbf{r}'}} d\Gamma(\mathbf{r'}) = \int\_{\Gamma} \ q(\mathbf{r'}) \, G(\mathbf{r}, \mathbf{r'}) \, d\Gamma(\mathbf{r'}) \tag{41}$$

$$\mathbf{H}\_{\text{(2)}}^{\text{(1)}} \mathbf{p}\_{\text{(2)}} = \mathbf{G}\_{\text{(1)}}^{\text{(1)}} \mathbf{q}\_{\text{(1)}} + \mathbf{G}\_{\text{(2)}}^{\text{(1)}} \mathbf{q}\_{\text{(2)}} \tag{42}$$

$$\begin{cases} \mathbf{q}\_{(1)} = \mathbf{D}\_{qq}^{(1)} \mathbf{q}\_{(2)}\\ \mathbf{p}\_{(2)} = \mathbf{D}\_{pq}^{(1)} \mathbf{q}\_{(2)} \end{cases} \tag{43}$$

$$\mathbf{H}\_{\{2\}}^{\{1\}} \mathbf{D}\_{pq}^{\{1\}} \mathbf{q}\_{\{2\}} = \mathbf{G}\_{\{1\}}^{\{1\}} \mathbf{D}\_{qq}^{\{1\}} \mathbf{q}\_{\{2\}} + \mathbf{G}\_{\{2\}}^{\{1\}} \mathbf{q}\_{\{2\}}.\tag{44}$$

$$
\begin{pmatrix} \mathbf{D}\_{qq}^{(1)} \\ \mathbf{D}\_{pq}^{(1)} \end{pmatrix} = \begin{bmatrix} -\mathbf{G}\_{(1)}^{(1)} & \mathbf{H}\_{(2)}^{(1)} \end{bmatrix}^{-1} \mathbf{G}\_{(2)}^{(1)} \tag{45}
$$

$$\begin{aligned} \mathbf{t} \cdot \mathbf{s} &= 0\\ \mathbf{t} \cdot \mathbf{n} &= -p\\ \mathbf{u} \cdot \mathbf{n} &= -u\_{\mathbf{f}} = -\frac{1}{\rho\_{\mathbf{f}} \omega^{2}} q' \end{aligned} \tag{46}$$

$$\mathbf{D}\_{tu}^{\text{(1)}} = \rho\_{\text{f}} \omega^2 \begin{bmatrix} \mathbf{D}\_{pq}^{\text{(1)}} \cdot \mathbf{N}\_1^{\text{e}} \cdot (\mathbf{N}\_1^{\text{e}})^\text{T} & \mathbf{D}\_{pq}^{\text{(1)}} \cdot \mathbf{N}\_2^{\text{e}} \cdot (\mathbf{N}\_1^{\text{e}})^\text{T} \\ \mathbf{D}\_{pq}^{\text{(1)}} \cdot \mathbf{N}\_1^{\text{e}} \cdot (\mathbf{N}\_2^{\text{e}})^\text{T} & \mathbf{D}\_{pq}^{\text{(1)}} \cdot \mathbf{N}\_2^{\text{e}} \cdot (\mathbf{N}\_2^{\text{e}})^\text{T} \end{bmatrix} \tag{47}$$

$$\mathbf{N}\_{1}^{\mathbf{e}} = \begin{bmatrix} n\_{1}^{\mathbf{e}\_{1}} & n\_{1}^{\mathbf{e}\_{2}} & \cdots & n\_{1}^{\mathbf{e}\_{N}} \\ n\_{1}^{\mathbf{e}\_{1}} & n\_{1}^{\mathbf{e}\_{2}} & \cdots & n\_{1}^{\mathbf{e}\_{N}} \\ \vdots & \vdots & \ddots & \vdots \\ n\_{1}^{\mathbf{e}\_{1}} & n\_{1}^{\mathbf{e}\_{2}} & \cdots & n\_{1}^{\mathbf{e}\_{N}} \end{bmatrix}, \mathbf{N}\_{2}^{\mathbf{e}} = \begin{bmatrix} n\_{2}^{\mathbf{e}\_{1}} & n\_{2}^{\mathbf{e}\_{2}} & \cdots & n\_{2}^{\mathbf{e}\_{N}} \\ n\_{2}^{\mathbf{e}\_{1}} & n\_{2}^{\mathbf{e}\_{2}} & \cdots & n\_{2}^{\mathbf{e}\_{N}} \\ \vdots & \vdots & \ddots & \vdots \\ n\_{2}^{\mathbf{e}\_{1}} & n\_{2}^{\mathbf{e}\_{2}} & \cdots & n\_{2}^{\mathbf{e}\_{N}} \end{bmatrix} \tag{48}$$

$$\mathbf{H}\_{\text{(s)}}^{\text{(s)}}\mathbf{u}\_{\text{(s)}} + \mathbf{H}\_{\text{(s+1)}}^{\text{(s)}}\mathbf{u}\_{\text{(s+1)}} = \mathbf{G}\_{\text{(s)}}^{\text{(s)}}\mathbf{t}\_{\text{(s)}} - \mathbf{G}\_{\text{(s+1)}}^{\text{(s)}}\mathbf{t}\_{\text{(s+1)}} + \mathbf{F}\_{\text{(s)}}\tag{49}$$

$$\mathbf{t}\_{\{s\}} = \mathbf{D}\_{tu}^{\{s-1\}} \mathbf{u}\_{\{s\}}, \mathbf{t}\_{\{s+1\}} = \mathbf{D}\_{tu}^{\{s+1\}} \mathbf{u}\_{\{s+1\}}.\tag{50}$$

$$\left[\mathbf{H}\_{\text{(s)}}^{\text{(s)}} - \mathbf{G}\_{\text{(s)}}^{\text{(s)}} \mathbf{D}\_{\text{tu}}^{\text{(s-1)}}\right] \mathbf{u}\_{\text{(s)}} + \left[\mathbf{H}\_{\text{(s+1)}}^{\text{(s)}} + \mathbf{G}\_{\text{(s+1)}}^{\text{(s)}} \mathbf{D}\_{\text{tu}}^{\text{(s+1)}}\right] \mathbf{u}\_{\text{(s+1)}} = \mathbf{F},\tag{51}$$

$$\mathbf{u}\_{\text{(2)}} = \mathbf{D}\_{uu}^{\{2\}} \mathbf{D}\_{uu}^{\{3\}} \cdots \mathbf{D}\_{uu}^{\{s-1\}} \mathbf{u}\_{\text{(s)}} \tag{52}$$


order of calculating the global matrix propagators in the fluid layer and the solid layer, the details of which are suppressed here.

An Efficient Approach for Seismoacoustic Scattering Simulation Based on Boundary Element Method 145

where the pulse duration *T*S is set to 4.2 s and *α* is set to 2.1. Eq. (53) is exactly the same as

**Figure 12.** The time function (a) and its Fourier spectra (b) used for the two models in Figure 11

F r e q u e n c y ( H z ) (a) (b)

Figure 13 shows the displacements of both horizontal and vertical components along the fluid-solid interface of *Model* 1. The one on the left is calculated by the reflection/transmission matrix method in [20]. And the one on the right is calculated by the present method. Inside the irregular part of the model, i.e., the basin structure, the largeamplitude multiple reflections in the water basin are well observed in the synthetic waveforms calculated by both our method and the reflection/transmission matrix method. Both the scattered *P* wave (P) and *Rayleigh* wave (R), which are indicated in the left plot, can

(a) (b)

<sup>10</sup>-3

<sup>10</sup>-2

<sup>10</sup>-1

<sup>10</sup><sup>0</sup>

0.01 0.1 1

Figure 14 shows the displacements of both horizontal and vertical components at the fluidsolid interface of *Model* 2. This model is a little bit more complicated than *Model* 1. The flat fluid layer on both sides of *Model* 2 will cause multiple reflections of the incident plane *P*wave easily. The multiple reflections will be interfered by the waves scattered by the irregular fluid-solid interface, which leads to the complexity of the wave field. In this case, the 'observation points' are all located at the fluid layer bottom. The one on the left is calculated by the reflection/transmission matrix method in [20]. And the one on the right is calculated by the present method. Again, the synthetic seismograms calculated by both methods agree with each other very well. The periodic multiple reflections at the flat portion of the fluid layer are clearly observed and interfered by the scattered waves from the inside

of the irregular part of the model, which are modeled very well by both methods.

*3.2.1. Model 1 subjected to a plane wave incidence* 


T i m e ( s )

0

0.2

0.4

0.6

0.8

be well and clearly observed in our calculated seismograms.

*3.2.2. Model 2 subjected to a plane wave incidence* 

the Eq. (88) in [20] and its shape and spectrum are shown in Figure 12.

Although the boundary element method is known as the best way to model wave propagation problems in unbounded media, it is still necessary to make some special treatment on the truncation edges of the model in order to avoid the appearance of some unphysical waves. For such purpose, many different techniques [33-40] applicable to the boundary element method have been developed.
