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  and the generalized ray method  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 , such as the semi-cylindrical canyon and alluvial valleys, semi-elliptical canyon, and hemispherical canyon. The methods based on high-frequency approximation (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 , finite element , pseudo-spectral method , and spectral element method ) 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 . 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  proposed an improved block Gaussian elimination scheme and Bouchon et al.  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 , 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 . 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)  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  were developed based on the discrete wavenumber representation of the wave field, i.e. the Aki-Larner method . 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, 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
Consider the boundary element method problem in the
Suppose the seismic source distribution consists of a simple point source at a position vector s, i.e.,
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
where r and r′ are the position vectors of “field point” and “source point” on the boundary Γ, respectively. The coefficient
and they satisfy the following equations
where i is the imaginary unit different from the italic
To solve Eq. (3), we need to do discretization on the boundary. For simplicity, we discretize the boundary Γ(
where the superscript (
The boundary and continuity conditions for the whole model under consideration are:
the traction-free condition on the free surface,
the continuity conditions of displacement and traction at the inner interfaces,
the radiation boundary conditions imposed on the far-field behavior at infinity,
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
where and are the element displacements separately corresponding to the upper interface and lower interface of the
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)
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 Γ(
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  with minor modification.
The global matrix propagators and are separately defined as 
for the layers above the source, and
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
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
Similarly, for the other layers above the source, i.e.,
In Eq. (18)2, the matrix propagators and are of the dimension
Similarly, for the other layers below the source, i.e.,
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
Here the subscript
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
Then the displacement at any observation points in the whole model can be obtained by the integral equation in the domain
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
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.
Figures 3 (a) and (b) separately show the antiplane responses of a canyon calculated by the present method due to incident
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
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
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.
The time responses at the surface of the basin by a vertically incident
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
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
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
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
where only one subscript index is needed to distinct element displacements and tractions at different interfaces, and the minus sign in front of is due to the defined positive direction for tractions.
The global matrix propagators in the layers above the source are defined as
Once the global matrix propagator 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
where is the global matrix propagatorin the bottom layer which can be obtained by using the radiation condition at infinity, i.e. u(
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 . Thus the acoustic wave equation in frequency domain in a homogeneous isotropic fluid layer is given by
where generally depends on the local geometry at a point
where p and q with subscripts (1) and (2) correspond separately to the element pressure and its normal derivative at Γ(1) and Γ(2) (p(1) 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
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
Since q(2) 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
where one point noted here is that the global matrix propagators obtained in Eq. (45) have the dimension
where and are the element displacements and tractions of the solid layer at the fluid-solid interface. And the and are the unit outward normal vector and the inplane tangential vector of the solid layer neighboring the fluid layer.
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. 2
where the matrices and are defined as follows:,
in which and are the two components of the unit outward normal vector of the neighboring solid layer at the element ‘e
3.1.3. Solution in source layer
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
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
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.
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
In order to make our comparison convincing, not only are the model parameters set to be the same as those used in  but also the time function of the input plane wave is selected as
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
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
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 
where the source function
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 . 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
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
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
The earlier arrivals on the horizontal component generated by the impact and subsequent reflection of the incident plane
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
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.
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
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, which is the other key step. We understand that the global matrix propagator 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 from as
where the matrices and are defined as Eq. (48). The symbolin Eq. (57) means the summation of the four quarters of, 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.
4.2.2. Numerical example
First consider the same model as
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.
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.