A Finite Element Model for Underwater Sound Propagation in 2-D Environment

: The ﬁnite element method is a popular numerical method in engineering applications. However, there is not enough research about the ﬁnite element method in underwater sound propagation. The ﬁnite element method can achieve high accuracy and great universality. We aim to develop a three-dimensional ﬁnite element model focusing on underwater sound propagation. As the foundation of this research, we put forward a ﬁnite element model in the Cartesian coordinate system for a sound ﬁeld in a two-dimensional environment. We ﬁrstly introduce the details of the implementation of the ﬁnite element model, as well as different methods to deal with boundary conditions and a comparison of these methods. Then, we use four-node quadrilateral elements to discretize the physical domain, and apply the perfectly matched layer approach to deal with the inﬁnite region. After that, we apply the model to underwater sound propagation problems including the wedge-shaped waveguide benchmark problem and the problem where the bathymetry consists of a sloping region and a ﬂat region. The results by the presented ﬁnite element model are in excellent agreement with analytical and benchmark numerical solutions, implying that the presented ﬁnite element model is able to solve complex two-dimensional underwater sound propagation problems accurately. In the end, we compare the ﬁnite element model with the popular normal mode model KRAKEN by calculating sound ﬁelds in Pekeris waveguides, and ﬁnd that the ﬁnite element model has better universality than KRAKEN.


Introduction
There are a variety of methods for calculating the underwater sound fields, such as the ray method, the normal mode method, and the parabolic equation method.However, the finite element method is used by only a few researchers to solve underwater sound propagation problems.
As we all know, the finite element method is a universal numerical calculation method in many engineering fields.The basic principle of the finite element method is to discretize the continuous area into limited number of elements and nodes.We can use elements of different shapes to deal with the boundaries of complex shapes.
The finite element method was first put forward by Clough [1] in 1960 for a plane elasticity problem.Finite element methods for time-harmonic acoustics governed by the Helmholtz equation have been an active research area over the past 50 years.Initial applications of finite element methods for time-harmonic acoustics focused on interior problems with complex geometries including direct and modal coupling of structural acoustic systems for forced vibration analysis, frequency response of acoustic enclosures, and waveguides [2][3][4].In recent years, tremendous progress has been made; for example, extending the finite element method to exterior problems in unbounded domains and higher frequency regimes, which incorporates knowledge of wave behavior into the algorithm, combined with parallel sparse interactive and domain decomposition solvers.
The exterior acoustics problem in unbounded domains presents a special challenge for finite element methods.In order to apply the FEM to exterior problems, the unbounded domain is usually truncated by an artificial boundary, yielding a bounded computational domain.The first complete finite element approach for modeling acoustic radiation and scattering in unbounded domains appears in the impedance matching technique presented by Hunt [5,6].Recent numerical treatments including infinite elements, absorbing layers, local absorbing boundary conditions, and exact nonlocal boundary conditions have proven to be effective in handling acoustic scattering problems in unbounded domains, especially for large-scale problems requiring interactive and parallel solution methods and for modeling inhomogeneities and acoustic-structure interaction.The choice of methods depends on the shapes and complexity of the scattering objects, inhomogeneities, frequency range, and resolution requirements, among other parameters.
The infinite elements are constructed with radial wave functions that automatically satisfy the radiation boundary condition.Bettess [7] pioneered the infinite element concept and selected the test function to be the same as the trial solution basis.However, the performance of both conjugated and unconjugated formulations deteriorates at larger wavenumbers and highly elongated artificial interfaces [8].The perfectly matched layer (PML) concept was originally introduced by Berenger [9,10] for electromagnetic waves.It is another option for modeling the far field for the exterior acoustics problem.The idea is to introduce an exterior layer of finite thickness at an artificial interface such that outgoing waves are absorbed prior to reaching the outer-layer truncation boundary.By splitting the scalar field into nonphysical components, satisfying equations that describe decaying waves and proper selection of PML coefficients, plane-wave reflection for an arbitrary angle of incidence is theoretically eliminated.There is also a great deal of research on absorbing boundary conditions [11] and the DtN nonreflecting boundary conditions [12,13].
The finite element method is able to solve problems in inhomogeneous media and allows for a natural coupling with complex structures.For exterior problems in unbounded domains, special techniques are required to reduce spurious reflection to a level below that of the discretization error.The numerical advantage of the FEM is that it leads to sparse matrices, which significantly speeds up computations and reduces memory requirements.Complexity estimates [14,15] show that domain-based methods such as the FEM are an effective alternative to the BEM for exterior acoustics problems, especially for large systems due to the sparse structure of the resulting system matrices.With the recent development in fast multiple methods, which accelerates the calculation of matrix-vector products in iterative integral methods [16], the method with the best efficiency is less clear, yet the FEM retains the advantages of robustness and natural integration with other discrete models in coupled problems.It is also possible to couple the finite method with boundary integral methods [17] and other domain-based methods, such as global Trefftz-based wave methods, which are effective at high frequencies on moderate geometrical complexity [18].
However, there is not such large amount of research on finite element methods in ocean acoustics.In the 1990s, Murphy and Chin-Bing [19][20][21] established a finite element model FOAM for shear waves in marine sediment.They put forward a wide-angle radiation condition, which is actually a kind of parabolic equation, to simulate the infinite area.In 1999, Kampanis and Dougalis [22] established a finite element model FENL in axially symmetric waveguides.They used the Dirichlet-to-Neumann boundary condition to simulate the infinite area.In 2006, Li Jun [23] established a finite element model FEMLJ for prediction of sound fields in shallow water, which used the same methods as those in FENL.In 2008, Athanassoulis [24] compared the finite element model FENL with coupled normal mode model CCMM.Then, in the next 10 years, more and more researchers began to focus on commercial finite element software for oceanographic phenomena and scattering problems [25,26].In 2013, Goldsberry and Isakson [27] mentioned in their paper that performing a fully three-dimensional finite element model is currently unfeasible due to the difficulties in implementation and limits in computational power.They compared two promising methods to represent the 3D field, the longitudinally invariant (LI) method and a 2D axial-symmetric reduction of the 3D Helmholtz equation.They used these two methods for ASA benchmark problem, but failed to agree with the PE solution.From 2016 to 2020, Chai Yingbing [28][29][30] focused on solving the numerical dispersion problem for the high-frequency problem by the traditional finite element method and put forward several methods to soften the stiffness matrix.
As can be seen, there is not enough research on finite element modeling in underwater sound propagation.Most of the existing finite element models rely on commercial software or focus on reverberation problems rather than propagation problems.In this paper, we will propose a finite element model for underwater sound propagation in a two-dimensional environment.We will also introduce the implementation in detail and show the results of this model.

