1 INTRODUCTION

When a ship or a submarine sails in a shallow sea, an elastic wave is usually generated on the seabed by the low-frequency noise of the ship. This type of elastic wave is commonly called 'ship seismic wave' and can be used to identify ships (Lu et al., 2004). Previous experimental studies investigated the transformation of underwater acoustic radiation into seismoacoustic waves (Averbakh et al., 2002; Dolgikh and Chupin, 2005). Theoretical modeling and analyses of seismic waves on the seafloor generated by lowfrequency sound waves were carried out using the normal mode model (Yan et al., 2006), fast field program (Lu et al., 2011) and parabolic equations (Zhang et al., 2010). These researchers focused on the propagating properties of seismic waves in the frequency domain by calculating transmission loss curves, frequency response curves, and so on. However, analysis in the frequency domain cannot distinguish accurately between the different components of the seismic waves. These components—P-waves, S-waves, normal modes, and interface waves—are important for our understanding of the propagating mechanism of ship seismic waves. A numerical simulation in the time domain of seismic waves generated by a low-frequency source is needed to analyze the wave components of ship seismic waves. The staggered-grid finite difference method was developed and is widely used to simulate seismic waves in the time domain (Zhang et al., 2008; Yang et al., 2012); however, little attention has been paid to seismic waves caused by low-frequency sound sources in shallow seas and in far-field cases. In this work, we first introduce an algorithm for simulating seismic waves at the seafloor in shallow water based on the staggered-grid finite difference method. Then we present numerical simulation of seismic waves caused by low-frequency sound waves from a point source in a typical shallow sea environment. The primary goal of this work is to establish a finitedifference scheme and demonstrate its effectiveness by analyzing the various types of ship seismic waves. The results of this work may serve as a guideline for other studies on similar wave phenomena.

2 STAGGERED-GRID FINITE DIFFERENCE METHOD

To model wave propagation in shallow water with a perfectly elastic, isotropic seabed, we follow the staggered-grid finite difference method proposed by Virieux (1986) and Graves (1996), which solves the following elastodynamic equations:

where v_i is the particle velocity component; τ_ij is the stress tensor component; λ and μ are Lame's constants; ρ represents the medium density; and δ_ij is the component of the Kronecker tensor; θ is the bulk strain.

In this study, we only carried out a 2-dimensional (2D) analysis which is simpler than a 3-dimensional approach and can provide the primary propagation characteristic of seismic waves in a shallow sea generated by a low-frequency sound source. 2D analysis gives a qualitative representation of the wave processes only and does not evaluate the level of acoustic and bottom seismic waves quantitatively (especially in the case of a sloping bottom). The staggered grid defines some of the velocity and stress components shifted from the locations of other components by half a grid length in space as shown in Fig. 1. By using a second-order approximation for the time derivatives and tenth-order approximation for the spatial derivatives (Dong et al., 2000; Pei and Mu, 2003), the solutions of Eqs.1 and 2 can be expressed in discrete forms by:

where U and V are the discretized particle velocity variables v_x and v_z, and P, Q, and R are the discretized stress variables τ_xx, τ_zz, and τ_xz, respectively. In these equations, the superscripts refer to the time indices and the subscripts refer to the spatial indices, for example, represents the x-component of the velocity at the point (i+1/2)Δx, j Δz at time (n+1/2)Δt. The central finite-difference operators, and represent the spatial derivatives in the x and z direction, respectively. o(△x^2N), N=5 indicates that there are 5 neighboring grids at each side of the grid to be evaluated, 2N=10 corresponds to the 10th-order approximation for the spatial derivatives, and C_Nⁿ is the coefficient of difference scheme. The superscripts of D_x⁺ and D_z^- refer to forward difference and backward difference, respectively.

In this article, the perfectly matched layer (PML) technology proposed by Collino and Sogka (2001) is applied to reduce the reflections from the artificial boundaries of the numerical model. The whole computational domain is divided into an interior region (white region) and PML regions (gray regions) as shown in Fig. 2.

