Development and application of a pre-corrected fast Fourier transform accelerated multi-layer boundary element method for the simulation of shallow water acoustic propagation

: Because of the complexities associated with the domain geometry and environments, accurate prediction of acoustics propagation and scattering in realistic shallow water environments by direct numerical simulation is challenging. Based on the pre-corrected Fast Fourier Transform (PFFT) method, we accelerated the classical boundary element method (BEM) to predict the acoustic propagation in a multi-layer shallow water environment. The classical boundary element method formulate the acoustics propagation problem as a linear equation system in the form of [A]{x} = {b}, where [A] is an N × N dense matrix composed of inﬂuence coe ﬃ cients. Solving such linear equation system requires O (N 2 / N 3 ) computational cost for iterative / direct methods. The developed method, PFFT-BEM, can e ﬀ ectively reduce the computational e ﬀ orts for direct numerical simulations from O (N 2~3 ) to O (Nlog N) , where N is the total number of boundary unknowns. To numerically simulate the sound propagation in a shallow water environment, we applied the ﬁrst-order non-reﬂecting boundary condition in the truncated numerical domain boundary to eliminate the errors due to reﬂected waves. Multi-layer coupled formulation was used to include the environment inhomogeneity in PFFT-BEM. multiple convergence tests on the number of layers and elements, we validated and quantiﬁed the accuracy of PFFT-BEM. To demonstrate the usefulness and capability of the developed PFFT-BEM, we simulated three-dimensional (3D) underwater sound propagation through 3D geometries to check the e ﬃ cacy of the established classical method: the 3D Parabolic equation model. Finally, PFFT-BEM was employed to simulate sound propagation through a complex multi-layer shallow water environment with internal waves. The “3D + T” results obtained by PFFT-BEM compared well with the physical test, thereby proving the capability and correctness of this method. with the PFFT-BEM showed strong variability as the internal wave train moved. The magnitude and trend of these numerical variations agree fairly well with the ﬁeld experiments. This indicates that the 3D + T PFFT-BEM is applicable for the direct simulation of sound propagation in a complex 3D shallow water environment. We concluded that the developed direct numerical method, PFFT-BEM, is very useful to quantify the accuracy of exiting numerical methods and simulate sound propagation in realistic multi-layer shallow water environments.


Introduction
The accurate forecasting of acoustics wave propagation and scattering in shallow water environments is important in naval operations and underwater communications. A large number of theoretical studies have been carried out to model the sound scattering effects in complex underwater environments in three dimensions (3D). The 3D coupled mode method could produce efficient and accurate prediction of the sound field for specified geometries (e.g., 3D axisymmetric geometry [1,2], longitudinal invariant geometry [3]). The 3D ray model has been proposed for high frequency sound propagation ( [4]). The 3D parabolic equation (PE) method is another powerful tool for the prediction of 3D sound propagation (e.g., [5]). Although several newly developed PE models, including two-way scattering, have been developed, generally for two-dimensional (2D) cases [6,7] and for 3D cases formulated in cylindrical coordinates [8], PE is mostly used for the prediction of one-way sound propagation as backward propagation is generally neglected in this method.
For direct numerical simulations (i.e., discretization methods without any assumptions), despite the continued development of computing capabilities, it is a challenge to perform 3D simulations of ocean acoustics in realistic environments due to the high computational costs. Volume discretization methods, such as the finite element method (FEM), finite volume method (FVM), and finite difference method (FDM), require fine meshes over the entire volume of the computational domain. Due to the large computational domain (~O(10) km) and the small acoustic wavelength (~O(1~10) m), the 3D mesh entails a large computational cost. By only discretizing the domain, the boundary discretization method or boundary element method (BEM) has a much lower computational cost than the volume discretization methods. By discretizing the domain boundaries and interfaces into boundary elements, BEM numerically formulate the acoustics propagation problem as a  (10 6~8 ). This induces a large computational cost and limits the use of BEM in real shallow water acoustics. Another major drawback of BEM is that classical BEM can only formulate the linear equation system for underwater medium with homogeneous properties (e.g., density and sound speed). However, complex shallow water environments with density and sound speed variations cannot be represented by a single homogeneous fluid layer.
In this study, to overcome the computational cost difficulty in BEM, we used the pre-corrected fast Fourier transform (PFFT) approach to accelerate the BEM solver for the calculation of acoustics propagation in a shallow water environment. The pre-corrected fast Fourier transform (PFFT) approach has been used to increase the speed of the BEM solver from (N 2 ) to (N log N) in different areas (e.g., [9][10][11]). By using PFFT, we can reduce the computational cost in the order of O(10 6~8 ). Compared with existing studies using PFFT for underwater acoustics problems (i.e., [11]), the present study made the following contributions. First, to overcome drawbacks due to homogenous assumptions, we simulated the medium inhomogeneity by coupling of multiple homogeneous sub-layers at their interfaces. Through a numerical test, we checked the convergence of the multilayer PFFT-BEM (Version 2019, Massachusetts Institute of Technology, Cambridge, MA, USA, 2019) with respect to the number of sub-layers to validate the accuracy of the multi-layer formulation in our method. Secondly, the accuracy of the BEM method will become progressively worse when moving closer to the domain boundary if an inaccurate far field radiation boundary condition is imposed at the truncated computational domain boundary. In the present method, we applied accurate first-order non-reflecting boundary conditions (NRBCs) in the truncated domain boundary. We utilized the 2D normal mode method to obtain the propagation wave modes near the domain boundary, which are needed for the application of the first-order NRBC. The accuracy of the truncated boundary condition was validated against existing benchmark solutions. Due to the full-field simulation capability, our method could serve as a benchmark solution for comparing and assessing the performance of the classical established model: the 3D Parabolic Equation model. To do that, we conduct 3D simulations to benchmark low-mid frequency acoustics over kilometer ranges for the prediction of 3D sound scattering in three complex underwater environments: a 3D Acoustical Society of America (ASA) wedge, a 3D seamount, and an underwater canyon. We also demonstrate the practical 3D+T simulation capability of PFFT-BEM for underwater acoustics scattering problems. Here, 3D+T refers to time-stepped, frozen-field modeling of underwater sound. Because sound moves at 1.5 km/s in the ocean, pulsed sound traveling hundreds of kilometers could, in principle, be modeled by updating the medium condition over the excursion time of hundreds of seconds as the sound moves away from the sound sources. Using 3D+T PFFT-BEM (Version 2019, Massachusetts Institute of Technology, Cambridge, MA, USA, 2019), we simulate sound scattering variations due to the change in medium properties caused by internal traveling waves, and we compare the numerical results with field tests.