Practical Aspects of the Implementation
The basic idea of the finite element method is transforming an infinite dimensional problem into a finite dimensional problem, through discretizing the physical domain into several elements with finite number of nodes.
In this section, the specific steps of our finite element model will be described.Firstly, the weak form of the Helmholtz equation is derived through the traditional Galerkin method.Then, four-node quadrilateral elements are used to discretize the physical domain.We assemble all the local stiffness matrices into a global stiffness matrix.In the end, the treatment of boundaries is discussed.

Weak Form of Helmholtz Equation
This finite element model mainly focuses on calculation of the two-dimensional sound field in the Cartesian coordinate system at a specific frequency.So, the following discussion is based on the Helmholtz equation.
The Helmholtz equation [31] in an inhomogeneous medium without sources is where Ω is a bounded domain in the (x, z) plane with the boundary ∂Ω; p is a function of (x, z) indicating the sound pressure; and ρ and k are functions of (x, z) indicating density and wave number, respectively.Multiply both sides of Equation (1) with a test function q(x, z) and obtain Then, integrate the equation over Ω, Finally, apply the Gauss theorem to obtain the weak formation of the Helmholtz equation, where the first term indicates the integral on the boundary of Ω, the second and third terms indicate the integral on Ω.We will discuss the discretization of Ω and ∂Ω separately.

The Finite Element Discretization
In order to calculate the last two terms of Equation ( 4), we use four-node quadrilateral elements to discretize Ω. Figure 1 is a schematic diagram of a four-node quadrilateral element and its local indices, midpoint coordinates (x 0 , z 0 ), length 2a, and width 2b.We define the element as Ω e .We use N n (x, z) for the shape functions or the basis functions, where n = 1, 2, 3, 4 indicates the local index of four nodes.N n (x, z) can be written as the following form using the Lagrange interpolation method, For the area outside Ω e , the values of shape functions are zero.

