Underwater Sound Propagation Modeling in a Complex Shallow Water Environment

Three-dimensional (3D) effects can profoundly influence underwater sound propagation in shallow-water environments, hence, affecting the underwater soundscape. Various geological features and coastal oceanographic processes can cause horizontal reflection, refraction, and diffraction of underwater sound. In this work, the ability of a parabolic equation (PE) model to simulate sound propagation in the extremely complicated shallow water environment of Long Island Sound (United States east coast) is investigated. First, the 2D and 3D versions of the PE model are compared with state-of-the-art normal mode and beam tracing models for two idealized cases representing the local environment in the Sound: (i) a 2D 50-m flat bottom and (ii) a 3D shallow water wedge. After that, the PE model is utilized to model sound propagation in three realistic local scenarios in the Sound. Frequencies of 500 and 1500 Hz are considered in all the simulations. In general, transmission loss (TL) results provided by the PE, normal mode and beam tracing models tend to agree with each other. Differences found emerge with (1) increasing the bathymetry complexity, (2) expanding the propagation range, and (3) approaching the limits of model applicability. The TL results from 3D PE simulations indicate that sound propagating along sand bars can experience significant 3D effects. Indeed, for the complex shallow bathymetry found in some areas of Long Island Sound, it is challenging for the models to track the interference effects in the sound pattern. Results emphasize that when choosing an underwater sound propagation model for practical applications in a complex shallow-water environment, a compromise will be made between the numerical model accuracy, computational time, and validity.

Three-dimensional (3D) effects can profoundly influence underwater sound propagation in shallow-water environments, hence, affecting the underwater soundscape. Various geological features and coastal oceanographic processes can cause horizontal reflection, refraction, and diffraction of underwater sound. In this work, the ability of a parabolic equation (PE) model to simulate sound propagation in the extremely complicated shallow water environment of Long Island Sound (United States east coast) is investigated. First, the 2D and 3D versions of the PE model are compared with stateof-the-art normal mode and beam tracing models for two idealized cases representing the local environment in the Sound: (i) a 2D 50-m flat bottom and (ii) a 3D shallow water wedge. After that, the PE model is utilized to model sound propagation in three realistic local scenarios in the Sound. Frequencies of 500 and 1500 Hz are considered in all the simulations. In general, transmission loss (TL) results provided by the PE, normal mode and beam tracing models tend to agree with each other. Differences found emerge with (1) increasing the bathymetry complexity, (2) expanding the propagation range, and (3) approaching the limits of model applicability. The TL results from 3D PE simulations indicate that sound propagating along sand bars can experience significant 3D effects. Indeed, for the complex shallow bathymetry found in some areas of Long Island Sound, it is challenging for the models to track the interference effects in the sound pattern. Results emphasize that when choosing an underwater sound propagation model for practical applications in a complex shallow-water environment, a compromise will be made between the numerical model accuracy, computational time, and validity.