Formulation of the Boundary Element Method
Consider the following 3D Helmholtz wave equation: where x = (x,y,z), p(x) is acoustics pressure, and ρ, V, k is the density, volume and wavenumber of the medium, respectively. At the interface between the medium and atmosphere, pressure release boundary condition is satisfied as: where S F is the free surface. The Sommerfeld radiation condition is defined at S ∞ : where n is the normal direction on the boundaries. For shallow water environments, a multilayer model is adopted to account for the inhomogeneous medium properties. In this model, we assume the density ρ and sound speed c in each layer to be constant at every instant of time. A sketch of the multilayer domain is shown in Figure 1.
the change in medium properties caused by internal traveling waves, and we compare the numerical results with field tests.

Formulation of the Boundary Element Method
Consider the following 3D Helmholtz wave equation: where = (x,y,z), ( ) is acoustics pressure, and , V, k is the density, volume and wavenumber of the medium, respectively. At the interface between the medium and atmosphere, pressure release boundary condition is satisfied as: where is the free surface. The Sommerfeld radiation condition is defined at : where is the normal direction on the boundaries. For shallow water environments, a multilayer model is adopted to account for the inhomogeneous medium properties. In this model, we assume the density and sound speed in each layer to be constant at every instant of time. A sketch of the multilayer domain is shown in Figure 1. Then, at every instant of time, the homogeneous Helmholtz equation governing the sound field in each layer can be expressed as where i is the number of layers with i = 1,....,I, and is the volume of layer i. Continuity of pressure and particle velocity are required across the interface ( ) of neighboring layers and (including at the penetrable bottoms). We define the bottom and top boundaries of the layer i as the Then, at every instant of time, the homogeneous Helmholtz equation governing the sound field in each layer can be expressed as where i is the number of layers with i = 1, ∆, I, and v i is the volume of layer i. Continuity of pressure and particle velocity are required across the interfaceŜ i+1 (Š i ) of neighboring layers v i+1 and v i (including at the penetrable bottoms). We define the bottom and top boundaries of the layer i as the symbols ' From Green's second theorem, the boundary integral equation (BIE) for the pressure ( ) in h layer is derived from Equation (1) as ', respectively. Continuity boundary condition across the interface can be expressed as where S i represents the top and bottom surface of each layer, S Ri is the radiation boundary in the far field, G(x;ξ) = 1 |x−ξ| e ik|x−ξ| and ξ is a dummy variable representing the point on the layer boundary S. Because the Green function G(x;ξ) satisfies the radiation condition Equation (3) itself, the bottom boundary S ∞ does not need to be meshed in order to satisfy the far-field radiation condition. Similarly, due to the attenuation effect in the bottom layer, the radiation condition Equation (3) is also automatically satisfied at the side boundary S RI in the bottom layer. However, due to the finite truncation of the side domain boundary, a proper radiation condition needs to be applied on S Ri in the fluid layer (i.e., i I). This is further discussed in Section 2.2.
By taking the limit of x → S i and using the far-field boundary conditions, pressure release boundary condition, and continuity interface boundary conditions, Equation (7) for each layer can be coupled into one large set of integral equations as where α is the solid angle and β i = ρ i /ρ i+1 . Using the constant panel method, we discretize the top and bottom surface of the ith layer into N i quadrilateral elements and the side boundary, where the radiation boundary condition is applied in layer i, into N Ri quadrilateral elements. On each element, quantities such as the acoustic pressure p and the normal gradient of the pressure p n are assumed to be constant. This transfers Equation (8) into a set of linear algebraic equations of the form where E j represents the jth constant element, and N is the total number of elements.

