Cite this Article
Tang Jun, Piao Sheng-Chun, Zhang Hai-Gang. Three-dimensional parabolic equation model for seismo-acoustic propagation: Theoretical development and preliminary numerical implementation . Chinese Physics B, 2017, 26(11): 114301  
Permissions
Three-dimensional parabolic equation model for seismo-acoustic propagation: Theoretical development and preliminary numerical implementation
Tang Jun1, 2, Piao Sheng-Chun1, 2, Zhang Hai-Gang1, 2, †
Acoustic Science and Technology Laboratory, Harbin Engineering University, Harbin 150001, China
College of Underwater Acoustic Engineering, Harbin Engineering University, Harbin 150001, China

 

† Corresponding author. E-mail: zhanghaigang@hrbeu.edu.cn

Project supported by the National Nature Science Foundation of China (Grant Nos. 11234002 and 11704337) and the National Key esearch Program of China (Grant No. 2016YFC1400100).

Abstract

A three-dimensional (3D) parabolic equation (PE) model for sound propagation in a seismo-acoustic waveguide is developed in Cartesian coordinates, with x, y, and z representing the marching direction, the longitudinal direction, and the depth direction, respectively. Two sets of 3D PEs for horizontally homogenous media are derived by rewriting the 3D elastic motion equations and simultaneously choosing proper dependent variables. The numerical scheme is for now restricted to the y-independent bathymetry. Accuracy of the numerical scheme is validated, and its azimuthal limitation is analyzed. In addition, effects of horizontal refraction in a wedge-shaped waveguide and another waveguide with a polyline bottom are illustrated. Great efforts should be made in future to provide this model with the ability to handle arbitrarily irregular fluid–elastic interfaces.

PACS: 43.20.+g;43.30.+m
Keyword:three-dimensional parabolic equation;sound propagation;seismo-acoustic waveguides
1. Introduction

The modelling of sound propagation is of great importance in underwater acoustics. Elasticity in the ocean bottom can induce remarkable effects on sound propagation in the oceanic water column.[1,2] Meanwhile, in many cases the effects of horizontal refraction are non-negligible, and thus three-dimensional (3D) propagation should be considered.[1,3] In predicting sound fields in shallow water, where the terrain is probably very complicated and sound waves interact frequently with the bottom, it is necessary to consider both bottom elasticity and 3D propagation. Significant efforts have been made to model the two-dimensional (2D) sound propagation in oceans overlying elastic bottom (referred to as seismo-acoustic waveguides) based on the parabolic equation (PE) method. Readers are referred to a review article by Xu et al. for a summary of the relevant studies.[4] Also, there have been a variety of 3D models, including but not limited to 3D PEs,[4–13] 3D coupled-mode methods,[14–16] hybrid methods,[17–20] and analytical methods.[21–26] Nevertheless, studies that incorporate both bottom elasticity and 3D propagation are very scarce. The state-of-the-art PE model for this purpose is the 3D seismo-acoustic PE formulated in Refs. [8]–[12] and numerically implemented in Ref. [13] (referred to as Lee’s model). However, in its present form this model is not capable of incorporating horizontal refraction caused by asymmetric bathymetry (bathymetric refraction), because cylindrically symmetric bathymetry is assumed in its numerical algorithm.[12,13] The source images method[24–26] is also suited for this purpose, but is only applicable in the case of wedge-shaped environments. Another noticeable work is the elastic mode PE in the case of low shear modulus, which employs the parabolic approximation in the horizontal plane and the mode representation in the vertical direction.[20] In an attempt to enrich the studies of the above issue, here we present a new model for 3D seismo-acoustic wave propagation problems. Following Lee et al.,[8–13] we establish our model on the basis of PE. The differences between our model and Lee’s model, which are also arguably advantages of our model, are mainly in three aspects as follows.

First, our model is established in a Cartesian coordinate system, while Lee’s model is in a cylindrical coordinate system. An advantage of Cartesian coordinate system is that, in contrast to cylindrical coordinate system, it is more natural to incorporate axial asymmetry. In addition, Cartesian coordinate system does not suffer the problem of growth in computational grid size.

Second, the approaches to factorizing the full wave equations into one-way equations (or PEs) are different. Needless to say, the starting point of any PE method is to obtain one-way equations. For purely acoustic problems, one can directly factorize the Helmholtz equation into a one-way equation; whereas, for seismo-acoustic problems, the factorization of the elastic motion equations is tricky because of the existence of the first-order range derivatives of particle displacements.[27] Thus, extra effort should be made to rewrite the elastic motion equations into a readily factorable form. In Lee’s model, the elastic motion equations are decomposed into one P-wave equation and three S-wave equations, which can be easily factorized after applying the far-field approximation. An alternative approach is to perform linear operations on the elastic motion equations and meanwhile choose proper dependent variables.[27–29] A comparison between derivations in Refs. [12] and [27] reveals that the latter approach brings about much simpler expressions for the interface conditions. In this paper, two sets of 3D PEs are derived by using the latter approach, which is shown in Section 2.

