**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*

148 Wave Processes in Classical and New Solids

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 *R*2 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

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.

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

the validity of the present method and the correctness of the fluid-solid formulation.

**4. Water effects and water reverberation modeling** 

show the calculated results of the two cases one by one in the following.

=6000 m/s

=3460 m/s =2500 kg/m3

vertically incident plane *P* wave: (a) with a fluid layer; (b) without a fluid layer.

Plane P-wave

Solid Vp

Vs

0 32 64

Distance (km)

water surface 

Fluid Vf

<sup>f</sup>

=1500 m/s

=1000 kg/m3

Depth (km)

(a)

**Figure 16.** The model used to show the effect of a fluid layer on the synthetic seismograms due to a

(a) (b)

(b)

Depth (km)

0 32 64

Plane P-wave

=3460 m/s =2500 kg/m3

Solid Vp=6000 m/s Vs

Air

Distance (km)

**4.1. Water effects on seismogram synthesis** 

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.

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

U x

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

**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.

0

(a)

2 5

3 2 D i s t a n c e ( k m )

6 4

U z

(a) (b)

0

5

) s( m e i T

1 0

1 5

2 0

2 5

0

5

1 0

T i m e (s)

1 5

2 0

0

(b)

2 5

3 2 D i s t a n c e ( k m )

6 4

U z

So far, we have given the formulation for seismic excitation in a multi-layered solid halfspace 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

ሺଵሻܙሺଵሻ ۵ሺଶሻ

where **p** and **q** with subscripts (2) are connected by Eq. (43). In order to solve Eq. (55), we

ሺଵሻܙሺଶሻ ۴, (55)

ሺଵሻ defined in Eq. (43) to change Eq. (55) into

**4.2. Water reverberation modeling** 

need to use the global matrix propagator ۲

*4.2.1. Formulation when source is located in fluid layer* 

۶ሺଶሻ

ሺଵሻܘሺଶሻ ൌ ۵ሺଵሻ

reverberation.

U x

0

5

1 0

T i m e (s)

1 5

2 0

2 5

0

5

1 0

T i m e (s)

1 5

2 0

#### *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**

150 Wave Processes in Classical and New Solids

*4.1.2. Explosive source* 

wave incidence.

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.

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

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

#### *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 halfspace 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}\_{\text{(2)}}^{\text{(1)}}\mathbf{p}\_{\text{(2)}} = \mathbf{G}\_{\text{(1)}}^{\text{(1)}}\mathbf{q}\_{\text{(1)}} + \mathbf{G}\_{\text{(2)}}^{\text{(1)}}\mathbf{q}\_{\text{(2)}} + \mathbf{F}\_{\text{'}} \tag{55}$$

where **p** and **q** with subscripts (2) are connected by Eq. (43). In order to solve Eq. (55), we need to use the global matrix propagator ۲ ሺଵሻ defined in Eq. (43) to change Eq. (55) into

$$
\begin{pmatrix}
\end{pmatrix}
\begin{pmatrix}
\mathbf{q}\_{\{1\}} \\
\mathbf{q}\_{\{2\}}
\end{pmatrix} = \mathbf{F}\_{\prime} \tag{56}
$$

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

U x

**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.

0

(a)

2 5

3 2 D i s t a n c e ( km )

6 4

U z

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

(a) (b)

0

5

1 0

) s ( T i m e

1 5

2 0

2 5

0

5

1 0

) s ( T i m e

1 5

2 0

0

(b)

2 5

3 2 D i s t a n c e ( km )

6 4

U z

**5. Summary** 

0

5

1 0

) s ( m e i T

1 5

U x

2 0

2 5

0

5

1 0

T i m e ( s )

1 5

2 0

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

$$\mathbf{D}\_{pq}^{(1)} = \frac{1}{\rho\_l \omega^2} \Sigma \mathbf{D}\_{tu}^{(1)} / [\mathbf{N}\_1^\mathbf{e} \cdot (\mathbf{N}\_1^\mathbf{e})^\mathbf{T} + \mathbf{N}\_2^\mathbf{e} \cdot (\mathbf{N}\_1^\mathbf{e})^\mathbf{T} + \mathbf{N}\_1^\mathbf{e} \cdot (\mathbf{N}\_2^\mathbf{e})^\mathbf{T} + \mathbf{N}\_2^\mathbf{e} \cdot (\mathbf{N}\_2^\mathbf{e})^\mathbf{T}],\tag{57}$$

where the matrices N� � and N� � are defined as Eq. (48). The symbol Σ in 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 *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**

152 Wave Processes in Classical and New Solids

��� (�) <sup>=</sup> �

where the matrices N�

���(�)

problem. What is left is to get the global matrix propagator ���

We understand that the global matrix propagator ���

� and N�

**Step 1.** Calculate the global matrix propagators ���

the fluid layer into the one in the fluid layer by Eq. (57). **Step 4.** Substitute the reduced global matrix propagator ���

**Step 6** to get the ground motion in time domain finally.

(�)/�**N**�

�f�� ���

summation of the four quarters of ���

summarized into the following.

layer by Eqs. (38-39).

*4.2.2. Numerical example* 

(�) �(�)

continuity conditions at the fluid-solid interface Eq. (43) are used to obtain ���

�)� + **N**�

� ∙ (**N**�

(�)���

(�) � �(�)

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

fluid layer can be obtained recursively from the lowermost layer by using Eqs. (38-39). Then the

� ∙ (**N**�

are operated in the sense of 'element to element'. The general solution steps can be

**Step 2.** Calculate the global matrix propagators in all the solid layers above the bottom

**Step 3.** Reduce the global matrix propagators obtained in the solid layer right neighboring

**Step 7.** Execute the Inverse Fourier Transform on the displacement solutions obtained in

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.

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 ���

(�)� ��(�) �(�)

�)� + **N**�

� ∙ (**N**�

� are defined as Eq. (48). The symbol Σ in Eq. (57) means the

(�)

(�), and the multiplication and division of the matrices

�)� + **N**�

(�) in the lowermost layer by Eq. (39).

�=�, (56)

(�) , which is the other key step.

(�) from ���

obtained in **Step 3** into the

(�).

(�) as

�)��, (57)

(�) in the solid layer right neighboring the

� ∙ (**N**�

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.

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

[2] Miklowitz J, Achenbach JD (1978) Modern Problems in Elastic Wave Propagation. New

[3] Giovine P, Oliveri F (1995) Dynamics and Wave Propagation in Dilatant Granular

[4] Godano C, Oliveri F (1999) Nonlinear seismic waves: a model for site effects. Int. J. Non-

[5] Zhao CG, Dong J, Gao FP, Jeng DS (2006) Seismic responses of a hemispherical alluvial valley to SV waves: a three-dimensional analytical approximation. Acta Mechanica

[6] Cerveny V (2001) Seismic Ray Theory. Cambridge: Cambridge University Press. 722p. [7] Chen XF (1990) Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method-I. Theory of two-

[8] Chen XF (1995) Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method–II. Applications for 2D

[9] Chen, XF (1996) Seismogram synthesis for multi-layered media with irregular interfaces by global generalized reflection/transmission matrices method–III. Theory of 2D P-SV

[10] Fu LY, Bouchon M (2004) Discrete wavenumber solutions to numerical wave propagation in piecewise heterogeneous media–I. Theory of two-dimensional SH case.

[11] Zhang W, Chen XF (2006) Traction image method for irregular free surface boundaries

[12] Koketsu K, Fujiwara H, Ikegami Y (2004) Finite-element simulation of seismic ground

[13] Wang YB, Takenaka H, Furumura T (2001) Modelling seismic wave propagation in a two-dimensional cylindrical whole-earth model using the pseudospectral method.

[14] Komititsch D, Tromp J (2002) Spectral-element simulations of global seismic wave

[15] Fu LY, Mu YG (1994) Boundary element method for elastic wave forward modeling.

[16] Fu LY (2002) Seismogram synthesis for piecewise heterogeneous media. Geophys. J. Int.

[17] Bouchon M, Schultz CA, Toksoz MN (1995) A fast implementation of boundary integral equation methods to calculate the propagation of seismic waves in laterally varying

[18] Ge ZX, Chen XF (2007) Wave propagation in irregularly layered elastic models: a boundary element approach with a global reflection/transmission matrix propagators.

[19] Ge ZX, Chen XF (2008) An efficient approach for simulating wave propagation with the boundary element method in multilayered media with irregular interfaces. Bull. Seism.

[20] Okamoto T, Takenaka H (1999) A reflection/transmission matrix formulation for seismoacoustic scattering by an irregular fluid-solid interface. Geophys. J. Int. 139: 531–546.

in finite difference seismic wave simulation. Geophys. J. Int. 167: 337–353.

motion with a voxel mesh. Pure Appl. Geophys. 161: 2183–2198.

propagation-I. Validation. Geophys. J. Int. 149: 390–412.

layered media. Bull. Seism. Soc. Am. 85: 1679–1687.

dimensional SH case. Bull. Seism. Soc. Am. 80: 1696–1724.

SH case. Bull. Seism. Soc. Am. 85: 1094–1106.

case. Bull. Seism. Soc. Am. 86: 389–405.

Geophys. J. Int. 157: 481–498.

Geophys. J. Int. 145: 689–708.

150: 800–808.

Acta Geophys. Sinica 37: 521–529.

Bull. Seism. Soc. Am. 97: 1025–1031.

Soc. Am. 98: 3007–3016.

York: John Wiley & Sons. 561p.

Materials. Meccanica 30: 341–357.

linear Mech. 34:457–468.

Sinica 22: 547–557.

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.