Calculation of Local Stiffness Matrices
We have already known the shape functions of each element.The sound pressure in the element is only related to the sound pressure of the four nodes and the shape functions, according to the features of these shape functions.As a result, we define p n as the sound pressure of each node, where n = 1, 2, 3, 4. So, the sound pressure in the element can be expressed as Then, we choose to set the test functions equal to the shape functions, We also set the values of test functions to zero on ∂Ω.We will consider boundary conditions later.
Supposing that ρ(x, y) is a constant in each single element and substituting Equations ( 6) and (7) into Equation (4), the finite element equations of each element without considering boundary conditions can be written as From Equation (5), we can rewrite Equation (8) as x 0 +a We define k mn as the coefficients of p n .They make up the elements of local stiffness matrix, where m, n = 1, 2, 3, 4 indicate the local indices of nodes.We can regard wave number k(x, z) = k e as a constant in one element.Then, k mn can be written as

Assembly of Global Stiffness Matrix
We have already obtained the local stiffness matrices using Equation (10).Then, we can assemble the global stiffness matrix with the local stiffness matrices.
There are four nodes in every element, and every node has a local index and a global index.We need to add the elements of local stiffness matrices to the corresponding elements of global stiffness matrix.For example, if the four local indices of nodes in an element are 1, 2, 3, and 4, the global indices of them are 20, 21, 45, and 46, respectively.If we define the global stiffness matrix as K, we can assemble the global stiffness matrix of this element by adding k 1,1 , k 1,2 , . . ., k 4,4 to K 20,20 , K 20,21 , . . ., K 46,46 .

Boundary Conditions
The boundary conditions can be calculated using the first term in Equation (4) as Below, we discuss different methods to deal with different kinds of boundary conditions.

Dirichlet Boundary
In underwater acoustics, the sea surface is always regarded as a Dirichlet boundary, which satisfies p = 0.
Besides the sea surface, we define the Dirichlet boundary as where l is the global indices of node (x l , z l ); then, we can treat the boundary conditions with the following steps: 1. Define the global stiffness matrix as K, the global load vector as f, and the sound pressure vector as p.They satisfy the linear equation system Kp = f.2. Subtract the product of the lth column of K and the lth element of p from f. 3. Set the elements of the lth row and the lth column of K to 0 and the lth diagonal element to 1. 4. Set the lth element of f to 0.
We must consider that the treatment of the Dirichlet boundary condition is always the last step when processing the global stiffness matrix.

Neumann Boundary
The most common Neumann boundary condition in underwater sound propagation is the totally reflected rigid bottom, which satisfies ∂p ∂ n = 0. Besides, if the environment is symmetric at a certain range, the axis of symmetry can be treated as the same Neumann boundary, too.
We define this Neumann boundary as Γ 2 ⊂ ∂Ω.When (x, z) ∈ Γ 2 , the first term of Equation ( 4) is This means that in this finite element model, there is no need for special treatment for this kind of boundary condition.

Free-Space Simulation
The finite element method can only be used in a bounded domain.Nevertheless, open problems involving theoretically boundless space extension can be solved when applying special conditions on the boundaries of the computational domain, in order to absorb the outgoing waves.This kind of outgoing radiation condition can be described as generalized Sommerfeld-Rellich radiation condition [32] to ensure the uniqueness and existence of solutions for the Helmholtz equation.
We will apply two methods to absorb the outgoing waves.One method is the traditional radiation boundary condition and the other is the perfectly matched layer approach.
The traditional radiation condition means that the normal derivative of sound pressure is equal to the product of sound pressure and a certain constant.It can be expressed as For example, we write the open boundary as Γ 3 = {x = X, z ∈ (0, D)}.The area at x > X denotes the infinite sea water.If the right side of the physical domain is an open boundary, where (x, z) ∈ Γ 3 , the first term of Equation ( 4) can be written as We substitute Equations ( 6) and ( 7) into Equation ( 12) so that we can obtain the local stiffness matrices of the elements on the open boundary, Finally, we obtain the global stiffness matrix by adding corresponding elements of Equation ( 13) to global stiffness matrix.
Perfectly matched layer (PML) approach is another method to deal with the open boundary.It was first put forward by Berenger in 1994 [9,10].The idea is to introduce an exterior layer of finite thickness at an artificial interface such that outgoing waves are absorbed prior to reaching the outer-layer truncation boundary.By splitting the scalar field into nonphysical components, satisfying equations that describe decaying waves and proper selection of PML coefficients, plane-wave reflection for an arbitrary angle of incidence is theoretically eliminated.In order to calculate the local stiffness matrices of the elements in PML, we make a coordinate transformation from (x, z) to ( x, z) = T(x, z).The sound pressure p in Ω PML satisfies where Ω PML = {T(x, z) : (x, z) ∈ Ω PML } is the mapping of physical domain after coordinate transformation and ∇ is the gradient operator after coordinate transformation.
The boundaries of PML can be selected arbitrarily without influencing the results.So, we set the boundary condition of PML as where ∂ Ω PML is the boundary of Ω PML after coordinate transformation.The weak form of the Helmholtz equation in Ω PML is According to the chain rule of derivation local stiffness matrices of elements in PML can be calculated by substituting Equations ( 6) and ( 7) into Equation ( 16).
Taking the previous problem as an example, we set an extra layer as perfectly matched layer , where δ x is the thickness of PML.We define the form of T(x) as where σ xm is a constant that decides the maximum loss in PML in the x-direction.From Equations ( 18) and ( 19), we can obtain k mn is the (m, n)th element of the local stiffness matrix in Ω PML .Then, we need to add k mn to the corresponding element of the global stiffness matrix K, just as above, so that we can obtain the global stiffness matrix.