Non-Reflecting Boundary Condition
The accuracy of PFFT-BEM will become progressively worse when moving closer to the domain boundary if the accurate far-field radiation boundary condition cannot be imposed exactly at the truncated computational domain boundary. The non-reflecting boundary condition (NRBC) [12] is a popular method used to eliminate these errors. This method first expresses the solutions of the Helmholtz equation in terms of a convergent expansion of (1/|r|), where r denotes the radial coordinates. The radiation condition at S Ri at a certain order in terms of 1/|r| can then be obtained in terms of the truncated domain size and the propagating wavenumber [13]. In the present study, the first-order Appl. Sci. 2020, 10, 2393 5 of 20 NRBC is applied at the boundary of the truncated computational domain in PFFT-BEM. In a radial coordinate, the first-order NRBC can be expressed as [13]: and p m is the acoustic pressure due to the mth propagating modes with where α m is the damping coefficient, which for 3D is defined as: where k rm is the horizontal wavenumber for the propagating mode m. Substituting Equation (13) into Equation (12) yields the following boundary condition at the truncated domain boundary: In this study, we compute the acoustic pressure for the propagating mode in the far-field p m and the corresponding wavenumber k rm (α m ) using the well-known normal mode method. A detailed derivation and good review of the normal mode method can be found in Chapter 5 of [14]. A similar method has been applied for 2D finite element analysis for shallow water acoustics problems [15]. However, for 3D boundary element method, the effect of the first-order NRBC is still unknown and need to be investigated. To demonstrate the efficiency of first-order NRBC, we compare the numerical results using Equation (15) at the truncated boundary with the results based on the radiation condition with k rm = k. The latter can be expressed as Equation (16) at the truncated domain boundary: Substituting Equation (15) into Equation (8) yields If we put Equation (17) in each layer into one single equation, the entire linear system of equations, Equation (17), can then be expressed as where [A] is an N × N dense and nonsymmetric matrix of influence coefficients, {x} is the vector of unknowns (p and p n ) on the sub-domain surfaces, and {b} is the known vector due to the point source.
An iterative solver such as the generalized minimum residual (GMRES) method could also be applied to solve the system shown in Equation (18). GMRES normally requires O(N 2 ) operations. The most expensive step in the iterative solver is to compute the product of ([A]{x}), which requires a minimum of O(N 2 ) operations. Physically, each multiplication F s (or F d ) in [A]{x} represents the influence from element j to element i on the domain boundary. As a result, N elements contain N 2 interactions. The O(N 2 ) effort of the iterative methods limits the applicability of the BEM to 3D underwater acoustic problems involving a large number (N) of unknowns (e.g., large domains or high frequency).

Evaluation of [A]{x} Pre-Corrected Fast Fourier Transform Method
An alternative approach is to use the pre-corrected fast Fourier transform (PFFT) accelerated algorithm, which can evaluate the product [A]{x} with O(NlogN) operations without explicit construction of the matrix [A]. The PFFT-BEM, which reduces the computational effort to O(NlogN), has been applied in various fields, such as the design of interaction circuits, the modeling of microelectromechanical systems, wave hydrodynamics problems, and acoustic scattering problems (e.g., [9][10][11]).
Here, we briefly introduce the PFFT technique. For illustration, we take a radiation sphere as an example in each iteration step of GMRES. As shown in Figure 2a, using the boundary element method, we discretize the surface of the sphere into N small triangular elements with surface sizes of S j , j = 1, . . . , N. We assume p n is unknown on the sphere surface and we try to evaluate Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 21

Validation of Multi-Layer PFFT-BEM Solutions
In this section, we first demonstrate the correctness and efficacy of PFFT-BEM by existing benchmark solutions by an ideal waveguide shown in Figure 3. Both the top and bottom boundary satisfy pressure release boundary condition with no bottom layer. The sound source is at a depth of 36 m with frequency at 20 Hz. The other water properties of the waveguide are shown in Figure 3. The computational domain size spans 3 km × 3 km × 100 m with uniform meshing on both the top and bottom surfaces using ∆ = /12.  In the step of grid definition, we define N g = N x × N y × N z grids in the 3D domain covering the entire surface of the sphere, as shown in Figure 2b. We refer to the cube made up by nearby grids as a cell. Based on a previous study using PFFT [10], the size of the cells h and the elements ∆l are assumed to be the same.
In the step of projection, the normal gradient of the pressure p n at element j is projected to the grids of the cell that contain element j. The formula for projection is where H m is the spatial interpolation function (see, e.g., [10]) for the m th grid of the cell covering element j, and q jm is the strength of the projected (p n ) j at the m th grid of the cell covering element.
In the step of convolution, we evaluate the total influence at grid point l due to all the projected (p n ) j at other grids as where G l−m = 1 |x l −ξ m| e ik|x l −ξ m | . We note that the summation in Equation (21) is in convolution form and can be computed using FFT with a computation operation of O(N g log N g ) and a memory requirement In the step of interpolation, we evaluate the total influence (IF s ) ij at element i (i = 1, . . . , N) due to all (p n ) j using the interpolation function H m and the grid values (IF s ) lm . The steps of projection convolution and interpolation are demonstrated in Figure 2c as steps (1) It is clear that the accuracy of the interpolation of the influence between elements j and i deteriorates when where h is the size of the cells. In order reduce the error in the near-field, the step of near-field correction is conducted by applying where represents the correction, and (IF s ) ij N and (IF s ) ij NF represent the total and near-field influences obtained by direct calculation and the convolution method. Based on earlier studies (e.g., [10]), the size of the near field used in the present study was chosen to be equal to 3h, which has been shown to be accurate. In summary, the evaluation of Other properties of our PFFT-BEM method include (i) a mesh-neighbor (MN) preconditioner for fast iterative convergence (2+ orders faster than the solver without a preconditioner) [16] and (ii) a high performance library PETSC [17] to optimize the implementation of the PFFT-BEM on massive parallel HPC platforms to achieve linear CPU scalability up to O(1000) CPUs. The details of these properties are omitted here and can be found in [16,17].

Validation of Multi-Layer PFFT-BEM Solutions
In this section, we first demonstrate the correctness and efficacy of PFFT-BEM by existing benchmark solutions by an ideal waveguide shown in