The third difference lies in the ability to incorporate bathymetric refraction. For the sake of simplicity, the numerical scheme in this paper solely focuses on y-independent cases. Although the variations along the y direction are ignored, irregularity of the seafloor along the x direction is still able to induce bathymetric refraction. Whereas, in Refs. [12] and [13], the bathymetry is assumed to be cylindrically symmetric, and thus the bathymetric refraction is eliminated.

The rest of this paper is organized as follows. The 3D PEs are derived in Section 2. Then, the interface conditions are transformed into PE-adaptable expressions in Section 3. After that, a numerical scheme combining the 2D PE method and the Fourier synthesis technique is outlined in Section 4. Finally, accuracies of the numerical scheme in both range-independent and range-dependent cases are tested, and the azimuthal limitation of our solution is analyzed. In addition, the effects of horizontal refraction in the range-dependent examples are discussed.

2 Motion equations

According to Newton’s second law, the motion equations for an elastic medium in the absence of external forces are where ρ is the density, ω is the circular frequency, Ui is the ith component of the displacement vector (u, v, w), σij is the ith component of the stress vector acting on the face normal to the jth direction and the subscript j following the comma denotes the partial derivative in the jth direction. Note that Einstein’s summation convention is adopted in Eq. (1). According to Hooke’s law, the stress tensor, in the case of an isotropic elastic medium, is related to the strain tensor by where δij is the Kronecker delta, Δ is the divergence of displacement given by εij is the strain tensor defined by and λ and μ are the Lamé constants.

Substituting Eqs. (2)–(4) into Eq. (1), we obtain the following motion equations in terms of displacements: 0 = ( λ + 2 μ ) ∂ 2 u ∂ x 2 + ρ ω 2 u + ∂ μ ∂ y ∂ u ∂ y + μ ∂ 2 u ∂ y 2 + ∂ μ ∂ z ∂ u ∂ z + μ ∂ 2 u ∂ z 2 + λ ∂ 2 v ∂ x ∂ y + ∂ μ ∂ y ∂ v ∂ x + μ ∂ v 2 ∂ x ∂ y + λ ∂ w 2 ∂ x ∂ z + ∂ μ ∂ z ∂ w ∂ x + μ ∂ 2 w ∂ x ∂ z , 0 = ( λ + 2 μ ) ∂ 2 v ∂ y 2 + ρ ω 2 v + ∂ λ ∂ y ∂ v ∂ y + 2 ∂ μ ∂ y ∂ v ∂ y + ∂ μ ∂ z ∂ v ∂ z + μ ∂ 2 v ∂ z 2 + μ ∂ v 2 ∂ x 2 + ∂ λ ∂ y ∂ w ∂ z + λ ∂ 2 w ∂ y ∂ z + ∂ μ ∂ z ∂ w ∂ y + μ ∂ 2 w ∂ y ∂ z + ∂ λ ∂ y ∂ u ∂ x + λ ∂ 2 u ∂ x ∂ y + μ ∂ 2 u ∂ x ∂ y , 0 = ( λ + 2 μ ) ∂ 2 w ∂ z 2 + ρ ω 2 w + ∂ λ ∂ z ∂ w ∂ z + 2 ∂ μ ∂ z ∂ w ∂ z + μ ∂ 2 w ∂ x 2 + ∂ μ ∂ y ∂ w ∂ y + μ ∂ 2 w ∂ y 2 + ∂ λ ∂ z ∂ u ∂ x + λ ∂ u 2 ∂ x ∂ z + μ ∂ 2 u ∂ x ∂ z + ∂ λ ∂ z ∂ v ∂ y + λ ∂ 2 v ∂ y ∂ z + ∂ μ ∂ y ∂ v ∂ z + μ ∂ 2 v ∂ y ∂ z . Note that in Eqs. (5)–(7) we have ignored all terms that contain ∂λ/∂x, ∂μ/∂x, or ∂ρ/∂x since the waveguides are divided into a series of x-independent regions in the framework of our model.

It is clear that in equations (5)–(7) it is difficult to factor in the x direction because of the existence of the ∂2/∂x∂y and ∂2/∂x∂z derivatives of u, v, and w. To be suitable for the PE method, equations (5)–(7) are suggested to be rewritten in the form[27] where Q is a vector of dependent variables, and K is a matrix composed of spatial (y and z) operators and medium parameters.

2.1 Dependent variables (Δ, v, w)