For modeling a semi-infinite space, the free boundary condition is the stress-free condition at the free surface, requiring (τ_iz)_z=0=0. We satisfy the free surface boundary condition by the zero-stress formulation of Graves (1996). If we assume symmetry for the stress components about z=0 and extend the grid several nodes above z=0, we can use the boundary condition to satisfy the free-surface condition. In Fig. 3, the zero-stress free-surface boundary is shown by the heavy solid line and coincides with the location of the normal stress nodes. Choosing the z axis as the vertical direction (positive downward) and setting the free surface at z=0, the following expression must be satisfied for the 2-D model:

Upon discretizing the model, certain values of the stress components need to be specified at and above the free-surface boundary to solve the interior system given by Eqs.3–7. Setting the free-surface boundary at iz=j (Fig. 3), the values of the stress field at and above the free surface are obtained using the antisymmetric property:

At the seafloor interface, arithmetic averaging must be used for the density values, and harmonic averaging is required for the shear modulus (Graves, 1996). In Fig. 3, if we replace the zero-stress free surface with the seafloor interface, these averages can be represented by

The wave equation of the fluid layer can be derived from the elastodynamic equations (Eqs.1 and 2). We simply set ρ as the seawater density, μ=0, and substitute λ with the bulk modulus K=ρc² of sea water (Shapiro et al., 2002).

For the numerical modeling, the time function of the source is given by,

where the variable f represents the vertical force or acoustic pressure in the following numerical examples; f_c is the central frequency and t_c describes the frequency bandwidth of the wavelet. A wavelet with f_c=10 Hz and t_c=0.2 s is shown in Fig. 4.

To minimize the numerical dispersion and guarantee the accuracy of the staggered-grid finitedifference (FD) scheme, we selected a spatial discretization scheme that satisfies the following inequality (Pitarka, 1999),

where f_max is the maximum frequency of the propagating signal, v_min is the lowest velocity, and Δ is the grid spacing. A theoretical analysis of the stability condition of the proposed scheme is difficult to obtain; therefore, we adopted the stability condition of the staggered-grid scheme simplified for a constant grid spacing (Dong et al., 2000),

where Δt is the time increment, and λ, μ are the Lame constants.

3 NUMERICAL EXAMPLES

In this study we performed two types of tests. Test 1 was designed to analyze the accuracy of the FD scheme with a vertical-force source in the seabed. In the second group of tests, sound sources located in the water were used to simulate the propagation of acoustic and seismic waves in a shallow sea.

3.1 Accuracy analysis

To assess the accuracy of the proposed staggeredgrid FD scheme, we modeled the wave propagation in a two-layered laterally homogeneous media (Fig. 2) (Test 1). The model is 1 250 m wide and 1 250 m deep. The spatial and temporal steps are Δx=Δz=5 m and Δt=1 ms. There are 10 cells of PML along the bottom and vertical sides of the computational domain. The source is a vertical force located at a depth of 575 m below the sea floor at the center of the model. The source function is a cosine envelope wavelet given by Eq.15 with a central frequency of 40 Hz and a bandwidth of 80 Hz. The upper layer is sea water with a depth of 100 m, v_p=1 500 m/s, and ρ=1 000 kg/m³. The lower layer represents the seabed with v_p=2 400 m/s, v_s=1 600 m/s and ρ=1 800 kg/m³. The receiver is 300 m from the source in both the x and z directions. In this model the waveform at the receiver is equivalent to the waveform recorded in an infinite space before any reflected waves have reached the receiver. Thus, the analytic solutions in the time domain for a 2D unbound space with a vertical-force source given by Pilant (1979) can be used to test the accuracy of the FD scheme above.

where u(x, z, t) and v(x, z, t) are the horizontal and vertical displacements at location (x, z) and time t; v_p and v_s are the compressional and shear wave speeds; H represents the Heaviside function; f is the source wavelet function; r and θ are the distance and the azimuthal angle between the receiver and the source.