Numerical Results
With the previous theoretical analysis, we are now able to build a finite element model for underwater sound propagation.The finite element model can be used to calculate sound fields in a complex 2-D environment.
In this section, some of the results will be displayed.Firstly, the finite element model is applied to a range-dependent benchmark problem of a totally reflected wedge-shaped waveguide.Then, it is used to calculate range-dependent problems with either a totally reflected bottom or a penetrable bottom, where the bathymetry consists of a sloping region and a flat region.

Range-Dependent Benchmark Problem
We consider the range-dependent benchmark problem first.In Figure 2, there is a wedge-shaped waveguide bounded above by a planar pressure-release sea surface and below by a planar rigid bottom.A line source is involved, parallel to the apex and located at (r s , θ s ).The wedge angle is denoted by θ 0 .The analytical solution of sound field in the waveguide is [33] where (r, θ) is the polar coordinate representation of the observation point, r s is the distance from the origin to the source, θ s is the angle between the x-axis and the line connecting the origin and the source, and θ 0 is the inclination of the wedge-shaped waveguide.Then, we will apply the presented finite element model to calculate the sound field in the wedge-shaped waveguide and compare its results with those by the above analytical solution.Figure 3    Then, we apply the presented finite element model to solve this problem.We discretize the wedge-shaped waveguide according to Figure 5 using four-node quadrilateral elements.If we use the traditional radiation boundary condition to simulate the boundless space, the elements from X to X + δ x are needless and only the physical domain of sea water is discretized and calculated.The size of the elements is 0.5 m × 0.5 m.Compared with Figure 4, obvious error caused by interference can be seen in Figure 6a.This is because the traditional radiation boundary condition cannot totally absorb the outgoing wave-that is, the reflection coefficient at the boundary is not zero.However, the result in Figure 6b is consistent with that in Figure 4.
Figure 7 is the transmission loss in the wedge-shaped waveguide at 250 Hz at 10 m receiving depth.The blue solid line indicates the analytical solution, the red dashed line indicates finite element solution with traditional radiation boundary condition, and the green dashed line indicates finite element solution with PML.It can be seen that the green dashed line matches the blue solid line better, indicating that PML is more appropriate for simulating the boundless area than the traditional radiation boundary condition, and hence, yields higher precision.

Slope Problem with Totally Reflected Bottom
Below, we consider the sound field of the wedge-shaped waveguide with a section of flat region, also called a slope problem. Figure 8 shows the schematic diagram of a slope-shaped waveguide.When 0 < x < X, there is a wedge-shaped waveguide.When x > X, there is a section of flat region.A line source is involved, which is parallel to the apex and located at (x s , z s ) in the wedge-shaped part.Since there is no analytical solution for this problem, the reference solution is provided by the coupled normal mode model DGMCM2D [34].Then, we use the finite element method to calculate the sound field in specific slopeshaped waveguide environment and compare it with the reference solutions.Figure 9    Then, we will use the finite element method to solve this problem.We discretize the slope-shaped waveguide according to Figure 11 using four-node quadrilateral elements.Only the physical domain of sea water is discretized and calculated.The area from X 2 to X 2 + δ x denotes the PML.The size of the elements is 0.5 m × 0.5 m.
Figure 12 is the transmission loss calculated by the finite element method with PML, where σ xm = 50, δ x = 12 m and the PML is divided into 15 layers vertically.It takes MATLAB almost 100 s to calculate on an i7-4790@3.60GHz CPU.