The dependent-variable formulation which was first adopted in 2D elastic PE[28] is (Δ, w). Guided by this, a corresponding dependent-variable formulation for 3D elastic PE is developed here, namely, (Δ, v, w).

Differentiating Eqs. (5)–(7) with respect to x, y, and z, respectively, and then summing the results, we obtain 0 = ( λ + 2 μ ) ∂ 2 ∂ x 2 Δ + ρ ω 2 Δ + ( λ + 2 μ ) ∂ 2 ∂ y 2 Δ + ( λ + 2 μ ) ∂ 2 ∂ z 2 Δ + 2 ∂ ( λ + μ ) ∂ z ∂ ∂ z Δ + ∂ 2 λ ∂ z 2 Δ + 2 ∂ μ ∂ z ∂ 2 w ∂ x 2 + 2 ∂ μ ∂ z ∂ 2 w ∂ y 2 + 2 ∂ μ ∂ z ∂ 2 w ∂ z 2 + 2 ∂ 2 μ ∂ z 2 ∂ w ∂ z + ω 2 ∂ ρ ∂ z w + ∂ μ ∂ y ∂ Δ ∂ y + ∂ μ ∂ y ∂ 2 v ∂ x 2 + ∂ μ ∂ y ∂ 2 v ∂ y 2 + ∂ μ ∂ y ∂ 2 v ∂ z 2 + ∂ λ ∂ y ∂ 2 v ∂ y 2 + ∂ λ ∂ y ∂ Δ ∂ x − ∂ λ ∂ y ∂ 2 v ∂ x ∂ y . Equation (9) is difficult to factor in the x direction because of the existence of the last two terms. Whereas, these terms vanish as long as ∂ λ/∂ y = 0. Hence, below we simply assume the medium parameters to be constant over the entire horizontal plane. This assumption is reasonable if we do not consider very long-distance propagation (say hundreds of kilometers). Then equation (9) can be reduced to 0 = ( λ + 2 μ ) ∂ 2 ∂ x 2 Δ + ρ ω 2 Δ + ( λ + 2 μ ) ∂ 2 ∂ y 2 Δ + ( λ + 2 μ ) ∂ 2 ∂ z 2 Δ + 2 ∂ ( λ + μ ) ∂ z ∂ ∂ z Δ + ∂ 2 λ ∂ z 2 Δ + 2 ∂ μ ∂ z ∂ 2 w ∂ x 2 + 2 ∂ μ ∂ z ∂ 2 w ∂ y 2 + 2 ∂ μ ∂ z ∂ 2 w ∂ z 2 + 2 ∂ 2 μ ∂ z 2 ∂ w ∂ z + ω 2 ∂ ρ ∂ z w . Substituting Eq. (3) into Eqs. (6) and (7), respectively, we obtain It is clear that equations (10)–(12) can be directly factored into forward and backward components in the x direction. Of these three equations, equations (10) and (12) involve terms of only Δ and w. Hence, instead of imposing all the three equations, we can impose just Eqs. (10) and (12) to obtain the solutions for Δ and w, which will greatly improve computational efficiency. The solution for v can be obtained by substituting the solution for Δ into Eq. (11), but will not be discussed in this paper.

Rewriting Eqs. (10) and (12) in the form of Eq. (8), we have where L and M are depth operators, and q = (Δ, w)T.

Although the above derivations are for an elastic medium, equations (10)–(13) are also applicable to a fluid medium, as long as the shear modulus μ is set to be 0.

2.2 Dependent variables (Λ, vy, w)

A second dependent-variable formulation that is suitable for the PE method is (Λ, vy, w), where Λ = ux + vy, ux is the x derivative of u, and vy is the y derivative of v.

Differentiating Eqs. (5) and (6) with respect to x and y, respectively, and summing the results, we obtain 0 = ( λ + 2 μ ) ∂ 2 Λ ∂ x 2 + ( λ + μ ) ∂ 3 w ∂ 2 x ∂ z + ∂ μ ∂ z ∂ 2 w ∂ 2 x + ( λ + 2 μ ) ∂ 2 Λ ∂ y 2 + ( λ + μ ) ∂ 3 w ∂ 2 y ∂ z + ∂ μ ∂ z ∂ 2 w ∂ y 2 + ∂ ∂ z ( μ ∂ Λ ∂ z ) + ρ ω 2 Λ . Differentiating Eq. (6) with respect to y, rewritingEqs. (6) and (7) in terms of (Λ, vy,w), we obtain 0 = μ ∂ 2 v y ∂ x 2 + ∂ μ ∂ z ∂ v y ∂ z + μ ∂ 2 v y ∂ y 2 + ρ ω 2 v y + μ ∂ 2 v y ∂ z 2 + ( λ + μ ) ∂ 2 Λ ∂ y 2 + ( λ + μ ) ∂ 3 w ∂ y 2 ∂ z , 0 = μ ∂ 2 w ∂ x 2 + μ ∂ 2 w ∂ y 2 + ∂ λ ∂ z Λ + ( λ + μ ) ∂ Λ ∂ z + ∂ ∂ z [ ( λ + 2 μ ) ∂ w ∂ z ] + ρ ω 2 w . It is clear that equations (14)–(16) can be factorized in the x direction readily. Like the (Δ, v, w) formulation, here Λ and w can be calculated by using just Eqs. (14) and (16). Equations (14) and (16) can also be written in the form of Eqs. (13) or (8), except with different q, L, and M. Note that the media parameters are assumed to be constant over the horizontal plane in the above derivations just as in Subsection 2.1.

