Characteristics of Very Low Frequency Sound Propagation in Full Waveguides of Shallow Water

This work is concerned with the characteristics of very low frequency sound propagation (VLF, ≤100 Hz) in the shallow marine environment. Under these conditions, the classical hypothesis of considering the sea bottom as a fluid environment is no longer appropriate, and the sound propagation characteristics at the sea bottom should be also considered. Hence, based on the finite element method (FEM), and setting the sea bottom as an elastic medium, a proposed model which unifies the sea water and sea bottom is established, and the propagation characteristics in full waveguides of shallow water can be synchronously discussed. Using this model, the effects of the sea bottom topography and the various geoacoustic parameters on VLF sound propagation and its corresponding mechanisms are investigated through numerical examples and acoustic theory. The simulation results demonstrate the adaptability of the proposed model to complex shallow water waveguides and the accuracy of the calculated acoustic field. For the sea bottom topography, the greater the inclination angle of an up-sloping sea bottom, the stronger the leak of acoustic energy to the sea bottom, and the more rapid the attenuation of the acoustic energy in sea water. The effect of a down-sloping sea bottom on acoustic energy is the opposite. Moreover, the greater the pressure wave (P-wave) speed in the sea bottom, the more acoustic energy remains in the water rather than leaking into the bottom; the influence laws of the density and the shear wave (S-wave) speed in the sea bottom are opposite.


Introduction
Sound propagation in shallow water has always been a hotspot in the field of underwater acoustics, and is the basis for understanding, predicting, and applying various shallow sea acoustic phenomena [1]. Due to the development of submarine stealth technology and the requirements for monitoring various physical ocean phenomena, the detection methods for targets in shallow water gradually turn to very low frequency (VLF, ≤100 Hz); this is closely related to the propagation characteristics of VLF acoustic signals in shallow water. Therefore, the research on low-frequency sound propagation in shallow water has received increasing attention. The oceans surrounding China are mostly shallow and have a typical shallow water waveguide environment [2,3]. As a result, in recent years VLF sound propagation in shallow water has received more considerable attention within the acoustic modeling community in China [4]. Nevertheless, the existing research on sound propagation in shallow water uses horizontal layered waveguides for the environment models and treats the sea bottom as a range-independent fluid medium [5]. Research on modeling of VLF acoustic fields in shallow water with complex terrain changes that treats the sea bottom as an elastic medium is scarce [6]. Actually, the sea bottom structures in shallow water are composed of a sedimentary layer of a near-porous medium and a basement layer of a near-elastic medium [7]. Because of their long wavelength and strong penetrating power, the VLF underwater acoustic signals would penetrate both the sedimentary layer and the basement layer when propagating in shallow water. The sea bottom should therefore be considered as a range-dependent elastic medium when studying the characteristics of VLF sound propagation in shallow water.
In this type of sea environment model, the strong penetrability of VLF acoustic signals could cause substantial acoustic energy to leak into the sea bottom, generating seismic acoustic signals that propagate through the sea bottom or bottom surface [8][9][10]. These seismic acoustic signals generally attenuate more slowly than acoustic signals propagating through sea water, thus they can be detected and identified at longer distances [9] and have broad applications in detection for any targets in shallow water [11][12][13]. But the existing research has mostly focused on the modeling and distribution of the propagation characteristics of VLF acoustic signals in only the sea water layer. Due to the lack of research findings, there are still few kinds of sensors or equipment that have been developed for any applications with the seismic acoustic signals. Therefore, when discussing the characteristics of VLF sound propagation in shallow water, it is important to comprehensively analyze the propagation characteristics in a full waveguide model which unifies the sea water and sea bottom in shallow water.
Existing research on sound propagation is based on reasonable acoustic field calculation methods. Over the years, a number of acoustic field calculation methods, such as the normal mode method, ray method, parabolic equation method, fast field method, and their derivatives, have been established to study sound propagation in shallow water. These methods are based on limiting assumptions and approximations of the wave equation and sea environment, restricting their universality, especially for the calculation of sound propagation at VLF in range-dependent shallow water with an elastic bottom [14]. As increasing attention is being paid to VLF sound propagation in range-dependent waveguides, it is critical to establish a shallow water acoustic field calculation method with greater universality. The finite element method (FEM) accurately describes the changes in an acoustic field by dividing the environment into discrete units [15]. In the past, beyond providing reference solutions, the FEM has seldom been used for sound propagation in the sea due to the large calculations involved. However, due to advances in computer technology, it is now possible to implement the FEM for sound propagation in shallow water.
To improve the existing research on the propagation characteristics of VLF acoustic signals in shallow water, in this paper we propose a finite element calculation method for VLF sound propagation. Taking the acoustic energy flux as the research object [16], we discuss the propagation characteristics of VLF acoustic signals in range-dependent full waveguides of shallow water with an elastic bottom. This paper is organized as follows: a calculation method for VLF underwater acoustic field based on FEM is described in Section 2.1, and the formulas for the acoustic energy flux in sea water and sea bottom of shallow water are given in Section 2.2. In Section 3, the effects of the sea bottom topography and the geoacoustic parameters on VLF sound propagation and its corresponding mechanisms are elaborated with numerical examples. Finally, our conclusions are presented in Section 4.