Validation of Multi-Layer PFFT-BEM Solutions
In this section, we first demonstrate the correctness and efficacy of PFFT-BEM by existing benchmark solutions by an ideal waveguide shown in Figure 3. Both the top and bottom boundary satisfy pressure release boundary condition with no bottom layer. The sound source is at a depth of 36 m with frequency at 20 Hz. The other water properties of the waveguide are shown in Figure 3. The computational domain size spans 3 km × 3 km × 100 m with uniform meshing on both the top and bottom surfaces using ∆ = /12.  In Figure 5, we show the convergence of the numerical results towards the wavenumber integration method results with respect to the boundary element size. We compare the acoustical pressure by PFFT-BEM with from wavenumber integration method at z = 36m. The errors shown in Figure 5 are defined in Equations (24) and (25). As shown in Figure 5, PFFT-BEM results converge clearly towards the wavenumber integration method solution as the element size decreases. We then check the efficiency of the first-order NRBC at the truncated numerical domain boundary. As shown in Figure 6, PFFT-BEM results clearly improve towards the wave integration method results after radiation condition at the numerical boundary is changed from krm = k radiation condition (Equation (16)) to first-order NRBC (Equation (12)) with m = 2. In addition, the errors of In Figure 5, we show the convergence of the numerical results towards the wavenumber integration method results with respect to the boundary element size. We compare the acoustical pressure p PFFT−BEM by PFFT-BEM with p wim from wavenumber integration method at z = 36 m. The errors shown in Figure 5 are defined in Equations (24) and (25). As shown in Figure 5, PFFT-BEM results converge clearly towards the wavenumber integration method solution as the element size decreases.
Max error = max p PFFT−BEM − p wim p theory .
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 21 In Figure 5, we show the convergence of the numerical results towards the wavenumber integration method results with respect to the boundary element size. We compare the acoustical pressure by PFFT-BEM with from wavenumber integration method at z = 36m. The errors shown in Figure 5 are defined in Equations (24) and (25). As shown in Figure 5, PFFT-BEM results converge clearly towards the wavenumber integration method solution as the element size decreases. We then check the efficiency of the first-order NRBC at the truncated numerical domain boundary. As shown in Figure 6, PFFT-BEM results clearly improve towards the wave integration method results after radiation condition at the numerical boundary is changed from krm = k radiation condition (Equation (16)) to first-order NRBC (Equation (12)) with m = 2. In addition, the errors of We then check the efficiency of the first-order NRBC at the truncated numerical domain boundary. As shown in Figure 6, PFFT-BEM results clearly improve towards the wave integration method results after radiation condition at the numerical boundary is changed from k rm = k radiation condition (Equation (16)) to first-order NRBC (Equation (12)) with m = 2. In addition, the errors of numerical results using the k rm = k radiation condition become progressively larger near the truncated domain boundary. When the first order NRBC is used, there are no errors due to the boundary reflection which demonstrate the effectiveness of the first-order NRBCs used by PFFT-BEM for shallow water environments.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 21 numerical results using the krm = k radiation condition become progressively larger near the truncated domain boundary. When the first order NRBC is used, there are no errors due to the boundary reflection which demonstrate the effectiveness of the first-order NRBCs used by PFFT-BEM for shallow water environments. To demonstrate the computational efficiency of the PFFT-BEM, we vary the boundary element size to increase the number of N. As shown in Figure 7, both the memory and operational count of PFFT-BEM are proportional to O(N) and much smaller than the ones of the classical BEM, which is proportional to ( ). This result shows the computational advantage of the PFFT-BEM over the classical BEM for shallow water acoustic simulations. For higher sound source frequency (i.e., smaller acoustics wavelength), the computational cost, which is larger due to the element size, must be much smaller than the acoustics wavelength throughout the entire domain surface. From Figure 8, we show that the PFFT-BEM memory and computational cost is proportional to the sound frequency while it is also much smaller than classical BEM.  To demonstrate the computational efficiency of the PFFT-BEM, we vary the boundary element size to increase the number of N. As shown in Figure 7, both the memory and operational count of PFFT-BEM are proportional to O(N) and much smaller than the ones of the classical BEM, which is proportional to O N 2 . This result shows the computational advantage of the PFFT-BEM over the classical BEM for shallow water acoustic simulations. For higher sound source frequency (i.e., smaller acoustics wavelength), the computational cost, which is larger due to the element size, must be much smaller than the acoustics wavelength throughout the entire domain surface. From Figure 8, we show that the PFFT-BEM memory and computational cost is proportional to the sound frequency while it is also much smaller than classical BEM.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 21 numerical results using the krm = k radiation condition become progressively larger near the truncated domain boundary. When the first order NRBC is used, there are no errors due to the boundary reflection which demonstrate the effectiveness of the first-order NRBCs used by PFFT-BEM for shallow water environments. To demonstrate the computational efficiency of the PFFT-BEM, we vary the boundary element size to increase the number of N. As shown in Figure 7, both the memory and operational count of PFFT-BEM are proportional to O(N) and much smaller than the ones of the classical BEM, which is proportional to ( ). This result shows the computational advantage of the PFFT-BEM over the classical BEM for shallow water acoustic simulations. For higher sound source frequency (i.e., smaller acoustics wavelength), the computational cost, which is larger due to the element size, must be much smaller than the acoustics wavelength throughout the entire domain surface. From Figure 8, we show that the PFFT-BEM memory and computational cost is proportional to the sound frequency while it is also much smaller than classical BEM.  operational counts.
In Figure 9, we show the convergence of the numerical results with respect to the number of layers used in the PFFT-BEM. The sound profiles with different gradients, which are used for the convergence tests, are shown in Figure 9a. We compare the averaged difference of the acoustic pressure at z = 36 m between with layers ( ) in the water column and + 1 layers ( ) in the water column. The difference used for comparison is defined in Equation (26). Figure 9b shows that the numerical results converge when the number of layers increases. This result demonstrates that the PFFT-BEM can be used for shallow water environments with inhomogeneous sound profiles: We then demonstrate the use of the PFFT-BEM to realistic range dependent underwater acoustics problems. The same waveguide geometry is used but with the summer sound profile from [18]. The bottom is homogeneous with sound speed 1700 m/s, density 1.5 g/cm 3 and medium attenuation = 0.5 dB/ . Figure 10a shows the spatial variations of sound speed in the water column. In Figure 9, we show the convergence of the numerical results with respect to the number of layers used in the PFFT-BEM. The sound profiles with different gradients, which are used for the convergence tests, are shown in Figure 9a. We compare the averaged difference of the acoustic pressure at z = 36 m between with I layers (p I ) in the water column and I + 1 layers (p I+1 ) in the water column. The difference used for comparison is defined in Equation (26). operational counts.
In Figure 9, we show the convergence of the numerical results with respect to the number of layers used in the PFFT-BEM. The sound profiles with different gradients, which are used for the convergence tests, are shown in Figure 9a. We compare the averaged difference of the acoustic pressure at z = 36 m between with layers ( ) in the water column and + 1 layers ( ) in the water column. The difference used for comparison is defined in Equation (26). Figure 9b shows that the numerical results converge when the number of layers increases. This result demonstrates that the PFFT-BEM can be used for shallow water environments with inhomogeneous sound profiles: We then demonstrate the use of the PFFT-BEM to realistic range dependent underwater acoustics problems. The same waveguide geometry is used but with the summer sound profile from [18]. The bottom is homogeneous with sound speed 1700 m/s, density 1.5 g/cm 3 and medium attenuation = 0.5 dB/ . Figure 10a shows the spatial variations of sound speed in the water column.  Figure 9b shows that the numerical results converge when the number of layers increases. This result demonstrates that the PFFT-BEM can be used for shallow water environments with inhomogeneous sound profiles: We then demonstrate the use of the PFFT-BEM to realistic range dependent underwater acoustics problems. The same waveguide geometry is used but with the summer sound profile from [18]. The bottom is homogeneous with sound speed 1700 m/s, density 1.5 g/cm 3 and medium attenuation α = 0.5 dB/λ. Figure 10a shows the spatial variations of sound speed in the water column. The sound frequency is 60Hz located at 70 m below the sea surface. The computational domain is 2 km × 2 km ×