In the fluid layer, we still use the (Δ, v, w) formulation since Δ can be related to sound pressure p by Δ = −λp.

3. Interface conditions

Using the stair-step approximation in the x direction, the waveguide is divided into a series of x-independent regions. For the sake of simplicity, we assume the geometry in each region to be y-independent as well. In this way, there are only horizontal interfaces in each region, while between adjacent regions only vertical interfaces.

The ocean surface is assumed to be soft boundary, i.e., p|z = 0 = 0.

The horizontal fluid–solid interface conditions are continuity of vertical displacement, continuity of normal stress, and vanishing of tangential stress, which can be written as follows: where the subscripts w and b indicate the water layer and bottom layer.

Considering Eq. (12) for the fluid case (μ = 0) and substituting it into Eq. (17), we obtain Substitution of Eq. (2) into Eq. (18) yields Substitution of Eqs. (2) and (19) into Eq. (1) for i = 3 yields Equations (20)–(22) constitute a set of interface conditions adapted for the (Δ, v, w) formulation.

Similarly, for the (Λ, vy, w) formulation, the interface conditions in Eqs. (17)–(19) can be rewritten as

Conditions for vertical interfaces between adjacent regions can never be strictly satisfied in a PE algorithm.[30–33] The conventional way to handle the vertical interfaces in PE is to fix the fields with an approximated method when marching them across the vertical interfaces. The approximated methods to do this includes coordinate rotation,[30] coordinate mapping,[31] energy conservation,[32] and single scattering.[33]

4. Outline of numerical scheme: PE combined with Fourier synthesis

Since the inhomogeneity of the medium parameters in the horizontal plane has been ignored in Section 2, and the irregularity of the seafloor along the y direction has been ignored in Section 3, it is now natural to employ the Fourier transform (the wavenumber integration method) along the y direction to simplify the numerical implementation. In this way, the 3D problem is divided into a series of 2D components at different values of the out-of-plane wavenumber,[15,16,23] each of which can be solved with the 2D PE method.

Performing the Fourier transform on Eq. (13), we obtain where I is the identity matrix, ky is the out-of-plane wavenumber, and is defined by Factoring Eq. (26) and retaining only the component for the forward waves, we obtain where k0 is the reference wavenumber, and the operator X is defined by A solution of Eq. (28) that marches the field in the x direction is q ¯ | x + Δ x = e i k 0 Δ x I + X q | ¯ x , where Δx is the marching step. Approximating the exponentiated square-root operator by a Padé series,[34] we obtain q ¯ | x + Δ x ≈ e i k 0 Δ x ∏ j = 1 n I + α j , n X I + β j , n X q ¯ | x , where αj,n and βj,n are the Padé coefficients and n is the order of the Padé approximation. Equation (31) can be discretized in depth using Galerkin’s method.[29] The discretized system can be solved by Gaussian elimination.

Note that discretizing Eq. (31) at grids right above or below the fluid–elastic interface brings about unphysical quantities. The interface conditions in Eqs. (20)–(22) or Eqs. (23)–(25) should be used to eliminate these unphysical quantities, and thus they should be Fourier transformed as well. Performing the Fourier transform on Eqs. (20)–(22), we obtain ∂ ( λ w Δ ¯ w ) ∂ z + ρ w ω 2 w ¯ b = 0 , where the overbars denote the Fourier-transformed quantities. In the same way, equations (23)–(25) are transformed into In using Eqs. (32)–(34) or Eqs. (35)–(37) to eliminate the unphysical quantities, we can adopt either the approach in Ref. [29] or that in Ref. [33]. In this paper, we adopt the former one.

All the above steps are for marching from the beginning of an x-independent region to the end of it. In marching across the vertical interfaces between two adjacent regions, one of the approximate methods in Refs. [30]–[33] should be used.

For the model as a whole, it is necessary to provide the starting fields for each 2D component, i.e., . If we consider a point source in a realistic waveguide, then the source in each 2D component is a line source. The self-starter method[35] is employed to produce the starting fields for each 2D component.

