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

1. Introduction

A seismogram synthesis is the theoretically calculated ground motion that would be recorded for a given crust structure and seismic source, which is an important tool that can be applied to the study of earthquake source from strong motion records. At the early stage of seismology, the reflectivity method [1] and the generalized ray method [2] were widely used for simulating seismic wave excitation and propagation, where a laterally homogeneous model is considered. The real Earth’s crust, however, is a laterally inhomogeneous structure indeed. To study the site effect of strong ground motion, a stratified structure with irregular interfaces is an eligible model for simulating seismic wave propagation. For the study of seismic waves propagating through a sediment filled basin in the case of rigid grains, the work [3-4] is noteworthy, where a model was proposed to reproduct a nonlinear effect experimentally observed for real seismic waves: site amplification decreases as the amplitude of the incident wave increases.

In the last three decades, the study of seismic wave excitation and propagation in laterally inhomogeneous media has become extensively important in seismology and geophysics, including both analytical solutions and non-analytical results. The analytical solutions, however, only exist for a few special cases [5], such as the semi-cylindrical canyon and alluvial valleys, semi-elliptical canyon, and hemispherical canyon. The methods based on high-frequency approximation [6](such as asymptotic ray theory, Gaussian beam method, and so on) can give a proper solution by including the contribution from the neighboring rays but are only appropriate for handling body wave propagation under relatively high frequency. The methods based on plane wave decomposition (such as discrete wave number method [7-10]) can provide a complete wave field for any finite frequency but usually involve the truncation of matrices and vectors having an infinite dimension. Nowadays, the 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. However, there are still difficuilties in using these methods to simulate wave propagation in models with highly irregular topography since these methods do not explicitly consider the boundary and continuity conditions between different formations. That disables the methods to sufficient accuracy for modeling the wave reflection/transmission across irregular interfaces. Furthermore, the two methods usually require very large computational resources (i.e. a large amount of memory and long calculation time).

To overcome the problem, the reflection/transmission matrices for the fluid-solid irregular interface [20] were developed based on the discrete wavenumber representation of the wave field, i.e. the Aki-Larner method [25]. But their method is only suitable for a long-wavelength irregularity and becomes unstable at very steep interfaces (e.g. a vertical interface), which is due to the intrinsic feature of the method itself. On the other hand, BEM is more suitable than the other methods to model the seismoacoustic scattering at very steep irregular interfaces[26], especially in the cases when better accuracy is required near boundaries and sources. However, the traditional BEM requires a tremendous increase in memory capacity and computational time when it's used to model the wave propagation in multilayered media. Motivated by the work in [19-20], we develop a method based on the combination of the traditional BEM and a global matrix propagator to simulate the seismoacoustic scattering due to an irregular fluid-solid interface. This method takes the advantage of the global matrix propagator to suppress the tremendous increase in computational resources required for simulating wave propagation in multilayered media, but also retains the merits of the boundary element method. This implies that higher-frequency seismoacoustic scattering in complex structures can be calculated with reduced computer resources.

For simplicity, the Chapter is organized as follows: the mathematical formulation for an SH-wave (Shear wave polarized in the horizontal plane) propagation in a multilayered solid is first presented in an easily-understood way in Section 2, where some simple examples are calculated to test the formulation in solid layers; then we gratually go deep into the full formulation for the fluid-solid scattering simulation in Section 3, where two irregular fluid-solid models are used to test the validity of the fluid-solid formulation; based on the formulation in Section3, one of the two fluid-solid models is calculated to show the effects of water layer and the water reverberation in water layer in Section 4; finally summary is drawn in Section 5.

2. SH-wave propagation modeling in mutilayered solids

2.1. Methodology statement

The problem to be studied is illustrated in Figure 1. In this model, there are L homogeneous layers Ω⁽ⁱ⁾ (i=1, 2, …, L) over a half-space, among which the ith layer is bounded by two irregular interfaces Γ_(i-1) and Γ_(i) (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 sth 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 μ⁽ⁱ⁾ and mass density ρ⁽ⁱ⁾. 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.

Figure 1.

A multilayered 2D model with irregular interfaces

Consider the boundary element method problem in the ith 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(r,\omega )+{\rho }^{(i)}{\omega }^{2}u(r,\omega )=-F(r,\omega ){\delta }_{si},$E1

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

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

$F(r,\omega )=F(\omega )\delta (r-s).$E2

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 Γ = Γ_(i)+Γ_(i+1)

$Cu(r)+{\displaystyle \underset{\Gamma }{\int }{T}^{*}(r,r^{\prime })u(r^{\prime })\text{d}\Gamma (r^{\prime })}={\displaystyle \underset{\Gamma }{\int }{U}^{*}(r,r^{\prime })t(r^{\prime })\text{d}\Gamma }(r^{\prime })+{\delta }_{si}{U}^{*}(r,s)F(r,\omega ),$E3

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{array}{l}{U}^{*}(r,r^{\prime })=-\frac{\text{i}}{4{\mu }^{(i)}}{H}_0^2\left(k^{(i)}\left|r-r^{\prime }\right|\right)\\ T^{*}(r,r^{\prime })=\frac{\text{i}{k}^{(i)}}{4}{H}_1^2\left(k^{(i)}\left|r-r^{\prime }\right|\right)(r^{\prime }-r)n,\end{array}$E4

and they satisfy the following equations

$\begin{array}{l}{\mu }^{(i)}{\nabla }^2{U}^{*}(r,r^{\prime })+{\rho }^{(i)}{\omega }^2{U}^{*}(r,r^{\prime })=-\Delta (r^{\prime }-r)\\ T^{*}(r,r^{\prime })={\mu }^{(i)}\frac{\partial }{\partial n}{U}^{*}(r,r^{\prime }),\end{array}$E5

where i is the imaginary unit different from the italic i being index, k⁽ⁱ⁾=ω /c⁽ⁱ⁾ with c⁽ⁱ⁾ being the shear wave velocity in the ith 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 Ω⁽ⁱ⁾ 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² 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 Γ_(i) and Γ_(i+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}{u}_j+{\displaystyle \sum _{k=1}^{2N}{H}_{jk}}{u}_k={\displaystyle \sum _{k=1}^{2N}{G}_{jk}}{t}_k+F_j{\delta }_{si},$E6

where H_jk is an integral of T(r^j, r^k) over field element k, G_jk is an integral of U(r^j, r^k) over field element k, and F_j is the source term. When the source element goes from j = 1 to 2N, one can obtain a system of equations which can be expressed in matrix form

$H{u}^{(i)}=G{t}^{(i)}+F{\delta }_{si},$E7

where the superscript (i) denotes the layer index. Obviously, the column vector u⁽ⁱ⁾ of the dimension 2N×1 consists of the element displacements on both the boundary Γ_(i) and the boundary Γ_(i+1), the same situation applies to the column vector t⁽ⁱ⁾. And the coefficient matrices H and G are of the dimension 2N×2N.

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

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

$t_{(1)}^{(1)}=0,$E8

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

$\begin{array}{l}{u}_{(i+1)}^{(i)}=u_{(i+1)}^{(i+1)}\\ t_{(i+1)}^{(i)}=-t_{(i+1)}^{(i+1)}\end{array}\},$E9

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

$\begin{array}{l}{u}_{(L+1)}^{(L)}=0\\ t_{(L+1)}^{(L)}=0\end{array}\},$E10

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)₂ 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 2N×1.

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