Slope Problem with Penetrable Bottom
As shown in the above sections, we have already seen the consistency between the finite element model and reference solution in a waveguide with a totally reflected bottom.Now, we will discuss the sound field in a waveguide with a penetrable bottom, as a slope bounded above by the sea surface and below by a penetrable bottom can be regarded as the simplest model of a shallow-sea waveguide near the coastline.
Figure 14 shows the schematic diagram of a slope with a penetrable bottom.X 1 < x < X 2 denotes a section of flat region of sea water.A line source is involved, which is parallel to the apex and located at (x s , z s ).For an environment with a penetrable bottom such as this problem, the physical domain is infinite not only in the horizontal direction but also in the vertical direction.As a result, we have to set a PML in the x-direction, a PML in the z-direction, and a PML in both the xand z-directions, as Figure 15 shows.In Ω PMLx = {x ∈ (X 2 , X 2 + δ x ), z ∈ (0, D 2 )}, the coordinate transformation is only set on the x-axis, according to Equations ( 18)- (22), where X = X 2 .Further, the coordinate transformation is only set on the z-axis in where σ zm is a constant that decides the maximum loss in PML in the z-direction.From Equations ( 24) and ( 25), we can obtain Here, k mn is the (m, n)th element of the local stiffness matrix of elements in Ω PMLz .As for , coordinate transformation needs to be set on both the xand z-axis-that is, both σ xm and σ zm are not zero in this area.We substitute Equations ( 19) and ( 25) into Equation ( 16) to obtain the elements of local stiffness matrices Figure 15.Schematic diagram of a slope with penetrable bottom and perfectly matched layers.
Then, we will use the finite element method to calculate the sound field in a specific slope-shaped waveguide environment with a penetrable bottom and compare it with the reference solution.Since there is no analytical solution for this problem, the reference solution is provided by DGMCM2D.It should be mentioned that for the wedge-shaped penetrable bottom problem, it is more appropriate to adopt the source images solution [35,36] as the reference solution.Figure 16 is the schematic diagram of specific environmental parameters.The sound speed and density of sea water are c 1 = 1500 m/s and ρ 1 = 1000 kg/m 3 ; the sound speed and density of sediment are c 2 = 2400 m/s and ρ 2 = 2200 kg/m 3 ; and the depth and range of source are z s = 20 m and x s = 50 m, respectively.We will consider sound fields at the frequencies of 30 Hz and 50 Hz.For the case of the source frequency of 30 Hz, Figure 17a is the reference solution of transmission loss calculated by DGMCM2D and Figure 17b is calculated by the finite element model.It can be seen that the two sets of results are in excellent agreement.For the case of the source frequency of 50 Hz, Figure 18a is the reference solution of transmission loss calculated by DGMCM2D and Figure 18b is calculated by finite element model.It is obvious that the two sets of results are also very consistent at 50 Hz.