Finite Element Form of Acoustic Field in Shallow Water
To explore the characteristics of VLF sound propagation in shallow water, a waveguide model that unifies the sea water and the sea bottom is established in three-dimensional cylindrical coordinates (r, θ, z), and the N × 2D hypothesis is adopted for this model. Under this hypothesis, the acoustic wave coupling between adjacent two-dimensional vertical planes (r, z) in the θ direction is neglected, and the original three-dimensional solution region is solved into a two-dimensional region in the (r, z) plane. Schematic diagrams of the full waveguide researched here are shown in Figure 1, with the model in three-dimensional cylindrical coordinates and the solution region in the (r, z) plane shown in Figure 1a,b, respectively.

Finite Element Form of Acoustic Field in Shallow Water
To explore the characteristics of VLF sound propagation in shallow water, a waveguide model that unifies the sea water and the sea bottom is established in three-dimensional cylindrical coordinates (r, θ, z), and the N × 2D hypothesis is adopted for this model. Under this hypothesis, the acoustic wave coupling between adjacent two-dimensional vertical planes (r, z) in the θ direction is neglected, and the original three-dimensional solution region is solved into a two-dimensional region in the (r, z) plane. Schematic diagrams of the full waveguide researched here are shown in Figure 1, with the model in three-dimensional cylindrical coordinates and the solution region in the (r, z) plane shown in Figure 1a and 1b, respectively. In the solution region, the waveguide is regarded as a range-dependent environment that presents a uniform fluid medium, Sw, above an infinite uniform elastic half space, Sb, and the sea depth H(r) depends on the horizontal range, r. The sea surface is defined by the plane z = 0, and the z-axis indicates the depth of this region. The r-axis represents the direction of horizontal propagation of the acoustic signals. ρ1 and c1 are the density and the sound speed, respectively, in sea water. ρb, cp, and cs are the density, pressure wave (Pwave) speed, and shear wave (S-wave) speed, respectively, in the sea bottom [6], and zs is the depth of the point acoustic source which is set on the symmetric axis of the cylindrical coordinate, L1. L2 is the air/fluid interaction boundary between the air and the sea water, while L3 is the fluid/elastic interaction boundary between the sea water and the sea bottom.
The FEM is based on the weak form equation [16,17], which discretizes the calculated models to finite elements and then calculates the solution of all the elements. The approximate solution of the model in each element can be obtained according to the variable value of the element node.
By combining the weighted integration of the Helmholtz equation with Gaussian theory [18], the finite element equation for Sw and Sb in Figure 1 can be expressed as In Equation (1), pi is the sound pressure at the ith node in the sea water layer, Mw, Kw, and Cw are respectively the mass matrix, the stiffness matrix, and the damping matrix in In the solution region, the waveguide is regarded as a range-dependent environment that presents a uniform fluid medium, S w , above an infinite uniform elastic half space, S b , and the sea depth H(r) depends on the horizontal range, r. The sea surface is defined by the plane z = 0, and the z-axis indicates the depth of this region. The r-axis represents the direction of horizontal propagation of the acoustic signals. ρ 1 and c 1 are the density and the sound speed, respectively, in sea water. ρ b , c p , and c s are the density, pressure wave (P-wave) speed, and shear wave (S-wave) speed, respectively, in the sea bottom [6], and z s is the depth of the point acoustic source which is set on the symmetric axis of the cylindrical coordinate, L 1 . L 2 is the air/fluid interaction boundary between the air and the sea water, while L 3 is the fluid/elastic interaction boundary between the sea water and the sea bottom.
The FEM is based on the weak form equation [16,17], which discretizes the calculated models to finite elements and then calculates the solution of all the elements. The approximate solution of the model in each element can be obtained according to the variable value of the element node.
By combining the weighted integration of the Helmholtz equation with Gaussian theory [18], the finite element equation for S w and S b in Figure 1 can be expressed as In Equation (1), p i is the sound pressure at the ith node in the sea water layer, M w , K w , and C w are respectively the mass matrix, the stiffness matrix, and the damping matrix in the fluid medium, {F wi } is the acoustic excitation at this node. In Equation (2), s i is the displacement at the ith node in the sea bottom, M b , K b , and C b are respectively the mass matrix, the stiffness matrix, and the damping matrix in elastic medium without boundary constraints, and {F bi } is the excitation load at this node.
When calculating the air/fluid interaction boundary, the L 2 is set as a Dirichlet boundary, as defined in Equation (3).
L 3 is set as a fluid/elastic coupling boundary, which should satisfy coupling conditions such as continuous normal displacement, continuous normal stress, and zero tangential stress. The definitions are as follows: where s n r and s n z (n = w, b) represent the displacement along the r-axis and z-axis in medium S w and medium S b . λ n and µ n (n = w, b) respectively represent the Lame constant in medium S w and medium S b . To simulate sound propagation in an infinite sea environment, a perfectly matched layer (PML) is used for the boundary treatment, as shown in Figure 1 [16].
Based on the FEM, by coupling the Helmholtz equation, fluid-elastic coupling boundary, and PML boundary, the acoustic pressure field p in fluid medium S w and the displacement field s in elastic half space S b can be calculated.