Finally, the 3D solution can be obtained by performing the inverse Fourier transform on .

5. Numerical results
5.1 Range-independent case

Example A is for testing the 3D PE model in a range-independent ocean acoustic environment overlying an elastic half-space. The geometry is shown in Fig. 1(a). A 25 Hz point source is located at a depth of 100 m in a 200 m water column with constant sound speed cw = 1500 m/s and density ρ = 1.0 g/cm3. In the elastic bottom, the sound speeds for P- and S-waves are cp = 1700 m/s and cs = 800 m/s, respectively, and the density is ρb = 1.5 g/cm3. In addition, the attenuations for P- and S-waves in the bottom are βp = βs = 0.5 dB/λ. A stable and accurate solution is obtained using Δx = 20 m and Δz = 2 m. The sound pressure transmission loss (TL) contour on the horizontal x–y plane at a depth of 30 m is illustrated in Fig. 1(b), where the left half is calculated with the (Λ, vy, w) formulation while the right half with the (Δ, v, w) formulation. The perfect bilateral symmetry of Fig. 1(b) proves the nice agreement between the results of the two formulations. Meanwhile, the TL contour clearly reveals the ring-like interference pattern except in a narrow slice near the y axis, which is due to the angular limitation of Padé approximation. This issue will be discussed in detail in Subsection 5.2. The TL curves along the cut-lines marked in Fig. 1 are plotted in Fig. 2 and compared with the reference solutions. The azimuth angles φ of the cut-lines are 0°, 45°, 76°, and 82.9°, respectively. The reference solutions are calculated with an axisymmetric 2D finite element method (FEM).[36] A fully 3D FEM solution for this test case is still difficult to reach due to computer memory limitation, but this task will become easier in the near future. It is shown in Fig. 2 that the TL curves from the 3D PE model are all in excellent agreement with the reference solutions, except the one along φ = 82.9°. Thus, the presented model offers high accuracy at up to φ = 76° in this example.

Fig. 1. (a) Geometry of the waveguide in Example A. (b) TL contour of sound pressure on the horizontal plane at a depth of z = 30 m. The solid lines are along the azimuthal directions of φ = 0°, 45°, 76°, and 82.9°, respectively.
Fig. 2. TL curves of sound pressure at a depth of 30 m and along azimuths of 0°, 45°, 76°, and 82.9° in Example A.
5.2 Wide-angle capability

Angular limitation of the PE model arises from the use of Padé approximation, and can be estimated by evaluating the phase error of Padé approximation.

We write the solution of the Fourier-transformed dependent variables as a sum of normal modes: where Ψm and kxm are respectively the eigenfunction vector and the horizontal wavenumber of the mth mode. Note that both Ψm and kxm are functions of ky.

From Eq. (38), we obtain the marching solution for each mode, On the other hand, substitution of Eq. (38) into Eq. (30) yields the marching solution in another form, Comparison between Eqs. (39) and (40) reveals that the depth operator X is related to modal horizontal wavenumber by Substituting Eqs. (41) and (38) into Eq. (31), we have Comparison between Eqs. (42) and (39) reveals that the use of Padé approximation is equivalent to taking an approximation for kxm, i.e., where is defined by For simplicity, below we assume that the sound speed is constant in the water layer and choose k0 = k, where k is the wavenumber in the water layer.

It is natural to relate the eigenvalues to propagation directions of the modes as follows: where φm and θm are the azimuth and pitch angles of the mth mode (see Fig. 3). Following Jensen et al.,[1] we define the phase error of Padé approximation as More generally, we can drop the subscript m in Eqs. (45) and (46). With the values of ε, kΔx, and θ fixed, the value of φ is determined by Eq. (46), which is named limitation angle here. The value of kΔx is proportional to Δx/d, where d represents the wavelength in the water layer. Normally, in PE algorithms Δx is usually less than d/2.

Fig. 3. (color online) Definitions of the azimuth and pitch angles.

Without loss of generality, we assume θ = 0 (in the far field of shallow water, θ is approximately 0) and give the curves of ε versus φ. The results are displayed in Fig. 4 for the cases n = 4, 6, and 8. Here Δx is set to be d/3. The horizontal curve in Fig. 4 represents the error tolerance chosen as 0.0002[1] and its intersections with the ε–φ curves represent the corresponding limitation angles. It is noticeable that the speed at which the limitation angle approaches to 90° becomes slower as the value of n increases. Note that for each case of n we consider two subcases with and without one stability constraint. The stability constraint is used to impose a proper attenuation on each of the evanescent modes to ensure the stability of the numerical implementation, at the sacrifice of introducing slight artificial attenuation to the trapped modes.[37] In most cases, a single stability constraint is enough to ensure stability; while for the case with extremely low shear modulus or with a very thin elastic layer, two stability constraints are required. In all the numerical examples in this paper adopted is one stability constraint in Padé approximation. An interesting phenomenon in Fig. 4 is that the larger the value of n, the less the accuracy loss by a stability constraint will be.