INTRODUCTION
Three-dimensional (3D) effects can profoundly influence underwater sound propagation and hence soundscape at different scales in the ocean (e.g., Duda et al., 2011;Ballard et al., 2015;Heaney and Campbell, 2016;Reilly et al., 2016;Oliveira and Lin, 2019;Reeder and Lin, 2019). In the particular case of coastal seas, a range of physical oceanographic and geological features can cause horizontal reflection, refraction, and diffraction of sound. Numerical models are often used to solve underwater acoustics related problems in realistic complex coastal environments. Examples include the assessment of underwater noise induced by offshore wind farms (Dahl and Dall'Osto, 2017;Lin et al., 2019) and the influence of estuarine salt wedges on sound propagation (Reeder and Lin, 2019). A number of 3D ocean acoustic propagation models have been developed over the past decades (e.g., Jones et al., 1986;Lee et al., 1992;Porter, 1992, Porter, 2016Bucker, 1994;Smith, 1999;Luo and Schmidt, 2009;Heaney and Campbell, 2016;Lin, 2019;DeCourcy and Duda, 2020). Still, simulating underwater sound propagation accurately for fully 3D environments involves significant scientific challenges and can demand high computational costs (e.g., Jensen et al., 2011;Lin et al., 2019). According to their governing equations and numerical schemes, 3D underwater acoustic models can be divided into three main groups: parabolic equation (PE) models (e.g., Lin and Duda, 2012;Heaney and Campbell, 2016), normal mode models (e.g., Porter, 1992;DeCourcy and Duda, 2020), and ray and beam tracing models (e.g., Porter, 2016;Calazan and Rodríguez, 2018).
Human activities in the ocean have increasingly added anthropogenic sounds to underwater environments (e.g., Duarte et al., 2021). Recent research has shown that underwater noise made by human activities, such as seismic airguns, ships, sonars, explosives, or pile drivers, has the potential to impact marine ecosystems. More specifically, the very loud noise of relatively short exposure can harm marine mammals, fish and marine invertebrates (e.g., Hastie et al., 2015;Shannon et al., 2016). As an important application of sound propagation modeling, the European Commission recommends that the Member States combine underwater sound measurements and models to ascertain levels and trends of underwater noise in the oceans and coastal areas (Dekeling et al., 2014). Since there is an increasing interest in using 3D underwater acoustic models for many other shallow water applications, there is a need to understand better the limits and performance of different models that are available. There are a few intermodel comparisons studies (e.g., Porter, 2019), but they are based on idealized benchmark problems. In this work, the ability of a PE model (Lin and Duda, 2012) to simulate sound propagation in extremely complicated shallow water environments as observed in Long Island Sound (United States east coast) is investigated. First, to benchmark the PE model, its performance is evaluated and compared to the normal mode (Porter, 1992) and beam tracing (Porter, 2016) models for two idealized shallow water cases of a 2D 50-m flat bottom and a 3D shallow water wedge representing the local environment in the Sound. After that, the PE model is utilized to simulate sound propagation in three local areas in the Sound. The inter-model comparison is also performed considering realistic environmental conditions, and the modeled source frequencies are 500 and 1500 Hz.
This article is composed of five sections. The background on 3D underwater acoustic models is discussed in section "Underwater Acoustic Models." The study area is introduced in section "Study Area and Simulation Scenarios." Next, the numerical results are presented in section "Numerical Results." Finally, concluding remarks are given in section "Conclusion."

Parabolic Equation
After the introduction of the PE method to underwater acoustics in the 1970s (Tappert, 1977), there has been a continuous development of this method for 3D ocean acoustic problems (e.g., Lee et al., 1992;Avilov, 1995;Smith, 1999;Sturm, 2005;Lin and Duda, 2012;Heaney and Campbell, 2016;Lin, 2019). Nowadays, 3D underwater acoustic PE models are considered one of the most efficient and accurate methods for modeling sound propagation in complex range-dependent environments. Recent applications of 3D PE models in ocean acoustics include, among others, the sound propagation effects of internal waves (Duda et al., 2011), estuarine salt-wedges (Reeder and Lin, 2019), global scale low-frequency acoustics (Heaney and Campbell, 2016;Oliveira and Lin, 2019), and offshore wind farm construction underwater noise .
Different solution schemes have been used to solve 3D PE in underwater acoustics. The finite difference (Lee et al., 1992), Split-Step Fourier (Tappert, 1977;Smith, 1999;Lin and Duda, 2012) and the Split-Step Padé (Collins, 1993;Lin, 2019) are among the most popular schemes. In this regard, the Split-Step Fourier scheme has been considered a good option for practical applications because of its computation speeds. On the other hand, the finite difference and Split-Step Padé schemes have been viewed as more suitable for idealized 3D shallow water benchmark problems because of better accuracy in treating the bottom interface (Jensen et al., 2011). Although these two schemes could undoubtedly provide the most accurate PE solutions for extremely complex and rangedependent bathymetries, associated computational costs could be too high for practical applications. To handle these complex bathymetry cases with a more efficient solution scheme, the Split-Step Fourier PE model has been improved to include higherorder approximation terms for sharper interface smoothing (Lin and Duda, 2012). Even though it cannot treat the exact interface condition, the Split-Step Fourier scheme still has the advantage in faster computation speeds. Thus, it is chosen over the finite difference and Split-Step Padé schemes in many practical 3D problems.
The Split-Step Fourier PE model used in this article is briefed here, and readers are referred to Lin and Duda (2012) for detailed discussion on the solution algorithm. The model solves the approximated 3D Helmholtz wave equation in a Cartesian system by taking a square root of the propagation operator and utilizing the Split-Step Fourier solution marching algorithm (Tappert, 1977). This model uses a density-reduced pressure variable to handle the density variation across the seafloor interface. To maintain the accuracy of the square root Helmholtz operator, the model utilizes the higher order square root approximation with cross terms (Lin and Duda, 2012) that improves the accuracy of the original wide-angle approximation proposed by Thomson and Chapman (1983). This higher order approximation with cross terms also better handles the bottom reflection with sharper interface smoothing.

Normal Mode
The application of normal mode theory in underwater acoustics goes back to the 1940s (Pekeris, 1948). Since the early 1990s several 3D underwater propagation models based on the normal mode theory have been proposed (e.g., Porter, 1992;Luo and Schmidt, 2009;Ballard et al., 2015;DeCourcy and Duda, 2020). In normal mode models, the 2D horizontal refraction equation can be handled by several different techniques, including rays (e.g., Weinberg and Burridge, 1974), Gaussian beams (Porter, 1992), and also PEs (e.g., Petrov et al., 2020). One of the most popular 3D normal mode models in the underwater acoustics community is Kraken3D (Porter, 1992), which is a combined software package of Kraken and Field3D (Porter, 1992). The former code solves the acoustic vertical normal modes, and the latter one computes the modal amplitude as a function of range and azimuth in 3D environments. Unlike another software Field used for 2D environments, Field3D does not have the modecoupling option, and the sound pressure field in Kraken3D is calculated with an adiabatic (uncoupled) mode assumption. Although the adiabatic approach has many advantages, it is not expected to always provide accurate solutions for some frequencies and environments as shown by DeCourcy and Duda (2020). Kraken3D is often run in the azimuth independent (Nx2D) mode, but can consider horizontal refraction using the beam tracing method (Porter, 1992). In this work, the 2D and 3D versions of Kraken are used.

Ray and Beam Tracing
The use of ray tracing in underwater acoustics has a long history. From its first applications in the first quarter of the 20th century until the 1970s, ray tracing was the predominant method for underwater sound propagation modeling. Since the 1980s, several 3D rays and Gaussian beams models have been developed (Červený et al., 1982;Jones et al., 1986;Bucker, 1994;Porter, 2016;Calazan and Rodríguez, 2018).
Due to its approximate nature, the ray tracing theory (in 2D or 3D) is typically better at high frequencies (Jensen et al., 2011), which could present some limitations on some practical applications. However, when dealing with broadband simulations and reverberation calculations, 3D ray models can offer some advantages over other full-waves approaches (Porter, 2019). Moreover, ray models have the advantage that they can trace rays backward if the bathymetry leads to such paths, and that may sometimes make a difference with respect to other models. Recent applications of 3D ray and beam tracing models in ocean acoustics include the study of high-frequency horizontal refraction on the continental shelf of the Florida Straits (Reilly et al., 2016), effects of complex oceanographic and bathymetric variations on sound propagation in the East China Sea (Porter, 2019), and long reverberation tails in the Norwegian coast (Jenserud and Ivansson, 2015).
The 3D beam tracing model Bellhop3D (Porter, 2016(Porter, , 2019) is among the most popular 3D models in underwater acoustics. This model is an extension of the 2D Bellhop model to 3D environments and takes into account ray reflection and refraction in the horizontal plane. In this work, both the 2D and 3D Bellhop models are utilized.

STUDY AREA AND SIMULATION SCENARIOS
The underwater acoustic models are applied in the shallow water environment of Long Island Sound (United States east coast). The reason for selecting Long Island Sound as a modeling example is due to its extremely complicated environment, especially the highly variable bathymetry. The bathymetric model used in the acoustic simulations was constructed using the Montauk, New York 1/3 arc-second (∼10 m) Digital Elevation Model (DEM) from the National Oceanic and Atmospheric Administration (NOAA). The Montauk DEM [National Geophysical Data Center [NGDC], 2007] is part of NOAA's NGDC effort to build high-resolution integrated bathymetric-topographic DEMs for select United States coastal regions.
Two idealized cases of a 2D 50-m flat bottom and a 3D shallow water wedge representing the local environment in the Sound are considered for inter-model comparison. In addition, three realistic scenarios (see Table 1 and   realistic scenarios are modeled in 2D. Finally, in the third realistic scenario, sound propagation along the sand bar considered in the previous scenario is investigated for 3D effects. This scenario along the sand bar is 4.23 km long (from 41.2106 • N 72.0848 • W to 41.2436 • N 72.0597 • W) and 3 km in width with the water depth ranging between roughly 120 to 8 m.
Three types of sound speed profiles (SSPs) are also considered in the simulations: (a) SSP constant c 0 = 1500 m/s, (b) SSP range-independent c (z), and (c) SSP range-dependent c (r, z), where r is the range and z depth. The latter two SSP scenarios are based on temperature and salinity from the Regional Ocean Modeling System (ROMS) Experimental System for Predicting Shelf and Slope Optics (ESPreSSO) model covering the Mid-Atlantic Bight (Rutgers Ocean Modeling Group, 2020). The ROMS ESPreSSO model has a 5-km horizontal resolution and 36-levels in the water column. Sound speed fields are derived from temperature and salinity based on the Mackenzie (1981) equation. For the sake of simplicity, only results from one ROMS-ESPreSSO time is used. Then, simulations represent environmental condition for September 10, 2018 at 00:00.

NUMERICAL RESULTS
Simulations for a 500 and 1500 Hz source were performed for all the idealized and realistic scenarios with a sound source placed at a 20-m depth. The properties of the bottom considered in all the simulations are sound speed c b = 1700m/s, density ρ = 1500kg/m 3 and attenuation β b = 0.5 dB/λ, where λ is the wavelength. Also, for all the simulations, water density is considered to be a constant of ρ w = 1000 kg/m 3 .
Numerical results are presented first for the idealized scenarios and then for the realistic scenarios in the following order: (i) 2D flat bottom, (ii) 3D shallow water wedge, (iii) 2D rough, weakly varying bottom, (iv) 2D propagation over a sand bar, and (v) 3D propagation along a sand bar. In the realistic scenarios, the bathymetry and SSPs considered in the 2D simulations correspond to the central part (y = 0, where y is the cross-range) of the propagation domains (see Figure 1).

Idealized 2D Flat Bottom
The performance of the 2D versions of PE, Kraken, and Bellhop is compared on a 50-m flat bottom with SSP constant (c 0 = 1500 m/s). This simple comparison aims to calibrate the models and ensure that potential differences on more complex simulations are not due to incorrect model setups.
In PE and normal mode models, different frequencies require different numerical discretizations. As the wavelength reduces with increasing frequency, smaller spatial distances between Frontiers in Marine Science | www.frontiersin.org calculation points must be used. Therefore, higher resolutions should be used in a PE model in the range ( x) and water depth ( z). Likewise, z should be smaller for higher frequencies in Kraken, leading to more computational points and increasing the computational cost. As a rule of thumb, to solve the normal mode equation accurately in Kraken, a minimum of 10 points per wavelength should be considered in the water column (z). A x smaller or equal to one wavelength and z approximately lower or equal than one-quarter of wavelength should be used in PE simulations.
For simplicity, the minimum required numerical discretizations were calculated for the most demanding w FIGURE 3 | Schematic diagram of the idealized shallow water wedge. Red dot presents the sound source. Frontiers in Marine Science | www.frontiersin.org frequency, 1500 Hz (nominal wavelength of 1 m), and also used for 500 Hz (nominal wavelength of 3 m). A grid with x = z = 0.1 m was used in Kraken and PE 2D simulations, which correspond to a very high resolution in range. It should be noticed that under standard conditions, the PE needs a much higher range resolution than Kraken-Field to run and provide a result. On the other hand, Kraken-Field can run with lower range resolution. However, if sharp bathymetric or sound speed variations are present, Kraken-Field cannot provide accurate results considering a low range resolution.
In the Bellhop simulations, 10000 Gaussian beams launched in an angular fan from −70 • and 70 • angle were considered. Bellhop settings were decided from a series of numerical convergence tests. Figure 2 compares transmission loss (TL) obtained by the three models for the flat bottom with the isovelocity waveguide. The agreement among the three models for this scenario is excellent for both 1500 and 500 Hz.

Idealized 3D Shallow Water Wedge
An idealized shallow water wedge case (see Figure 3) is simulated with the 3D versions of PE, Kraken, and Bellhop for 500 and 1500 Hz. This case is a shallow water higher frequency adaptation of the wedge problem with a penetrable lossy bottom proposed by Jensen and Ferla (1990) for the Acoustical Society of America (ASA) and has been widely used by the underwater acoustic modeling community for 3D model validation and to study 3D horizontal reflection (e.g., Lin, 2019;Porter, 2019). Similar to the Jensen and Ferla's wedge problem, the slope angle of the shallow water wedge is π/63 rad (1/20 slope), the sound speed in the water is 1500 m/s and the density is 1000 kg/m 3 ; the sound speed in the bottom is c b = 1700 m/s, and the density is ρ = 1500 kg/m 3 . The bottom also incorporates a volume attenuation of β b = 0.5 dB per wavelength. However, in the shallow water wedge, a point source of unit intensity is placed 1 km (horizontal distance) away from the wedge at 20 m depth below the sea surface. The water depth at the source is 50 m, intending to represent the typical characteristics of the propagation domains chosen for Long Island Sound.
In the PE simulations, the marching (range) and transverse (cross-range) grid size for 500 Hz was set to be 0.5 and 0.25 m, respectively. These values decreased to 0.25 and 0.125 m, respectively, for 1500 Hz. In the water depth, the grid size was set to 0.2 and 0.1 m for 500 and 1500 Hz, respectively. Parabolic and Bellhop results for the shallow water wedge scenario agree very well (see Figure 4). The Kraken solution tends to fall after 1.0 km for 1500 Hz and after 0.2 km for 500 Hz. The errors found in 3D Kraken solution are due to neglecting the mode coupling in the model. Sound propagation over a bathymetry like a wedge can induce strong horizontal refraction and mode coupling, which makes it challenging for a model that neglects mode coupling to provide an accurate solution.

Realistic 2D Rough, Weakly Varying Bottom
The rough, weakly varying bottom scenario results are presented in Figure 5 for the 2D versions of PE, Kraken, and Bellhop. The range-independent SSP c(z) is used. The same model parameters considered in the flat bottom scenario (see section "Idealized 2D Flat Bottom") are used here. Overall, the three models agree very well near the source. However, after about 3 km in distance, some minor differences between the three model solutions can be observed for both 500 and 1500 Hz. Although the three solutions are still comparable, the agreement is not as well as the 2D flat bottom case shown in section "Idealized 2D Flat Bottom." In this realistic case, the PE solution tends to be closer to Kraken than the Bellhop solution. Differences can be due to the interpolation of environment (bathymetry and SSP) variables over depth and range. In this regard, PE and Kraken use stair-step profiles to approximate the rough bottom interface condition. The simulation of this scenario highlights the challenge of considering the same waveguide parameters in the three models for the same environmental data.

Realistic 2D Propagation Over a Sand Bar
Transmission loss results from the PE, normal mode and beam tracing models for sound propagation over a sand bar are presented with the range-dependent SSP c(r, z) in Figure 6. These results are all done on a 2D slice, neglecting horizontal refraction. This particular case cuts through the bottom because the receiver is at 20 m (plots on the right of Figure 6), and the seafloor shallows to nearly 8 m deep at the sand bar. The bathymetry of this scenario presents a steep slope at the beginning, with bathymetry changing in 3 km from approximately 80 m depth  to the top of the sand bar. After the top of the sand bar, the depth increases again to 30 m, and then bathymetry presents slight changes for 2 km. As it is observed in Figure 6, PE and Kraken TL comparison looks very good, with 1-2 dB differences in the peak levels after the sand bar. Bellhop agrees with PE and Kraken before the sound reaches the sand bar. However, the quality of Bellhop results tends to fall after beams cross the top of the sand bar. It should be noted that going down to 1500 Hz is pushing a ray model a bit for these water depths. With the complicated bathymetry of the propagation over a sand bar scenario and diffraction happening at every break in the profile, it can be considered hard for a ray model to track the interference pattern. Also, the low-frequency case of 500 Hz at such shallow water is below the threshold of a ray/beam code's validity.

Realistic 3D Propagation Along a Sand Bar
With the PE model being tested in simpler environmental scenarios, simulations were carried out for the propagation along a sand bar with the range-dependent SSP c(r, z) to investigate 3D effects. The same model parameters values used in the shallow water wedge scenario are used here. Apart from the 3D simulations, the azimuth independent Nx2D simulations were also performed. This type of simulation used a 2D model with the same PE technique but neglected the transverse coupling of sound energy and only considered its variation in the radial direction. A comparison between 3D and Nx2D results identifies 3D propagation effects that 2D models cannot detect.
Focusing and horizontal refraction can be seen in the 3D model at different ranges. However, these effects are not reproduced by the Nx2D model (see Figure 7). These 3D effects are clearly induced by the bathymetric changes along the propagation paths. For instance, between 1.5 and 2.8 km range, where the sound bar reaches the lowest depths (from 20 to 10 m depth) at cross-ranges between 0.0 and 0.5 km, strong horizontal reflection can be observed on 3D results. Differences between 3D and Nx2D depth-integrated intensity are higher for 500 Hz than 1500 Hz.

CONCLUSION
A 3D PE underwater sound propagation model was applied to the complex shallow water environment in Long Island Sound, off the United States east coast. First, the 2D and 3D versions of the PE model were compared with state-of-the-art normal mode and beam tracing models for two idealized cases representative of the local environment: (i) a 2D 50-m flat bottom, and (ii) a 3D shallow water wedge. After verifying its performance by comparing to the normal mode and beam tracing models, the PE model is utilized to simulate sound propagation in three realistic local scenarios. Frequencies of 500 and 1500 Hz are considered in all the simulations. Bathymetric models were constructed based on a high-resolution bathymetry (∼10 m horizontal resolution) dataset. Three different realistic sound speed fields were considered in the simulations (constant value, constant profile over the range, and range-dependent).
Overall, TL results provided by the PE, normal mode, and beam tracing models tend to agree with each other. However, differences can emerge with increasing the bathymetry complexity and expanding the range of propagation. For the complex shallow bathymetry found in some areas of Long Island Sound, it can be indeed considered challenging for the models to track the interference sound pattern. Also, the low-frequency cases of 500 Hz at such shallow water are below the threshold of a ray/beam model's validity. In this regard, the model results indicate that when choosing an underwater sound propagation model for practical applications in a complex shallow water environment, a compromise will be made between numerical model accuracy, computational time, and validity.
Results obtained with 3D PE suggest that in the shallow water environment of Long Island Sound bathymetric features such as sand bars can induce significant 3D effects on sound propagation. These effects were confirmed to happen more often for 500 Hz than 1500 Hz sound and can occur, for instance, on sound propagating along local sand bars. This can potentially affect the underwater soundscape dominated by low frequency marine traffic sound in the area.
The methodology used in this study can also be used for other 3D underwater sound propagation studies. More specifically, before using the 3D model, its 2D version can be compared with other state-of-the-art models. In this way, a better understanding of model performance and limitations can be obtained. Moreover, this inter-model comparison can be used as model calibration when field acoustic data is not available for validating the model.
Although this investigation is focused on Long Island Sound, it can provide helpful information for underwater acoustic applications in other shallow areas. This fact assumes special significance given the increasing interest in using underwater acoustic modeling for environmental impacts assessments. Future work also includes inter-model comparison in shallow water environments considering more physical processes known to influence sound propagation, such as scattering from the sea surface. Passive acoustic monitoring of the underwater soundscape with distributed hydrophone arrays in Long Island Sound is also suggested to investigate the 3D propagation effects shown in the modeling work reported in this article.

DATA AVAILABILITY STATEMENT
Kraken and Bellhop input files are available upon request.

AUTHOR CONTRIBUTIONS
TO: numerical simulations, methodology, formal analysis, and writing -original draft. Y-TL and MP: methodology, formal analysis, and writing -review and editing. All authors contributed to the article and approved the submitted version.

FUNDING
TO thanks FCT/MCTES for the financial support to CESAM (UIDP/50017/2020 + UIDB/50017/2020), through national funds. The funding support from the Office of Naval Research for Y-TL via the grant N00014-21-1-2416 was also acknowledged. MP was supported by the Office of Naval Research under contracts N68335-17-C-0553 and N00014-18-C-7007.