The numerical results based on the staggered-grid FD scheme agree well with the analytical solutions (Fig. 5), thus validating the staggered-grid FD method. The slight differences between the two solutions were caused by the numerical dispersion that occurred because the maximal frequency of the source was limited to 80 Hz. As most of the frequency of sound source in this study is very low and the wavelength is much wider than the grid size, the numerical dispersion in Test 1 will be relatively weak in other tests. Figure 6 shows a snapshot of the velocity component in the x-z plane at time 0.315 s. The source excites two direct waves that propagate through the seabed, a compressional wave P₂ and a shear wave S₂, where the subscript 2 indicates the sequence number of the layer in the half-space (1-water layer, 2-seabed layer). The horizontal velocities of the P-and S-waves in Fig. 6 are zero in the horizontal and vertical directions. This shape is coincident with Eq.19 where u(x, z, t) is antisymmetric about the horizontal and vertical directions. Similarly, the shape of the P-wave and the S-wave in Fig. 7 is symmetric about the horizontal and vertical directions according to the trigonometric functions in Eq.20.

The PML layers on the bottom and sides of the computational domain absorb the outgoing waves very well and spurious reflection is almost invisible. At t=0.315 s the upgoing direct wave P₂ had interacted with the seafloor interface and produced a downgoing reflected compressional wave P₂P₂, a reflected shear wave P₂S₂ in the seabed, and an upgoing transmitted wave P₂P₁ in the sea water. The downgoing reflected wave P₂P₁P₁ is produced from the transmitted wave P₂P₁ at the water surface. While both compressional and shear waves are generated in the sea bed layer, only compressional waves are generated in the water layer. The above propagation characteristics demonstrate that the numerical model presented here achieves good simulation of the artificial boundary, seafloor interface, and free surface.

3.2 Numerical simulation in the near-field case

Based on numerical Test 1, the forward numerical simulation of the seismic wave generated by a lowfrequency sound source in a shallow sea can be carried out if we replace the vertical force located in the seabed with a sound source located in the water and adjust the size of the computational domain accordingly (with the same grid size as in Test 1). The far-field case is defined as the situation where the source and the receiver or the interest area are separated by a distance that is more than ten times the depth of the water. In this case, most of the acoustic rays arriving at the receiver were reflected by the sea floor and sea surface many times, with the water layer acting as a waveguide. The near-field case occurs when the receiver or the interest area are separated from the source by a distance that is only a few times the depth of the water.

First, we calculated the seismic wave generated by a low-frequency sound source in a shallow sea for the near-field case (Test 2). Snapshots of the vertical component of the particle velocity at time 0.1 s and 0.125 s are shown in Figs. 8 and 9, respectively. The "+" symbol in Fig. 8 represents the sound source which is placed 20 m above the seafloor and has the same source function as in Test 1. Placing the source at a water depth of 80 m aims to simulate a submarine moving in shallow sea. At 0.1 s (Fig. 8), the source pulse labeled P₁ has interacted with the seafloor and the sea surface. The upgoing reflected wave P₁P_1u and the downgoing reflected wave P₁P_1d were produced from the lower and upper boundaries of the water layer. In Fig. 8 we also observe two wavefronts in the seabed layer, namely, the transmitted compressional wave P₁P₂ and the transmitted shear wave P₁S₂ which are produced from the source pulse P₁ when it hits the seafloor. At t=0.125 s (Fig. 9) the interface waves have propagated from the sea floor and their amplitudes decay rapidly away from the interface. The first seismic wave in the seabed is the low-amplitude compressional wave P₁P₂, followed by the transmitted shear wave P₁S₂ and the interface waves along the x direction. In the near-field case, the seismic waves in the seabed are composed mostly of the transmitted shear wave P₁S₂ and the interface waves.

3.3 Numerical simulation in the far-field case