Acoustic Energy Flux in Full Waveguide of Shallow Water
In previous studies on the characteristics of acoustic intensity in a shallow water waveguide, the study object was the pressure field p [3,5]. However, p only reflects the propagation characteristics of the P-wave in a fluid medium. Because the sea bottom is considered as an elastic medium in this paper, it is impossible to correctly describe the propagation characteristics of acoustic intensity in a full waveguide of shallow water with p. Unlike p, the acoustic energy flux I exists in both the fluid medium and the elastic medium, making it more suitable for revealing the propagation characteristics of acoustic intensity in the full waveguide. Therefore, in this paper, the propagation characteristics of VLF acoustic intensity in a shallow water waveguide is discussed in terms of the acoustic energy flux.
In the time domain, the acoustic energy flux I is the time average of the instantaneous acoustic intensity I(t). For an isotropic elastic medium, its definition is [19,20]: where T and v are respectively the stress matrix and particle velocity field in the medium under three-dimensional cylindrical coordinates; T and v can be obtained through the corresponding relationship between p, s [21]. As the acoustic field is independent to the azimuthal angle θ in the cylindrical coordinate, Equation (7) can be simplified as: In the time domain, the acoustic energy flux I can be obtained by computing the time average of Equation (8). According to Parseval's theorem, the average acoustic energy flux I(ω) in the frequency domain can be expressed as: where the asterisk (*) denotes the complex conjugate.
As the elements in the vector I(ω) are complex, the real part Re[I(ω)] represents the acoustic energy that can propagate over long distances, while the imaginary part Im[I(ω)] represents the acoustic energy that does not propagate. When sound waves propagate through an elastic medium, the amplitude and direction of the complex vector I(ω) are defined as follows: Since µ = 0, T rz = T zr = 0 and T rr = T zz = λ∇ 2 ϕ p = p in the fluid, the equation for the acoustic energy flux in the fluid layer can be written as: Using Equations (9) and (12), the acoustic energy flux I(ω) in each layer of the full waveguide in shallow water can be calculated. In the following sections, we consider I(ω) as the research object for investigation of the propagation characteristics of VLF acoustic signals in shallow water.

Simulation and Discussion
In this section, using numerical examples, we first verify the accuracy of the FEM simulation results for VLF sound propagation in shallow water, and then elaborate the effects of the sea bottom topography and the geoacoustic parameters on the characteristics of VLF sound propagation and its corresponding mechanisms.