Model Comparisons of the PFFT-BEM with the 3D Parabolic Equation Method
To demonstrate the usefulness of PFFT-BEM, PFFT-BEM is compared with 3D Parabolic equation (PE) models to check the correctness and performance of different PE schemes. We first briefly review the formulation and assumptions of the 3D PE theory. Then, we present three numerical examples for comparison: ASA wedge, underwater seamount, and Gaussian canyon.

Formulation and Assumptions of the PE Method
In order to demonstrate the comparison between PFFT-BEM and 3D PE, we briefly describe the formulation of the PE method following previous studies (e.g., [14,19]). We begin with the Helmholtz equation in the frequency domain to determine the acoustic wave pressure p(r) in a cylindrical coordinate system:

Model Comparisons of the PFFT-BEM with the 3D Parabolic Equation Method
To demonstrate the usefulness of PFFT-BEM, PFFT-BEM is compared with 3D Parabolic equation (PE) models to check the correctness and performance of different PE schemes. We first briefly review the formulation and assumptions of the 3D PE theory. Then, we present three numerical examples for comparison: ASA wedge, underwater seamount, and Gaussian canyon.

Formulation and Assumptions of the PE Method
In order to demonstrate the comparison between PFFT-BEM and 3D PE, we briefly describe the formulation of the PE method following previous studies (e.g., [14,19]). We begin with the Helmholtz equation in the frequency domain to determine the acoustic wave pressure p(r) in a cylindrical coordinate system: where k 0 = ω/c 0 is the reference wavenumber, c 0 is the reference sound speed (defined based on the sound profiles [14]), ω is the sound frequency, and n = c 0 /c = k/ k 0 is the index of refraction. The derivation of PE methods in a Cartesian coordinate system is similar to that one in a cylindrical coordinate system and is therefore omitted here. By applying the separation of variables and accounting for only outgoing propagation sound waves, we can approximate p(r) by where H 0 is the zeroth-order Hankel function of the first kind. If we assume that there is only outgoing propagation in the range, the Helmholtz equation reduces to an initial boundary value problem governing by a parabolic equation as with ϕ(r = r 0 , θ, z) = ϕ 0 (θ, z), where, ϕ 0 denotes the initial outgoing field at r = r 0 , I is the identity operator, X is the 2D depth operator in the r-z plane, and Y is the azimuthal operator, defined as Solving Equation (29) usually has a large computational cost. Several methods approximate the nonlinear term (square root operator in Equation (29)) to reduce the computational cost. The 3D PE proposed by Tappert [5] expands the square root operator and only retain the linear terms as It is known that this formulation allows only for narrow angle propagation along θ. For wide angle approximation, Sturm [20] approximates the square root operator using the Padé approximation along z and a linear approximation long θ as √ I + X + Y = I + n p k=1 a k,n p X where n p is the number of terms, and a k,n p , b k,n p , 1 ≤ k ≤ n p are the real coefficients given by a k,n p = 2 2n p + 1 Another approximation form, which includes part of the nonlinear term, was proposed by Feit and Fleck [21] in a cylindrical coordinate 3D PE method as The numerical results calculated based on these approximations usually have a phase error in the far field [20]. These errors are quantified in the present study by comparison between direct simulation results from PFFT-BEM and results from 3D PE.
Since PE only has one range derivative in the marching direction, it is not possible to enforce both boundary conditions across the vertical boundaries of the marching directions [14]. For most 3D PE methods, only continuity of pressure is enforced across these vertical boundaries. To conserve energy at each marching step, Collins [22] introduced a variable p = p/ √ ρc, where c is the sound speed and ρ is the medium density. Then, X and Y shown in Equation (25) can be expressed as However, as indicated by Collins [22], this method tends to overcorrect the solution, and for problems involving a large slope ocean bottom (> 10 • ), the energy conservation PE error becomes significant. In summary, three main assumptions are made during the formulation of the 3D PE method: • Higher-order term effects: the higher-order terms of the square root operator are small and negligible; • Backscattering effects: only forward scattering is included (no backscattering effects); • Large slope effects: only pressure continuity is satisfied across the marching vertical boundaries.
As a result, we further demonstrate the differences between 3D direct numerical simulation results by the PFFT-BEM method and 3D PE results due to the assumptions used in 3D PE. To achieve this, we examine three classical shallow water examples: the ASA wedge, 3D underwater seamount, and a Gaussian canyon problem.