Fig. 4. (color online) Plots of phase error ε versus azimuth φ for Padé approximations for n = 4, 6, and 8. The quantity ns denotes the number of stability constraints. Intersections between the curves and the horizontal line of ε = 0.0002, which represents the error tolerance, shows the limitation angle in each case.

For clarity, listed in Table 1 are the limitation angles for the subcases with a single constraint condition in Fig. 4, as well as those for the scenarios of Δx = d/6 and d (with all other settings unchanged). The results in Table 1 imply that the use of smaller marching steps brings about better wide-angle capability. Note that the settings of Δx = d/3 and n = 8 have been used in Example A in Subsection 5.1. According to Table 1, the limitation angle for this case is 75.2°, which is in accord with the results shown in Fig. 2.

Finally, it is pertinent to admit that the convergence of the procedure for calculating Padé coefficients, which is described in Refs. [34] and [37], becomes problematic as n increases.[38] This is why we do not involve the results for n > 10 in Table 1.

Table 1.

Limitation angles of Padé approximation with different values of n and Δx/d.

.
5.3. Range-dependent case and horizontal refraction

Example B is for considering the geometry as shown in Fig. 5(a). The depth is 200 m at x = 0 m and decreases linearly to 0 m at x = 4 km. A 25 Hz source is located at a depth of 100 m. This scenario can be regarded as a 3D version of the ASA problem.[39] The parameters of the water layer are the same as those in Example A. At the bottom, we consider two cases of P- and S-wave speeds: (i) cp = 1700 m/s, cs = 300 m/s; (ii) cp = 1700 m/s, cs = 800 m/s. The attenuations for P- and S-waves are both 0.5 dB/λ, and density of the bottom medium is 1.5 g/cm3. In both cases, Δx = 20 m, Δz = 1 m, the order of Padé approximation is 10, and the formulation of (Δ, v, w) is used. The range dependence is handled by the mapping approach.[31] The mapping approach is of high accuracy when the slope is sufficiently small.[31,40] The advantage of the mapping approach is its briefness. A more accurate solution might be obtained by using the single-scattering correction, and this may be implemented and tested in future. Since there is no horizontal refraction in the along-slope direction (φ = 0°), we can validate the 3D PE solution for Example B with a 2D benchmark solution along this specific direction. We display in Fig. 5(b) the 3D PE solutions at a depth of 30 m along φ = 0° together with the corresponding FEM solutions. Solutions from the two methods agree well with each other, which verifies that the 3D PE solutions of this range-dependent case are of high accuracy. This also indicates that backscattering along the track of the TL curve is negligible and that the mapping approach is of sufficient accuracy for handling the sloping interfaces here. In order to check the 3D PE solutions in all azimuths, we compare them with the source images solutions.[24–26] Slices of the 3D PE solutions along the cut-plane bisecting the wedge are plotted in Fig. 6 and compared with the corresponding source images solutions. We observe that except for the areas beyond the limitation angle, the 3D PE solutions agree well with the source images solutions.

Fig. 5. (color online) (a) Geometry of the waveguide in Example C. (b) TL curves of sound pressure at a depth of 30 m along the φ = 0° direction for Cases 1 (upper) and 2 (lower) in Example B. The lower two curves have been shifted downward by 30 dB to avoid overlapping.
Fig. 6. TL contours of sound pressure along the cut-plane bisecting the wedge for Cases 1 (upper) and 2 (lower) in Example B. The panels in the left column ((a), (c)) shows the results calculated using the 3D PE model, while the panels in the right column ((b), (d)) refer to the results calculated using the source images method. The solid lines denote the limitation angle of the 3D PE solutions, i.e., 78.5°.

Now we analyze the effects of horizontal refraction in Example B. The TL curves at a depth of 30 m along φ = 45° and 63.4° for both Cases 1 and 2 are plotted in Fig. 7. Note that in these plots the horizontal axis is x instead of the horizontal distance (x2 + y2)−1/2. It can be seen that deviations between the 3D solutions and the 2D solutions along φ = 45° are not remarkable, implying that horizontal refraction effects are weak along φ = 45° in these two cases. It is for sure that in areas where φ < 45° horizontal refraction effects are even weaker. Whereas, as shown in Figs. 7(b) and 7(d), horizontal refraction effects along φ = 63.4° are much stronger. With increasing azimuth φ, the importance of horizontal refraction is believed to keep growing. Interestingly, all the 3D solutions displayed in Fig. 7 decay faster than the corresponding 2D solutions. This is in accordance with the hyperbola-like edges (which should be nearly straight lines in results of N × 2D calculations) of shadow zones as shown in Fig. 6. The hyperbola-like pattern can be explained in a perspective of the ray theory.[3]