Comparison of Simulation Results
We verify the FEM simulation of the VLF underwater acoustic field for three types of sea bottom topography, horizontal [20], wedge-shaped uphill [14], and wedge-shaped downhill [22], by comparison with existing research results. The two-dimensional solution region and environmental parameters for the three types of shallow water environments are shown in cylindrical coordinates in Figure 2.
ASA model [23,24], in which the sea water depth is decreased linearly from 200 m at 0 km to 0 m at 4 km, so the uphill slope angle is α1 ≈ 2.86°, the density and P-wave speed in this model are ρb = 1500 kg/m 3 , cp = 1700 m/s (attenuation for P-wave is αp = 0.5 dB·λ −1 ), the depth of the 100 Hz point acoustic source is 100 m. For the downhill sea bottom model in Figure  2c, the sea water depth increases linearly from 200 m at 0 km to 400 m at 4 km at a slope angle of α2 ≈ 2.86°; the density, P-wave speed, and acoustic source are set the same as in Figure 2b. For the three simulation conditions, the sound pressure Transmission Loss (TL) curves calculated by FEM and other sound field simulation methods are shown for comparison in Figure 3. The receiving position is set at a depth of z = 30 m for all simulations. The sound pressure TL calculation formula is given by Equation (13)  In each environment illustrated in Figure 2, the sound speed and density of the sea water are respectively set to c 1 = 1500 m/s and ρ 1 = 1000 kg/m 3 . Figure 2a depicts a horizontally layered model of a shallow water waveguide that treats the sea bottom as an elastic medium, in which the density, P-wave speed, and S-wave speed are respectively set to ρ b = 1500 kg/m 3 , c p = 2000 m/s (attenuation for P-wave is α p = 0.1 dB·λ −1 ), and c s = 1000 m/s (attenuation for S-wave is α s = 0.1 dB·λ −1 ), the depth of the 100 Hz point acoustic source is 20 m. In Figure 2b, the sea bottom for the shallow water is wedge-shaped uphill, as the ASA model [23,24], in which the sea water depth is decreased linearly from 200 m at 0 km to 0 m at 4 km, so the uphill slope angle is α 1 ≈ 2.86 • , the density and P-wave speed in this model are ρ b = 1500 kg/m 3 , c p = 1700 m/s (attenuation for P-wave is α p = 0.5 dB·λ −1 ), the depth of the 100 Hz point acoustic source is 100 m. For the downhill sea bottom model in Figure 2c, the sea water depth increases linearly from 200 m at 0 km to 400 m at 4 km at a slope angle of α 2 ≈ 2.86 • ; the density, P-wave speed, and acoustic source are set the same as in Figure 2b.
For the three simulation conditions, the sound pressure Transmission Loss (TL) curves calculated by FEM and other sound field simulation methods are shown for comparison in Figure 3. The receiving position is set at a depth of z = 30 m for all simulations. The sound pressure TL calculation formula is given by Equation (13):  In Figure 3, the dotted lines are calculated results by FEM, and the solid lines are the results by three comparison acoustic field simulation methods. Based on their applicability, the fast field method (Scotter program) [17,25,26] and the parabolic equation method (RAM program) [26,27] are used to examine the error of FEM for the models in Figure 2ac, respectively. The comparison results in Figure 3 show that the TL curves that were simulated by FEM are highly consistent with the results obtained by the existing simulation programs. This comparison verifies the accuracy of FEM for the simulation of underwater acoustic fields in shallow water. Considering the FEM is more suitable for the acoustic field simulation in complex ocean environments and can be used to solve the acoustic field in full waveguides of shallow water, the FEM is used for simulation of the characteristics of VLF sound propagation.

The Effects of Sea Bottom Topography on the Characteristics of VLF Sound Propagation
In this section, the effects of sea bottom topography on VLF sound propagation in full waveguides of shallow water are discussed. The two most common types of sea bottom topography, wedge-shaped uphill and wedge-shaped downhill, are discussed in Sections 3.2.1 and 3.2.2, respectively. In Figure 3, the dotted lines are calculated results by FEM, and the solid lines are the results by three comparison acoustic field simulation methods. Based on their applicability, the fast field method (Scotter program) [17,25,26] and the parabolic equation method (RAM program) [26,27] are used to examine the error of FEM for the models in Figure 2a-c, respectively. The comparison results in Figure 3 show that the TL curves that were simulated by FEM are highly consistent with the results obtained by the existing simulation programs. This comparison verifies the accuracy of FEM for the simulation of underwater acoustic fields in shallow water. Considering the FEM is more suitable for the acoustic field simulation in complex ocean environments and can be used to solve the acoustic field in full waveguides of shallow water, the FEM is used for simulation of the characteristics of VLF sound propagation.