Higher-Order Term Effects: Acoustical Society of America (ASA) Wedge
The first benchmark example is the ASA wedge, as shown in Figure 11. The sound source is fixed at (0, 0, 100 m) with a frequency at 75 Hz. The wedge meets free surface (z = 0) at y = -2000 m. The slope angle of the wedge is π/36. The computational domain size covers the range of 28 km × 4 km with a maximum depth at 400 m. The medium properties are defined in Figure 11, with sound speed c 1 = 1500 m/s and density ρ 1 = 1 g/cm 3 in the water and c 2 = 1700 m/s, density ρ 2 = 1.5 g/cm 3 and medium attenuation α = 0.5 dB/λ in the bottom. The boundary element model uses quadrilateral element with ∆l = λ/12. We compared the 3D PFFT-BEM with 3D cylindrical coordinate PE results based on the nonlinear operator Equation (32). The far-field solution of 3D PFFT-BEM and a 3D PE method are compared with benchmark solutions. From Figure 12a, quantitative validation of large-scale 3D sound field prediction by direct numerical computation is obtained. In contrast, phase error is observed in the far field of 3D PE results. The approximation of the square root operator in Equation (21) is the main reason for this phase error. Figure 12b shows the phase differences between PFFT-BEM and the 3D PE method. As expected, the phase errors of 3D PE increase as a function of distance. It should be noted that higher-order terms can be included to make the 3D PE solutions more accurate [23,24]. However, the small marching steps of this higher-order 3D PE method will introduce large computational costs, especially in the cylindrical coordinates. We show the 3D feature of the TL contour for the ASA wedge in Figure 13.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 21 order 3D PE method will introduce large computational costs, especially in the cylindrical coordinates. We show the 3D feature of the TL contour for the ASA wedge in Figure 13.   Appl. Sci. 2020, 10, x FOR PEER REVIEW 15 of 21 order 3D PE method will introduce large computational costs, especially in the cylindrical coordinates. We show the 3D feature of the TL contour for the ASA wedge in Figure 13.   order 3D PE method will introduce large computational costs, especially in the cylindrical coordinates. We show the 3D feature of the TL contour for the ASA wedge in Figure 13.

Backscattering Effects: 3D Seamount
In the second numerical example, we show a comparison of the 3D PFFT-BEM and a 3D PE method for an idealized seamount. The water column is homogeneous with a sound speed of 1500 m/s, a density of 1.0 g/cm 3 , and no medium absorption for sound waves. The bottom is also homogeneous with a sound speed of 1700 m/s, a density of 1.5 g/cm 3 , and medium attenuation at 0.5 dB/λ. The origin of the coordinate system is located above the seamount at the sea surface with z pointing down towards the ocean. The 100 Hz sound source in this example is placed away from the tip of the seamount at x = −2000 m, y = 850 m, and z = 250 m. The summit of the seamount is at a depth of 200 m, and its slope is 7 • . The radius of the seamount is 2500 m, and the water depth outside the seamount area is 512.5 m. Figure 14a compares the 3D PFFT-BEM numerical results of TL to the 3D PE results by [19] along the x−axis at z = 250 m in front of the seamount (2 km away from the source). First, as shown in Figure 14a, convergence is observed with an increase in the element numbers, which demonstrates the accuracy of PFFT-BEM. The difference between 3D PE and PFFT-BEM is more significant in front of the seamount. This difference is mainly due to the backscattering effects, which are not included in the 3D PE. We zoom in on the range from 750 to 1100 m, which is in front of the seamount. As shown in Figure 14b, significant differences between the PFFT-BEM and 3D PE results (3 ∼ 5 dB) are observed. In Figure 14b, we compare the TL along the x-axis for two cases with different slopes (1/6 and 1/8). As expected, the TL is smaller (i.e., larger acoustic pressure) in the larger slope case (1/6), which indicates stronger backscattering effects. Under such conditions, we expect that the errors due to neglecting the backscattering effects are more significant, and the use of a 3D full field method such as PFFT-BEM is more necessary. To demonstrate the 3D scattering effects, we plot Figure 15 to show the TL contour computed by 3D PFFT-BEM.