In Test 3, we carried out forward numerical simulation of the seismic wave in the far-field case. The source pulse has a central frequency of 20 Hz and a bandwidth of 40 Hz. The computational domain is extended to 2 500 m in the horizontal direction. Figures 10 and 11 are snapshots of the horizontal and vertical velocity, respectively, in the x-z plane at 1.5 s. The source position in Test 3 is at x=0 m, z=80 m. We also added 30 cells on the left side of the sound source (to the left of x=0 m). This area is not included in the computational domain because it is used only to separate between the sound source and the left PML boundary. There are two distinct wave packets along the propagation direction in the snapshots. The wave packet at a range of about 1 800 m corresponds to the interface wave because its peak amplitude is highest along the seafloor with an exponentially decaying amplitude away from the guiding interface. The wave packet at the radial range of 2 000–2 300 m has a different pattern; the horizontal velocity reaches its maximum while the vertical velocity diminishes at the mid-depth of the water layer. Moreover, the horizontal velocity signal is in phase throughout the water column while the vertical velocity in the middle of the water column is in antiphase to the velocity in the lower and upper parts of the column. This corresponds exactly to the depth function of the first mode according to normal mode theory (Jensen et al., 1994; Li et al., 2000). Hence, the wave components of the seismic wave in the far-field case are mainly the normal mode and the interface wave caused by the waveguide effect of the sea water and the interface effect of the sea floor.

To analyze the influence of the source frequency on the seismic wave in the far field, in Test 4 we reduced the central frequency of the source pulse to 10 Hz with a bandwidth of 20 Hz. Snapshots of the horizontal and vertical velocity are shown in Figs. 12 and 13, respectively. Compared with Figs. 10 and 11, the peak amplitude of the interface wave is higher while that of the normal mode is lower. This effect is also evident in the time series of the vertical velocity in Fig. 14 (f_c=20 Hz) and Fig. 15 (f_c=10 Hz). Thus, a lowfrequency source is favorable for generating interface waves in the far field. This propagation characteristic was analyzed by Lu et al. (2013) who derived the phase velocity dispersion curves for waves in a shallow sea. As shown in Fig. 16 (Lu et al., 2013), there are usually an interface wave and several normal modes displayed in the dispersion curves when the frequency range is below 100 Hz. Each normal mode has a well-defined low-frequency cutoff while the interface wave does not. Thus, the time series of the seismic wave is dominated by the interface wave when the source frequency is less than the minimal cut-off frequency of the normal modes ( < 20 Hz). Higher-order normal modes appear when the source frequency increases and gradually passes their cut-off frequencies (>20 Hz). The normal modes gradually become the dominant carriers of the source energy, and the interface wave weakens accordingly.

4 CONCLUSION

To analyze the wave components and the propagation properties of ship seismic waves in a shallow sea, we developed an algorithm for simulating seismic waves at the seafloor based on the staggeredgrid finite difference method. Forward numerical simulations of seismic waves generated by a lowfrequency sound source in a shallow-water environment were carried out. Because the frequency of the sound-wave source in this work is very low and the bottom attenuation is relatively small compared with that of a high-frequency sound source, we ignored the dissipation losses in this simulation. Based on the results of the numerical simulations without bottom attenuation we can draw the following conclusions.

(1) The seismic waves in the near-field case are mostly composed of transmitted S-waves and interface waves while the transmitted P-waves are relatively weak near the seafloor.

(2) In the far-field case, the wave components of the seismic waves are mainly the normal modes and the interface waves.

(3) As the frequency of the source decreases, the peak amplitude of the interface wave increases while that of the normal mode decreases. Thus, a lowfrequency source is favorable for generating interface waves in the far field.

The wave components and the propagation properties of ship seismic waves in shallow water are discussed to some extent in this paper. The hydrogeology conditions in a shallow sea are complex; these include shallow water, negative sound-speed gradients, sloping sea bottom and poroelastic sediments. In this work we did not account for dissipation losses of the seismic wave propagation in the bottom medium. Wave absorption by propagation in the bottom medium is higher than in the water medium; this can affect the simulation results. Therefore, this effect should be incorporated in future studies on the propagation properties of ship seismic waves in a shallow sea.