The Effects of Sea Bottom Topography on the Characteristics of VLF Sound Propagation
In this section, the effects of sea bottom topography on VLF sound propagation in full waveguides of shallow water are discussed. The two most common types of sea bottom topography, wedge-shaped uphill and wedge-shaped downhill, are discussed in Sections 3. acoustic energy flux in the (r, z) region are compared for a VLF acoustic source at 100 Hz. The sound intensity level (SIL) for acoustic energy flux is defined by Equation (14): For the simulation, aside from the sea bottom topography, the environmental parameters have the same values as in Figure 2a. ∆H 1 represents the vertical distance between the sea surface and the sea bottom at r = 4 km. The amplitude and direction of the acoustic energy flux at each grid point are represented by the length and direction of the arrows. Figure 5 depicts the distribution of the 100 Hz acoustic energy flux in a full waveguide of shallow water for a wedge-shaped uphill sea bottom. Aside from ∆H 1 , the environmental parameters are the same as in Figure 2a. Figure 5a-c respectively depict the simulation results at ∆H 1 = 80, 50, and 0 m with corresponding uphill slope angles of α 1 = 0.29 • , 0.72 • , and 1.43 • . Figure 5d shows a comparison of the propagation characteristics of the 100 Hz acoustic energy flux for four sea bottom simulation conditions: ∆H 1 = 100, 80, 50, and 0 m. The receiving depth in the sea water is z = 50 m. Figure 6 depicts the distribution of the ray tracing when acoustic signals propagate in a shallow water environment with wedgeshaped uphill sea bottom, and the uphill slope angle of the wedge-shaped bottom is α 1 .
According to normal mode method, under the simulation conditions in Figure 2a, nine orders of normal modes would be excited in the shallow water waveguide at 100 Hz, and the grazing angles θ of these normal modes on the horizontal sea bottom interface are in the range [4 • , 40 • ]. The higher the normal mode order is, the larger the grazing angle [17]. By comparing the simulation results in Figures 4 and 5, it can be seen that under the influence of the uphill slope angle of the sea bottom α 1 , after its nth reflection on the sea bottom interface, the grazing angle of each excited normal mode is increased by (2n − 1)α 1 degrees, as shown in Figure 6. For VLF sound propagation in shallow water with a wedge-shaped uphill sea bottom, the propagation characteristics of the low-order normal modes, which correspond to small grazing angles, are continuously coupled to the high-order normal modes, which correspond to large grazing angles, significantly strengthening the leakage effect of the acoustic energy to the sea bottom and relieving the fluctuation. The larger the slope angle α 1 is, the greater the number of low-order normal modes that are coupled to higher orders, the more the acoustic energy leaks into the sea bottom, the faster the attenuation of the acoustic energy in the sea water, the fewer normal mode orders, and the simpler the normal mode interference in far-field regions. For these reasons, under the simulation conditions in Figure 6, as ∆H 1 changes from 100 to 0 m, the 100 Hz acoustic energy flux reveals a loss of more than 20 dB.
For the simulation, aside from the sea bottom topography, the environmental parameters have the same values as in Figure 2a. △H1 represents the vertical distance between the sea surface and the sea bottom at r = 4 km. The amplitude and direction of the acoustic energy flux at each grid point are represented by the length and direction of the arrows. Figure 5 depicts the distribution of the 100 Hz acoustic energy flux in a full waveguide of shallow water for a wedge-shaped uphill sea bottom. Aside from △H1, the environmental parameters are the same as in Figure 2a. Figure 5a-c respectively depict the simulation results at △H1 = 80, 50, and 0 m with corresponding uphill slope angles of α1 = 0.29°, 0.72°, and 1.43°. Figure 5d shows a comparison of the propagation characteristics of the 100 Hz acoustic energy flux for four sea bottom simulation conditions: △H1 = 100, 80, 50, and 0 m. The receiving depth in the sea water is z = 50 m. Figure 6 depicts the distribution of the ray tracing when acoustic signals propagate in a shallow water environment with wedge-shaped uphill sea bottom, and the uphill slope angle of the wedge-shaped bottom is α1.   According to normal mode method, under the simulation conditions in Figure 2a, nine orders of normal modes would be excited in the shallow water waveguide at 100 Hz,  According to normal mode method, under the simulation conditions in Figure 2a, nine orders of normal modes would be excited in the shallow water waveguide at 100 Hz, Figure 6. Distribution of the ray tracing when acoustic signals propagate in a shallow water environment with wedge-shaped uphill sea bottom, and the uphill slope angle of the wedge-shaped bottom is α 1 .