Backscattering Effects: 3D Seamount
In the second numerical example, we show a comparison of the 3D PFFT-BEM and a 3D PE method for an idealized seamount. The water column is homogeneous with a sound speed of 1500 m/s, a density of 1.0 g/cm 3 , and no medium absorption for sound waves. The bottom is also homogeneous with a sound speed of 1700 m/s, a density of 1.5 g/cm 3 , and medium attenuation at 0.5 dB/ . The origin of the coordinate system is located above the seamount at the sea surface with z pointing down towards the ocean. The 100 Hz sound source in this example is placed away from the tip of the seamount at x = -2000 m, y = 850 m, and z = 250 m. The summit of the seamount is at a depth of 200 m, and its slope is 7 ∘ . The radius of the seamount is 2500 m, and the water depth outside the seamount area is 512.5 m. Figure 14a compares the 3D PFFT-BEM numerical results of TL to the 3D PE results by [19] along the −axis at = 250 in front of the seamount (2 km away from the source). First, as shown in Figure 14a, convergence is observed with an increase in the element numbers, which demonstrates the accuracy of PFFT-BEM. The difference between 3D PE and PFFT-BEM is more significant in front of the seamount. This difference is mainly due to the backscattering effects, which are not included in the 3D PE. We zoom in on the range from 750 to 1100 m, which is in front of the seamount. As shown in Figure 14b, significant differences between the PFFT-BEM and 3D PE results (3 ∼ 5 dB) are observed. In Figure 14b, we compare the TL along the x-axis for two cases with different slopes (1/6 and 1/8). As expected, the TL is smaller (i.e., larger acoustic pressure) in the larger slope case (1/6), which indicates stronger backscattering effects. Under such conditions, we expect that the errors due to neglecting the backscattering effects are more significant, and the use of a 3D full field method such as PFFT-BEM is more necessary. To demonstrate the 3D scattering effects, we plot Figure 15 to show the TL contour computed by 3D PFFT-BEM.  For problems involving an ocean bottom with a very large slope (> 10°), the validity of 3D PE needs to be examined by a direct numerical method. In the next benchmark problem, we consider a waveguide problem with a Gaussian canyon. The water column has a depth of 700 m with c1 = 1500

Large Slope Effects: Gaussian Canyon
For problems involving an ocean bottom with a very large slope (> 10 • ), the validity of 3D PE needs to be examined by a direct numerical method. In the next benchmark problem, we consider a waveguide problem with a Gaussian canyon. The water column has a depth of 700 m with c 1 = 1500 m/s and ρ 1 = 1.0 g/cm 3 . The sound source has frequency of 25 Hz and is fixed at (0, 0, 30 m). The bottom has the properties of ρ 2 = 1.5 g/cm 3 , c 2 = 1700 m/s and medium attenuation α = 0.5 dB/λ. The geometry of the Gaussian canyon bottom is defined as where the parameters are h 0 = 20 m, h 1 = 180 m, and σ = 720 m. The largest slope in the direction across the canyon is approximately 0.21 (angle of 12.1 • ). The geometry of the canyon is plotted in Figure 16.
We compare the 3D PFFT-BEM solution along the canyon with the 3D cylindrical coordinate PE results from [20]. As shown in Figure 17a, numerical convergence of PFFT-BEM is observed with increasing element numbers. This proves the accuracy of the PFFT-BEM. In Figure 17a, differences (∼8 dB) between the 3D PE and 3D PFFT-BEM can be observed in energy-focusing regions along the canyon axis where horizontal refractions by the sidewalls of the canyon are the strongest (shown in Figure 17b). These differences are mainly due to the fact that, for 3D PE, only one of the two boundary continuity conditions can be satisfied in the marching direction. As a result, for environment geometries with large slopes in the marching directions, such as this Gaussian canyon, the 3D PE method becomes inaccurate and deviates from the direct numerical solutions by PFFT-BEM. We also show the contours of the transmission loss in Figure 17b where strong 3D scattering effects from the canyon are observed.

Prediction of Sound Scattering by Traveling Internal Waves
Tidal interactions with bottom features often generate internal solitary waves (ISWs) in a shallow water environment. Previous experiments [26][27][28] found that acoustic signals fluctuate when sound waves are transmitted through ISWs. The amplitudes and phases of these fluctuations are

Prediction of Sound Scattering by Traveling Internal Waves
Tidal interactions with bottom features often generate internal solitary waves (ISWs) in a shallow water environment. Previous experiments [26][27][28] found that acoustic signals fluctuate when sound waves are transmitted through ISWs. The amplitudes and phases of these fluctuations are highly correlated with those of ISWs. Quantifying the interactions between the acoustic waves and ISWs is important for shallow acoustics. Here, we compare our numerical simulation of sound

