A three-dimensional, longitudinally-invariant finite element model for acoustic propagation in shallow water waveguides

High fidelity models of acoustic propagation in shallow water waveguides are necessary for sonar performance predictions and acoustic communication. However, the waveguides are often very complex including interface roughness, sediment patchiness, sound speed profiles, and complex bathymetry. The finite element method (FEM) provides an excellent tool to model these environments because it is fully customizable and provides a full wave solution. One drawback of the method is its intense computational load. There are currently no fully three-dimensional finite element shallow water waveguide acoustic propagation models although two-dimensional models do exist.¹ Because two-dimensional models cannot be compared with experimental data, a longitudinally-invariant three-dimensional model is proposed to bridge the gap. In this scheme, a cosine transform is performed on the three-dimensional Helmholtz equation, producing a sum of two-dimensional components that can be solved using current methods. Although this produces a fully three-dimensional model, one constraint is that the geometry must be invariant in one direction. Many environments including continental slopes, ridges, and canyons can be approximated as longitudinally invariant. Within this constraint, the method can still include range dependent effects such as sediment and sound speed variation.

In this FEM, the variational form of the Helmholtz equation is solved over small subdomains or elements using polynomial basis sets. The solution is obtained by considering the boundary conditions between the elements that lead to a linear system of equations. The collection of interconnected elements, know as a mesh, represents the total field. A rigorous mathematical description of the finite element method is given in Ref. 2, while a description of its application to underwater acoustics is given in Ref. 3. The FEM makes no approximations to the Helmholtz equation, and the system geometry can be expressed accurately to the size of the element. This makes the model extremely versatile. For the calculations in this paper, a commercially available finite element program, comsol, is used for meshing and solving.⁴ 

Full three-dimensional FEMs for ocean waveguides are still beyond the computation capabilities available. However, the physics of three-dimensional propagation can be captured if the geometry modeled is invariant in one direction. In this scheme, a cosine transform is performed on the three-dimensional Helmholtz equations along one spatial axis. The resulting two-dimensional pressure components can be calculated using two-dimensional FEM.

Following Ref. 5, the derivation begins with the three-dimensional Helmholtz equation, appropriate for a time harmonic acoustic propagation problem,

Here $P¯$ is the pressure field, k is the acoustic wavenumber, k = ω/c(x, y, z) where ω is the radial acoustic frequency and c is the sound speed that can vary with position. The source frequency spectrum is s(ω). In this study, only single frequencies are considered. Last, (x_o, y_o, z_o) is the point source location.

If the sound speed and other environmental factors are only functions of range, x, and depth, z, the three-dimensional Helmholtz equation can be reduced to

where

The integral can be written discretely as

Here k_y is the out-of-plane wavenumber and $P¯¯(ω,x,ky,z)$ is the transformed field. Using Eq. (2), the two-dimensional (2D) pressure fields are computed for a series of out-of-plane wavenumbers using 2D FE and summed to determine the 3D field.

Two parameters in this derivation must be addressed: The limits of the sum in Eq. (4) and the discretization of the sum. To determine these parameters, it is helpful to consider the canonical problem of a point source in an unbounded medium at 500 Hz. Shown in Fig. 1(a) is the integrand of Eq. (3) as a function of the normalized out-of-plane wavenumber and range from the source. Note that the value of the integrand goes to zero smoothly past k_y = k where the effective wavenumber in Eq. (2) is complex. This evanescent part of the integral describes the curving wavefronts near the source and is only significant at short ranges. Also, note that the integrand is smooth as k_y → 0. This suggests that the integrand should be sampled more finely near k_y = k. Two different sampling schemes are shown in Fig. 1(b). The straight line has constant discretization in k_y. The curved line uses a slightly offset gamma cumulative distribution function (CDF) coupled with a linear function just after k_y = k. This results in a much finer sampling near k_y = k for the same number of samples. For k_y > k the samples are discretized much more coarsely because the integrand is smooth. The sampling scheme is shown for the integral at range 10 m from the source in Fig. 1(c). Note how the variable discretization samples the integrand more finely near k_y = k where the function is more variable. Both the variable and the constant discretization schemes used 300 points.