VLF Sound Propagation in Shallow Water with Wedge-shaped Downhill Sea Bottom
Like Figure 5, Figure 7 illustrates the distribution of the 100 Hz acoustic energy flux in a full waveguide of shallow water with a wedge-shaped downhill sea bottom. Aside from ∆H 1 , the environmental parameters are the same as in Figure 2a. larger the slope angle α1 is, the greater the number of low-order normal modes that are coupled to higher orders, the more the acoustic energy leaks into the sea bottom, the faster the attenuation of the acoustic energy in the sea water, the fewer normal mode orders, and the simpler the normal mode interference in far-field regions. For these reasons, under the simulation conditions in Figure 6, as △H1 changes from 100 to 0 m, the 100 Hz acoustic energy flux reveals a loss of more than 20 dB.

VLF Sound Propagation in Shallow Water with Wedge-shaped Downhill Sea Bottom
Like Figure 5, Figure 7 illustrates the distribution of the 100 Hz acoustic energy flux in a full waveguide of shallow water with a wedge-shaped downhill sea bottom. Aside from △H1, the environmental parameters are the same as in Figure 2a. By comparing the simulation results in Figures 4 and 7, it can be seen that under the influence of the downhill slope angle of the sea bottom α2, after its nth reflection on the sea bottom interface, the grazing angle of each excited normal mode is decreased by (2n-1)α2 degrees, as shown in Figure 8. This is the opposite influence of that observed for the wedge-shaped uphill sea bottom. For VLF sound propagation in shallow water with a wedge-shaped downhill sea bottom, the propagation characteristics of the high-order normal modes, which correspond to large grazing angles, are continuously coupled to the low-order normal modes, which correspond to small grazing angles, significantly weakening the leakage effect of the acoustic energy to the sea bottom. The larger the grazing angle α2 is, the more high-order normal modes are coupled to the low orders, and the less acoustic energy leaks into the sea bottom. The wedge-shaped downhill sea bottom would By comparing the simulation results in Figures 4 and 7, it can be seen that under the influence of the downhill slope angle of the sea bottom α 2 , after its nth reflection on the sea bottom interface, the grazing angle of each excited normal mode is decreased by (2n-1)α 2 degrees, as shown in Figure 8. This is the opposite influence of that observed for the wedge-shaped uphill sea bottom. For VLF sound propagation in shallow water with a wedge-shaped downhill sea bottom, the propagation characteristics of the high-order normal modes, which correspond to large grazing angles, are continuously coupled to the low-order normal modes, which correspond to small grazing angles, significantly weakening the leakage effect of the acoustic energy to the sea bottom. The larger the grazing angle α 2 is, the more high-order normal modes are coupled to the low orders, and the less acoustic energy leaks into the sea bottom. The wedge-shaped downhill sea bottom would have an enhancement effect on the propagation of VLF signals in shallow water.
200 m, with a receiving depth of z = 50 m.
By comparing the simulation results in Figures 4 and 7, it can be seen that under the influence of the downhill slope angle of the sea bottom α2, after its nth reflection on the sea bottom interface, the grazing angle of each excited normal mode is decreased by (2n-1)α2 degrees, as shown in Figure 8. This is the opposite influence of that observed for the wedge-shaped uphill sea bottom. For VLF sound propagation in shallow water with a wedge-shaped downhill sea bottom, the propagation characteristics of the high-order normal modes, which correspond to large grazing angles, are continuously coupled to the low-order normal modes, which correspond to small grazing angles, significantly weakening the leakage effect of the acoustic energy to the sea bottom. The larger the grazing angle α2 is, the more high-order normal modes are coupled to the low orders, and the less acoustic energy leaks into the sea bottom. The wedge-shaped downhill sea bottom would have an enhancement effect on the propagation of VLF signals in shallow water.

