Three-dimensional boundary fitted parabolic-equation model of underwater sound propagation

The parabolic-equation (PE) method, initially introduced by Hardin and Tappert (1973) and Tappert (1974) to underwater acoustics, has been considered to be one of the most efficient and accurate methods for modeling underwater sound propagation in complex environments. The method has been applied to many two-dimensional (2-D) problems (e.g., Tappert, 1977; Thomson and Chapman, 1983; Collins 1989, 1993; Jensen et al., 1994), and there are also a number of three-dimensional (3-D) sound propagation numerical models employing the PE method (Siegmann et al., 1985; Lee and Schultz, 1995; Smith, 1999; Sturm, 2005; Lin et al., 2012; Lin et al., 2013; Heaney and Campbell, 2016). In general, PE models utilize some type of split-step algorithms to either split the free field propagator and the medium sound speed anomalies (separating the physical processes), or split the directional derivatives in the acoustics governing equation (separating the mathematical operators). In this paper, the latter algorithm is adapted for establishing a 3-D PE method with a boundary fitted grid to solve problems with non-planar or tilted boundaries.

The principle of the PE method is to approximate the Helmholtz wave equation of elliptic type to a one-way wave equation of parabolic type by first factorizing the Helmholtz operator to forward and backward propagation components of square-root power and then neglecting the backward component for backscattering (Hardin and Tappert 1973; Tappert, 1974). This enables marching solution schemes to solve the problem as a transport equation with one less directional derivative at each marching step. Thus, in 3-D applications, only 2-D derivatives are solved at a time (Lin et al., 2012). However, this can still be a computational burden because approximations of 2-D derivatives do not consist of banded matrices that can be easily handled in numerical computation. The alternative is to further split the 2-D derivatives into one-dimensional (1-D) operators (Siegmann et al., 1985; Lee and Schultz, 1995; Sturm, 2005; Lin et al., 2012; Heaney and Campbell, 2016), following the method of alternating directions in computational fluid dynamics (CFD), which is also referred to the alternating direction implicit (ADI) method (Peaceman and Rachford, 1955; Douglas and Rachford, 1956). Although there are similarities, splitting a PE operator can be very different from splitting a CFD operator. The reason is that the PE method involves the square-root Helmholtz operator, which needs an additional iterative scheme to account for higher-order operator splitting (Lin et al., 2012). The details will be addressed in Sec. II.

A boundary fitted grid requires variable spacing between grid points, and such a solution mesh design has been applied to 2-D PE methods for predicting sound scattering from a rough sea surface (Rosenberg, 1999) and for efficiently handling thin sediment layers in the seabed (Sanders and Collins, 2013). Both of the applications were made possible by generalizing the Galerkin approximation of the depth operator in the 2-D PE (Collins, 1989) onto a non-uniform grid. In the application of modeling scattering from a rough sea surface, the boundary grid points were fixed exactly onto the rough surface, and the pressure-release boundary condition, which does not have directivity, was enforced at those grid points. In another application of handling thin sediment layers, the computational efficiency was improved by only reducing grid spacing in the layers instead of the whole model domain. This non-uniform Galerkin method for the 2-D PE is extended in this paper for 3-D modeling.

The rest of the paper is organized as follows. The 3-D PE approximation using a higher-order operator splitting with cross-derivative terms is reviewed in Sec. II, and the 3-D PE method with a boundary fitted grid is introduced in Sec. III. Two benchmark problems to test the method are shown in Sec. IV. Last, the paper is concluded in Sec. V with remarks for future research directions.

This section reviews the 3-D PE approximation method and a higher-order split-step Padé algorithm. First consider the following Helmholtz wave equation:

where p is the sound pressure, $∇=∂/∂x ex+∂/∂y ey+∂/∂z ez$ is the vector differential operator in a 3-D Cartesian coordinate system (x, y, z), and e_x, e_y, and e_z are unit position vectors on each of the dimensions, respectively. Also in Eq. (1), ρ is the medium density, k₀ = ω/c₀ is the reference wavenumber, c₀ is the reference sound speed, ω = 2πf is the angular frequency of sound, and n is the index of refraction n = c₀/c = k/k₀. The index of refraction n and the medium density ρ are considered to be 3-D spatial functions. To enable marching solution algorithms which conserve energy, one can change the pressure variable to $u=p exp(−ik0x)/ρc$⁠, which is the so-called impedance reduced pressure variable (Collins and Westwood, 1991), and the baseline phase is removed according to the reference wavenumber k₀.

Mathematically, the 3-D Helmholtz equation is classified as an elliptic partial differential equation yielding a boundary value problem. To solve such an equation, all boundary values or conditions around the 3-D problem domain must be given, and the whole sound pressure field has to be solved at once, as depicted in Fig. 1(a). This will require a tremendous amount of computation resources when the problem domain is large. To enable more efficient solution schemes, one can apply the PE approximation to neglect backscattering (Hardin and Tappert, 1973; Tappert, 1974) and transfer the Helmholtz equation of elliptic type to the one-way parabolic wave equation that can be solved with marching solution algorithms as shown in Fig. 1(b). This PE approximation involves a factorization that takes a square root of the Helmholtz operator and keeps only the outgoing propagation component, and the one-way parabolic wave equation in the Cartesian coordinate system can be written as the following:

where x is the solution marching direction, and $Y$ and $Z$ are two separated operators containing partial derivatives with respect to either y or z:

Here, a differential operator splitting is implemented, and the coefficient w is a weighting factor between 0 and 1 to control the splitting of medium sound speed anomalies into $Y$ and $Z$⁠. This weighting factor can be $w=Δz/(Δy+Δz)$⁠, which is related to the ratio of the model grid intervals. The one-way wave equation resulting from using cylindrical coordinates will also have a similar form (Lee and Schultz, 1995; Sturm, 2005; Heaney and Campbell, 2016). Note that this operator splitting can be implemented according to a more physical intuition, i.e., to split the free field propagator and the medium sound speed anomalies, which will result in the so-call split-step Fourier PE originally introduced by Tappert (1974 and 1977). In this paper, the operator splitting according to the mathematical operators is chosen, and it eventually leads to a split-step Padé algorithm that enables the ADI method and improves the computation efficiency by converting a 2-D differential operation into a sequence of 1-D operations, as depicted in Fig. 2.

Lin et al. (2012) introduced a higher-order split-step Padé algorithm to solve Eq. (2a) based on the following approximation to the square root of the Helmholtz operator:

where the non-commutative property of $Y$ and $Z$ is kept, and the last two cross terms are essential to improve the approximation accuracy. The derivation of this approximation is reviewed in  Appendix A, and an error analysis is provided in  Appendix B. The higher-order PE solution to Eq. (2a) can be formally written as the next marching solution from x to x + Δx:

where the two-step splitting on $−1+1+Z$ is the so-called Strang splitting (Strang, 1968), making the error to be of order (Δx)², namely, $Err=(Δx)2⟦(−1+1+Y),(−1+1+Y)(−1+1+Z)+(−1+1+Z)(−1+1+Y)⟧$⁠. The square brackets denote $⟦A,B⟧=AB−BA$⁠, which is a measure of non-commutativity between $A$ and $B$⁠. Other details on error analysis of split-step PE solutions are provided by Jensen et al. (1994). In this higher-order solution, the exponentiated cross term is acting as a solution corrector for wider propagation angles, and it can be carried out with the following iteration scheme proposed by Lin et al. (2012):

Combining Eqs. (4) and (5) yields an efficient numerical scheme to implement the 2-D square-root Helmholtz operator $eik0Δx(−1+1+Y+Z)$ as a sequence of 1-D differential operations by alternating derivatives between ∂y and ∂z at each PE marching step from x to x+Δx. This alternating scheme indeed follows the principle of the ADI method (Peaceman and Rachford, 1955; Douglas and Rachford, 1956). Furthermore, the 1-D differential operators in Eqs. (4) and (5) can be computed through the following Padé approximations:

and

where L is the order of the approximation, and $G$ can be either $Y$ or $Z$⁠. Equation (6a) was proposed by Collins (1993) for higher-order 2-D PE models, and Eq. (6b) was proposed by Milinazzo et al. (1997) to improve the accuracy of the PE approximation in the evanescent spectrum by rotating the complex branch line. A 60° rotation is chosen for the examples shown in this paper. The Padé coefficients $αl,L$⁠, $βl,L$⁠, $C0,L$⁠, $Al,L,$ and $Bl,L$ in Eqs. (6a) and (6b) can be computed using formulas provided by Collins (1993) and Milinazzo et al. (1997). These two approximations along with Eq. (3) are the key components of the higher-order split-step Padé algorithm. A brief error analysis for the square-root operator splitting shown in Eq. (3) was provided by Lin et al. (2012), and the analysis is extended in  Appendix B to facilitate discussions given later in this paper.

The key to form a boundary fitted grid is to allow variable grid spacing. Specifically, as shown in Fig. 3, the grid points at the edge of the model domain (denoted as open circles) are fixed onto the curved boundary with variable intervals to the nearest grid points. The advantage of the boundary fitted grid is most profound in the case of fitting to a pressure release surface because the boundary condition p = 0 is a scalar and has no direction. This boundary fitting principally requires 2-D non-uniform mesh generation at each PE marching step. However, with the split-step Padé algorithm, it only requires 1-D non-uniform mesh generation for the unidirectional derivative in the alternating direction procedure.

Exploiting the advantage of the split-step Padé algorithm on operator splitting, one can simply adapt the 1-D non-uniform Galerkin discretization proposed by Rosenberg [1999, Eqs. (14) and (15) therein] and Sanders and Collins [2013, Eqs. (15) and (16) therein] to approximate either $(−1+1+Y)$ or $(−1+1+Z)$⁠. This leads to the following two essential finite-difference formulas for the operators:

where a, b, and c denote general functions of a single variable η, which can be either y or z, and they match with the coefficient functions in Eqs. (2b) and (2c). In addition, $hi=ηi−ηi−1$ is the grid interval between the (i−1)th and ith grid points, i.e., Δy_i = y_i − y_i−1, or Δz_i = z_i − z_i−1, and $γi=(hi+hi+1)/2$⁠.

As shown by Sanders and Collins (2013), the non-uniform Glaerkin discretization can also be employed to model the sharpness of medium interfaces in sound propagation models. The implementation is straightforward by fitting a set of grid points onto the interface and placing another set of grid points very close to the interface to preserve the sharpness. Such an application is also applied in the boundary fitted PE model presented in this paper.

The first benchmark problem to test the boundary fitted 3-D PE method for curved boundaries is a semi-circular waveguide problem, depicted in Fig. 4. This is a canonical problem, and the PE solution can be rigorously examined by comparing to an analytical solution. Recall the 3-D wave equation with cylindrical coordinates:

where x is the waveguide axis, and r and ϕ are the radius and orientation angle of the semi-circular plane perpendicular to the x axis. The variable p is sound pressure, and k is the medium wavenumber parameter. The right hand side of the equation is a point source function with unit intensity at (x, r, ϕ) = (0, r₀, ϕ₀), and a pressure release boundary is enforced around the boundary, i.e., $p(x,a,ϕ)=p(x,r,0)=p(x,r,π)=0$⁠, where a is the radius of the waveguide.

This semi-circular waveguide problem is in fact a proper Sturm–Liouville problem (Arfken and Weber, 2005) with two sets of orthogonal normal modes along ϕ and r:

where m is the angular mode number, n is the radial mode number, and together they give a unique identification to each of the radial modes. Besides that, the radial modes are composed of the Bessel function of the first kind, J_m(κ_mnr), where the angular mode number is the order of the Bessel function, and the coefficient κ_mn in the argument is determined by the Bessel function zeros divided by the waveguide radius, i.e., κ_mn = j_m,n/a. With these two sets of orthogonal normal modes, which are known as cylindrical harmonics, the sound pressure field inside the waveguide can be determined as the following:

and the propagation angle of each Fourier−Bessel mode (m, n) is

The upper panels of Fig. 5 show the first nine normal modes of the semi-circular waveguide for m = 1 − 3 and n = 1–3 with the following waveguide parameters: the source frequency 200 Hz, the medium sound speed 1500 m/s, and the waveguide radius 50 m. One can see that the mode numbers (m, n) correspond to the number of antinodes along ϕ and r. The lower panel of Fig. 5 shows the modal propagation angle for m = 1–9 and n = 1–5, and we will examine the dependency of PE solution errors on the propagation angle later.

The pressure field solution, Eq. (10), expresses that each normal mode propagates in the waveguide at its own phase speed $2πf/k2−κmn2$ without any amplitude attenuation. This behavior readily provides an excellent means to check 3-D PE solutions. One can simply apply a single mode function at x = 0 to initiate the PE solution marching scheme and examine its error as the solution marches outwards along the x axis. Such computations were made for initial mode number m = 1–9 and n = 1–5 and a propagation distance 5 km in this example.

Two different PE solution schemes were implemented; the first scheme employed the boundary fitted grid, and the second scheme employed a rectangular grid which did not fit to the semi-circular boundary. Both of the PE methods still utilized Cartesian coordinates, and the grid intervals were dy = 0.25 m and dz = 0.25 m with a solution marching step of 1.875 m (quarter wavelength). The reference sound speed was given to be 1500 m/s. In addition, the order of the Padé approximation was chosen to be 2 for the exponentiated PE propagator, Eq. (6a), and 3 for the square-root Helmholtz operator, Eq. (6b). Figure 6 shows the amplitude errors of these two PE solutions compared to the Fourier−Bessel mode solution in terms of the error rate dB/km computed over a 5-km distance along the waveguide. It is clear that the boundary fitted PE outperformed the rectangular grid PE, especially around the boundary and the nodal curves of modes, and its error can be even smaller than 0.05 dB/km for lower order modes which have narrower propagation angles.

The phase errors of the PE solutions for the semi-circular waveguide problem are shown in Fig. 7. As seen in the amplitude error plot Fig. 6, the boundary fitted PE can also have less phase errors than the rectangular grid PE. For low-order modes with propagation angles less than 20°, which in ocean acoustics applications are particularly important as they have less transmission loss and can propagate for long distance, the relative phase errors, defined as $Δk/k2−κmn2$⁠, were as low as 0.001 per mille (‰) in the boundary fitted solution, while they were about 0.05‰ (50 times worse) in the rectangular grid solution. The error dependency on the modal propagation angle is obvious. Note that, even though the relative phase error of the boundary fitted PE solution can increase to 1.5‰ for higher-order modes with propagation angles greater than 35°, it is only slightly greater than the inherent error (0.001) due to the square-root Helmholtz operator splitting (see the discussion in  Appendix B). This example indeed demonstrates the accuracy of the boundary fitted PE method.

The second benchmark problem depicted in Fig. 8 is the idealized wedge problem proposed by Jensen and Ferla (1990) for the Acoustical Society of America (ASA) to test range-dependent numerical models of underwater sound propagation. The slope angle is π/63 rad (∼2.86°), yielding a 1/20 slope, and the medium properties, including sound speed, density, and attenuation coefficient, are 1500 m/s, 1000 kg/m³, and 0 dB per wavelength in the water and 1700 m/s, 1500 kg/m³, and 0.5 dB per wavelength in the bottom. A point source of unit intensity is located 4 km (horizontal distance) away from the wedge apex at 100 m depth below the sea surface, and the source frequency is 25 Hz. A theoretical 3-D sound field solution obtained from the method of images by Deane and Buckingham (1993) is used to examine 3-D PE solutions, where the Padé approximations had the same orders as in the preceding example.