Fig. 7. TL curves of sound pressure at a depth of 30 m along φ = 45° and 63.4° for Case 1 ((a), (b)) and Case 2 ((c), (d)) in Example B. The panels in the left column ((a), (c)) are for φ = 45°, while the panels in the right column ((b), (d)) are for φ = 63.4°.

In Example C, we analyze the horizontal refraction effects in the waveguide shown in Fig. 8(a), where the bottom interface is a polyline along the x direction. The depth is 300 m at x = 0 m and decreases linearly to 100 m at x = 4 km. In the area where x > 4 km, the depth is fixed at 100 m. In the elastic bottom, the sound speeds are cp = 3800 m/s, cs = 1800 m/s. Other parameters or settings are the same as those in Example B. This waveguide extends to infinity in the x direction, and thus offers the possibility to investigate horizontal refraction effects in relatively long distances. The 3D PE solutions at z = 30 m along the directions of φ = 0°, 45°, and 63.4° are plotted in Figs. 9(b)–9(d) and compared with the corresponding 2D PE solutions. It is shown that like the results in Example B, there is no horizontal refraction in the direction of φ = 0°. In fact, this can be easily proved as follows.

Fig. 8. (color online) (a) Geometry of the waveguide in Example C. (b)–(d) TL curves of sound pressure at a depth of 30 m along φ = 0°, 45°, and 63.4° in Example C.

Since our problem is symmetric with respect to the vertical plane of y = 0 (or φ = 0°), we have ∂p/∂y = 0 in this plane. Since ∂/∂φ = (x2 + y2)−1/2 (cos φ·∂/∂y − sin φ∂/∂x), we also have ∂p/∂φ = 0 in this plane. Then, the 3D Helmholtz equation for the water column in this plane can be reduced to a form that is identical to that for the cylindrically symmetric 2D case, and so can the boundary conditions. Therefore, the 3D solution in the plane of φ = 0° is equivalent to a corresponding 2D solution, i.e., there is no horizontal refraction in the plane of φ = 0°.

Now we continue to investigate horizontal refraction effects in other azimuths. Once again, we observe that the horizontal refraction effects in the direction of φ = 63.4° are much stronger than in the direction of φ = 45°. In addition, it can be seen that the discrepancies between the 3D and 2D solutions, basically lags of TL maxima, accumulate gradually with distance increasing in the sloping bottom area, which is also true in the previous example. The reason is that in our examples the horizontal refraction arises during reflection of sound waves at the sloping interface, and the number of reflections increases with distance increasing. Interestingly, lags of TL maxima still keep growing in the flat bottom area. This is because the rays, which have been refracted into larger azimuths in the sloping area,[16] continue to propagate to even larger azimuths after leaving the sloping area.

6. Conclusions

A 3D PE model for seismo-acoustic propagation is presented. By carefully choosing dependent variables and performing linear operations on the elastic motion equations, two sets of seismo-acoustic PEs are derived, i.e., the (Δ, v, w) formulation and the (Λ, vy, w) formulation. Numerical results demonstrate the feasibility of both of the two sets of PEs. In addition, bathymetric refraction due to the x-dependent sea bottom interface is well depicted by this model. It is illustrated that bathymetric refraction may induce remarkable effects on sound propagation in our simulation cases.

We also analyze the wide-angle capability of our model in detail. It is shown that the combination of 2D PE and Fourier transform along the y-direction gives rise to azimuthal limitation on the horizontal plane for the resulting 3D solution. Pertinently, even for a fully Cartesian 3D PE model,[6,7] accuracy of the solution along or near the longitudinal direction (the horizontal direction that crosses the marching direction) is not assured because it is inevitable to employ asymptotic approximation on the square root operator. A possible solution to this problem is to conduct multiple calculations with the corresponding marching directions along different azimuths, retaining only the results that are assured accurate in each calculation, and then assembling them. However, for now this approach is not possible to implement in our model, since we have not incorporated longitudinal variation of bathymetry.

Great efforts are still required in order to provide this model with the ability to handle arbitrarily irregular fluid–elastic interfaces. For example, by using artificial absorption layers or perfectly matched layers on the two sides of the longitudinal direction, we can probably solve the 3D PE in Eq. (13) directly without employing Fourier transform along the y-axis, and thus relieve the limitation of y-invariant bathymetry. In addition, placing the fluid–elastic interface in the midway of two grid points[33] would also bring convenience to the extension to arbitrary bathymetry.