The Effects of Geoacoustic Parameters on the Characteristics of VLF Sound Propagation
In this section, we discuss the effects of three geoacoustic parameters, the density ρ b of the sea bottom, P-wave speed c p , and S-wave speed c s in the sea bottom, on VLF sound propagation in full shallow water waveguides. The effect laws of the above three parameters are respectively discussed in Sections 3.3.1-3.3.  Figure 9a-c are ρ b = 1.2 × 10 3 , 1.8 × 10 3 , and ρ b = 2.0 × 10 3 kg/m 3 , respectively. The remaining parameters of the shallow water environment are the same as in Figure 2a.
By comparison ofFigures 2a and 9, it can be seen that, as the density of the sea bottom increases, the acoustic energy flux in the sea water decreases and more acoustic energy is transmitted into the sea bottom. These phenomena become more pronounced as the propagation range of the simulation increases, and based on the reflection rules for fluid/elastic interfaces [28], they can be further analyzed as follows.
The reflection coefficient R of the fluid/solid interface is given by Equation (15): where ρ 1 c 1 is the acoustic wave impedance in sea water, ρ b c p and ρ b c s are respectively the P-wave impedance and S-wave impedance at the sea bottom, θ is the grazing angle when an acoustic wave is incident on the fluid/elastic interface, and θ pt and θ st are respectively the refraction angles of the P-wave and S-wave in the elastic medium below the interface. The three sets of parameters (c 1 , θ), (c p , θ pt ), and (c s , θ st ) satisfy Snell's law [28] By comparison of Figures 2a and 9, it can be seen that, as the density of the sea bottom increases, the acoustic energy flux in the sea water decreases and more acoustic energy is Under the simulation conditions shown in Figures 2a and 9, the sound speed on sides of the water/bottom interface satisfies the relationship c p > c 1 > c s . Under this premise, there is a critical grazing angle θ that satisfies the relationship sinθ pt = c p /c 1 ·cosθ = 1. When an acoustic wave is incident on the interface at a grazing angle larger than θ , at the sea bottom the P-wave, which is formed by the acoustic wave refraction, propagates as a nonuniform wave, and it cannot carry acoustic energy into the sea bottom. However, at the sea bottom the S-wave, which is formed by the acoustic wave, always propagates normally, and it still carries acoustic energy into the sea bottom. As the S-wave impedance increases, more acoustic energy is carried from the sea water to the sea bottom, and the modulus of the reflection coefficient R decreases. Therefore, when the sea bottom density increases, acoustic energy is more easily transmitted into the sea bottom via propagation, leading to the attenuation of the SIL in propagation. A comparison of Figures 2a and 9 reveals that they are consistent with this analysis. Figure 10a-c demonstrates the propagation characteristics of acoustic energy flux for P-wave speeds of c p = 2500, 2800, and 3000 m/s, respectively, under the simulation conditions in Figure 2a. The remaining parameters of the shallow water environment are the same as in Figure 2a. By comparing Figures 2a and 10, it can be seen that, as the P-wave speed in the sea bottom increases, the acoustic energy flux in the sea water increases and the interference structure of the SIL becomes more complex. Using Equation (15) and Snell's Law, these simulation results can be further analyzed as follows.