To test the boundary fitted 3-D PE, the Cartesian coordinates are rotated so that the x-y plane is aligned with the sloping bottom, and the sea surface becomes a tilted surface in the rotated coordinate system, as shown in Fig. 8. In the calculation, the nominal grid intervals were dy = 5 m and dz = 1 m with a solution marching step of 10 m. Besides fitting to the sea surface with variable grid intervals, the boundary fitted PE also preserved the sharpness of the bottom interface by placing a series of grid points just below and parallel to the sea floor, which enables a better approximation of the interface conditions using stair-step profiles. The sea floor grid interval was 0.2 m (a 340th of the acoustic wavelength in the bottom). To demonstrate the performance of the boundary fitted PE, a rectangular grid with the same nominal intervals was also utilized without either fitting to the sea surface or adding extract grid points with small intervals to preserve the sharpness of the bottom interface. The PE solutions with the reference sound speed c₀ = 1500 m/s are shown in Fig. 9, along with comparisons to the theoretical solution obtained from the method of images (Deane and Buckingham, 1993).

A 3-D perspective plot showing the boundary fitted PE solution for the ASA wedge problem is provided in Fig. 9(a), and one can clearly observe the modal interference patterns in both horizontal and vertical planes. The change of the vertical interference pattern along the x axis, as seen in the cross sections every 5 km, suggests that each of the vertical modes had different horizontal refraction, and a lower order mode could travel a longer horizontal distance before deflecting away from the wedge apex. This makes physical sense because a lower order mode has a smaller vertical propagation angle, hence, yielding a smaller out-of-plane reflection from the sloping bottom. The horizontal refraction of modes eventually formed acoustic caustics along the hyperbolic envelops bounding modal propagation paths. The caustic loci of the first three modes determined by a modal theory (Buckingham, 1987) are shown in Fig. 9(b), and the outermost border of the PE transmission loss (TL) indeed followed the first modal caustic. One can also observe in Fig. 9(c) that the vertical modal interference pattern resulted from the boundary fitted PE method underwent transitions right at the caustic locations determined by the normal mode theory. This confirms the accuracy of the boundary fitted PE solution. Furthermore, Fig. 9(d) shows comparisons of the TL solutions obtained from three different methods, and the boundary fitted 3-D PE solution had an excellent agreement with the reference solution from the method of images by Deane and Buckingham (1993). Although it is not shown here, the PE solution without cross-derivative terms in either cylindrical coordinates (Sturm, 2005) or Cartesian coordinates (Lin et al., 2012) can have significant errors around the modal caustics. Both of the boundary fitted and rectangular grid PE methods in this example were implemented with cross-derivative terms, and they indeed captured the TL transition at the modal caustics. However, the rectangular grid PE solution had evident phase shift, indicating non-negligible phase errors. On the other hand, the phase of the boundary fitted PE solution almost perfectly matched with the reference solution even over a propagation distance of more than 400 wavelengths.

A boundary fitted PE method to solve the 3-D Helmholtz wave equation with non-planar or tilted boundaries has been developed. This method exploits the advantage of a higher-order split-step Padé algorithm to apply the non-uniform 1-D Galerkin discretization on the 2-D solution mesh at each PE marching step, and it has been tested with both semi-circular waveguide and tilted wedge problems. The benchmark results have confirmed the accuracy of this boundary fitted PE method.

There are a number of potential applications, including modeling 3-D acoustic scattering from a rough sea surface by fitting model grid points to the rough surface. Underwater sound reflection from a rough sea surface can cause acoustic field fluctuations due to focusing, defocusing, and scattering effects, and the presented boundary fitted PE method is readily applicable to investigate these effects with a forward scattering assumption. When also considering backward scattering, an extension of the 2-D methods proposed by Collins and Evans (1992) and Lingevitch and LePage (2010) for two-way scattering needs to be made, and it is proposed for future research.

