**2.1. Methodology statement**

126 Wave Processes in Classical and New Solids

major approaches to study a laterally inhomogeneous model mainly rely on numerical methods which include domain methods and boundary methods. The domain numerical methods (such as finite difference [11], finite element [12], pseudo-spectral method [13], and spectral element method [14]) have become the standard tools for seismic wave modeling, however, these methods do not explicitly consider the boundary continuity conditions between different formations. That disables the methods to sufficient accuracy for modeling the reflection/transmission across irregular interfaces. Furthermore, they have difficulties in dealing with a large-scaled model because they require the subdivision of the domain into elements, especially the case when the domain needs to be extended to infinity. Boundary numerical methods have emerged as good alternatives to the domain numerical methods in dealing with an infinite continuum model because the radiation condition at infinity is easily fulfilled. They also have an obvious advantage over the domain numerical methods: the dimensionality of the problem under consideration is reduced by one order, because only a boundary instead of a domain discretization is required. Among this type of methods, the boundary element method (BEM) has been extensively used in simulating seismic wave excitation and propagation [15]. However, for a stratified model with irregular interfaces, the dimension of the final system of simultaneous equations in the traditional boundary element metod is proportional to the production of the element number and the interface number, which leads to an exponential increase in computing time and memory requirement. In order to overcome that problem, different approaches have been tried, for example, Fu [16] proposed an improved block Gaussian elimination scheme and Bouchon et al. [17] introduced a sparse method. Recently, global matrix propagators were introduced to improve the efficiency of the traditional BEM for modeling seismic wave excitation and propagation in multilayered solids

[18-19], which expects great savings in computing time and memory requirement.

On the other hand, the seismoacoustic scattering due to an irregular fluid-solid interface must be considered when we model the seismic wave propagation in oceanic regions or gulf areas. This kind of problem is related to a wide range of seismic research conducted at or close to oceanic regions or gulf areas [20], such as deep ocean acoustic experiments, ocean bottom seismic observations, or the interpretation of the effects of water layer on the observed surface wave traveling through gulf areas. Many metropolitan cities all over the world are located at or close to seaside, so it is important to take into account the effect of the sea water when conducting the earthquake resistant design for high-rise buildings in such big gulf cities. Also, the underground structures beneath sea bottom are known to be strongly inhomogeneous. Therefore, it is necessary to simulate the seismoacoustic scattering in irregularly multilayered elastic media overlain by a fluid layer due to some scenario earthquake event. Furthermore, ocean bottom observations are a key feature in the modern world-wide standardized seismograph network [21]. Correct seismic wave modeling can help to obtain valid interpretation of the observed waveforms. In the past three decades, many techniques, classified into different categories [7, 18], have been developed for calculating seismic wave excitation and propagation in laterally inhomogeneous media. Among those methods, the finite difference method (FDM) [22-23] and the finite element method (FEM) [24] become very popular in modelling seismoacoustic scattering due to irregular fluid-solid interface because of their flexibility in modeling complex media.

The problem to be studied is illustrated in Figure 1. 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 superscript with round brackets indicates layer index and the subscript with round brackets interface index). The uppermost interface is a free surface, and an arbitrary source is embedded in the *s*th layer. Assume each individual layer to be isotropic, linearly elastic material, so the tensor of elastic constants is determined by two independent Lame constants because of the rotational and translational symmetry. For the two-dimensional *SH* case under consideration, the physical properties of each layer can only be described by the shear modulus (*i*) and mass density (*i*) . This problem will be solved by a boundary element method with a global matrix propagator. In this section, we will present the formulations of the method in the solid layers.

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

(, ) (, ) ( )

**r r rr r r**

*i*

<sup>i</sup> (, ) ( ), <sup>4</sup>

**r r r r r rn**

shear wave velocity in the *i*th layer, represents a Dirac delta function which goes to infinity at the point **r r** and is equal to zero elsewhere, and **n** denotes the outward normal vector of the surface on which the traction is acting. Please note that the position vectors **r** and **r** are

generated by a concentrated unit source acting at the "source point" **r**. The symbol *H*2 means

To solve Eq. (3), we need to do discretization on the boundary. For simplicity, we discretize the boundary (*<sup>i</sup>*) and (*<sup>i</sup>*+1) each into *N* constant elements (the results in the case of linear element can be deduced by analog), so that the coefficient *C* is taken as 1/2 always. For a

*i*

(**r**, **r**) are the fundamental solutions for