Analysis of Universality
The accuracy of the presented finite element model has been demonstrated previously.Here, we discuss the universality of this finite element model compared with the normal mode model KRAKEN.
Consider the problem of sound propagation excited by a point source of strength S ω in a Pekeris waveguide involving a homogeneous water column of depth D 1 , sound speed c 1 , and density ρ 1 , overlying a homogeneous fluid half-space of sound speed c 2 and density ρ 2 .A schematic of this problem is shown in Figure 19a.A Cartesian coordinate system is chosen, with the vertical z axis passing through the source and the x axis being parallel with the sea surface.The location of the source is denoted by (0, z s ).For this problem, the two-dimensional Helmholtz equation for the displacement potential ψ(x, z) takes the form With the following Fourier transform pair, the depth-dependent wave equation can be shown to be where x .The depth-dependent Green's function in the water (0 where x , with k 1 denoting the water wavenumber at frequency ω.At the bottom (z > D 1 ), the depth-dependent Green's function is represented by where k z,2 is the vertical wavenumber at the bottom, defined as follows, to satisfy the radiation condition for z → ∞, By imposing the boundary conditions of vanishing pressure at the sea surface, and the continuity of the normal partial displacement and pressure at the bottom, we reach the following linear system of equations for the amplitudes of the solution to Equation (33).
Once the amplitudes A + 1 , A − 1 , and A + 2 are obtained by solving the above linear system of equations, the depth-dependent Green's function can then be evaluated by Equations (34) and (35).Then, the displacement potential ψ(x, z) can be computed by evaluating the inverse Fourier transform given in Equation (31).
The transmission loss for the line source problem is defined by TL(x, z) = −20 log 10 p(x, z) p 0 (r = 1) , where and the reference pressure takes the form (assuming the line source is located in water) (1) The solution calculated by the wavenumber integration method can serve as the benchmark solution.
We now consider a specific environment, as schematically shown in Figure 20.The source depth is 20 m.The sound speed and density in water are 1500 m/s and 1000 kg/m 3 , respectively; the sound speed and density at the bottom are 1650 m/s and 1500 kg/m 3 , respectively.We firstly calculate the transmission loss with frequency from 105 Hz to 125 Hz. Figure 21 shows the results of transmission loss versus range and frequency at 40 m receiver depth using different methods.The horizontal axis denotes range and the longitudinal axis denotes frequency.Figure 21a is the transmission loss calculated by the popular normal mode model, KRAKEN.Figure 21b is calculated by the wavenumber integration method given above, which is regarded as the standard solution.Figure 21c is calculated by the finite element model.It can be seen from Figure 21 that the finite element solution is consistent with the standard solution.However, the result calculated by KRAKEN shows obvious errors at certain frequencies.It is evident that the result by KRAKEN is irregular around 112 Hz, whereas the results by the finite element model and by the wavenumber integration method are in excellent agreement.Below, we will consider two specific frequencies, 112 Hz and 118 Hz.At frequency 112 Hz, the bottom wavenumber is approximately 0.426495 m −1 , while the horizontal wavenumbers of the first five modes by KRAKEN are 0.464455, 0.449802, 0.381011 + 0.916079E-2i , 0.309800 + 0.178554E-1i, 0.188504 + 0.392254E-1i.At frequency 118 Hz the bottom wavenumber is approximately 0.449343 m −1 , while the horizontal wavenumbers of the first five modes by KRAKEN are 0.489759, 0.475685, 0.451877, 0.411541 + 0.790571E-2i, and 0.346650 + 0.154633E-1i.From Figure 22a, we can see that at frequency 112 Hz, the horizontal wavenumber of mode 3 is so close to the bottom wavenumber, which is the branch point for this problem, that KRAKEN fails to find this mode.As a result, as shown in Figures 21 and 22c, the field results of KRAKEN are inaccurate.However, at frequency 118 Hz, Figure 22b shows that the horizontal wavenumber of mode 3 is sufficiently far away from the bottom wavenumber such that KRAKEN succeeds in finding this mode.As a result, as shown in Figures 21 and 22d, the field results by KRAKEN are in excellent agreement with the standard solution calculated by the wavenumber integration method.In addition, it is evident from Figures 21 and 22 that the results by the finite element model and the wavenumber integration method are in excellent agreement throughout the whole frequency band of interest.
The above numerical results indicate that the finite element model can provide accurate field results for general Pekeris problems and is well-suited for the problems of sound propagation.However, in the special case where a mode lies very close to the branch cut, KRAKEN might fail to provide accurate field results due to the neglect of the branch line integral.For a single-frequency calculation, the probability of the occurrence of a mode lying near the branch cut is not great.For broadband calculations, however, the probability is high.Another scenario with a high probability encountering this phenomenon for single-frequency problems is when the bathymetry varies with range.However, the finite element method calculates the sound field by forming the stiffness matrix and the load vector and is stable wherever the branch cut is.As a result, the finite element model has more universality than KRAKEN.

Conclusions
In this paper, we introduce a finite element model in Cartesian coordinates for underwater sound field calculation in a two-dimensional environment.
At the beginning, the details of implementation were given.We firstly derived the weak form of the Helmholtz equation using the Galerkin method, then discretized the physical domain using four-node quadrilateral elements, and introduced the calculation of local stiffness matrices and assembly of global stiffness matrix.
Then, different methods to deal with boundary conditions, especially simulating the boundless area, were discussed and compared.The Dirichlet boundary and Neumann boundary were easy to deal with.As for the boundless area-that is, regarding how to set an artificial boundary condition to simulate media in an infinite area, we compared the traditional radiation boundary condition and perfectly matched layer approach and found that the traditional radiation boundary may cause noticeable errors because of the nonzero reflection coefficient at the artificial boundary.
After that, range-dependent problems were calculated using the finite element model, such as the wedge-shaped benchmark problem and the slope problem, with either a totally reflected bottom or a penetrable bottom.All of the problems were solved with satisfying results.
In the end, we took the Pekeris waveguide problem as an example to illustrate the universality of the presented finite element model.We used the wavenumber integration method to provide the benchmark solution and compared the results calculated by KRAKEN and the finite element model.The results show that KRAKEN may be inaccurate at the frequencies where a mode lies too close to the branch cut, whereas the results of our finite element model are in excellent agreement with the benchmark solution.
In summary, our finite element model has the following advantages: • High accuracy.Our finite element model achieves high accuracy in both rangedependent problems and inhomogeneous media problems.• Great universality.Compared with the popular KRAKEN model, our finite element model has better universality and is well-suited for more complex problems.• Focused on range-dependent underwater sound propagation problems.