Another application of the boundary fitted PE method is to model 3-D sound propagation along the coast in beach environments, which was originally investigated by Collins et al. (1995) for sound propagation perpendicular to the coast in a 2-D configuration. Because the pressure-release boundary condition does not have direction, one can rotate the model coordinates to align with the sloping beach for better numerical accuracy at the bottom interface meanwhile enforcing the pressure-release boundary condition precisely on the sea surface with a boundary fitted grid, as implemented in the ASA wedge example presented in this paper. This application is straightforward and can provide physical insights into the 3-D coupling between acoustic propagation in water and seismic propagation on land.

This work was sponsored by the Office of Naval Research, USA, under Grant Nos. N00014-11-1-0701 and N00014-17-1-2692. The author also thanks Timothy F. Duda and Arthur E. Newhall for useful discussions for developing the method.

This appendix reviews the derivation of the higher-order square-root Helmholtz operator splitting with cross-derivative terms proposed by Lin and Duda (2012). With the two partial differential operators $Y$ and $Z$ defined in Eqs. (2b) and (2c), the square-root Helmholtz operator can be first written as $1+Y+Z$⁠. Now, let $ξ1=−1+1+Y$⁠, and $ξ2=−1+1+Z$⁠, the square-root operator can then be re-written in terms of $ξ1$ and $ξ2$ as follows:

The operator splitting is done by expanding $F(ξ1, ξ2)$ with a two-variable Taylor series around ξ₁ = 0 and ξ₂ = 0. This leads to the following expansion:

Because

for j = 1–2, k = 1–2, and $ℓ$ = 1–2, but j ≠ k ≠ $ℓ$⁠, the expansion can be explicitly expressed as

where the non-commutative property of operators $ξ1$ and $ξ2$ is kept. Finally, keeping terms up to the second order and substituting back $Y$ and $Z$ yield the second order operation splitting formula shown in Eq. (3).

A brief error analysis for the square-root operator splitting shown in Eq. (3) was provided by Lin et al. (2012), and the analysis is extended here to facilitate discussions provided in the computation examples shown in this paper. The analysis follows closely a modal approach proposed by McDaniel (1975) to examine the phase error resulted from the approximation. Considering an idealized case in a free space, as opposed to McDaniel's single modal waveguide, the square-root operator splitting essentially approximates the following dispersion (wavenumber) relation of a plane wave:

where k_x, k_y, and k_z are the wavenumber components on the x, y, and z axis, respectively, and k₀ is the medium reference wavenumber in the free space. To relate the error with the propagation angle, we can denote the wavenumber vector as $k→0=k0 (cos θ e→x+sin θ cos ϕ e→y+sin θ sin ϕ e→z)$⁠, where θ is the inclination angle between $k→0$ and its projection on the x axis (the solution marching direction) and ϕ is the orientation angle between the projection of $k→0$ on the y-z plane and its component in the y axis as shown in Fig. 10(a). Substituting the plane-wave dispersion in Eq. (B1) with the higher-order square-root operator splitting in Eq. (3) yields the following approximation for k_x:

The approximation shown in Eq. (B2) has a relative error compared to the reference wavenumber:

where $εy=1−1−sin2θ cos2ϕ$⁠, $εz=1−1−sin2θ sin2ϕ$ and both are positive. For 2-D propagation, when the wavenumber vector aligns with either the y or z axis, i.e., ϕ = 0, π, or ±π/2, the approximation error vanishes, yielding a perfect operator splitting. The maximum error occurs when ϕ = ±π/4 or ±3π/4 (when the propagation vector directs to the diagonal).

Figures 10(b) and 10(c) show the error contours as a function of the inclination and orientation angles for a range of phase error tolerances. Comparisons between splitting errors with and without the cross-derivative terms are also presented. One can clearly see that the operator splitting with the cross terms has greater approximation accuracy, and the effective inclination angle is in general 10° better. For example, if considering the relative phase error tolerance to be 0.001, which means the phase error will be less than π before the PE solution marches to 500 wavelengths in range, the effective inclination angle can reach up to 35° compared to 20° without cross terms included in the worst case when the orientation angle equals ±π/4 or ±3π/4, see Fig. 10(c).