The effect of the different discretization schemes is shown in Fig. 1(d) where the results of the two schemes are compared with the analytic solution. The constant scheme did not resolve the integrand near k_y = k and results in an unstable solution. The variable discretization scheme, on the other hand, results in a near perfect solution that deviates less than the width of the analytic solution line. This example illustrates how a variable sampling scheme can provide a stable solution with fewer samples.

The structure of the integrand for an ocean waveguide is much more complex. However, it retains many of the same qualities as the free space solution. It approaches zero smoothly for values of k_y > k and only has significant values for the evanescent wavenumbers at short ranges. It is much more complex near k_y = k and becomes smooth as k_y → 0. Therefore a similar scheme in which the discretization is determined by the gamma CDF can be used. For propagation in a range dependent waveguide described in the following text, the variable discretization method resulted in a 30-fold decrease in the number of evaluations of the 2D model relative to the equivalent constant method. To put this in perspective, running the fully parallelized model on the Texas Advanced Computing Center LoneStar Cluster took 3 days with the variable discretization method compared with a projected 90 days using the constant method.⁶ 

At first glance, it may seem that the integral could be evaluated with quadrature schemes. However, the integrand suffers from many of the same problems discussed by Jensen, Kuperman, Porter, and Schmidt in Ref. 7, Sec. 4.5.2, with respect to the evaluation of the inverse Hankel transform for axially symmetric problems. In particular, the oscillatory nature of the integral is highly dependent on range and environment requiring different quadrature points for each instance. The authors of Ref. 7 rejected quadrature scheme in favor of the simple equidistant discretization. However, as discussed previously, this is impractical for this problem. Through trial and error, the discretization based on the cumulative CDF has been found to be the most robust and efficient for this problem.

To validate the model, range independent and dependent test cases were selected to compare with existing models. The model parameters are given in Table I. The range independent case consisted of a 200 m deep waveguide with a sandy bottom. This problem required 2500 evaluations of the 2D field components. The solution was compared with a wavenumber integration solution provided by oases (Ocean Acoustic and Seismic Exploration Synthesis).⁸ Results are shown in Fig. 2. The models agree with excellent precision.

The range dependent case is the ASA benchmark wedge, case III, with a lossy penetrable bottom.⁹ This case required 6100 evaluations of the 2D field components. The pressure field along the source axis up the wedge is shown in Fig. 3. This can be directly compared with Fig. 8 in Ref. 9. It should be noted that unlike the calculations provided by Jensen and Ferla in Ref. 9, the FE model includes all orders of coupling and scattering making it an excellent benchmark solution. To display the out-of-plane features of the model, a cut was taken along the gray line in Fig. 3 and shown in Fig. 4. Note the modal cut-offs as the field approaches the wedge apex. This can be compared qualitatively with Fig. 7 of Ref. 10 although in Ref. 10, the source was considerably closer to the wedge apex. Many of the same physical features such as shadow zones and mode cut-offs are evident.

A 2D finite element propagation model has been extended to three dimensions through a cosine transform that is appropriate for scenarios in which the geometry along one spatial coordinate is invariant. The model requires that a 2D field be calculated for each out-of-plane wavenumber. An inverse cosine transform yields the final solution. It was determined that a constant sampling interval scheme to compute the transform was computationally prohibitive. A non-uniform spacing based on an offset gamma cumulative distribution function was found to be robust and efficient.

The model was validated by comparison with range independent and range dependent benchmark problems. The range independent case was a shallow water waveguide of water over sand, appropriate for a continental shelf environment. The range dependent case was a wedge environment appropriate for propagation near shore. In both instances, the finite element model agreed well with available solutions.

Because the FEM solves the Helmholtz equation exactly to the degree of the discretization and is fully customizable, it can be applied to complex range dependent environments such as continental shelf regions where interface roughness, range dependent sediment parameters, and water column sound speed profiles are important. Unlike many other models, it provides the entire pressure field, both forward and backward propagation. Therefore it can be used for both propagation and reverberation modeling.

This work was supported by ONR, Ocean Acoustics under the direction of Robert Headrick.