=/*c*(*i*)

(**r**, **r**) represents the displacement field at the "field point" **r**

) over field element *k*, *Gjk* is an integral of *U*(**r***<sup>j</sup>*

(6)

(4)

being the

bounded

(5)

with *c*(*i*)

, **r***<sup>k</sup>*

. And the coefficient

(7)

(1) **t** 0, (8)

) over field

of the

(**r**, **r**) and *T*\*

\* 2 () ( ) 0

1

*i i*

**r r r r**

( ) \* 2 ()

() 2 \* () 2 \*

*U U*

 

(, ) ( , ),

**r r r r n**

now separately defined as "field point" and "source point" in the whole domain (*i*)

given element *j* (taken as the "source point"), Eq. (3) can be discretized as follows

2 2

*k k u Hu Gt F*

1 1 <sup>1</sup> , <sup>2</sup> *N N*

where the superscript (*i*) denotes the layer index. Obviously, the column vector **u**(*i*)

The boundary and continuity conditions for the whole model under consideration are:

(1)

*j jk k jk k j si*

element *k*, and *Fj* is the source term. When the source element goes from *j* = 1 to 2*N*, one can

() () , *i i si* **Hu Gt F**

dimension 2*N*1 consists of the element displacements on both the boundary (*<sup>i</sup>*) and the

*i*

\* () \*

the Hankel function of the second kind and its subscript indicates the order.

, **r***<sup>k</sup>*

obtain a system of equations which can be expressed in matrix form

boundary (*<sup>i</sup>*+1), the same situation applies to the column vector **t**(*i*)

matrices **H** and **G** are of the dimension 2*N*2*N*.

1. the traction-free condition on the free surface,

*T U*

where i is the imaginary unit different from the italic *i* being index, *k*(*i*)

*i i*

*<sup>k</sup> T Hk*

*U H k*

<sup>i</sup> (, ) <sup>4</sup>

traction components, and *U*\*

displacements and tractions, respectively,

and they satisfy the following equations

by the boundary Here, *U*\*

where *Hjk* is an integral of *T*(**r***<sup>j</sup>*

**Figure 1.** A multilayered 2D model with irregular interfaces

Consider the boundary element method problem in the *i*th layer, the scalar wave equation in frequency domain is given by (the details on the derivation of the SH-wave equation (1) and PSv-wave equation (29), which are suppressed here for brevity, refer to [27])

$$
\mu^{(i)}\nabla^2 u(\mathbf{r}, \alpha) + \rho^{(i)}\alpha^2 u(\mathbf{r}, \alpha) = -F(\mathbf{r}, \alpha)\delta\_{si'} \tag{1}
$$

where *u*(**r**, ) is the seismic response at a position vector **r** in frequency domain, *si* = 1 for *i* = *s* and vanishes else, and *F*(**r**, ) is the seismic source in the *s*th layer.

Suppose the seismic source distribution consists of a simple point source at a position vector **s**, i.e.,

$$F(\mathbf{r}, \alpha) = F(\alpha) \delta(\mathbf{r} - \mathbf{s}). \tag{2}$$

With the aid of the free-space Green's function, Eq. (1) can thus be transformed into the following boundary integral equation for the displacement *u* on the boundary (*<sup>i</sup>*)+(*i*+1)

$$\frac{1}{\Gamma} \mathrm{Cu}(\mathbf{r}) + \int \mathrm{T}^\*(\mathbf{r}, \mathbf{r}') \mu(\mathbf{r}') \mathrm{d}\Gamma(\mathbf{r}') = \int \mathrm{L}^\*(\mathbf{r}, \mathbf{r}') t(\mathbf{r}') \mathrm{d}\Gamma(\mathbf{r}') + \delta\_{\mathrm{si}} \mathrm{L}^\*(\mathbf{r}, \mathbf{s}) F(\mathbf{r}, \alpha) \,\tag{3}$$

where **r** and **r** are the position vectors of "field point" and "source point" on the boundary , respectively. The coefficient *C* generally depends on the local geometry at **r**, *t*(**r**) are the traction components, and *U*\* (**r**, **r**) and *T*\* (**r**, **r**) are the fundamental solutions for displacements and tractions, respectively,

$$\begin{aligned} \text{iL}^\*\left(\mathbf{r}, \mathbf{r}'\right) &= -\frac{\dot{\mathbf{i}}}{4\mu^{(i)}} H\_0^2 \left(k^{(i)} \left|\mathbf{r} - \mathbf{r}'\right|\right) \\ T \stackrel{\text{i}}{=} \frac{\text{i} k^{(i)}}{4} H\_1^2 \left(k^{(i)} \left|\mathbf{r} - \mathbf{r}'\right|\right) (\mathbf{r}' - \mathbf{r}) \mathbf{n}, \end{aligned} \tag{4}$$

and they satisfy the following equations

128 Wave Processes in Classical and New Solids

each layer can only be described by the shear modulus

**Figure 1.** A multilayered 2D model with irregular interfaces

where *u*(**r**,

**s**, i.e.,

*s* and vanishes else, and *F*(**r**,

we will present the formulations of the method in the solid layers.

determined by two independent Lame constants because of the rotational and translational symmetry. For the two-dimensional *SH* case under consideration, the physical properties of

will be solved by a boundary element method with a global matrix propagator. In this section,

Consider the boundary element method problem in the *i*th layer, the scalar wave equation in frequency domain is given by (the details on the derivation of the SH-wave equation (1)

() 2 () 2 (, ) (, ) (, ) , *i i*

) is the seismic response at a position vector **r** in frequency domain,

Suppose the seismic source distribution consists of a simple point source at a position vector

*F F* ( , ) ( ) ( ). **r rs**

 

With the aid of the free-space Green's function, Eq. (1) can thus be transformed into the following boundary integral equation for the displacement *u* on the boundary (*<sup>i</sup>*)+(*i*+1)

\*\* \* ( ) ( , ) ( )d ( ) ( , ) ( )d ( ) ( , ) ( , ), *Cu T u U t*

where **r** and **r** are the position vectors of "field point" and "source point" on the boundary , respectively. The coefficient *C* generally depends on the local geometry at **r**, *t*(**r**) are the

 

) is the seismic source in the *s*th layer.

*u uF* **r rr**

*si*

(2)

*siU F*

**r rr r r rr r r rs r** (3)

(1)

*si* = 1 for *i* =

 

and PSv-wave equation (29), which are suppressed here for brevity, refer to [27])

  (*i*)

and mass density

(*i*)

. This problem

$$\begin{aligned} \mu^{(i)} \nabla^2 \mathcal{U}^\*(\mathbf{r}, \mathbf{r}') + \rho^{(i)} \alpha^2 \mathcal{U}^\*(\mathbf{r}, \mathbf{r}') &= -\Delta(\mathbf{r}' - \mathbf{r}) \\ \mathcal{T}^\*(\mathbf{r}, \mathbf{r}') = \mu^{(i)} \frac{\partial}{\partial \mathbf{n}} \mathcal{U}^\*(\mathbf{r}, \mathbf{r}'), \end{aligned} \tag{5}$$

where i is the imaginary unit different from the italic *i* being index, *k*(*i*) =/*c*(*i*) with *c*(*i*) being the shear wave velocity in the *i*th layer, represents a Dirac delta function which goes to infinity at the point **r r** and is equal to zero elsewhere, and **n** denotes the outward normal vector of the surface on which the traction is acting. Please note that the position vectors **r** and **r** are now separately defined as "field point" and "source point" in the whole domain (*i*) bounded by the boundary Here, *U*\* (**r**, **r**) represents the displacement field at the "field point" **r** generated by a concentrated unit source acting at the "source point" **r**. The symbol *H*2 means the Hankel function of the second kind and its subscript indicates the order.

To solve Eq. (3), we need to do discretization on the boundary. For simplicity, we discretize the boundary (*<sup>i</sup>*) and (*<sup>i</sup>*+1) each into *N* constant elements (the results in the case of linear element can be deduced by analog), so that the coefficient *C* is taken as 1/2 always. For a given element *j* (taken as the "source point"), Eq. (3) can be discretized as follows

$$\frac{1}{2}\mu\_j + \sum\_{k=1}^{2N} H\_{jk}\mu\_k = \sum\_{k=1}^{2N} G\_{jk}t\_k + F\_j\delta\_{si'} \tag{6}$$

where *Hjk* is an integral of *T*(**r***<sup>j</sup>* , **r***<sup>k</sup>* ) over field element *k*, *Gjk* is an integral of *U*(**r***<sup>j</sup>* , **r***<sup>k</sup>* ) over field element *k*, and *Fj* is the source term. When the source element goes from *j* = 1 to 2*N*, one can obtain a system of equations which can be expressed in matrix form

$$\mathbf{H}\mathbf{u}^{(i)} = \mathbf{G}\mathbf{t}^{(i)} + \mathbf{F}\boldsymbol{\delta}\_{si'} \tag{7}$$

where the superscript (*i*) denotes the layer index. Obviously, the column vector **u**(*i*) of the dimension 2*N*1 consists of the element displacements on both the boundary (*<sup>i</sup>*) and the boundary (*<sup>i</sup>*+1), the same situation applies to the column vector **t**(*i*) . And the coefficient matrices **H** and **G** are of the dimension 2*N*2*N*.

The boundary and continuity conditions for the whole model under consideration are:

1. the traction-free condition on the free surface,

$$\mathbf{t}\_{(1)}^{(1)} = \mathbf{0},\tag{8}$$

2. the continuity conditions of displacement and traction at the inner interfaces,

$$\begin{aligned} \mathbf{u}\_{(i+1)}^{(i)} &= \mathbf{u}\_{(i+1)}^{(i+1)} \\ \mathbf{t}\_{(i+1)}^{(i)} &= -\mathbf{t}\_{(i+1)}^{(i+1)} \end{aligned} \tag{9}$$

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

, 1,2, , 1,

, 1, 2, , ,

**t Du** (13)

**t Du** (14)

(1) (2) (2) (2) (2) (2) 0. *uu tu* **HDu Hu GDu** (15)

**D** (16)

*i s*

is due to the above-defined positive direction. Equation (12) gives the relationship between the element tractions and the element displacements at the interfaces (*<sup>i</sup>*) and (*<sup>i</sup>*+1) which

From now on, we introduce the global matrix propagators which is the key to the introduced method under discussion. For completeness, we cite some formula directly from

*is s L* 

In Eqs. (13) and (14), the superscripts in the global matrix propagators denote the layer index, and the subscripts in the element displacements and tractions denote the interface index. Each layer has two global matrix propagators which function in different ways for the layers above the source and the layers below the source. For one arbitrary layer above the source, **D***uu* transfers information of the element displacements from the upper interface to the lower interface in the layer whilst **D***tu* bridges the element tractions and the element displacements at the lower interface of the layer. However, for one arbitrary layer below the source, **D***uu* transfers information of the element displacements from the lower interface to the upper interface in the layer whilst **D***tu* bridges the element tractions and the element displacements at the upper interface of the layer. That's the main idea of the global matrix propagators which propagate information from "top and bottom" where the boundary conditions are known to "middle" where the source is located by the displacement and traction continuity conditions at

In the following part, how to obtain the global matrix propagators recursively will be briefly

Because **u**(2) could be the displacement at the second interface excited by an arbitrary source, the necessary condition of Eq. (15) is that its coefficient matrix equals zero, which gives

> (1) <sup>1</sup> (1) (1) (1) (1) (2) (2) (1) . *uu*

**HG H**

reviewed. Applying boundary condition (8) and Eq. (13) with *i* 1 to Eq. (12) gives

(1) (1) (1) (1) (1)

 

enclose the *i*th layer.

the work by Ge and Chen [18] with minor modification.

for the layers above the source, and

for the layers below the source.

The global matrix propagators ( )*<sup>i</sup>* **<sup>D</sup>***uu* and ( )*<sup>i</sup>* **<sup>D</sup>***tu* are separately defined as [18]

**u Du**

**u Du**

inner interfaces. Both of the two matrices are of the dimension *NN*.

*tu*

**D**

( ) ( ) ( 1) ( ) ( 1) ( 1)

*i i uu i i i tu i*

( ) ( 1) ( ) ( ) () ()

*i i uu i i i tu i*

3. the radiation boundary conditions imposed on the far-field behavior at infinity,

$$\begin{aligned} \mathbf{u}\_{(L+1)}^{(L)} &= \mathbf{0} \\ \mathbf{t}\_{(L+1)}^{(L)} &= \mathbf{0} \end{aligned} \tag{10}$$

where the superscript denotes the layer index and the subscript denotes the interface index, which is followed as the same hereafter. And the minus sign in equation (9)2 is due to the definition of outward normal vector in boundary element method. Please note that the displacements and tractions in equations (8-10) are different from those in equation (7). The former are of the dimension *N*1 while the latter 2*N*1.

In order to apply those boundary conditions to equation (7), we need to rewrite equation (7) into

$$
\begin{bmatrix}
\mathbf{H}\_{\{i\}}^{(i)} & \mathbf{H}\_{\{i+1\}}^{(i)}
\end{bmatrix}
\begin{bmatrix}
\mathbf{u}\_{\{i\}}^{(i)} \\
\mathbf{u}\_{\{i+1\}}^{(i)}
\end{bmatrix} = \begin{bmatrix}
\mathbf{G}\_{\{i\}}^{(i)} & \mathbf{G}\_{\{i+1\}}^{(i)}
\end{bmatrix}
\begin{bmatrix}
\mathbf{t}\_{\{i\}}^{(i)} \\
\mathbf{t}\_{\{i+1\}}^{(i)}
\end{bmatrix} + \mathbf{F} \boldsymbol{\delta}\_{si'} \tag{11}
$$

where ��(�) (�) and ��(���) (�) are the element displacements separately corresponding to the upper interface and lower interface of the *i*th layer, similarly ��(�) (�) and ��(���) (�) are the element tractions corresponding to the upper interface and lower interface of the *i*th layer. Obviously, they are column vectors of the dimension *N*1. It should be pointed out that for one inner interface we should number the element displacements and tractions in the same order for the upper neighboring layer domain as that for the lower neighboring layer domain. Although this is a simple numbering, it is a common bug in computer programs. Look out for it! More attention should be paid to the coefficient matrices ��(�) (�) and ��(���) (�) , ��(�) (�) and ��(���) (�) which are all of the dimension 2*NN* (except the case when *i L*) and called matrix propagators. They will be used to form the global matrix propagators for the hybrid boundary element method being discussed in this paper.

For the sake of convenient derivation, we first define the upward direction as the positive direction of the tractions in the whole model. Following this, we just need to make a minor change to Eq. (11)

$$
\begin{bmatrix} \mathbf{H}\_{\text{(i)}}^{(i)} & \mathbf{H}\_{\text{(i+1)}}^{(i)} \end{bmatrix} \begin{bmatrix} \mathbf{u}\_{\text{(i)}} \\ \mathbf{u}\_{\text{(i+1)}} \end{bmatrix} = \begin{bmatrix} \mathbf{G}\_{\text{(i)}}^{(i)} & -\mathbf{G}\_{\text{(i+1)}}^{(i)} \end{bmatrix} \begin{bmatrix} \mathbf{t}\_{\text{(i)}} \\ \mathbf{t}\_{\text{(i+1)}} \end{bmatrix} + \mathbf{F} \delta\_{si} \tag{12}
$$

where only one subscript index is needed to indicate the element displacement and traction vector at different interfaces, and the minus sign in the front of the matrix propagator ��(���) (�)

is due to the above-defined positive direction. Equation (12) gives the relationship between the element tractions and the element displacements at the interfaces (*<sup>i</sup>*) and (*<sup>i</sup>*+1) which enclose the *i*th layer.

From now on, we introduce the global matrix propagators which is the key to the introduced method under discussion. For completeness, we cite some formula directly from the work by Ge and Chen [18] with minor modification.

The global matrix propagators ( )*<sup>i</sup>* **<sup>D</sup>***uu* and ( )*<sup>i</sup>* **<sup>D</sup>***tu* are separately defined as [18]

$$\begin{cases} \mathbf{u}\_{(i)} = \mathbf{D}\_{\mu u}^{(i)} \mathbf{u}\_{(i+1)}\\ \mathbf{t}\_{(i+1)} = \mathbf{D}\_{tu}^{(i)} \mathbf{u}\_{(i+1)} \end{cases}, i = 1, 2, \cdots, s - 1,\tag{13}$$

for the layers above the source, and

130 Wave Processes in Classical and New Solids

into

where ��(�) (�)

��(�) (�)

and ��(���)

change to Eq. (11)

and ��(���)

2. the continuity conditions of displacement and traction at the inner interfaces,

( ) ( 1) ( 1) ( 1) ( ) ( 1) ( 1) ( 1)

 

0 , 0

 

where the superscript denotes the layer index and the subscript denotes the interface index, which is followed as the same hereafter. And the minus sign in equation (9)2 is due to the definition of outward normal vector in boundary element method. Please note that the displacements and tractions in equations (8-10) are different from those in equation (7). The

In order to apply those boundary conditions to equation (7), we need to rewrite equation (7)

tractions corresponding to the upper interface and lower interface of the *i*th layer. Obviously, they are column vectors of the dimension *N*1. It should be pointed out that for one inner interface we should number the element displacements and tractions in the same order for the upper neighboring layer domain as that for the lower neighboring layer domain. Although this is a simple numbering, it is a common bug in computer programs.

matrix propagators. They will be used to form the global matrix propagators for the hybrid

For the sake of convenient derivation, we first define the upward direction as the positive direction of the tractions in the whole model. Following this, we just need to make a minor

> ( 1) ( 1) , *i i i i i i i i i i si i i*

 

**u t HH G G F**

where only one subscript index is needed to indicate the element displacement and traction vector at different interfaces, and the minus sign in the front of the matrix propagator ��(���)

( ) ( ) () () () () ( ) ( 1) ( ) ( 1)

Look out for it! More attention should be paid to the coefficient matrices ��(�)

( ) ( ) () () () () ( ) ( 1) ( ) ( 1) ( ) ( )

*i i i i i i i i i i si i i i i*

( ) ( )

*i i*

( 1) ( 1)

(�) which are all of the dimension 2*NN* (except the case when *i L*) and called

(�) are the element displacements separately corresponding to the upper

(�)

 

**u t H H GG F u t**

 

*i i i i i i i i*

**u u**

3. the radiation boundary conditions imposed on the far-field behavior at infinity,

**u**

former are of the dimension *N*1 while the latter 2*N*1.

interface and lower interface of the *i*th layer, similarly ��(�)

boundary element method being discussed in this paper.

( ) ( 1) ( ) ( 1)

*L L L L* 

,

**t t** (9)

**t** (10)

,

(11)

(�) are the element

(�)

 and ��(���) (�) ,

(�)

and ��(���)

**u t** (12)

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

for the layers below the source.

In Eqs. (13) and (14), the superscripts in the global matrix propagators denote the layer index, and the subscripts in the element displacements and tractions denote the interface index. Each layer has two global matrix propagators which function in different ways for the layers above the source and the layers below the source. For one arbitrary layer above the source, **D***uu* transfers information of the element displacements from the upper interface to the lower interface in the layer whilst **D***tu* bridges the element tractions and the element displacements at the lower interface of the layer. However, for one arbitrary layer below the source, **D***uu* transfers information of the element displacements from the lower interface to the upper interface in the layer whilst **D***tu* bridges the element tractions and the element displacements at the upper interface of the layer. That's the main idea of the global matrix propagators which propagate information from "top and bottom" where the boundary conditions are known to "middle" where the source is located by the displacement and traction continuity conditions at inner interfaces. Both of the two matrices are of the dimension *NN*.

In the following part, how to obtain the global matrix propagators recursively will be briefly reviewed. Applying boundary condition (8) and Eq. (13) with *i* 1 to Eq. (12) gives

$$\mathbf{H}\_{\{1\}}^{(1)}\mathbf{D}\_{\mu\mu}^{(1)}\mathbf{u}\_{\{2\}} + \mathbf{H}\_{\{2\}}^{(1)}\mathbf{u}\_{\{2\}} + \mathbf{G}\_{\{2\}}^{(1)}\mathbf{D}\_{\mu}^{(1)}\mathbf{u}\_{\{2\}} = 0. \tag{15}$$

Because **u**(2) could be the displacement at the second interface excited by an arbitrary source, the necessary condition of Eq. (15) is that its coefficient matrix equals zero, which gives

$$
\begin{bmatrix} \mathbf{D}\_{uu}^{(1)} \\ \mathbf{D}\_{tu}^{(1)} \end{bmatrix} = - \begin{bmatrix} \mathbf{H}\_{\{1\}}^{(1)} & \mathbf{G}\_{\{2\}}^{(1)} \end{bmatrix}^{-1} \mathbf{H}\_{\{2\}}^{(1)}.\tag{16}
$$

Similarly, for the other layers above the source, i.e., *i* 2, , *s*1, we have

$$
\begin{bmatrix} \mathbf{D}\_{uu}^{(i)} \\ \mathbf{D}\_{tu}^{(i)} \end{bmatrix} = - \begin{bmatrix} \mathbf{H}\_{\{i\}}^{(i)} - \mathbf{G}\_{\{i\}}^{(i)} \mathbf{D}\_{tu}^{(i-1)} & \mathbf{G}\_{\{i+1\}}^{(i)} \end{bmatrix}^{-1} \mathbf{H}\_{\{i+1\}}^{(i)}.\tag{17}
$$

Substituting boundary condition (10) and Eq. (14) with *i L* into Eq. (12) and considering the same reason as that giving rise to Eq. (16), we get

$$\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{18}$$

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

, for 1, 2,..., 1,

*i s*

, for 1, 2,..., .

**r rr r r rr r r rs r r** (25)

*is s L*

(24)

 

*f t* (26)

(23)

( ) ( 1) ( 1) ( ) ( )

Then the displacement at any observation points in the whole model can be obtained by the

\* \* \* ( ) ( ) ( , ) ( )d ( ) ( , ) ( )d ( ) ( , ) ( , ), , *<sup>i</sup>*

which can be transformed into a summation of the element displacements and tractions obtained just now without extra effort. Usually, the displacements at the free surface is of the most interest. By taking the inverse Fourier transform on the above frequency domain

In the following, we present the calculated synthetic seismograms for some typical models for which exiting results are available. Throughout the following examples, the source time

> 2 22 2 22 <sup>0</sup> ( ) 2 1 exp , *<sup>c</sup> ut f t*

where *f*0 is the characteristic frequency of the wavelet and *t* is the time. Observation points

The first model is a semi-circular canyon in an elastic homogenous half-space, as shown in Figure 2. Although the shape is too simple for a realistic canyon, the wave reflection and diffraction due to the topographical irregularity in time domain are quite complex, which can serve for an eligible test to the present method and our computation program. Here, we calculate the time-domain responses of the surface along the semi-circular canyon subject to incident SH waves. The shear velocity and mass density are taken as 1 km/s and 1 g/cm3

 

solution, we can finally get the time domain solution, i.e., the synthetic seismogram.

*ii s i uu uu uu s*

( ) ( 1) ( 1) ( 1) ( 1)

*si uT u U t UF*

*ii s i uu uu uu s*

function of the synthetic seismograms is a Ricker wavelet defined as

**u DD D u**

**u DD D u**

( ) ( 1) ( 1)

*i i tu i*

**t Du**

( ) () ()

*i i tu i*

**t Du**

are set along the free surface of all the models.

*2.2.1. Semi-circular canyon model* 

respectively in our calculation.

**Figure 2.** Semi-circular canyon half-space model

integral equation in the domain

**2.2. Numerical examples** 

In Eq. (18)2, the matrix propagators ��(�) (�) and ��(�) (�) are of the dimension *NN*, which is different from the cases when *i* < *L*.

Similarly, for the other layers below the source, i.e., *i L*1, *L*2,, *s*1, we have

$$
\begin{bmatrix} \mathbf{D}\_{uu}^{(i)} \\ \mathbf{D}\_{tu}^{(i)} \end{bmatrix} = - \begin{bmatrix} \mathbf{H}\_{\{i+1\}}^{(i)} + \mathbf{G}\_{\{i+1\}}^{(i)} \mathbf{D}\_{tu}^{(i+1)} & -\mathbf{G}\_{\{i\}}^{(i)} \end{bmatrix}^{-1} \mathbf{H}\_{\{i\}}^{(i)}.\tag{19}
$$

So far, the global matrix propagators are obtained for all the layers except for the source layer. But the problem has not been solved yet. In order to get the final solution of the problem, we need to form the system of equations in the source layer, which can be easily done by setting *i s* in Eq. (12)

$$
\mathbb{E}\left[\mathbf{H}\_{\text{(s)}}^{\text{(s)}} \quad \mathbf{H}\_{\text{(s+1)}}^{\text{(s)}}\right] \begin{bmatrix} \mathbf{u}\_{\text{(s)}} \\ \mathbf{u}\_{\text{(s+1)}} \end{bmatrix} = \begin{bmatrix} \mathbf{G}\_{\text{(s)}}^{\text{(s)}} & -\mathbf{G}\_{\text{(s+1)}}^{\text{(s)}} \\ \end{bmatrix} \begin{bmatrix} \mathbf{t}\_{\text{(s)}} \\ \mathbf{t}\_{\text{(s+1)}} \end{bmatrix} + \mathbf{F}\_{\text{s}}.\tag{20}
$$

Here the subscript *s* in the **F***s* vector just indicates the source layer index. One more step needed to solve the above equations is to substitute *i s*1 and *s*1 separately into Eqs. (17) and (19), which gives

$$\begin{cases} \mathbf{t}\_{(s)} = \mathbf{D}\_{tu}^{(s-1)} \mathbf{u}\_{(s)} \\ \mathbf{t}\_{(s+1)} = \mathbf{D}\_{tu}^{(s+1)} \mathbf{u}\_{(s+1)} \end{cases} . \tag{21}$$

Then Eq. (20) can be rewritten into by Eq. (21)

$$\left[\mathbf{H}\_{\{s\}}^{(s)} - \mathbf{G}\_{\{s\}}^{(s)}\mathbf{D}\_{tu}^{(s-1)} \quad \mathbf{H}\_{\{s+1\}}^{(s)} + \mathbf{G}\_{\{s+1\}}^{(s)}\mathbf{D}\_{tu}^{(s+1)}\right] \bigg|\_{\mathbf{u}} \begin{matrix} \mathbf{u}\_{\{s\}} \\ \mathbf{u}\_{\{s+1\}} \end{matrix} = \mathbf{F}\_{\mathbf{s}'} \tag{22}$$

which has the same number of unknowns and equations. Solving Eq. (22) by Gaussian elimination scheme, we can get the element displacements at the boundary of the source layer so that the element displacements and tractions at the other interfaces can be obtained by the following recursive equations

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

$$\begin{cases} \mathbf{u}\_{(i)} = \mathbf{D}\_{uu}^{(i)} \mathbf{D}\_{uu}^{(i+1)} \cdots \mathbf{D}\_{uu}^{(s-1)} \mathbf{u}\_{(s)}\\ \mathbf{t}\_{(i+1)} = \mathbf{D}\_{tu}^{(i)} \mathbf{u}\_{(i+1)} \end{cases}, \text{ for } i = 1, 2, \dots, s - 1,\tag{23}$$

$$\begin{cases} \mathbf{u}\_{\{i+1\}} = \mathbf{D}\_{uu}^{(i)} \mathbf{D}\_{uu}^{(i+1)} \cdots \mathbf{D}\_{uu}^{(s+1)} \mathbf{u}\_{\{s+1\}}\\ \mathbf{t}\_{\{i\}} = \mathbf{D}\_{tu}^{(i)} \mathbf{u}\_{\{i\}} \end{cases}, \text{ for } i = s+1, \ s+2, \ldots, L. \tag{24}$$

Then the displacement at any observation points in the whole model can be obtained by the integral equation in the domain

$$\mu(\mathbf{r}) + \int\_{\Gamma} T^\*(\mathbf{r}, \mathbf{r}') \mu(\mathbf{r}') \mathrm{d}\Gamma(\mathbf{r}') = \int\_{\Gamma} \mathrm{d}^\*(\mathbf{r}, \mathbf{r}') t(\mathbf{r}') \mathrm{d}\Gamma(\mathbf{r}') + \delta\_{\mathrm{sl}} \mathrm{l} \mathrm{l}^\*(\mathbf{r}, \mathbf{s}) \mathrm{F}(\mathbf{r}, \boldsymbol{\alpha}), \quad \mathbf{r} \in \Omega^{(i)}, \tag{25}$$

which can be transformed into a summation of the element displacements and tractions obtained just now without extra effort. Usually, the displacements at the free surface is of the most interest. By taking the inverse Fourier transform on the above frequency domain solution, we can finally get the time domain solution, i.e., the synthetic seismogram.

#### **2.2. Numerical examples**

132 Wave Processes in Classical and New Solids

Similarly, for the other layers above the source, i.e., *i* 2, , *s*1, we have

*i*

*tu*

 **D**

*i*

*tu*

 **D**

Then Eq. (20) can be rewritten into by Eq. (21)

by the following recursive equations

same reason as that giving rise to Eq. (16), we get

In Eq. (18)2, the matrix propagators ��(�)

different from the cases when *i* < *L*.

done by setting *i s* in Eq. (12)

and (19), which gives

( ) <sup>1</sup> ( ) ( ) ( 1) ( ) ( ) () () ( 1) ( 1) ( ) .

Substituting boundary condition (10) and Eq. (14) with *i L* into Eq. (12) and considering the

<sup>1</sup> () () () () ()

*LLL tu L L* 

> and ��(�) (�)

( ) <sup>1</sup> ( ) ( ) ( 1) ( ) ( ) ( 1) ( 1) () () ( ) .

**H GD G H**

( 1) ( 1) . *s s s s s s s s s s s s s*

.

**u t H H GG F**

*uu i ii i i i i tu i i i*

So far, the global matrix propagators are obtained for all the layers except for the source layer. But the problem has not been solved yet. In order to get the final solution of the problem, we need to form the system of equations in the source layer, which can be easily

> ( ) ( ) () () () () ( ) ( 1) ( ) ( 1)

Here the subscript *s* in the **F***s* vector just indicates the source layer index. One more step needed to solve the above equations is to substitute *i s*1 and *s*1 separately into Eqs. (17)

> ( 1) ( ) ( ) ( 1) ( 1) ( 1)

 

*s s tu s s s tu s*

**t Du**

( ) ( ) ( ) ( 1) ( ) ( ) ( 1)

 

**H GD H G D F**

which has the same number of unknowns and equations. Solving Eq. (22) by Gaussian elimination scheme, we can get the element displacements at the boundary of the source layer so that the element displacements and tractions at the other interfaces can be obtained

, *<sup>s</sup> s ss s s s s s tu s s tu s*

 

() () ( 1) ( 1)

0

**DGH**

 

(�)

Similarly, for the other layers below the source, i.e., *i L*1, *L*2,, *s*1, we have

**H GD G H**

.

**D** (19)

**u t** (20)

**t Du** (21)

( 1)

(22)

*s*

**u**

**u**

(18)

are of the dimension *NN*, which is

**D** (17)

*uu i ii i i i i tu i i i*

( )

**D**

*L uu*

> In the following, we present the calculated synthetic seismograms for some typical models for which exiting results are available. Throughout the following examples, the source time function of the synthetic seismograms is a Ricker wavelet defined as

$$u(t) = \left(2\pi^2 f\_0^2 t^2 - 1\right) \exp\left(-\pi^2 f\_c^2 t^2\right),\tag{26}$$

where *f*0 is the characteristic frequency of the wavelet and *t* is the time. Observation points are set along the free surface of all the models.

#### *2.2.1. Semi-circular canyon model*

The first model is a semi-circular canyon in an elastic homogenous half-space, as shown in Figure 2. Although the shape is too simple for a realistic canyon, the wave reflection and diffraction due to the topographical irregularity in time domain are quite complex, which can serve for an eligible test to the present method and our computation program. Here, we calculate the time-domain responses of the surface along the semi-circular canyon subject to incident SH waves. The shear velocity and mass density are taken as 1 km/s and 1 g/cm3 respectively in our calculation.

**Figure 2.** Semi-circular canyon half-space model

Figures 3 (a) and (b) separately show the antiplane responses of a canyon calculated by the present method due to incident *SH* waves with the angle of 0° and 30°. In the case of vertical incidence, the direct waves keep the same amplitude along the surface of the canyon except near the edges, where the destructive interference occurs. In the case of inclined incidence, the amplitudes of the direct waves inside the canyon decrease toward the rear-side edge. At the horizontal surface near the front-side edge, the peak amplitude becomes larger due to the constructive interference between the direct wave and the wave reflected at the canyon. The two seismograms were calculated by the discrete wavenumber method in [28]. The nice agreement can be seen between our results and those in [28].

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

**Figure 4.** Mountain topography and the surrounding homogeneous half-space

matrix method. Our calculated results show a good agreement with those in [8].

**Figure 5.** Synthetic seismograms for the mountain model in which the point source is placed at (a) *x* = 0,

(a) (b)


OFFSET (km)

*z* = 2 km, (b) *x* = 1, *z* = 2 km. The characteristic frequency is 2.0 Hz.


OFFSET (km)

0

1

2

3

4

5

TIM E (s)

6

7

8

9

10

Figure 5(a) shows the time responses along the mountain topography shown in Fig. 4 when the point source is plated right below the geometric center of the mountain; thus the synthetic seismogram has a symmetric feature about the center point. It can be seen from Figure 5(a) that the first arrival of seismic waves comes late at the locations on the mountain due to the higher altitude. We can also see strong diffracted waves produced by the mountain, i.e. the secondary phase. Figure 5(b) shows the time responses along the same topography as in Figure 5(a) but for the case when the point source is located at *x* = 1, *z* = 2 km. It can be seen from Figure 5(b) that although the first arrival of seismic waves has a nonsymmetric feature as expected, the diffracted waves are symmetric as that for the case shown in Figure 5(a). That indicates that the diffracted waves for both cases propagate along the free surface with the shear wave velocity of the media, and are produced by the mountain topography. The two seismograms were first calculated in [8] by using a global generalized reflection/transmission

0

1

2

3

4

5

TIM E (s)

6

7

8

9

10

**Figure 3.** Time responses of the canyon model due to an *SH* wave incidence: (a) vertical incidence, (b) incidence with an angle of 30°. The characteristic frequency is 1.0 Hz.

#### *2.2.2. Mountain model*

Next, a mountain model with a point source is considered, as shown in Figure 4. The mountain shape topography is described by the equation

$$\mathbf{z}(\mathbf{x}) = \begin{cases} -h\left[\mathbf{1} + \cos\left(\pi \left|\mathbf{x}\right|\right)\right] / 2, \text{for } \left|\mathbf{x}\right| \le 1\\ 0, & \text{for } \left|\mathbf{x}\right| \ge 1 \end{cases} \tag{27}$$

where *h* = 1 km and the elastic constant and the mass density are 1 km/s and 1 g/cm3, respectively.

**Figure 4.** Mountain topography and the surrounding homogeneous half-space

0

1

2

3

4

TIM E (s)

5

6

7

8

Figures 3 (a) and (b) separately show the antiplane responses of a canyon calculated by the present method due to incident *SH* waves with the angle of 0° and 30°. In the case of vertical incidence, the direct waves keep the same amplitude along the surface of the canyon except near the edges, where the destructive interference occurs. In the case of inclined incidence, the amplitudes of the direct waves inside the canyon decrease toward the rear-side edge. At the horizontal surface near the front-side edge, the peak amplitude becomes larger due to the constructive interference between the direct wave and the wave reflected at the canyon. The two seismograms were calculated by the discrete wavenumber method in [28]. The nice

0

1

2

3

4

TIM E (s)

5

6

7

8

**Figure 3.** Time responses of the canyon model due to an *SH* wave incidence: (a) vertical incidence,

(a) (b)

Next, a mountain model with a point source is considered, as shown in Figure 4. The

1 cos / 2,for 1 ( ) , 0, for 1 *h xx*

where *h* = 1 km and the elastic constant and the mass density are 1 km/s and 1 g/cm3,

 

*x*


OFFSET (km)

(27)

(b) incidence with an angle of 30°. The characteristic frequency is 1.0 Hz.

mountain shape topography is described by the equation


OFFSET (km)

*z x*

*2.2.2. Mountain model* 

respectively.

agreement can be seen between our results and those in [28].

Figure 5(a) shows the time responses along the mountain topography shown in Fig. 4 when the point source is plated right below the geometric center of the mountain; thus the synthetic seismogram has a symmetric feature about the center point. It can be seen from Figure 5(a) that the first arrival of seismic waves comes late at the locations on the mountain due to the higher altitude. We can also see strong diffracted waves produced by the mountain, i.e. the secondary phase. Figure 5(b) shows the time responses along the same topography as in Figure 5(a) but for the case when the point source is located at *x* = 1, *z* = 2 km. It can be seen from Figure 5(b) that although the first arrival of seismic waves has a nonsymmetric feature as expected, the diffracted waves are symmetric as that for the case shown in Figure 5(a). That indicates that the diffracted waves for both cases propagate along the free surface with the shear wave velocity of the media, and are produced by the mountain topography. The two seismograms were first calculated in [8] by using a global generalized reflection/transmission matrix method. Our calculated results show a good agreement with those in [8].

**Figure 5.** Synthetic seismograms for the mountain model in which the point source is placed at (a) *x* = 0, *z* = 2 km, (b) *x* = 1, *z* = 2 km. The characteristic frequency is 2.0 Hz.