Prediction of Sound Scattering by Traveling Internal Waves
Tidal interactions with bottom features often generate internal solitary waves (ISWs) in a shallow water environment. Previous experiments [26][27][28] found that acoustic signals fluctuate when sound waves are transmitted through ISWs. The amplitudes and phases of these fluctuations are highly correlated with those of ISWs. Quantifying the interactions between the acoustic waves and ISWs is important for shallow acoustics. Here, we compare our numerical simulation of sound propagate through ISW with the experiment data from SWRAM'95 (shallow-water acoustics in a random medium) which is conducted in the New Jersey continental shelf in 1995. The experiment setup and data are presented in detail by previous studies [26][27][28][29][30]. The purpose of including this example is to demonstrate the capability of using large-scale direct simulations in the prediction of acoustic scattering by internal waves, which is critically needed in shallow water acoustics.
We plot the positions of the source, receiver and ISWs fronts during experiment in Figure 18a. During the interval 18:00 to 20:00 GMT on 4 August 1995, the source is fixed at z s = 12 m below the sea surface at the position approximately 15 km from the receiver, as shown in Figure 18b [28]. The receiver is placed at z R = 45 m below the sea surface. The broadband source frequency f was 30~160 Hz with three peak frequencies at 30, 60, and 90 Hz.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 19 of 21 Using the PFFT-BEM, we simulate the propagation of the acoustic waves through the internal waves. As shown in Figure 18b, a soliton wave packet model based on Duda and Preisig [31] is used in the 3D PFFT-BEM direct numerical simulation. The properties of the water column and sea bottom are shown in Figure 18b. To consider the change in the sound speed, we use four layers in the gradient region where the sound speed changes from 1530 to 1480 m/s. The properties of the bottom are cb = 1750 m/s, = 1.8 g/cm 3 with an attenuation parameter of = 0.4375 dB/λ. The velocity of the ISWs (in the y direction) is 0.65 m/s. Here, we ignore the differences between the speeds of separate solitary waves. As a result, the entire internal wave train is assumed to move at the same speed. The sound source frequencies used in the 3D PFFT-BEM are 30, 60 and 90 Hz. To demonstrate the accuracy of PFFT-BEM in simulating the acoustical scattering effects by the internal waves, we compare the acoustic pressure calculated using the 3D PFFT-BEM at the receiver with the field measurement during a 30 min period [28]. As shown in Figure 19, the numerical prediction by the PFFT-BEM agrees fairly well with the field experimental data for the time variation trend and the magnitude of the fluctuation of the received sound pressure. The difference is mainly because the internal waves are irregular during the experiment, while the regular wave model is used in our PFFT-BEM numerical model. The root mean square between the numerical prediction and experimental data is 3 %. For underwater sound measurements, the uncertainty of the receiver system in field conditions is approximately 0.5 dB [32]. As a result, a relative error of 3 % (i.e., 0.13 dB) would be sufficient for the underwater acoustic numerical simulation. Using the PFFT-BEM, we simulate the propagation of the acoustic waves through the internal waves. As shown in Figure 18b, a soliton wave packet model based on Duda and Preisig [31] is used in the 3D PFFT-BEM direct numerical simulation. The properties of the water column and sea bottom are shown in Figure 18b. To consider the change in the sound speed, we use four layers in the gradient region where the sound speed changes from 1530 to 1480 m/s. The properties of the bottom are c b = 1750 m/s, ρ b = 1.8 g/cm 3 with an attenuation parameter of α = 0.4375 dB/λ. The velocity of the ISWs (in the y direction) is 0.65 m/s. Here, we ignore the differences between the speeds of separate solitary waves. As a result, the entire internal wave train is assumed to move at the same speed. The sound source frequencies used in the 3D PFFT-BEM are 30, 60 and 90 Hz.
To demonstrate the accuracy of PFFT-BEM in simulating the acoustical scattering effects by the internal waves, we compare the acoustic pressure calculated using the 3D PFFT-BEM at the receiver with the field measurement during a 30 min period [28]. As shown in Figure 19, the numerical prediction by the PFFT-BEM agrees fairly well with the field experimental data for the time variation trend and the magnitude of the fluctuation of the received sound pressure. The difference is mainly because the internal waves are irregular during the experiment, while the regular wave model is used in our PFFT-BEM numerical model. The root mean square between the numerical prediction and experimental data is 3 %. For underwater sound measurements, the uncertainty of the receiver system in field conditions is approximately 0.5 dB [32]. As a result, a relative error of 3 % (i.e., 0.13 dB) would be sufficient for the underwater acoustic numerical simulation.
internal waves, we compare the acoustic pressure calculated using the 3D PFFT-BEM at the receiver with the field measurement during a 30 min period [28]. As shown in Figure 19, the numerical prediction by the PFFT-BEM agrees fairly well with the field experimental data for the time variation trend and the magnitude of the fluctuation of the received sound pressure. The difference is mainly because the internal waves are irregular during the experiment, while the regular wave model is used in our PFFT-BEM numerical model. The root mean square between the numerical prediction and experimental data is 3 %. For underwater sound measurements, the uncertainty of the receiver system in field conditions is approximately 0.5 dB [32]. As a result, a relative error of 3 % (i.e., 0.13 dB) would be sufficient for the underwater acoustic numerical simulation. Figure 19. Comparison of field normalized acoustic pressure data and normalized acoustic pressure predicted by the PFFT-BEM. Figure 19. Comparison of field normalized acoustic pressure data and normalized acoustic pressure predicted by the PFFT-BEM.

Conclusions
We have presented an efficient boundary element method with a multi-layer formulation and first-order non-reflecting boundary condition for the simulation of sound propagation in a shallow water environment. By using the pre-corrected fast Fourier transform method, the developed method, PFFT-BEM, which reduces the computational cost from O(N 2 ) to O(Nlog N), can be applied to the 3D full-field direct numerical simulation of sound propagation in environments with a sound frequency of up to O(100) Hz and a domain size of up to O(10) km. We found that the developed multi-layer numerical method is capable to simulate the sound propagation in inhomogeneous waveguide accurately. We also concluded that the first-order NRBC used in PFFT-BEM can efficiently eliminate the errors due to reflection from artificial numerical boundary. We demonstrated the usefulness of the developed PFFT-BEM for shallow water acoustic applications. First, the 3D direct numerical simulation capability of PFFT-BEM enable us to check the validity of classical approximation models. We show the evaluation of the 3D PE models for three examples: the ASA wedge, the 3D seamount, and the Gaussian canyon. Through these comparisons, we demonstrated and quantified the differences between 3D direct numerical simulation results by PFFT-BEM and 3D PE results due to the main assumptions used in the 3D PE. To show the usefulness of the PFFT-BEM in real applications, we compared the numerical results obtained by the 3D+T PFFT-BEM and the field measurement data from the SWARM' 95 internal wave experiment. The transmission loss obtained with the PFFT-BEM showed strong variability as the internal wave train moved. The magnitude and trend of these numerical variations agree fairly well with the field experiments. This indicates that the 3D+T PFFT-BEM is applicable for the direct simulation of sound propagation in a complex 3D shallow water environment. We concluded that the developed direct numerical method, PFFT-BEM, is very useful to quantify the accuracy of exiting numerical methods and simulate sound propagation in realistic multi-layer shallow water environments. Funding: This research was partially supported by the Office of Naval Research (N00014-16-1-2200) and was granted access to the HPC re-sources provided by the DOD HPC Modernization Program.