•
Flexible.Since our model is entirely developed by ourselves instead of relying on commercial software, it is easier to modify and optimize.
However, there are still shortcomings of the model, as listed below; we will improve this model in the future.The work in this paper forms a solid foundation for future research.

•
The integral method.While calculating the elements of stiffness matrices, we used the trapezoid method to evaluate the integral, which is not efficient enough.We will try other methods, such as the Gauss integral method, to improve the accuracy and efficiency.

•
The selection of meshes.Only four-node quadrilateral elements are used in this model.As a result, when calculating the sound field in a complex environment especially with complex shape of boundaries, the meshes cannot match the boundaries very well.We plan to optimize the selection of meshes, combining triangular elements and quadrilateral elements to discretize the physical domain.
In conclusion, in this paper, a finite element model has been presented to calculate the underwater sound field in complex environments with high precision and great universality.

Figure 1 .
Figure 1.Schematic diagram of local indices, midpoint coordinates, length, and width of the fournode quadrilateral element.

Figure 2 .
Figure 2. Schematic diagram of a wedge-shaped waveguide bounded above by a planar pressurerelease sea surface and below by a planar rigid bottom.A line source is involved, which is parallel to the apex and located at (r s , θ s ).The inclination of the waveguide is denoted by θ 0 .
is the schematic diagram of specific environmental parameters of the wedge-shaped waveguide.The sound speed and density of sea water are c = 1500 m/s and ρ = 1000 kg/m 3 , respectively.The source frequency is f = 250 Hz.The depth and range of source are z s = 20 m and x s = 50 m, respectively.The inclination of the wedge is 45 • .

Figure 3 .
Figure 3. Schematic diagram of specific environmental parameters of a wedge-shaped waveguide.

Figure 4 Figure 4 .
Figure4is the sound field calculated by the analytical solution.An obvious interference pattern between the source and wedge vertex can be observed.

Figure 5 .
Figure 5. Schematic diagram of discretization of physical domain and PML in a wedge-shaped waveguide using four-node quadrilateral elements.

Figure
Figure 6a is the sound field calculated by the finite element method with the traditional radiation boundary condition.It takes MATLAB almost 4 s to calculate on an i7-4790@3.60GHz CPU. Figure 6b is obtained with the use of PML, where σ xm = 50, δ x = 6 m, and the PML is divided into 15 layers vertically.It takes almost 10 s.Compared with Figure4, obvious error caused by interference can be seen in Figure6a.This is because the traditional radiation boundary condition cannot totally absorb the outgoing wave-that is, the reflection coefficient at the boundary is not zero.However, the result in Figure6bis consistent with that in Figure4.Figure7is the transmission loss in the wedge-shaped waveguide at 250 Hz at 10 m receiving depth.The blue solid line indicates the analytical solution, the red dashed line indicates finite element solution with traditional radiation boundary condition, and the green dashed line indicates finite element solution with PML.It can be seen that the green dashed line matches the blue solid line better, indicating that PML is more appropriate for simulating the boundless area than the traditional radiation boundary condition, and hence, yields higher precision.

Figure 6 .Figure 7 .
Figure 6.Transmission loss in a wedge-shaped waveguide at 250 Hz calculated by the finite element method (a) with traditional radiation condition and (b) with perfectly matched layer.

Figure 8 .
Figure 8. Schematic diagram of a slope-shaped waveguide bounded above by a planar pressurerelease sea surface and below by a rigid bottom.A line source is involved, which is parallel to the apex and located at (x s , z s ).
is the schematic diagram of specific environmental parameters of a slope-shaped waveguide.The sound speed and density of sea water are c = 1500 m/s and ρ = 1000 kg/m 3 , respectively.The source frequency is f = 250 Hz.The depth and range of source are z s = 20 m and x s = 50 m, respectively.

Figure 9 .
Figure 9. Schematic diagram of specific environmental parameters of a slope-shaped waveguide.

Figure 10 Figure 10 .
Figure 10 is the reference solution of transmission loss calculated by DGMCM2D.An evident interference pattern is observed between the source and wedge vertex and also in the flat region, compared with Figure 4.

Figure 11 .Figure 12 .
Figure 11.Schematic diagram of discretization of physical domain in a slope-shaped waveguide using four-node quadrilateral elements.