Effect of the P-Wave Speed in the Sea Bottom on the Characteristics of VLF Sound Propagation
As the P-wave speed c p in the sea bottom increases, the modulus of R increases and the critical grazing angle θ increases, which shortens the horizontal span of the acoustic energy in the sea water. Therefore, more acoustic energy is reflected into the sea water rather than transmitted to the sea bottom as c p increases, leading to an increase of the SIL in propagation. The comparison between Figures 2a and 10 is consistent with this analysis.
an acoustic wave is incident on the fluid/elastic interface, and θpt and θst are respectively the refraction angles of the P-wave and S-wave in the elastic medium below the interface. The three sets of parameters (c1, θ), (cp, θpt), and (cs, θst) satisfy Snell's law [28].
Under the simulation conditions shown in Figures 2a and 9, the sound speed on sides of the water/bottom interface satisfies the relationship cp > c1 > cs. Under this premise, there is a critical grazing angle θ′ that satisfies the relationship sinθptʹ = cp/c1·cosθ′ = 1. When an acoustic wave is incident on the interface at a grazing angle larger than θ′, at the sea bottom the P-wave, which is formed by the acoustic wave refraction, propagates as a nonuniform wave, and it cannot carry acoustic energy into the sea bottom. However, at the sea bottom the S-wave, which is formed by the acoustic wave, always propagates normally, and it still carries acoustic energy into the sea bottom. As the S-wave impedance increases, more acoustic energy is carried from the sea water to the sea bottom, and the modulus of the reflection coefficient R decreases. Therefore, when the sea bottom density increases, acoustic energy is more easily transmitted into the sea bottom via propagation, leading to the attenuation of the SIL in propagation. A comparison of Figures 2a and 9 reveals that they are consistent with this analysis.
3.3.2. Effect of the P-wave Speed in the Sea Bottom on the Characteristics of VLF Sound Propagation Figure 10a-c demonstrates the propagation characteristics of acoustic energy flux for P-wave speeds of cp = 2500, 2800, and 3000 m/s, respectively, under the simulation conditions in Figure 2a. The remaining parameters of the shallow water environment are the same as in Figure 2a. By comparing Figures 2a and 10, it can be seen that, as the P-wave speed in the sea bottom increases, the acoustic energy flux in the sea water increases and the interference structure of the SIL becomes more complex. Using Equation (15) and Snell's Law, these simulation results can be further analyzed as follows. As the P-wave speed cp in the sea bottom increases, the modulus of R increases and the critical grazing angle θʹ increases, which shortens the horizontal span of the acoustic energy in the sea water. Therefore, more acoustic energy is reflected into the sea water rather than transmitted to the sea bottom as cp increases, leading to an increase of the SIL in propagation. The comparison between Figures 2a and 10 is consistent with this analysis.   As the P-wave speed cp in the sea bottom increases, the modulus of R increases and the critical grazing angle θʹ increases, which shortens the horizontal span of the acoustic energy in the sea water. Therefore, more acoustic energy is reflected into the sea water rather than transmitted to the sea bottom as cp increases, leading to an increase of the SIL in propagation. The comparison between Figures 2a and 10 is consistent with this analysis. Figure 11a-c demonstrates the propagation characteristics of the acoustic energy flux for S-wave speeds of cs = 400, 600, and 800 m/s, respectively, under the simulation condition in Figure 2a. The remaining parameters of the shallow water environment are the same as in Figure 2a. Under the premise that c1 > cs, the S-wave formed by the acoustic wave refraction can propagate to the sea bottom and can carry some acoustic energy to the bottom. The greater the speed of the S-wave is, the more the acoustic energy in the sea water is able to penetrate into the sea bottom and the greater the transmission loss of the acoustic energy in the seawater.

Conclusions
This paper demonstrates the power of FEM for handling VLF sound propagation in a full waveguide of shallow water. By taking the acoustic energy flux as the research object and implementing FEM, the effects of the sea bottom topography and the geoacoustic parameters on the characteristics of VLF sound propagation and its corresponding mechanisms are investigated through numerical simulation and acoustic theory. The research results provide a theoretical reference for the research, development, and application of underwater acoustic engineering equipment for VLF sound in shallow water, and also can be used to develop new sensors and equipment which unify underwater acoustic signals and the seismic acoustic signals in underwater detection.
From this study, we draw the following specific conclusions: (1) According to the simulation results, FEM is highly applicable for the calculation of sound propagation in complex sea environments, especially for sound propagation in range-dependent full shallow water waveguides. Our results indicate that the underwater acoustic field prediction simulated by FEM is consistent with the predictions obtained by other shallow water simulation methods. (2) Both wedge-shaped uphill and downhill sea bottoms can significantly influence VLF sound propagation in shallow water. Compared with a horizontal sea bottom, under the influence of the uphill slope angle α1, the low-order normal modes, which correspond to small grazing angles, are continuously coupled to the high-order normal modes, which correspond to large grazing angles, significantly strengthening the leakage effect of the acoustic energy to the sea bottom and relieving the fluctuation. The larger the slope angle α1 is, the more the low-order normal modes are coupled to the higher order ones, the more the acoustic energy leaks into the sea bottom, the Under the premise that c 1 > c s , the S-wave formed by the acoustic wave refraction can propagate to the sea bottom and can carry some acoustic energy to the bottom. The greater the speed of the S-wave is, the more the acoustic energy in the sea water is able to penetrate into the sea bottom and the greater the transmission loss of the acoustic energy in the seawater.

Conclusions
This paper demonstrates the power of FEM for handling VLF sound propagation in a full waveguide of shallow water. By taking the acoustic energy flux as the research object and implementing FEM, the effects of the sea bottom topography and the geoacoustic parameters on the characteristics of VLF sound propagation and its corresponding mechanisms are investigated through numerical simulation and acoustic theory. The research results provide a theoretical reference for the research, development, and application of underwater acoustic engineering equipment for VLF sound in shallow water, and also can be used to develop new sensors and equipment which unify underwater acoustic signals and the seismic acoustic signals in underwater detection.
From this study, we draw the following specific conclusions: (1) According to the simulation results, FEM is highly applicable for the calculation of sound propagation in complex sea environments, especially for sound propagation in