PE models have been actively used as a tool to predict sound fields in range-dependent oceans. Apparently, it will continue to play this role due to its high efficiency compared with the direct methods, e.g., FEM. Once the capability of 3D seismo-acoustic PEs is well expanded, the predicted results for sound fields will be closer to truth than purely acoustic PEs.

Reference
[1] Jensen F B Kuperman W A Porter M B Schmidt H 2012 Computational Ocean Acoustics New York Springer 457 530
[2] Zhang H G 2010 Research on Modeling and Rule of Infrasound Propagation in Shallow Water Ph. D. Dissertation Harbin Harbin Engineering University in Chinese
[3] Katsnelson B G Petnikov V G Lynch J 2012 Fundamentals of Shallow Water Acoustics New York Springer 124 138
[4] Xu C X Tang J Piao S C Liu J Q Zhang S Z 2016 Chin. Phys. 25 124315
[5] Sturm F 2016 J. Acoust. Soc. Am. 139 263
[6] Lin Y T Duda T F Newhall A E 2013 J. Comput. Acoust. 21 1250018
[7] Xu C X Piao S C Yang S E Zhang H G Tang J 2016 Acta Acustica 41 477 in Chinese
[8] Lee D Nagem R J Resasco D C 1997 J. Comput. Acoust. 5 157
[9] Lee D Nagem R J Resasco D C Chen C F 1998 Appl. Anal. 68 147
[10] Lee D Nagem R J Resasco D C Chen C F 1999 Appl. Anal. 75 183
[11] Tony W H Cheu S C Chen C F Chiang T P Lee D 1999 J. Comput. Acoust. 7 185
[12] Nagem R J Lee D 2002 J. Comput. Acoust. 9 421
[13] Hsieh L W 2005 Modeling 3D Wave Propagation in the Ocean Coupled with Elastic Bottom and Irregular Interface Ph. D. Dissertation Taiwan National Taiwan University
[14] Luo W Y Zhang R H Schmidt H 2011 Sci. China: Phys. Mech. Astron. 54 1561
[15] Luo W Y Yu X L Yang X F Zhang Z Z Zhang R H 2016 Chin. Phys. 25 124309
[16] Mo Y X Piao S C Zhang H G Li L 2014 Acta Phys. Sin. 63 214302 in Chinese
[17] Peng Z H Li F H 2005 Acta Acustica 30 97 in Chinese
[18] Qin J X Katsnelson B Peng Z H Li Z L Zhang R H Luo W Y 2016 Acta Phys. Sin. 65 034301 in Chinese
[19] Trofimov M Yu Kozitskiy S B Zakharenko A D 2015 Wave Motion 58 42
[20] Trofimov M Yu Kozitskiy S B Zakharenko A D 2015 Proc. Meetings Acoust. 24 070025
[21] Petrov P S Sturm F 2016 J. Acoust. Soc. Am. 139 1343
[22] Petrov P S Prants S V Petrova T N 2017 Phys. Lett. 381 1921
[23] Luo W Y Zhang R H 2015 Sci. China: Phys. Mech. Astron. 58 594301
[24] Deane G B Buckingham M J 1993 J. Acoust. Soc. Am. 93 1319
[25] Petrov P S 2017 MATLAB code http://oalib.hlsresearch.com/ThreeD/ (Last viewed 31 July 2017)
[26] Tang J 2017 MATLAB code http://oalib.hlsresearch.com/ThreeD/(Last viewed 31 July 2017)
[27] Jerzak W Siegmann W L Collins M D 2005 J. Acoust. Soc. Am. 117 3497
[28] Greene R R 1985 J. Acoust. Soc. Am. 77 1991
[29] Collins M D 1989 J. Acoust. Soc. Am. 86 1459
[30] Outing O A Siegmann W L Collins M D Westwood E K 1990 J. Acoust. Soc. Am. 87 1035
[31] Collins M D Dacol D K 2000 J. Acoust. Soc. Am. 107 1937
[32] Collins M D 1993 J. Acoust. Soc. Am. 94 975
[33] Collins M D Siegmann W L 2015 J. Acoust. Soc. Am. 137 492
[34] Collins M D 1993 J. Acoust. Soc. Am. 93 1736
[35] Cederberg R J Collins M D 1997 IEEE J. Ocean Eng. 22 102
[36] COMSOL Multiphysics Reference Manual 5.2 (COMSOL AB, Stockholm, Sweden)
[37] Collins M D 1991 J. Acoust. Soc. Am. 89 1050
[38] Milinazzo F A ZalaGary C A Brooke H 1997 J. Acoust. Soc. Am. 101 760
[39] Goh J T Schmidt H Gerstoft P Seong W 1997 IEEE J. Ocean Eng. 22 226
[40] Outing D A Siegmann W L Collins M D 2007 IEEE J. Ocean Eng. 32 620