Figure 13 Figure 13 .
Figure 13 is the transmission loss in the slope-shaped waveguide at 250 Hz at 10 m receiving depth.The blue solid line indicates the reference solution calculated by DGMCM2D and the red dashed line indicates the finite element solution with PML.It can be seen that the red dashed line matches very well with the blue solid line, indicating that the finite element model has high precision.

Figure 14 .
Figure 14.Schematic diagram of a slope with penetrable bottom and a planar pressure-release sea surface.A line source is involved, which is parallel to the apex and located at (x s , z s ).

Figure 16 .
Figure 16.Schematic diagram of specific environmental parameters of a slope-shaped waveguide with a penetrable sea bottom.

Figure 17 .
Figure 17.Transmission loss at 30 Hz in a slope-shaped waveguide with a penetrable bottom.(a) Calculated by DGMCM2D.(b) Calculated by the finite element model.

Figure 18 .
Figure 18.Transmission loss at 50 Hz in a slope-shaped waveguide with a penetrable bottom.(a) Calculated by DGMCM2D.(b) Calculated by the finite element model .

with k 2 Figure 19 .
Figure 19.Pekeris waveguide with unlimited depth.(a) Schematic diagram of the environment.(b) Schematic diagram of upgoing and downgoing waves in two-layer medium.

Figure 20 .
Figure 20.Schematic diagram of specific environmental parameters of two-layer medium with unlimited depth.

Figure 21 .
Figure 21.Transmission loss at depth of 40 m from 105 Hz to 125 Hz using different methods.(a) Solution calculated by normal mode model KRAKEN.(b) Benchmark solution calculated by wavenumber integration method.(c) Finite element solution.The depth-dependent Green's function at depth 40 m calculated by Equation (34) at frequencies 112 Hz and 118 Hz are shown in Figure 22a,b, respectively, where the bottom wavenumbers are indicated by the red dashed lines.Figure 22c,d show results of transmission loss versus range at a depth of 40 m computed by the wavenumber integration method (blue solid lines), KRAKEN (green dashed lines), and the finite element model (red dashed lines) at the two specific frequencies-112 Hz and 118 Hz.At frequency 112 Hz, the bottom wavenumber is approximately 0.426495 m −1 , while the horizontal wavenumbers of the first five modes by KRAKEN are 0.464455, 0.449802, 0.381011 + 0.916079E-2i , 0.309800 + 0.178554E-1i, 0.188504 + 0.392254E-1i.At frequency 118 Hz the bottom wavenumber is approximately 0.449343 m −1 , while the horizontal wavenumbers of the first five modes by KRAKEN are 0.489759, 0.475685, 0.451877, 0.411541 + 0.790571E-2i, and 0.346650 + 0.154633E-1i.From Figure22a, we can see that at frequency 112 Hz, the horizontal wavenumber of mode 3 is so close to the bottom wavenumber, which is the branch point for this problem, that KRAKEN fails to find this mode.As a result, as shown in Figures21 and 22c, the field results of KRAKEN are inaccurate.However, at frequency 118 Hz, Figure22bshows that the horizontal wavenumber of mode 3 is sufficiently far away from the bottom wavenumber such that KRAKEN succeeds in finding this mode.As a result, as shown in Figures21 and 22d, the field results by KRAKEN are in excellent agreement with the standard solution calculated by the wavenumber integration method.In addition, it is evident from Figures21 and 22that the results by the finite element model and the wavenumber integration method are in excellent agreement throughout the whole frequency band of interest.The above numerical results indicate that the finite element model can provide accurate field results for general Pekeris problems and is well-suited for the problems of sound propagation.However, in the special case where a mode lies very close to the branch cut, KRAKEN might fail to provide accurate field results due to the neglect of the branch line integral.For a single-frequency calculation, the probability of the occurrence of a mode lying near the branch cut is not great.For broadband calculations, however, the probability is high.Another scenario with a high probability encountering this phenomenon for single-frequency problems is when the bathymetry varies with range.However, the finite

Figure 22 .
Figure 22.(a) Magnitude of the depth-dependent Green's function at a depth of 40 m at 112 Hz; (b) transmission loss at a depth of 40 m at 112 Hz; (c) magnitude of the depth-dependent Green's function at a depth of 40 m at 118 Hz; (d) transmission loss at a depth of 40 m at 118 Hz.In panels (a,c), the red dashed lines indicate bottom wavenumbers.In panels (b,d), the blue solid lines indicate benchmark solution calculated by wavenumber integration method, the red dashed lines indicate finite element solution, and the green dashed lines indicate transmission loss calculated by KRAKEN model.