$$\left[\begin{array}{cc}{H}_{(i)}^{(i)}& H_{(i+1)}^{(i)}\end{array}\right]\left[\begin{array}{c}{u}_{(i)}^{(i)}\\ u_{(i+1)}^{(i)}\end{array}\right]=\left[\begin{array}{cc}{G}_{(i)}^{(i)}& G_{(i+1)}^{(i)}\end{array}\right]\left[\begin{array}{c}{t}_{(i)}^{(i)}\\ t_{(i+1)}^{(i)}\end{array}\right]+F{\delta }_{si},$$E11

where ${\mathbf{ }\mathbf{u}}_{\left(i\right)}^{\left(i\right)}$ and ${\mathbf{ }\mathbf{u}}_{\left(i+1\right)}^{\left(i\right)}$ are the element displacements separately corresponding to the upper interface and lower interface of the ith layer, similarly ${\mathbf{ }\mathbf{t}}_{\left(i\right)}^{\left(i\right)}$ and ${\mathbf{ }\mathbf{t}}_{\left(i+1\right)}^{\left(i\right)}$ are the element tractions corresponding to the upper interface and lower interface of the ith 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 ${\mathbf{ }\mathbf{H}}_{\left(i\right)}^{\left(i\right)}$ and${\mathbf{ }\mathbf{H}}_{\left(i+1\right)}^{\left(i\right)}$, ${\mathbf{ }\mathbf{G}}_{\left(i\right)}^{\left(i\right)}$and ${\mathbf{ }\mathbf{G}}_{\left(i+1\right)}^{\left(i\right)}$ which are all of the dimension 2N×N (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)

$$\left[\begin{array}{cc}{H}_{(i)}^{(i)}& H_{(i+1)}^{(i)}\end{array}\right]\left[\begin{array}{c}{u}_{(i)}\\ u_{(i+1)}\end{array}\right]=\left[\begin{array}{cc}{G}_{(i)}^{(i)}& -G_{(i+1)}^{(i)}\end{array}\right]\left[\begin{array}{c}{t}_{(i)}\\ t_{(i+1)}\end{array}\right]+F{\delta }_{si},$$E12

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 ${\mathbf{ }\mathbf{G}}_{\left(i+1\right)}^{\left(i\right)}$ is due to the above-defined positive direction. Equation (12) gives the relationship between the element tractions and the element displacements at the interfaces Γ_(i) and Γ_(i+1) which enclose the ith 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 $D_{uu}^{(i)}$ and $D_{tu}^{(i)}$ are separately defined as ^[18^]

$\{\begin{array}{l}{u}_{(i)}=D_{uu}^{(i)}{u}_{(i+1)}\\ t_{(i+1)}=D_{tu}^{(i)}{u}_{(i+1)}\end{array},i=1,2,\cdots ,s-1,$E13

for the layers above the source, and

$\{\begin{array}{l}{u}_{(i+1)}=D_{uu}^{(i)}{u}_{(i)}\\ t_{(i)}=D_{tu}^{(i)}{u}_{(i)}\end{array},i=s+1,s+2,\cdots ,L,$E14

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

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

$$H_{(1)}^{(1)}{D}_{uu}^{(1)}{u}_{(2)}+H_{(2)}^{(1)}{u}_{(2)}+G_{(2)}^{(1)}{D}_{tu}^{(1)}{u}_{(2)}=0.$$E15

Because u₍₂₎ 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

$\left[\begin{array}{c}{D}_{uu}^{(1)}\\ D_{tu}^{(1)}\end{array}\right]=-{\left[\begin{array}{cc}{H}_{(1)}^{(1)}& G_{(2)}^{(1)}\end{array}\right]}^{-1}{H}_{(2)}^{(1)}.$E16

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

$\left[\begin{array}{c}{D}_{uu}^{(i)}\\ D_{tu}^{(i)}\end{array}\right]=-{\left[\begin{array}{cc}{H}_{(i)}^{(i)}-G_{(i)}^{(i)}{D}_{tu}^{(i-1)}& G_{(i+1)}^{(i)}\end{array}\right]}^{-1}{H}_{(i+1)}^{(i)}.$E17

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{array}{l}{D}_{uu}^{(L)}=0\\ D_{tu}^{(L)}={\left[G_{(L)}^{(L)}\right]}^{-1}{H}_{(L)}^{(L)}\end{array}.$$E18

In Eq. (18)₂, the matrix propagators ${\mathbf{ }\mathbf{H}}_{\left(L\right)}^{\left(L\right)}$ and ${\mathbf{ }\mathbf{G}}_{\left(L\right)}^{\left(L\right)}$ are of the dimension N×N, 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

$\left[\begin{array}{c}{D}_{uu}^{(i)}\\ D_{tu}^{(i)}\end{array}\right]=-{\left[\begin{array}{cc}{H}_{(i+1)}^{(i)}+G_{(i+1)}^{(i)}{D}_{tu}^{(i+1)}& -G_{(i)}^{(i)}\end{array}\right]}^{-1}{H}_{(i)}^{(i)}.$E19

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)

$$\left[\begin{array}{cc}{H}_{(s)}^{(s)}& H_{(s+1)}^{(s)}\end{array}\right]\left[\begin{array}{c}{u}_{(s)}\\ u_{(s+1)}\end{array}\right]=\left[\begin{array}{cc}{G}_{(s)}^{(s)}& -G_{(s+1)}^{(s)}\end{array}\right]\left[\begin{array}{c}{t}_{(s)}\\ t_{(s+1)}\end{array}\right]+F_s.$$E20

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{array}{l}{t}_{(s)}=D_{tu}^{(s-1)}{u}_{(s)}\\ t_{(s+1)}=D_{tu}^{(s+1)}{u}_{(s+1)}\end{array}.$E21

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

$$\left[\begin{array}{cc}{H}_{(s)}^{(s)}-G_{(s)}^{(s)}{D}_{tu}^{(s-1)}& H_{(s+1)}^{(s)}+G_{(s+1)}^{(s)}{D}_{tu}^{(s+1)}\end{array}\right]\left[\begin{array}{c}{u}_{(s)}\\ u_{(s+1)}\end{array}\right]=F_s,$$E22

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

$$\{\begin{array}{l}{u}_{(i)}=D_{uu}^{(i)}{D}_{uu}^{(i+1)}\cdots D_{uu}^{(s-1)}{u}_{(s)}\\ t_{(i+1)}=D_{tu}^{(i)}{u}_{(i+1)}\end{array},\text{ for }i=\text{1},\text{ 2},\mathrm{...},s-\text{1},$$E23

$$\{\begin{array}{l}{u}_{(i+1)}=D_{uu}^{(i)}{D}_{uu}^{(i+1)}\cdots D_{uu}^{(s+1)}{u}_{(s+1)}\\ t_{(i)}=D_{tu}^{(i)}{u}_{(i)}\end{array},\text{ for }i=s+\text{1},s+\text{2},\mathrm{...},L.$$E24

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

$u(r)+{\displaystyle \underset{\Gamma }{\int }{T}^{*}(r,r^{\prime })u(r^{\prime })\text{d}\Gamma (r^{\prime })}={\displaystyle \underset{\Gamma }{\int }{U}^{*}(r,r^{\prime })t(r^{\prime })\text{d}\Gamma }(r^{\prime })+{\delta }_{si}{U}^{*}(r,s)F(r,\omega ),r\in {\Omega }^{(i)},$E25

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

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),$$E26

where f₀ 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/cm³ 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].

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

$z(x)=\{\begin{array}{c}-h\left[1+\cos \left(\pi \left|x\right|\right)\right]/2,\text{for}\left|x\right|\le 1\\ 0,\text{ for}\left|x\right|\ge 1\end{array},$E27

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

Figure 4.

Mountain topography and the surrounding homogeneous half-space

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.

2.2.3. Deep basin model

Consider the model configuration consisting of a deep basin structure, as shown in Figure 6. The basin has a trapezoidal shape with 1 km depth and 10 km width at the surface. The shear wave velocity of the basin-filling sediments and the half-space are 2.5 km/s and 1 km/s, respectively, and the ratio of mass density is set to 1:1.

Figure 6.

Deep basin topography and the surrounding homogeneous half-space

Figure 7.

Time responses along the free surface of the basin in Figure 6 due to a vertically incident SH-wave: (a) characteristic frequency f_c is 0.25 Hz (4 sec), (b) f_c is 0.5 Hz (2 sec)

The time responses at the surface of the basin by a vertically incident SH wave are plotted in Figure 7(a). The characteristic frequency f_c is 0.25 Hz (4 s). This figure shows the horizontally propagating waves generated by the edges of the basin. The amplitude of Love waves (horizontally polarized shear surface waves) is smaller than that of the direct wave. Although these Love waves make the total duration longer, the time interval between the direct wave arrival and the Love wave arrival is less than 10 sec near the edges. Each arrival of reflected Love waves is well separated. Figure 7(b) shows the time responses for the same basin model as in Figure 7(a) but for a Ricker wavelet of characteristic frequency f_c = 0.5 Hz (2 s). Very similar phenomena to that shown in Figure 7(a) can be seen, but the amplitude of the Love wave is larger in this case. The X-shaped pattern of the surface-wave arrivals is more clearly shown for high-frequency incident waves, but the maximum amplitude of the Love wave decreases because energy splits into the fundamental mode and the other higher mode. The two seismograms were calculated first by the discrete wavenumber method in [29]. And nice agreement between our calculated and those in [29] can be observed.

2.2.4. Two-layer model with irregular interface

Finally, a two-layer model with an irregular interface is considered to have a point source, as shown in Figure 8. The irregular interface is described by

$$z(x)=h^{(1)}+\{\begin{array}{c}h,\text{ for }\left|x\right|\le d\\ \left\{1+\cos \left[\pi \left(\left|x\right|-d\right)/2d\right]\right\}h/2,\text{for }d\le \left|x\right|\le 3d\\ 0,\text{ for }\left|x\right|\ge 3d\end{array},$$E28

where h⁽¹⁾ = 2.0 km, h = 0.5 km, and d = 0.5 km. The material parameters are 1 km/s and 1 g/cm³ for the top layer, and 1.5 km/s and 1 g/cm³ for the bottom layer, respectively.

The synthetic seismograms for the two-layer model are shown in Figure 9. It can be seen from Figure 9 that all the seismic phases can be divided into two groups of hyperbolas, namely, a group of hyperbolas with a center point of x = x_s and another group of hyperbolas with a center point of x = 0. The latter ones are due to the reflections by the flat free surface and the flat part of the interface. The former ones, however, are due to the diffractions by the irregular part of the interface. This model was first calculated by the global generalized reflection/transmission matrices method in [8]. Our result shows a nice agreement with that in [8].

Figure 8.

A two-layer model with an irregular interface: the source is at (2 km, 1 km)

Figure 9.

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

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=1, 2, …, L) over a half-space, among which the ith layer is bounded by two irregular interfaces Γ_(i-1) and Γ_(i) (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 sth layer. The properties of an arbitrary solid layer Ω⁽ⁱ⁾ are described by the elastic constants λ⁽ⁱ⁾, μ⁽ⁱ⁾ and mass density ρ⁽ⁱ⁾. 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

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

$$\mu {\nabla }^2\mathbf{u}+\left(\lambda +\mu \right)\nabla \nabla \bullet \mathbf{u}+\rho {\omega }^2\mathbf{u}=-\mathbf{f}$$E29

,

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 u_j(r) on the boundary Γ of the region Ω [30]

$$∆\left(\mathbf{r}\right)u_k\left(\mathbf{r}\right)+{\int }_{\mathrm{\Gamma }}^{}{u}_j\left(\mathbf{r}\mathbf{\text{'}}\right)T_{jk}\left(\mathbf{r},\mathbf{ }\mathbf{r}\mathbf{\text{'}}\right)d\mathrm{\Gamma }\left(\mathbf{r}\mathbf{\text{'}}\right)={\int }_{\mathrm{\Gamma }}^{}{t}_j\left(\mathbf{r}\mathbf{\text{'}}\right)U_{jk}\left(\mathbf{r},\mathbf{ }\mathbf{r}\mathbf{\text{'}}\right)d\mathrm{\Gamma }\left(\mathbf{r}\mathbf{\text{'}}\right)+f_j{U}_{jk}\left(\mathbf{r},\mathbf{ }\mathbf{s}\right)$$E30

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, t_j(r′) are the traction components, and U_jk(r, r′) and T_jk(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 f_j. Application of Eq. (30) in domain Ω⁽ⁱ⁾ and discretization of interfaces Γ_(i) and Γ_(i+1) each into N elements (for simplicity, we use constant element) give rise to

$${\mathbf{H}}^{\left(i\right)}{\mathbf{u}}^{\left(i\right)}={\mathbf{G}}^{\left(i\right)}{\mathbf{t}}^{\left(i\right)}+\mathbf{F}{\delta }_{is}$$E31

,

where H⁽ⁱ⁾ and G⁽ⁱ⁾ are the coefficient matrices obtained by integrating the fundamental solutions T_jk and U_jk over elements at both interfaces Γ_(i) and Γ_(i+1), and u⁽ⁱ⁾ and t⁽ⁱ⁾ are vectors containing element displacements and tractions at both interfaces Γ_(i) and Γ_(i+1), respectively. Obviously, H⁽ⁱ⁾ and G⁽ⁱ⁾ are of the dimension 4N×4N, and u⁽ⁱ⁾ and t⁽ⁱ⁾ are of the dimension 4N×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

$$\left[\begin{array}{cc}{\mathbf{H}}_{\left(i\right)}^{\left(i\right)}& {\mathbf{H}}_{(i+1)}^{\left(i\right)}\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_{\left(i\right)}^{\left(i\right)}\\ {\mathbf{u}}_{(i+1)}^{\left(i\right)}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{G}}_{\left(i\right)}^{\left(i\right)}& {\mathbf{G}}_{(i+1)}^{\left(i\right)}\end{array}\right]\left[\begin{array}{c}{\mathbf{t}}_{\left(i\right)}^{\left(i\right)}\\ {\mathbf{t}}_{(i+1)}^{\left(i\right)}\end{array}\right]+\mathbf{F}{\delta }_{si}$$E32

,

where u⁽ⁱ⁾ and t⁽ⁱ⁾ with subscripts (i) and (i+1) corresponds separately to the element displacements and tractions at Γ_(i) and Γ_(i+1). Similarly, the subscripts (i) and (i+1) in H⁽ⁱ⁾ and G⁽ⁱ⁾ indicate the integration over elements at Γ_(i) and Γ_(i+1), respectively. Obviously, the matrices H⁽ⁱ⁾ and G⁽ⁱ⁾ with subscripts (i) or (i+1) are of the dimension 4N×2N.

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}}_{\left(i\right)}^{\left(i\right)}{\mathbf{u}}_{\left(i\right)}^{}+{\mathbf{H}}_{(i+1)}^{\left(i\right)}{\mathbf{u}}_{(i+1)}^{}={\mathbf{G}}_{\left(i\right)}^{\left(i\right)}{\mathbf{t}}_{\left(i\right)}^{}-{\mathbf{G}}_{(i+1)}^{\left(i\right)}{\mathbf{t}}_{(i+1)}^{}+\mathbf{F}{\delta }_{si}$$E33

,

where only one subscript index is needed to distinct element displacements and tractions at different interfaces, and the minus sign in front of ${\mathbf{G}}_{(i+1)}^{\left(i\right)}$ is due to the defined positive direction for tractions.

The global matrix propagators in the layers above the source are defined as

$$\left\{\begin{array}{c}{\mathbf{u}}_{\left(i\right)}^{}={\mathbf{D}}_{uu}^{\left(i\right)}{\mathbf{u}}_{(i+1)}^{}\\ {\mathbf{t}}_{(i+1)}^{}={\mathbf{D}}_{tu}^{\left(i\right)}{\mathbf{u}}_{(i+1)}^{}\end{array},\mathrm{ }\mathrm{ }\mathrm{ }i=\mathrm{2,3},\cdots ,s-1,\right.$$E34

in which one should note that the layer index i does not start from 1 but 2 since the first layer is a fluid layer in our model. Substitution of Eq. (34) into Eq. (33) yields

$${\mathbf{H}}_{\left(i\right)}^{\left(i\right)}{{\mathbf{D}}_{uu}^{\left(i\right)}\mathbf{u}}_{(i+1)}^{}+{\mathbf{H}}_{(i+1)}^{\left(i\right)}{\mathbf{u}}_{(i+1)}^{}={\mathbf{G}}_{\left(i\right)}^{\left(i\right)}{{{\mathbf{D}}_{tu}^{(i-1)}\mathbf{D}}_{uu}^{\left(i\right)}\mathbf{u}}_{(i+1)}^{}-{\mathbf{G}}_{(i+1)}^{\left(i\right)}{\mathbf{D}}_{tu}^{\left(i\right)}{\mathbf{u}}_{(i+1)}^{}$$E35

.

Since u_(i+1) in Eq. (35) could be the displacement along the lower interface excited by an arbitrary source, the necessary condition of Eq. (35) is that its coefficient matrix equals zero, which leads to

$$\left(\begin{array}{c}{\mathbf{D}}_{uu}^{\left(i\right)}\\ {\mathbf{D}}_{tu}^{\left(i\right)}\end{array}\right)=-{\left[\begin{array}{cc}{\mathbf{H}}_{\left(i\right)}^{\left(i\right)}-{\mathbf{G}}_{\left(i\right)}^{\left(i\right)}{\mathbf{D}}_{tu}^{\left(i-1\right)},& {\mathbf{G}}_{\left(i+1\right)}^{\left(i\right)}\end{array}\right]}^{-1}{\mathbf{H}}_{\left(i+1\right)}^{\left(i\right)},\mathrm{ }\mathrm{ }\mathrm{ }i=\mathrm{2,3},\cdots ,s-1.\mathrm{ }$$E36

Once the global matrix propagator ${\mathbf{D}}_{tu}^{\left(1\right)}$in the fluid layer is obtained, those for the layers above the source can be obtained recursively from Eq. (36). For the moment, we have not known the expression of the global matrix propagator in the fluid layer yet.

Similarly, the global matrix propagators in the layers below the source are defined as

$$\left\{\begin{array}{c}{\mathbf{u}}_{(i+1)}^{}={\mathbf{D}}_{uu}^{\left(i\right)}{\mathbf{u}}_{\left(i\right)}^{}\\ {\mathbf{t}}_{\left(i\right)}^{}={\mathbf{D}}_{tu}^{\left(i\right)}{\mathbf{u}}_{\left(i\right)}^{}\end{array},\mathrm{ }\mathrm{ }\mathrm{ }i=s+1,s+2,\cdots ,L.\right.$$E37

Substituting Eq. (37) into Eq. (33) and considering the same reason as in Eq. (36), we have

$$\left(\begin{array}{c}{\mathbf{D}}_{uu}^{\left(i\right)}\\ {\mathbf{D}}_{tu}^{\left(i\right)}\end{array}\right)=-{\left[\begin{array}{cc}{\mathbf{H}}_{\left(i+1\right)}^{\left(i\right)}+{\mathbf{G}}_{\left(i+1\right)}^{\left(i\right)}{\mathbf{D}}_{tu}^{\left(i+1\right)},& {-\mathbf{G}}_{\left(i\right)}^{\left(i\right)}\end{array}\right]}^{-1}{\mathbf{H}}_{\left(i\right)}^{\left(i\right)},\mathrm{ }\mathrm{ }\mathrm{ }i=L-1,L-2,\cdots ,s+1,$$E38

where ${\mathbf{D}}_{tu}^{\left(L\right)}$ is the global matrix propagator$\mathrm{ }$in the bottom layer which can be obtained by using the radiation condition at infinity, i.e. u^(L+1)=0 and t^(L+1)=0, as follows

$$\left\{\begin{array}{c}{\mathbf{D}}_{uu}^{\left(L\right)}=0\\ {\mathbf{D}}_{tu}^{\left(L\right)}={\left[{\mathbf{G}}_{\left(L\right)}^{\left(L\right)}\right]}^{-1}{\mathbf{H}}_{\left(L\right)}^{\left(L\right)}\end{array}\right.$$E39

.

So far, we have only defined the global matrix propagators for all the solid layers. The global matrix propagator in the fluid layer will be given next.

3.1.2. Formulations in fluid layer

The pressure (i.e. the infinitesimal pressure perturbation) is used as the variable to treat the interaction between the elastic waves and the acoustic waves in the fluid layer [20]. Thus the acoustic wave equation in frequency domain in a homogeneous isotropic fluid layer is given by

$${\nabla }^{2}p+k^{2}p=0$$E40

,

where p is the fluid pressure and k = ω/V is the wavenumber with ω and V being the angular frequency and the acoustic wave velocity, respectively. Using the fundamental solution of Helmholtz equation, we can replace Eq. (40) by the following boundary integral equation for the fluid pressure p(r) on the boundary Γ surrounding the fluid area

$$\mathrm{c}\left(\mathbf{r}\right)p\left(\mathbf{r}\right)+{\int }_{\mathrm{\Gamma }}^{}p\left(\mathbf{r}\mathbf{\text{'}}\right)\frac{\partial G\left(\mathbf{r},\mathbf{ }\mathbf{r}\mathbf{\text{'}}\right)}{\partial n_{\mathbf{r}\mathbf{\text{'}}}}d\mathrm{\Gamma }\left(\mathbf{r}\mathbf{\text{'}}\right)={\int }_{\mathrm{\Gamma }}^{}q\left(\mathbf{r}\mathbf{\text{'}}\right)G\left(\mathbf{r},\mathbf{ }\mathbf{r}\mathbf{\text{'}}\right)d\mathrm{\Gamma }\left(\mathbf{r}\mathbf{\text{'}}\right)$$E41

where $\mathrm{c}\left(\mathbf{r}\right)$generally depends on the local geometry at a point r, and $G\left(\mathbf{r},\mathbf{ }\mathbf{r}\mathbf{\text{'}}\right)$is known as the Green’s function of acoustics which physically represents the effect observed at the point $\mathbf{r}$of a unit source at the point $\mathbf{r}\mathbf{\text{'}}$[31].$q\left(\mathbf{r}\mathbf{\text{'}}\right)=\frac{\partial p\left(\mathbf{r}\mathbf{\text{'}}\right)}{\partial n_{\mathbf{r}\mathbf{\text{'}}}}$, where the terminology $\frac{\partial *}{\partial n_{\mathbf{r}\mathbf{\text{'}}}}$ represents the partial derivative of the function * with respect to the unit outward normal at the point $\mathbf{r}\mathbf{\text{'}}$on the boundary Γ. Applying Eq. (41) in the uppermost layer Ω⁽¹⁾ and discretizing the interfaces Γ₍₁₎ and Γ₍₂₎ each into N elements (still constant element), we yield the matrix equation of the form

$${\mathbf{H}}_{\left(2\right)}^{\left(1\right)}{\mathbf{p}}_{\left(2\right)}^{}={\mathbf{G}}_{\left(1\right)}^{\left(1\right)}{\mathbf{q}}_{\left(1\right)}^{}+{\mathbf{G}}_{\left(2\right)}^{\left(1\right)}{\mathbf{q}}_{\left(2\right)}^{}$$E42

,

where p and q with subscripts ₍₁₎ and ₍₂₎ correspond separately to the element pressure and its normal derivative at Γ₍₁₎ and Γ₍₂₎ (p₍₁₎ in Eq. (42) is suppressed since it’s always equal to zero due to the free surface condition). Obviously, they are the vectors of the dimension N×1. Follow the rule in the solid layers, we first define the global matrix propagators in the fluid layer as

$$\left\{\begin{array}{c}{\mathbf{q}}_{\left(1\right)}^{}={{\mathbf{D}}_{qq}^{\left(1\right)}\mathbf{q}}_{\left(2\right)}^{}\\ {\mathbf{p}}_{\left(2\right)}^{}={\mathbf{D}}_{pq}^{\left(1\right)}{\mathbf{q}}_{\left(2\right)}^{}\end{array}\right.$$E43

,

where the first definition is actually unnecessary since we focus on the solution at the fluid-solid interface and ignore the normal derivative of the pressure at the free surface. Substitution of Eq. (43) into Eq. (42) yields

$${\mathbf{H}}_{\left(2\right)}^{\left(1\right)}{\mathbf{D}}_{pq}^{\left(1\right)}{\mathbf{q}}_{\left(2\right)}^{}={\mathbf{G}}_{\left(1\right)}^{\left(1\right)}{{\mathbf{D}}_{qq}^{\left(1\right)}\mathbf{q}}_{\left(2\right)}^{}+{\mathbf{G}}_{\left(2\right)}^{\left(1\right)}{\mathbf{q}}_{\left(2\right)}^{}$$E44

.

Since q₍₂₎ in Eq. (44) could be the normal derivative of the pressure at the fluid-solid interface excited by an arbitrary source, the necessary condition of Eq. (44) is that its coefficient matrix equals zero, which leads to

$$\left(\begin{array}{c}{\mathbf{D}}_{qq}^{\left(1\right)}\\ {\mathbf{D}}_{pq}^{\left(1\right)}\end{array}\right)={\left[\begin{array}{cc}{-\mathbf{G}}_{\left(1\right)}^{\left(1\right)},& {\mathbf{ }\mathbf{H}}_{\left(2\right)}^{\left(1\right)}\end{array}\right]}^{-1}{\text{G}}_{\left(2\right)}^{\left(1\right)}$$E45

,

where one point noted here is that the global matrix propagators obtained in Eq. (45) have the dimension N×N which is different from those in Eqs. (34) and (37). Therefore, we need to extend the global matrix propagators defined for the fluid layer into the same dimension as those defined for the solid layers, i.e., from N×N to 2N×2N, by using the continuity conditions at the fluid-solid interface. They are : 1) the tangential traction (not stress, be careful with the terminology) is zero; 2) the normal traction is continuous; and 3) the normal displacement is continuous, which can be expressed in frequency domain as follows [20, 32]

$$\begin{array}{c}\mathbf{t}\bullet \mathbf{s}=0\\ \mathbf{t}\bullet \mathbf{n}=-p\\ \mathbf{u}\bullet \mathbf{n}=-u_{\mathrm{f}}=-\frac{1}{{\rho }_{\mathrm{f}}{\omega }^2}q\end{array}$$E46

,

where $\mathbf{u}$and $\mathbf{t}$are the element displacements and tractions of the solid layer at the fluid-solid interface. And the $\mathbf{n}$and $\mathbf{s}$are the unit outward normal vector and the inplane tangential vector of the solid layer neighboring the fluid layer. ρ_f is the mass density of the fluid layer. The displacement of the fluid u_f is along the normal direction but points to the neighboring solid layer (i.e., the opposite direction of the normal vector n), which leads to the minus sign in front of u_f.

Now, the global matrix propagators obtained in Eq. (45) can be extended into the same dimension as those defined in the solid layers, i.e. 2N×2N, by

$${\mathbf{D}}_{tu}^{\left(1\right)}={\rho }_{\text{f}}{\omega }^2\left[\begin{array}{cc}{\mathbf{D}}_{pq}^{\left(1\right)}\bullet {\text{N}}_1^{\mathrm{e}}{\bullet \left({\text{N}}_1^{\mathrm{e}}\right)}^{\mathrm{T}}& {\mathbf{D}}_{pq}^{\left(1\right)}\bullet {\text{N}}_2^{\mathrm{e}}{\bullet \left({\text{N}}_1^{\mathrm{e}}\right)}^{\mathrm{T}}\\ {\mathbf{D}}_{pq}^{\left(1\right)}\bullet {\text{N}}_1^{\mathrm{e}}{\bullet \left({\text{N}}_2^{\mathrm{e}}\right)}^{\mathrm{T}}& {\mathbf{D}}_{pq}^{\left(1\right)}\bullet {\text{N}}_2^{\mathrm{e}}{\bullet \left({\text{N}}_2^{\mathrm{e}}\right)}^{\mathrm{T}}\end{array}\right]$$E47

,

where the matrices ${\text{N}}_1^{\mathrm{e}}$ and ${\text{N}}_2^{\mathrm{e}}$are defined as follows:

${\text{N}}_1^{\mathrm{e}}=\left[\begin{array}{c}\begin{array}{cc}{n}_1^{{\text{e}}_1}& n_1^{{\text{e}}_2}\\ n_1^{{\text{e}}_1}& n_1^{{\text{e}}_2}\end{array} \begin{array}{cc}\cdots & n_1^{{\text{e}}_N}\\ \cdots & n_1^{{\text{e}}_N}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ n_1^{{\text{e}}_1}& n_1^{{\text{e}}_2}\end{array} \begin{array}{cc}\ddots & ⋮\\ \cdots & n_1^{{\text{e}}_N}\end{array}\end{array}\right]$,

$${\text{N}}_2^{\mathrm{e}}=\left[\begin{array}{c}\begin{array}{cc}{n}_2^{{\text{e}}_1}& n_2^{{\text{e}}_2}\\ n_2^{{\text{e}}_1}& n_2^{{\text{e}}_2}\end{array} \begin{array}{cc}\cdots & n_2^{{\text{e}}_N}\\ \cdots & n_2^{{\text{e}}_N}\end{array}\\ \begin{array}{cc}⋮& ⋮\\ n_2^{{\text{e}}_1}& n_2^{{\text{e}}_2}\end{array} \begin{array}{cc}\ddots & ⋮\\ \cdots & n_2^{{\text{e}}_N}\end{array}\end{array}\right]$$E48

,

in which $n_1^{{\mathrm{e}}_i}$ and $n_2^{{\mathrm{e}}_i}$are the two components of the unit outward normal vector of the neighboring solid layer at the element ‘e_i’, i = 1, 2, 3,…, N. Note that the multiplication between matrices in Eq. (47) is operated in the sense of ‘element to element’.

3.1.3. Solution in source layer

We have clearly obtained the global matrix propagators for all the layers except for the source layer which is the place to finally solve the problem. Setting i = s in Eq. (33), we have [19]

$${\mathbf{H}}_{\left(s\right)}^{\left(s\right)}{\mathbf{u}}_{\left(s\right)}^{}+{\mathbf{H}}_{(s+1)}^{\left(s\right)}{\mathbf{u}}_{(s+1)}^{}={\mathbf{G}}_{\left(s\right)}^{\left(s\right)}{\mathbf{t}}_{\left(s\right)}^{}-{\mathbf{G}}_{\left(s+1\right)}^{\left(s\right)}{\mathbf{t}}_{\left(s+1\right)}^{}+\mathbf{F}$$E49

,

in which both the element displacements and tractions are unknowns, however, they cannot be solved from the present form of Eq. (49) yet. From the aforementioned definitions of the global matrix propagators, i.e., take i = s−1 in Eq. (34) and i = s+1 in Eq. (37), we can get the relationship between the traction and displacement at the sth and (s+1)th interfaces as [19]

$${\mathbf{t}}_{\left(s\right)}^{}={\mathbf{D}}_{tu}^{(s-1)}{\mathbf{u}}_{\left(s\right)}^{},\mathrm{ }{\mathbf{t}}_{(s+1)}^{}={\mathbf{D}}_{tu}^{(s+1)}{\mathbf{u}}_{(s+1)}^{}$$E50

.

Substitution of Eq. (50) into Eq. (49) gives rise to [19]

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

,

which has the same number of unknowns and equations and can be solved without question.

Once the displacements at the interfaces enclosing the source layer are solved, the displacements at the fluid-solid interface can be obtained by

$${\mathbf{u}}_{\left(2\right)}={\mathbf{D}}_{uu}^{\left(2\right)}{\mathbf{D}}_{uu}^{\left(3\right)}\cdots {\mathbf{D}}_{uu}^{(s-1)}{\mathbf{u}}_{\left(s\right)}$$E52

,

which is the displacement solution in frequency domain at the fluid-solid interface. By applying the inverse Fourier transform to Eq. (52), one can get the solution in time domain, i.e., the synthetic seismograms. A general procedure of applying the introduced approach to calculate a seismic ground motion in an irregularly multilayered media lain by a fluid layer is summarized as follows.

Step 1. Calculate the global matrix propagators in the uppermost layer, i.e. the fluid layer, by Eq. (45).

Step 2. Extend the global matrix propagators obtained in Step 1 by Eq. (47).

Step 3. Use the extended global matrix propagator in Step 2 to calculate the global matrix propagators in the solid layers above the source layer by Eq. (36).

Step 4. In case of a plane wave incidence problem or a seismic source located at the bottom layer, skip this step and jump to Step 5. Otherwise, calculate the global matrix propagators in the solid layers below the seismic source layer by Eqs. (38-39).

Step 5. Substitute all the global matrix propagators obtained in the above steps into the simultaneous matrix equation in the seismic source layer and solve the problem by Eq. (51).

Step 6. Obtain the displacement solutions at the fluid-solid interface by Eq. (52).

Step 7. Execute the Inverse Fourier Transform on the displacement solutions obtained in Step 6 to get the seismic ground motion in time domain finally.

The above-mentioned procedure is for the case when the fluid layer is the top layer. In case that the fluid layer is not the uppermost layer of the model under consideration, such as a fluid layer surrounded by two solid layers, the treatment similar to Eq. (47) can be also applied at the two fluid-solid interfaces. The only difference of the solution procedure is the order of calculating the global matrix propagators in the fluid layer and the solid layer, the details of which are suppressed here.

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.

3.2. Numerical examples

In this subsection, we conduct the numerical implementation of the fluid-solid formulation of the introduced method. Three examples are calculated for validity testing, considering both plane wave incidence and SV-wave source excitation. Although the models used in this subsection seem very simple, they are often used as examples for validity testing of some seismogram synthetic method.

We first calculate the case of a plane P-wave (Primary or pressure wave) vertically incident onto irregularly layered fluid-solid structures in the two models, as shown in Figure 11. The reason for the selection of the two models is because the synthetic results of seismoacoustic scattering due to the irregular fluid-solid interface are available in [20]. The existing results for the two models calculated by the reflection/transmission matrix method in [20] were validated in detail by the finite difference method [20], which concluded that the reflection/transmission matrix method could give more accurate results. Therefore, we will compare our results with those calculated by the reflection/ transmission matrix method.

Figure 11.

The two models used to test the fluid-solid formulation: (a) Model 1; (b) Model 2 [20]

In order to make our comparison convincing, not only are the model parameters set to be the same as those used in [20] but also the time function of the input plane wave is selected as

$$f\left(t\right)=\left\{\begin{array}{c}\sqrt{\alpha /\pi }\mathrm{e}\mathrm{x}\mathrm{p}\left(-\alpha t^2\right),\mathrm{ }-T_s/2\le t\le T_s/2\\ \mathrm{ }\mathrm{ }0,\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }t\ge T_s/2\mathrm{ },\mathrm{ }\mathrm{ }t\le -T_s/2\mathrm{ }\end{array}\right.$$E53

,

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 the Eq. (88) in [20] and its shape and spectrum are shown in Figure 12.

Figure 12.

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

3.2.1. Model 1 subjected to a plane wave incidence

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 large-amplitude 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 be well and clearly observed in our calculated seismograms.

3.2.2. Model 2 subjected to a plane wave incidence

Figure 14 shows the displacements of both horizontal and vertical components at the fluid-solid 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.

Figure 13.

Time responses along the fluid-solid interface of Model 1: the one on the left calculated in [20], and the one on the right calculated by the present method

Figure 14.

The same as Figure 13, but for Model 2

3.2.3. Model 1 subjected to an SV-wave point source

The above two testing calculations show the validity of the present method and the correctness of our computational program, respectively. However, the two numerical examples are only plane wave incidence cases. To validate the present method further, we show one more example using an SV-wave (Shear wave polarized in the vertical plane) excitation source as [23]

$$M_{xz}=-M_{zx}=-Mf\left(t\right)\delta \left(x\right)\delta \left(z\right)$$E54

,

where the source function f(t) is a Ricker wavelet with central frequency 0.25 Hz. The model used for the calculation is similar to Model 1 but with a width of 80 km in total, and the basin-like fluid-solid interface having a width of 60 km and thickness of 2.1 km. The material properties of the model is little different from Model 1, details of which can refer to [23]. The SV-wave source is located at (distance, depth) = (−0.1 km, 4.1 km). As for this case, the results calculated by the discrete wavenumber method are available in [23] so that we can make a full waveform comparison.

Figure 15.

Time responses along the basin-like fluid-solid interface: the one on the left calculated in [23] with a time offset of 8 seconds, and the one on the right calculated by the present method.

Figure 15 shows the displacements of both horizontal and vertical components along the basin-like fluid-solid interface due to the SV-wave source presented in Eq. (54). The one on the left is calculated by the discrete wavenumber method in [23]. And the one on the right is calculated by the present method. It can be clearly seen that the large-amplitude multiple reflections in the fluid-solid basin part are well observed in the synthetic waveforms calculated by both the present method and the discrete wavenumber method. Furthermore, the first scattered P-waves, the first scattered Rayleigh waves, and the secondary scattered Rayleigh waves, which are separately indicated as P, R and R2 in the plot on the left, can be well and clearly observed in our calculated seismograms. Considering the time offset used in [23], the very good agreement between the results by the two methods further confirms the validity of the present method and the correctness of the fluid-solid formulation.

Up to now, we have validated the fluid-solid formulation by three examples, which results in the conclusion that the introduced approach can accurately cover the seismoacoustic scattering due to an irregular fluid-solid interface. In the next Section, we will show how to use the introduced approach to simulate the effects of a fluid layer on the synthetic seismograms and the water reverberation by three preliminary examples, respectively.

4. Water effects and water reverberation modeling

4.1. Water effects on seismogram synthesis

In this Subsection, we show the effect of a fluid layer on synthetic seismogram by using one of the testing models. In practice, seismologists or geophysicists are interested in the seismic ground motion at land in gulf areas where the fluid layer plays an important role in the recorded seismograms. Therefore, it is necessary and significant to simulate the seismoacoustic scattering in irregularly multilayered elastic media lain by a fluid layer due to some scenario earthquake event. Here, we select the first model used in Subsection 3.2 to show water layer effects due to the two reasons: a) this model is very simple and its calculated results are relatively easy to interpret; b) this model is kind of close to a practical gulf area although the solid layer is considered as an elastic half space. For the same model, we take into account two cases: the first case of a vertically incident plane P wave, as shown in Figure 16, the second case of an explosive source located in the solid half-space. We will show the calculated results of the two cases one by one in the following.

Figure 16.

The model used to show the effect of a fluid layer on the synthetic seismograms due to a vertically incident plane P wave: (a) with a fluid layer; (b) without a fluid layer.

4.1.1. Plane wave incidence

Figure 17 shows the time responses of the horizontal and vertical motions along the irregular solid interface due to a vertically incident plane P wave for the models with and without a fluid layer, respectively. It can be seen from Figure 17 that the direct waves appearing on the vertical component keep the same amplitude along the solid interface, regardless of the existence of the uppermost fluid layer. Outside the irregular part of the interface, later arrivals on the vertical component are mainly due to the wave reflection at the upper part of the irregular interface but they seem to contain some contribution of the diffracted waves since their amplitude changes very slowly. Judging from their particle motions and apparent velocities, we can say that they appear to become Rayleigh waves soon after the departure from the edges of the irregular part of the interface. Comparison between the results for the models with and without the fluid layer can be clearly seen in Figure 17. For the model shown in Figure 16(a), although the later arrivals on the vertical component inside the irregular part of the interface seem complicated, the multiple later arrivals outside the irregular part of the interface due to the multiple reflections caused by the fluid layer are clearly observed. In the case of the absence of the fluid layer, i.e. the model shown in Figure 16(b), there is only one later arrival on the vertical component, even the previous complicated later arrivals appearing inside the irregular part of the interface in the case of the presence of the fluid layer completely disappear, which appears to be very interesting results.

Figure 17.

Time responses along the irregular solid interface in the models shown in Figure 16 due to a vertically incident plane P wave: (a) with a fluid layer; (b) without a fluid layer.

The earlier arrivals on the horizontal component generated by the impact and subsequent reflection of the incident plane P wave at the irregular part of the interface grows as the incident wave propagates upward and reaches its maximum at the edges, regardless of the existence of the uppermost fluid layer. Those diffracted waves gradually separate into P waves and Rayleigh waves with the increasing distance away from the irregular part edges. The difference between the results on the horizontal component for the models with and without the fluid layer can also be clearly seen, no matter inside the irregular part or outside the irregular part of the interface. Inside the irregular part of the interface the diffracted waves for the model with the fluid layer have much more complicated features than those for the model without the fluid layer. Different from the vertical component, multiple diffracted waves inside the irregular part of the model with the fluid layer can be clearly observed as well. For the model without the fluid layer, however, no multiple diffracted waves can be observed inside the irregular part of the surface. Outside the irregular part of the interface, the multiple arrivals of P waves and Rayleigh waves on the horizontal component can be seen for the model with the fluid layer, although it is difficult to distinguish them clearly. While for the model without the fluid layer, no more later arrivals present on the horizontal component, which is similar to the observation on the vertical component. But the amplitude of the first diffracted waves outside the irregular part of the surface of the model without the fluid layer (i.e. P waves) is a little bit larger than that of the model with the fluid layer, which is due to the energy dissipation by the fluid layer.

4.1.2. Explosive source

Besides the plane wave incidence, let us see what happens if we have an explosive source located in the solid layer. The source time function is a Ricker wavelet with central frequency 0.25 Hz. Figure 18 shows the time responses of the horizontal and vertical motions along the irregular solid interface due to an explosive source located at (distance, depth) = (32 km, 4 km) in the solid layer for Model 1 with and without a fluid layer, respectively. It can be seen from Figure 18 that the direct arrivals appearing along the irregular interface come from the explosive source directly and separate into two waves with the distance away from the center, regardless of the existence of the uppermost fluid layer. Those can be recognized as P and Rayleigh waves. The main difference of the results for the model with and without the fluid layer exists in the later multiple arrivals, which can be clearly observed from the comparison between Figure 18 (a) and (b). For the model with the fluid layer, the multiple reflections inside the basin-like fluid-solid interface can be clearly seen to generate the later multiple Rayleigh wave arrivals along the solid surface outside the irregular interface. While for the model without the fluid layer, no such multiple arrivals appear since it’s just a half-space. Those phenomena resemble the case of plane wave incidence.

The above two examples clearly show the difference between the synthetic seismograms for the model with and without a fluid layer in the cases of plane wave incidence and explosive source excitation, from which we can say that the method presented in this Chapter can be used to study the effect of fluid layers on the seismic ground motions at land and ocean bottom seismic observations. Next, we will consider the case when we put the explosive source in a fluid layer and see how the introduced method can be used to simulate the water reverberation.

Figure 18.

Time responses along the irregular solid interface in Model 1 due to an explosive source located at (distance, depth) = (32 km, 4 km): (a) with a fluid layer; (b) without a fluid layer.

4.2. Water reverberation modeling

4.2.1. Formulation when source is located in fluid layer

So far, we have given the formulation for seismic excitation in a multi-layered solid half-space overlain with a fluid layer. Readers may ask if the method can be applied to obtain strong ground motion at land when we have the source located in the fluid layer. The answer is definitely positive. For doing so, there are two key steps. The first one is to form the solution matrix equations in the fluid layer, which can be done from Eq. (42) as follows

$${\mathbf{H}}_{\left(2\right)}^{\left(1\right)}{\mathbf{p}}_{\left(2\right)}^{}={\mathbf{G}}_{\left(1\right)}^{\left(1\right)}{\mathbf{q}}_{\left(1\right)}^{}+{\mathbf{G}}_{\left(2\right)}^{\left(1\right)}{\mathbf{q}}_{\left(2\right)}^{}+\mathbf{F}$$E55

,

where p and q with subscripts ₍₂₎ are connected by Eq. (43). In order to solve Eq. (55), we need to use the global matrix propagator ${\mathbf{D}}_{pq}^{\left(1\right)}$ defined in Eq. (43) to change Eq. (55) into

$$\left(-\begin{array}{cc}{\mathbf{G}}_{\left(1\right)}^{\left(1\right)}& {\mathbf{H}}_{\left(2\right)}^{\left(1\right)}{\mathbf{D}}_{pq}^{\left(1\right)}-{\mathbf{G}}_{\left(2\right)}^{\left(1\right)}\end{array}\right)\left(\begin{array}{c}{\mathbf{q}}_{\left(1\right)}^{}\\ {\mathbf{q}}_{\left(2\right)}^{}\end{array}\right)=\mathbf{F}$$E56

,

which has the same number of unknowns and equations and can be solved without problem. What is left is to get the global matrix propagator${\mathbf{D}}_{pq}^{\left(1\right)}$, which is the other key step. We understand that the global matrix propagator ${\mathbf{D}}_{tu}^{\left(1\right)}$ in the solid layer right neighboring the fluid layer can be obtained recursively from the lowermost layer by using Eqs. (38-39). Then the continuity conditions at the fluid-solid interface Eq. (43) are used to obtain ${\mathbf{D}}_{pq}^{\left(1\right)}$ from ${\mathbf{D}}_{tu}^{\left(1\right)}$ as

$${\mathbf{D}}_{pq}^{\left(1\right)}=\frac{1}{{\rho }_{\text{f}}{\omega }^2}\mathrm{\Sigma }{\mathbf{D}}_{tu}^{\left(1\right)}/\left[{\text{N}}_1^{\mathrm{e}}{\bullet \left({\text{N}}_1^{\mathrm{e}}\right)}^{\mathrm{T}}+{\text{N}}_2^{\mathrm{e}}{\bullet \left({\text{N}}_1^{\mathrm{e}}\right)}^{\mathrm{T}}+{\text{N}}_1^{\mathrm{e}}{\bullet \left({\text{N}}_2^{\mathrm{e}}\right)}^{\mathrm{T}}+{\text{N}}_2^{\mathrm{e}}{\bullet \left({\text{N}}_2^{\mathrm{e}}\right)}^{\mathrm{T}}\right]$$E57

,

where the matrices ${\text{N}}_1^{\mathrm{e}}$ and ${\text{N}}_2^{\mathrm{e}}$are defined as Eq. (48). The symbol$\mathrm{ }\mathrm{\Sigma }$ in Eq. (57) means the summation of the four quarters of${\mathbf{D}}_{tu}^{\left(1\right)}$, and the multiplication and division of the matrices are operated in the sense of ‘element to element’. The general solution steps can be summarized into the following.

Step 1. Calculate the global matrix propagators ${\mathbf{D}}_{tu}^{\left(1\right)}$ in the lowermost layer by Eq. (39).

Step 2. Calculate the global matrix propagators in all the solid layers above the bottom layer by Eqs. (38-39).

Step 3. Reduce the global matrix propagators obtained in the solid layer right neighboring the fluid layer into the one in the fluid layer by Eq. (57).

Step 4. Substitute the reduced global matrix propagator ${\mathbf{D}}_{pq}^{\left(1\right)}$ obtained in Step 3 into the solution matrix equation (56) and solve for the normal derivative of pressure.

Step 5. Obtain the pressure solutions at the fluid-solid interface by Eq. (43). Then obtain the traction solutions at the fluid-solid interface by the first two in Eq. (46).

Step 6. Obtain the displacement solutions at the fluid-solid interface by${\mathbf{D}}_{tu}^{\left(1\right)}$.

Step 7. Execute the Inverse Fourier Transform on the displacement solutions obtained in Step 6 to get the ground motion in time domain finally.

4.2.2. Numerical example

First consider the same model as Model 1, with an explosive source located at (distance, depth) = (32 km, 0.5 km) and receivers along the fluid-solid interface. The time function of the source is a Ricker wavelet with central frequency 0.25 Hz. The calculated time responses of displacements are shown in Figure 19(a), from which we can clearly see the multiple reflections inside the water basin and the generated multiple arrivals along the solid surface ourside the basin part. The second calculation taken into account is also for Model 1 but with the water width enlarged to 30 km, results of which are shown in Figure 19(b) and the time responses of the displacements along the fluid-solid interface ressembles those in Figure 19(a). The main difference between the two calculated results exists in the time difference between the multiple arrivals, which is due to the enlarged width of the water basin area. This example implies that the present method can be used to simulate the water revibaration in the sea, which is of importantly practical significance to deep ocean acoustic experiments and so on.

Figure 19.

Time responses along the fluid-solid interface due to an explosive source located at (distance, depth) = (32 km, 0.5 km) in: (a) Model 1 and (b) Model 1 with a water width of 30 km.

5. Summary

In this Chapter, we presented an efficient approach based on the combination of the traditional boundary element method and the global matrix propagators for seismoacoutic scattering simulation in multilayered fluid-solid media. To make readers understand the introduced method easily, we first gave the mathematical formulation for SH-wave propagation simulation in multilayered solids, followed by some examples to test the validity of the formulation in solids. Then, we gradually went deep into the seismoacoustic scattering simulation due to an irregular fluid-solid interface, where the fluid pressure was used as the variable inside the fluid layer and the global matrix propagators defined in the fluid layer were extended successfully to connect the pressure wave with the elastic wave by using the standard continuity conditions at the fluid-solid interface (i.e. zero tangential traction, continuity of normal traction, and continuity of normal displacement). The synthetic waveforms in time domain for some selected models calculated by the present method agree well with those calculated by the reflection/transmission matrix method and those calculated by the discrete wavenumber method. The effects of a fluid layer on the synthetic seismograms can be exactly covered by the present method, and the water reverberation in the sea can also be simulated as well. The introduced method is especially suitable to simulate seismoacoustic scattering in a multilayered elastic structure overlain by a fluid layer. Since the global matrix propagators can be calculated recursively, the computer memory required for a multilayered model is the same as that for a two-layer model. In case of the application of dynamic allocation of matrices and saving the global matrix propagators on hard drive, the array size assigned for calculating a two-layer model is sufficient for a multilayered model. In that case, the advantages of the boundary element method can be preserved, which implies that seismoacoustic scattering synthesis due to a high-frequency excitation can be modelled with reduced computer resources.

Furthermore, if a calculated model is only partially modified, such as increasing the number of layers below the uppermost fluid layer in the case of a plane wave incidence or changing the free surface profile in the case of a point source excitation, not all the matrices need recalculating. That is due to the two merits of the present method: the global matrix propagators for the layers above the source are calculated downwards; and the global matrix propagators for the layers below the source are independent of the rest of the model.

Acknowledgement

The support from the State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, P. R. China is greatly appreciated. The financial support from the Global COE Program entitled “International Urban Earthquake Engineering Center for Mitigating Seismic Mega Risk”, Tokyo Institute of Technology, Japan is gratefully acknowledged. The author (Qian) also appreciates very much the help of Prof. Chen X.F. from USTC and Prof. Ge Z. from Peking University for providing us with their program designed for multilayered solid media. And the author (Qian) thanks very much Prof. Fu L.Y. and Dr. Yu G.X. from the Institute of Geology and Geophysics, Chinese Academy of Science, P. R. China for their very much useful discussion and help on the application of BEM as well. The author (Qian) would also like to thank Prof. H. Takenaka from Kyushu University, Japan for giving the very constructive suggestions and comments on the research of seismoacoustic scattering simulation which help improve the introduced approach very much. The author (Qian) would like to appreciate the support from Prof. K. Kishimoto at Tokyo Institute of Technology, Japan. Finally, the author (Qian) would like to sincerely appreciate Prof. Sohichi Hirose from Tokyo Institute of Technolody, Japan for his so much inspiring discussion and comments on writing the programs of boundary element method.