Multiscale multiphysics data-informed modeling for three-dimensional ocean acoustic simulation and prediction

The propagation of sound in the sea is often complex in the sense that intricate spatial and temporal variations of energy and energy flux can arise due to heterogeneous sound speed, interaction with a structured seabed, and the motion of the heterogeneous water. In particular, acoustic propagation in shallow areas of the ocean can exhibit strong temporal and spatial variability because the properties of the water column and the bathymetry often have strong gradients. Because of this complexity, simulation of sound propagation is a commonly used tool for examining the effects (Jensen et al., 2011). To represent oceanic conditions, simulations can use sound-speed fields from volumetric ocean models that follow the dynamical rules (i.e., conservation of momentum, conservation of mass, and advection/diffusion of scalar properties). The ocean models can be observationally constrained to move their state toward nature (e.g., Robinson and Lermusiaux, 2004; Moore et al., 2011; Moore et al., 2019), or they can run without data to simulate conditions that are dynamically possible but not matching conditions at any specific time. Simulations can also involve idealized conditions like analytically shaped seabed features, rectilinear or meandering current jets, a single tidal constituent, and so on. What is lacking at this time is a procedure to rapidly simulate three-dimensional (3D) acoustically-relevant nonhydrostatic nonlinear internal waves (NNIWs) that exhibit higher acoustic property gradients than other features, and thus have strong acoustic effects, but which are difficult to predict. Here, we describe a numerical simulation method that allows data constraints to be implemented in a manner often used in regional oceanographic mesoscale simulation and prediction, then inserts NNIWs that are consistent with the data-informed model fields, then builds 3D sound-speed fields for 3D acoustic simulations, enabling study of important acoustic effects of NNIWs.

Generally speaking, aside from the long-recognized effects on sound of seafloor shape, seabed structure, entrained air bubbles, and ice cover, there are three technical hurdles to obtaining complete simulation of sound in the sea: (i) enumerating the state of the time-dependent and spatially-varying ocean, (ii) handling the large dynamic range of the scales of ocean variability, and (iii) handling the fact that the water and sea surface are moving. In virtually all studies so far, the fluid is approximated as motionless (i.e., the low Mach number of the flow is set to zero), decoupling the particle motion of the acoustic pressure wave from the water velocity; this approximation is also used here, so number three is not addressed. A notable exception to the motionless approximation is simulation of sound reflecting from moving boundaries, for example the water surface (Siderius and Porter, 2008). Addressing the other hurdles is the topic of this paper, namely, the accurate generation of data-informed temporal snapshots of ocean condition for 3D acoustic study and prediction in shallow water, and how to compute useable ocean fields in the working area that include all relevant scales of variability, from tens of meters in the horizontal to the mesoscale feature size of roughly 100 km.

Many papers have shown that NNIWs on continental shelves have strong acoustic effects. Propagation along the wave crests leads to strong refraction (e.g., Katsnelson and Pereselkov, 2000; Badiey et al., 2005; Collis et al., 2008; Lin et al., 2009; Duda et al., 2012). Propagation normal to crests gives strong mode coupling that can be examined as a deterministic process (Zhou et al., 1991; Preisig and Duda, 1997; Duda and Preisig, 1999; Duda, 2004) or as a stochastic process (Raghukumar and Colosi, 2014a,b). In fact, refraction and coupling can coexist (Shmelev et al., 2014). The effects of NNIW crest curvature, a common NNIW characteristic, have also been examined (Lynch et al., 2010; Duda et al., 2011; Lin et al., 2013b). The effects are strong where NNIWs are present, with Duda et al. (2011) presenting a formula for the NNIW amplitude required to induce significant effects 3D, but NNIWs are intermittent and often appear in groups separated by areas of considerably weaker wave activity, so the effects on acoustics will also be intermittent.

A framework to predict intermittent NNIW, and thus the intermittent acoustic effects they cause, is described in this paper. A key capability is to predict where and when NNIW packets occur, what direction they are going, and at what speed. Secondary would be prediction of NNIW wave amplitudes and other parameters, but details of nonlinear phenomena are known to be difficult to predict, so predictions of them would be relatively unreliable in comparison to the macroscopic wave packet properties. The framework consists of a sequential set of operations that result in an acoustic field prediction capability for a given location and time. Briefly, the output of a regional ocean model that includes long-wavelength internal waves and tides, but does not include short-wavelength waves which can have nonhydrostatic pressure, is used to build a field of small-scale waves that is added in to make a more inclusive 3D ocean state prediction, which is in turn passed to a 3D acoustic propagation simulator. When ocean data are available, data-driven ocean modeling (Robinson et al., 1998; Lermusiaux et al., 2006; Lermusiaux et al., 2010; Edwards et al., 2015; Wilkin et al., 2011) anchors the output to what is known about the ocean state. When they are unavailable, the framework can still be employed using idealized environmental conditions. We call the framework the Integrated Ocean Dynamics and Acoustics composite model (Duda et al., 2014a; Duda et al., 2014b). The first version of the entire composite model, denoted as the IODA-A model, is presented here.

The six families of operations in the composite model framework are each examined in sequence, shaping the structure of the paper. These are shown in flowchart format in Fig. 1, described as follows: (1) A regional model with surface tides and internal tides is the starting point. Specific aspects of this type of modeling that are important for acoustics are addressed. (2) Estimation of a background state (with no internal waves) in a region of interest for modeling of internal-wave propagation is required. This required separation of internal waves from a background state suffers from entangled time and space scales. An isopycnal surface tracking method is adopted. (3) Internal tides in the regional model must be characterized. For this, internal-tide signals must be extracted from the full field of isopycnal displacements (position differences from the background state) at critical locations, and their properties tabulated to form input for internal tide and NNIW propagation analysis. In this operation, internal tide propagation trajectories are computed. (4) Next, the extended rotation-modified Korteweg-de Vries nonhydrostatic wave model (eKdVf or reKdV) is run. It is initialized and constrained by both the background state and the characterized internal tides. (5) The regional model and eKdVf fields are merged into a set of volumetric ocean state snapshots, with NNIW across-crest resolution of order 10 m and along-crest resolution of order hundreds of meters to a kilometer. (6) Finally, 3D acoustic simulations are run; the parabolic equation method is used here (Lin et al., 2013a).

The need for the sequential framework is illustrated by considering the scale of the 3D acoustic simulations shown later in this paper. Each has a scale of 3 by 6 km, and given the 1-km grid of a regional model with good resolution, there would be at most 28 model profiles in each domain. Therefore, the acoustic results would depict propagation in a gradually varying medium; depicting this in the medium having additional features of smaller scale requires insertion of statistical small-scale sound-speed anomalies (Lermusiaux, 2006), deterministic NNIW anomalies (as we do here), or theoretically generated statistically described sound-field structures from small-scale features. The sequential nature of the system invokes questions of scale separation between the fields of the various model elements. Each model element describes a specific range of scales but can also act as a bridge between scales. For example, the eKdVf model links mesoscale variability, which it does not describe to dissipation-scale activity and which it also does not describe, by showing how large scale conditions control the formation of high-shear NNIW. Scale separation is not formally treated here, but rather, the complementary capabilities and the included processes of the three main dynamical components (regional model, internal tide ray model, eKdVf tide, and NNIW model) are emphasized.

The regional modeling is covered in Sec. II. Section III describes the background state estimation. Section IV covers analysis of internal tides in the regional model. The extended rotation-modified Korteweg-de Vries modeling is explained in Sec. V. Section VI shows how the volumetric sound-speed fields are formed. The 3D acoustic simulation method is discussed in Sec. VII. Section VIII shows some example calculation results. The entire framework is examined in Sec. IX, including issues that arise when the steps are interfaced because of dynamical or other inconsistencies. Section X is a summary.

The first step in the modeling framework is the generation of the best available estimate of the four-dimensional ocean state. For the intended shallow-water focus of this work, contemporary state-of-the-art regional models that incorporate constraints from observations (data assimilation) and extract boundary conditions from the fields of larger-domain models are ideal to use. This use of observations to create “data driven” model fields is desirable for prediction purposes, but non-data-driven studies with idealized or archetypal conditions can also be valuable for research, so data assimilation is an option. Two families of methods are used to build the constrained model fields, including those using sequential field adjustment at the times of observations, and four-dimensional variational methods that incorporate dynamics to allow observations spaced over time to be optimally and collectively incorporated (Robinson et al., 1998; Lermusiaux et al., 2006; Edwards et al., 2015; Wilkin et al., 2011). Regional models that have been used for IODA-A purposes are the Massachusetts Institution of Technology Multidisciplinary Simulation, Estimation, and Assimilation Systems (MIT-MSEAS) primitive equation model (Haley and Lermusiaux, 2010; Haley et al., 2015), the MIT General Circulation Model (MITgcm) (Marshall et al., 1997), the global Hybrid Coordinate Ocean Model (HYCOM; Bleck, 2002; Arbic et al., 2010), and the Regional Ocean Modeling System (ROMS, Shchepetkin and McWilliams, 2005).

The regional modeling setup requirements to fill the needs of the IODA-A framework are (1) proper surface tides; (2) a fine enough horizontal grid to resolve internal tidal waves; and (3) a transition from an offshore region of few NNIW (or none) to a shallow-water region of plentiful NNIW. Data-assimilative 3D regional models universally use the hydrostatic pressure approximation. Models with nonhydrostatic pressure (i.e., with the full vertical equation of motion) will contain NNIW if properly set up with a very fine mesh (Venayagamoorthy and Fringer, 2006) and will not need the nested eKdVf model, but these models are seldom run in 3D because of the computational expense for resolving tidally-forced NNIWs in shallow water, with an example calculation appearing in Duda et al. (2016). Surface tides in the model domain can result from tidal forcing intrinsic to the model, as with global HYCOM, or they can be introduced by including tidal forcing at open boundaries (e.g., Ponte and Cornuelle, 2013).

All surface forcing options (heat, fresh water, and momentum fluxes) for the regional model are allowable in these simulations; the internal-wave modeling is agnostic toward these, but predictive simulations would benefit from adoption of the highest quality regional model fields. Options include (1) data-informed forcing with surface fluxes from weather model reanalysis, real-time products, or forecast products, (2) climatic flux conditions, or (3) no surface fluxes.

Specific modeling output fields used to date for IODA-A NNIW simulations are (1) the MSEAS Primitive Equation “Shallow Water 2006 (SW06)” product for the Middle Atlantic Bight and Gulf Stream (Haley and Lermusiaux, 2010; Kelly and Lermusiaux, 2016), (2) an idealized canyon internal tide MITgcm configuration similar to that in Zhang et al. (2014), but with the stratification increased from the published conditions to foster NNIW development, and (3) the ONR Inner Shelf program 200-m grid and 600-m grid ROMS configurations (Suanda et al., 2017). HYCOM fields have been used only for internal tide modeling (Duda et al., 2018).

Data-assimilative MSEAS output is used extensively to generate the results shown in this paper. MIT-MSEAS has already been utilized extensively for coupled ocean-acoustics predictions. Dynamical effects of the ocean environment on underwater sound propagation were forecast in real time in several ocean regions, e.g., Dabob Bay (Xu et al., 2008), Mediterranean Sea (Lam et al., 2009), and Middle Atlantic Bight (Colin et al., 2013). In Lermusiaux et al. (2010), the complex tidal-to-large-scale dynamics of the northeastern Taiwan ocean region with strong internal tides and their effects of Nx2-D sound propagation were studied and successfully compared to oceanographic and acoustic transmission loss data. The results showed that with a realistic ensemble data assimilation scheme, the coupled ocean-acoustic modeling had predictive skill for both the ocean physics and acoustic fields and their uncertainties.

For the present study, MSEAS was set up for tidal-to-mesoscale dynamical studies in the Middle Atlantic Bight region, using its nonlinear free-surface hydrostatic primitive-equation model, in a configuration with generalized-level vertical-coordinates (100 levels) and implicit two-way nested domains (1 and 3 km resolutions). Figure 2 shows the MSEAS 3-km grid domain, the long-wave eigenspeed field for mode-one baroclinic waves (similar to the speed of mode-one internal tides, see Sec. III B), and extracted amplitudes of mode-one and mode-two semidiurnal internal tides. Features to note are the slow wave speeds and shorter wavelengths in shallow water compared to deep water, the ratio of mode-two to mode-one energy in the shallow areas (much less than one), and the effect of the Gulf Stream at the south end on the wave speeds, caused by a movement of the main thermocline down towards the center of the water column. The eigenspeeds in the figure are mode-one eigenvalues of the equation for long-wave normal modes with no rotation (mode one has the highest speed and the fewest zero crossings of the vertical mode shape ϕ)

which are useful for normal-mode decomposition of baroclinic ocean wavefields. Here, N is the buoyancy frequency, $N=[−g(ρ−1∂ρ/∂z)]1/2$ where ρ is potential density and g is the acceleration of gravity. Figure 3 shows a surface temperature snapshot for the 1-km domain that we use for results in this paper, which is nested into the 3-km domain of Fig. 2.

To test the IODA-A system in a shelf break canyon region where the basic pattern of inhomogeneous internal wave field is known, we carried out an idealized canyon simulation using hydrostatic MITgcm. The model spans 720 and 875 km in the along- and across-shelf directions, respectively. The boundaries are periodic in the along-shelf direction and open in the across-shelf direction. The model horizontal resolution changes from 5 km on the boundaries to 250 m in the central region of the domain. It has 170 vertical layers with grid spacing varying from 1.5 m on the surface to 40 m on the bottom. There are 600 × 720 × 170 (∼73 × 10⁶) grid cells. Similar to the ROMS model in Zhang et al. (2014), this MITgcm model has a Gaussian-shaped canyon incising into the outer-shelf and is forced by M₂ tides on the offshore boundary. The initial stratification of the model differs from that in Zhang et al. (2014) in the uppermost 85 m, where the climatological profile is replaced by a summertime profile measured on the continental shelf. A density snapshot at 30-m below surface from the MITgcm simulation [Fig. 4(a)] shows a patch of negative density anomaly ∼20 km to the right of the canyon axis (facing coast), indicating downward displaced density interface. The density anomaly patch reflects the spatial focus of the internal tides generated at the sources along the 220-m isobath contour line. Due to the canyon topography, the phases of the internal tide sources differ and cause a phased-array-like behavior of the internal tides, forming a horizontal beam of internal tide radiation, resulting in the focus of the internal tide energy (Zhang et al., 2014). Because of the low resolution and the hydrostatic approximation, this simulation does not produce any NNIW. For comparison, we also carried out a 3D nonhydrostatic MITgcm simulation with the same model domain, bathymetry, and forcing, but a much higher resolution (∼25 m in the canyon region) to examine the NNIWs. The total number of grid cells is 1440 × 2160 × 170. This computation is exceedingly expensive: The simulation of 6 M₂ tidal cycles takes four weeks with 720 CPUs. This cost effectively motivates the IODA-A composite model effort. The produced density field at the same time and depth [Fig. 4(b)] shows a packet of NNIWs with wavelength of 100–500 m propagating onshore to the right of the canyon, corresponding to the patch of negative density anomaly in the hydrostatic simulation.

The ocean background state is the temperature, salinity, and current structure that includes features of motion at subtidal frequencies. The background state, without tidal and internal tidal features, is needed to calculate time- and space-dependent internal wave parameters, such as mode shape and modal phase velocity, for our internal wave modeling. This section explains how the background state and the mode parameters are computed.

The next step after obtaining tidally-forced regional model fields is to separate the internal waves features in the model fields, including internal tides, from the background state that they propagate in. This is necessary because the eKdVf model solves for concurrent propagation, in a non-tidal background condition, of internal tides and higher frequency NNIW. A low-pass filtering procedure explained in the next paragraph is adequate to form the background stratification. Tidal band currents are removed by time averaging over an integral number of semidiurnal periods. The eKdVf will recreate a version of the removed internal tides, but advection of the waves by barotropic tidal currents is not included in results shown here. An optional initial step to later enable barotropic tidal advection of the internal waves is to compute and subtract barotropic tidal currents from the model current field using harmonic analysis, or subtract tidal current estimated from another source (the difference will be small for a properly set up model, Logutov and Lermusiaux, 2008) and correct later for their presence by advecting all water column features along tidal ellipse trajectories. This step is not taken here, so horizontal tidal displacements of order 500 m root-mean-square (rms) are absent from maps shown here.

The internal-tide isolation is implemented by finding mean isopycnal heights in moving time-windows (length W); typically W = 72 h. The heights (z′) of isopycnals (surfaces of uniform density ρ) on a dense vertical grid at all positions are found via interpolation from typically hourly model fields, yielding z′(x, y, ρ, t). The time mean of this is the mean height of the isopycnals. Temperatures (T) and salinities (S) at these heights are then tabulated, also via interpolation, creating T(x, y, ρ, t) and S(x, y, ρ, t). The time mean of the heights, T and S give the mean fields Z(x, y, ρ), T_ρ (x, y, ρ), S_ρ(x, y, ρ), and the functions T_B[x, y, Z(ρ)], and S_B[x, y, Z(ρ)]. Using these, T, S, and ρ can be interpolated onto a set of standard depths (z) everywhere in the model x–y grid. The residual heights η = z′-Z are the isopycnal displacements. The inverse function of Z(x, y, ρ) is the background stratification, ρ_B(x, y, z). Figure 5 shows results of this procedure for one specific computation. Finally, the background sound-speed c_B is computed from T_B and S_B.

The background current field is needed for internal tide and eKdVf analysis. Mean geostrophic current is computed using the thermal wind equations with a level of no motion at the seabed. A second background current estimate is computed by simply averaging currents point-by-point over the selected time window.

Once the background state is established, the parameters governing internal wave propagation can be computed and then passed to propagation algorithms. Two sets of parameters are needed. The first set is computed on the x, y plane and governs propagation of mode-one waves of tidal frequency (mode-one internal tides), and the second set contains the coefficients of the mode-one eKdVf. The first set is needed first and their calculation covered here, while their use to determine internal-tide behavior is covered in Sec. IV.

The internal-tide modal parameters within a water column with sheared background currents in geostrophic balance are found using the generalized Jones equation. A reduced version of this equation valid for parallel sheared geostrophic background flow was derived by Jones (1967), and the full version is examined in Duda et al. (2018). Written for horizontally progressive harmonic waves of the form $ψ(x,y,z,t)=ŵ(z)exp (−iωt+ikx+ily)$⁠, the equation is

Here, f is the Coriolis parameter, ρ is the potential density, and $ω′$ is the depth-dependent intrinsic wave frequency in a Lagrangian frame following a geostrophic background current (u₀, v₀), i.e., ω′(z) = ω – ku₀(z) – lv₀(z). Here, solutions of Eq. (2) are called the “full modes” for internal waves propagating in baroclinic currents with arbitrary stable density gradient N². Examination of this equation and of internal wave modes in shear flows appears in Duda et al. (2018). Solution of this equation with boundary conditions of zero displacement at the surface (an approximation) and at the seafloor yields the modal shapes $ŵn(z)$ and modal wavenumbers k_n and l_n, with subscript n denoting eigenvalue index.

The full mode Eq. (2) is a fifth-order equation in wavenumber and is a polynomial eigenvalue problem. The solution method that has been chosen is to find modal phase speed eigenvalues c_n(θ) as a function of geographic direction. The equation is first converted to a fifth-order polynomial function of c(θ)=ω/|k|, where θ is the angle of the wave vector k = (k, l). In general, the set of solutions c_n(θ) at each n is anisotropic, with fastest wave speed found in the downstream direction of the dominating current. For a given location, the N² and current profiles U = [u₀(z), v₀(z)] need to be provided along with a selection of θ values, yielding mode speed c_n (eigenvalue with index n) and $ŵn(z)$ as functions of direction. Mode one and c₁ are primarily used in this paper. Figure 6 shows fields from the 1-km MSEAS model output (also in Fig. 3), which is examined throughout this paper, extracted during one calculation. The figure also shows c₁(x, y) for four internal-wave wavenumber directions θ, all at the M₂ tidal frequency. Along the southwestern boundary of the subgrid area within which c₁ is computed it can be seen that waves moving to the southwest in the direction of the current [Fig. 6(e)] have significantly faster speeds than waves moving to the northeast [Fig. 6(b)].

Note that the modeling procedures (1)–(6) listed in Sec. I can be completed for any area of any tide-inclusive regional model at any modeled time. This paper shows results from an idealized model of a canyon area, and from one rectangular box of the MSEAS model where the SW06 field program took place.

Tracing internal tide rays from the shelf edge area on to the shelf in the background state is one way to estimate the behavior of internal tides in the model. Tidal line filtering or band filtering of the model fields is another method (Kelly and Lermusiaux, 2016), but this does not move us directly to our goal of modeling NNIW, using the eKdVf model, that would be consistent with the model state. The behaviors of internal tides that we require are direction of propagation on the shelf, plus tidal amplitude and tidal waveform at an offshore location prior to propagation onto the shelf. The direction of propagation is found by tracing internal tide rays from offshore locations using the c₁(x, y, θ) fields explained in Sec. III B. Ray initial directions are found by tidal analysis of the displacements η(x, y, z, t) estimated from the model.

This section will show the ray tracing method, methods for determining ray initiation locations, and ray initiation directions. Mode one, with a vertical displacement profile all of the same sign, will be given the majority of the attention because field studies show this to be dominant (e.g., Rippeth and Inall, 2002; Duda and Rainville, 2008).

A complete description of a method for tracing internal wave modal rays in ocean regions having baroclinic background currents appears in Duda et al. (2018). The method is briefly described here. The ray equations are written in terms of the angle α(x, y) = arctan(p_y/p_x), which is the direction of a vector normal to the phase front, with p_x= ∂ϕ/∂x and p_y = ∂ϕ/∂y being phase derivatives of the wave moving along the ray. For rays of mode-one internal waves of a defined frequency, the slowness field S(x, y, α) = c₁⁻¹ evaluated at the phase-normal angle is needed, as well as its derivative with respect to azimuth. Spatial variations of S will arise from x,y variations in the profiles N²(z), u₀(z), and v₀(z). The equations to be integrated to yield ray locus [x(s), y(s)] are

with normalization factor

Here, the coordinates (x, y) define the ray trajectory with arclength increment ds. An auxiliary quantity is the slope M of the ray path given by

At any position on the trajectory the ray angle is β(x, y) = arctan(M). Implementing the equations requires specification of the starting phase-front normal angle α. An initial β can be given and the consistent initial α can be found by solving Eq. (5) with a standard iterative method. For the case of isotropic S (i.e., ∂S/∂θ = 0) these equations simplify to established equations. In the isotropic case α equals β (i.e., skew angle α – β = 0).

Here, the starting locations for mode-one internal tide rays moving onshore are taken to be at or near the locus of the critical bathymetric slope. The critical bathymetric slope exists when the slope passes from subcritical to supercritical, with criticality defined using the normalized bottom slope

where S_B is the bathymetric slope and S_C is the (depth-dependent) internal wave energy propagation characteristic slope satisfying

Here, ζ is the angle of the characteristics with respect to horizontal, and the stratification strength (N²) is evaluated close to the seabed. The slope is critical when γ equals one. At criticality, singular situations can occur in inviscid linear theory such as unbounded energy upon reflection, but reflection is predicted to behave better in practice (Hall et al., 2013).

In many continental slope regions, the bathymetry is steep and γ > 1 on the main slope (the slope is supercritical), γ < 1 on the shelf (shelf is subcritical), and the locus of critical slope lies along the edge of the shelf at the top of the slope (the shelf break). Analysis of internal tide generation in such geometries by Zhang and Duda (2013) and Zhang et al. (2014) shows that internal tides develop at a critical shelf break above a supercritical slope, in agreement with many previous studies, including the illustrative laboratory experiment of Zhang et al. (2008). Furthermore, those papers quantify the energy flux as a function of position and find the maximum flux divergence to lie at the shelf break (the critical locus), with energy flux offshore in deep water and onshore in water shallower than the critical locus. In addition, good representations of internal tides on the shelf can be made by summing mode-one internal tides from point sources at the critical locus (Zhang et al., 2014), further evidence that the critical area makes a good starting point for propagation of mode-one internal tides.

Because of the flux divergence at the critical locus, we set ray initial locations at or near the critical locus. This is found by evaluating Eq. (7) then Eq. (6) everywhere, using the background N(z) (Sec. III A). In the example calculations shown here, for grid line in the along-shelf direction, the cross-shelf grid location that is immediately inshore of the supercritical to subcritical γ transition at the shelf edge is chosen as the critical locus. This makes a set of grid points, one point for each along-shelf grid index, denoted (x_C, y_C). Figure 7 shows an example line of ray starting points (x_S, y_S) in relation to (x_C, y_C) for one study. The (x_C, y_C) locations can be chosen as ray starting points, but other options can be chosen at the user's discretion, for example a line in the vicinity of the critical loci. Typically, the start is adjusted inshore by 1 km because the internal tides seem better suited for analysis. It is easiest to choose ray start points that are on the model grid, where processed data are available, without interpolation, for the further analysis needed for ray starting and eKdVf initialization.

The initial ray phase-front normal directions α, or the initial trajectory angles β, at the starting points are found by tidal analysis of the displacements η(x_S, y_S, z, t) at the ray starting points. The relative phases of the internal tides at the sites determine the ray starting angles. Internal tides created at a critical shelf break are known to have bottom-concentrated energy, with the waves in the form of beams that are consistent with the energy residing in many modes (Gerkema et al., 2004; Carter et al., 2008; Kelly et al., 2010; Zhang et al., 2008; Zhang and Duda, 2013). Therefore, the wave energy near the seabed is analyzed to determine the wave propagation direction. Tidal analysis is performed in time windows equal to the background state averaging duration W, which can be overlapping if high-resolution analysis is desired. Analysis is restricted to the deepest 60% of the water column. At each ray starting point, the displacement time series at the single depth of greatest displacement variance over the time window (typically 72 h) is processed with harmonic tidal analysis to yield the best-fit tidal signals (Pawlowicz et al., 2002), called η_IT(x_S, y_S, t). In the MSEAS-SW06 area simulations, between 10% and 35% of the analyzed displacement energy is fitted to tidal oscillations, varying from site to site. Figure 8(a) shows example displacement time series η(x_S, y_S, z, t) at one (x_S, y_S) ray start position. The computed displacements near the seabed in Fig. 8(a) are more clearly oscillatory and tidal than the displacements closer to the surface. In Fig. 8(b), phase shifts of the internal tidal displacements η_IT(x_S, y_S, t) computed at different ray starting positions are evident, consistent with wave crest bending, or equivalently, with variation of α along the line of starting points.

Using the distances between triplets of (x_S, y_S) points (arranged along a line in the example of Fig. 7), and using the wavelength of the local internal tide with no current assumed, a correlation analysis of the tidal waveforms is used to determine the phase lag and thus the phase propagation direction. For each (x_S, y_S), the time lag Δt of correlation maximum with the wave-forms to the right and left are converted to propagation distance, Δx = ω_M2Δtk_M2⁻¹ using the M₂ semidiurnal frequency and local baroclinic mode-one wavelength, and then to α = tan⁻¹(Δx/Δy).

Trajectories of rays computed with Eqs. (3) and (4) in the background environment are relatively straight, despite the complexity of the current and density variations, compared to behavior of deep-water mode-one waves crossing extreme features such as the Gulf Stream (Duda et al., 2018). Figure 9(c) shows one set of shelf rays, also computed for the conditions in the SW06 experiment modeled with MSEAS. The crossing of rays is governed largely by the initial angles, although adjacent rays (1-km initial spacing) do exhibit different curvatures to the extent that they will cross. In general, rays on the shelf curve much less than they do in deep-water areas with energetic flow such as the Gulf Stream (Duda et al., 2018).

A wave evolution model that includes wave motions with nonhydrostatic pressure as well as the hydrostatic internal tides can be implemented along the internal-tide rays to estimate NNIW behavior. A few models governing internal wave modal evolution are available from the Korteweg-de Vries (KdV) family, with the simplest original KdV equation having no rotation, first-order nonlinear expansion, and applying to plane waves or cylindrically symmetric modal waves such as mode-one baroclinic waves or ocean surface waves (Lee and Beardsley, 1974). The equation that is chosen here is the rotationally modified extended Korteweg-de Vries equation (eKdVf, Holloway et al., 1999; Grimshaw, 1985), written

where c is linear phase speed of the wave mode being analyzed (baroclinic mode one in our case, thus, c₁), and η(s) is the mode amplitude along the path with length increment ds. Variable coefficients that depend on the varying density and current profiles are A_n that appear in the quadratic and cubic nonlinearity terms, the coefficient of nonhydrostatic dispersion B, and finally G of the rotational term (equal to f²/2c in the absence of a background flow). Q is an amplification quantity related to conservation of wave action. When the equation is solved, the output for mode n is labeled η_n(t, s). Dissipation has been added to the system in a manner similar to that of Liu et al. (2004) (but not equivalent) because shelf NNIW have been observed to decay (Sandstrom and Oakey, 1995; Shroyer et al., 2010), and the waves that appear with no drag tend to be larger than observed in nature. The dissipation term takes the form of quadratic drag, F = DC_D|η|η, added to Eq. (8), where D is also a coefficient that depends on the density and current profiles. A feature not yet implemented is a horizontal convergence term to account for amplification of energy by adjacent ray convergence, similar to the Q term for action conservation under depth changes. Also, note that the quadratic drag applies to the internal waves and does account for any barotropic flows (either tides or mean flows). The drag term is an approximate treatment of all neglected dissipation processes.

Equation (8) is transformed for solution. With change in variable W = ηQ^1/2 and coordinate $ξ=∫c−1(s)ds−t$ such that ds = cdξ, the equation becomes

Coefficients: These are computed using properties along the track of the internal-wave modes for which Eq. (8) is valid. The equation is valid for modes Φ(z) consistent with the long-wave limit eigenvalue equation for vertical displacement

with boundary conditions Φ(0) = Φ(-H) = 0, where U_P(z) is the background current U(z) projected onto the internal-tide ray trajectory (the NNIW direction). Note that Eq. (10) reduces to Eq. (1) when U = 0. Note that the modes Φ(z) used for NNIW studies along lines differ from the modes $ŵn(z)$ used to study internal tidal propagation in two dimensions. The expressions for the coefficients are given in Grimshaw et al. (2004); see Alias et al. (2014) for G with background shear. For example

It is evident from the definition that the z unit vector is directed upward.

Initial conditions: Adjusting the mode-one amplitude initial condition η₁(t, s₀) at ray origin s₀ by comparing against data, such as SW06 experiment data, is critical to obtaining useful results from Eq. (9). The most objective way to initialize the mode-one eKdVf at the s₀ initiation points would be to use mode decomposition to extract the mode-one coefficient from MSEAS-derived η(x_S, y_S, z, t), either prior to or after a tidal harmonic analysis. However, the extraction of mode-one energy in this way gives small η₁(t, s₀) values that usually give no NNIW formation at locations where NNIW are common in nature, such as the SW06 site. The lack of NNIW is likely because the eKdVf does not allow for mode coupling which must be occurring and strengthening the mode-one amplitude η₁ as the mode-one wave move onto the shelf in nature. Essentially, mode coupling must occur at the outer shelf, converting most of the high-mode internal tide energy that is not dissipated to mode-one energy.

The conversion of multi-mode internal tidal beams at the shelf break [beam-like in the vertical; Zhang and Duda (2013)] to the dominantly mode-one internal tides and mode-one NNIW observed on shelves (Rippeth and Inall, 2002) is a process that is not well understood or modeled, so initializing Eq. (8) must be ad hoc at this time. Coupled linear equations have been examined (Kelly et al., 2013; Griffiths and Grimshaw, 2007) and weakly nonlinear (Kelly et al., 2016) but not fully nonlinear equations. Taking as the initial condition the displacement time series η_IT(x_S, y_S, t) at the depth of maximum displacement variance at the ray origin (Sec. IV C) resulted in NNIW that were unrealistically large. The ad hoc method used now is to multiply the tidally-analyzed result η_IT(x_S, y_S, t) at the depth of maximum variance [Sec. IV C, Fig. 7(b)] by an attenuation factor and using the result as η₁(t, s₀). Attenuation factors of from 0.5 to 1 have been used.

Displacement field solutions along rays are plotted in Fig. 9(b). Figure 9(d) also shows a displacement timeseries. Initialization attenuation factor 0.7 and drag coefficient zero are used for this example. These waves compare reasonably with waves measured in this area at the time for which the MSEAS reanalysis regional model simulation was built (Kelly and Lermusiaux, 2016), Fig. 10, although the modeled waves have shorter wavelength. Importantly, the ad hoc treatment for the initial conditions, required because KdV-type models do not include the mode coupling energy exchange that is known to occur on the outer shelf, gives uncertainty to the waves amplitudes and the number of waves in the modeled wave packets (or freedom to adjust to match observations).

The mode-one displacement fields with internal tides and NNIW appearing together which are produced by solving the eKdVf along internal tide rays are valid only along the rays, and these assume weak environmental variation normal to the rays. Without regard to whether variation normal to the rays invalidates the eKdVf solutions, we can interpolate features across the gaps between the rays to build 3D sound-speed fields (see ray gaps in Fig. 9). The procedure to build the 3D mode-one NNIW-perturbed sound speed field c_a1(x, y, z, t) is another step where incompatible assumptions in the component models of the system force an ad hoc nature onto the procedure.

The primary ad hoc treatment of the eKdVf output along rays is an imposed grouping of rays that have similar trajectories, and similar NNIW, that are likely to depict waves in nature that would have connected crests. Chosen groups usually have three to six rays, with the 1-km ray origin spacing used here. Groups may cross, and the waves independently connected. When waves cross, their sum will yield an approximation to the field that disregards the fact that sums of nonlinear waves do not yield the true solution.

Two recipes for interpolation after ray group selection are given here. In the first, a computation of the sound-speed field along the rays c_an from mode-n internal waves is made by first calculating the full-depth displacement field ζ₁(s, z, t) = η₁(t, s) Φ(s, z). The mode-one wave-displaced sound speed field is then computed by remapping, c_a1(s, z+ζ₁, t) = c_B(s, z, t). If desired, the field can then be interpolated back onto a uniform z-grid. The field c_a1 can then be found by interpolation between the internal-tidal lines. Because the NNIW have long crests that are relatively coherent (Jackson, 2004), a sensible way to do this is to identify a benchmark location for identifiable events at a specific time, such as the location of the leading NNIW in a packet, and fit piecewise cubic (x, y) curves through space onto which linear interpolations will be mapped. The method of mapping interpolated fields to a curve is explained in Duda et al. (2011).

The second recipe differs by interpolating the mode amplitudes η₁(s, t) onto curves between the lines, then using loosely gridded background conditions to compute ζ₁, c_a1, and c_B onto a finer grid. The mode-amplitude interpolation method was used in Duda et al. (2011) to interpolate into gaps between a line of sensor moorings as NNIW passed by, also for 3D modeling purposes. The interpolation is again best done along curves fitted to (assumed) long-crested NNIW events, as in the previous recipe. With background state known throughout the (x, y) plane, and thus the modes, the mode-one displacement can be computed, ζ₁(x, y, z, t) = η₁(x, y, t) Φ(x, y, z), and from this c_a1(x, y, z, t) can be computed.

Interpolated two-dimensional (2D) η₁ NNIW fields from the family of IODA model runs shown in Figs. 2, 3, and 5–9 are shown in Figs. 11 and 12. The second recipe is used. Figure 11 shows two time snapshots of the mode amplitude field, one 59 min after the other, showing the development of NNIW at the north end of the modeled domain in the 59-min time window. The waves move at approximately 0.75 m/s so the developing packet has moved about 3 km in this time.

The sound propagation component of the composite model is a time-stepped 3D parabolic equation (PE) simulation method (3.5 D method) used by numerous investigators (Oba and Finette, 2002; Finette and Oba, 2003; Heaney and Campbell, 2016). The specific method used in IODA-A is the Cartesian grid 3D split-step Fourier method described in Lin et al. (2013a). The code computes the 3D sound field from a harmonic point source embedded in one vertical boundary wall (x = 0) of a gridded domain, given by complex pressure coefficient as a function of x, y, and z, where the x and y unit vectors are in the horizontal direction and the z unit vector is vertical. Broadband results are computed by aggregating the complex pressures from many frequencies in the chosen frequency band. Many studies utilizing this code have been published (Lin et al., 2009; Lynch et al., 2010; Duda et al., 2011; Lin et al., 2013a).

First of all, the compute domain must be selected. This is typically a rectangle within the larger rectangle of internal tidal ray and NNIW modeling. The source location and frequency (or band) are then selected. The 3D time snapshot of water column sound-speed in the domain is taken from the interpolation process explained in Sec. VI. Before the code can operate, some other inputs and parameters must be provided. The marching step length Δx in the solution algorithm is typically one acoustic wavelength, the transverse grid interval Δy is typically one-half wavelength, and the vertical grid interval Δz is typically one-tenth wavelength. Two other parameters are the update rates for the volumetric sound-speed information, and for the bathymetry; these are updated after a specific marching distance in the “zero-angle” x direction. At larger distance settings, the 3D simulation domain becomes a sequence (in x) of ocean strips with y dependence that is interpolated onto the acoustic grid. An important input field is the water depth in the domain (the bathymetry), which for best results should be on a much finer grid than the grid for the regional dynamical model. US coastal relief model bathymetry (NOAA, 2011), for example, has three-arc-second resolution, about 92 meters. Seabed geoacoustic parameters are also needed to describe the layered medium beneath the water. Here, for illustration, a simple sandy type half-space seabed specification is used: sound speed 1650 m/s, density 1500 kg m⁻³, sound attenuation 0.5 dB/wavelength.

The 3D acoustic simulation in the area with curved internal waves gives results showing complicated acoustic fields that have been seen in similar environments (Duda et al., 2011, Lynch et al., 2010). Figure 13 shows the domain for a 500-Hz simulation of 3 km in the area of curved waves in the upper half of Fig. 12, with x axis directed to the northeast. The source depth is 60 m. Figure 14 shows the waves and the acoustic field level in the final plane at x = 3006 m. The acoustic field inside the NNIW packet, or after transiting the packet, has an entirely different spatial structure than it does away from the packet. This difference, relevant for array signal processing, which can be quantified by either the horizontal acoustic coherence scale or the phase structure function (Duda et al., 2012), is the motivation for this study of NNIW predictability methods.

In this section, the results from two example computations (case studies) are shown. The first case study is idealized in that climatology is used to define the environment, so it is only weakly guided by data. It focuses on the NNIW computation part of the system. The second case study is a full implementation of the system, with the regional model assimilating field data to constrain all downstream aspects of the modeling.

In the first study, the IODA-A model NNIW calculation is tested in the idealized shelf break canyon environment that was discussed in Sec. II. In this situation, the hydrostatic MITgcm canyon model cannot properly model the NNIW that are generated to the right of the canyon. Instead, it produces a smoothed large patch of isopycnal depression at the NNIW location [Fig. 3(a)]. On the contrary, the high-resolution nonhydrostatic model captures the small-wavelength NNIWs nicely [Fig. 3(b)]. Unfortunately, the high required resolution and the slow pace of the solver make this type of direct real-time nonhydrostatic predictions impossible at this time. With the IODA-A model, the NNIWs in the canyon vicinity can be produced in a timely manner. Figure 15 shows the results obtained using the IODA-A composite model. A concentration of NNIWs are reproduced at approximately the same location as in the nonhydrostatic MITgcm model. Figure 15(a) shows the internal-tide ray initial angles, as well as the highly variable internal-tide displacement amplitudes at the ray origins, depicted by the arrow lengths which are scaled by the rate of energy conversion from surface to internal tides (Zhang et al., 2014). Figure 15(b) shows the eKdVf model output along the rays. NNIW appear to the right of the canyon (facing into shallow water), consistent with the full nonhydrostatic modeling (Fig. 4). Note that some outlier trajectories are initiated on the left side of the canyon where internal tides are weak and difficult to isolate from other flow features.

In the second example, signals arriving at synthetic horizontal arrays in the New Jersey shelf NNIW simulation of Figs. 11 and 12 are examined. This example illustrates how modeling like this has a much different sonar performance prediction capability than modeling such as in Colin et al. (2013) would have. That work used an older coarser version of the MSEAS fields than used here and excluded internal tidal bores and NNIW. Figure 16(a) shows the computation domain for a 1000-Hz 3D parabolic equation simulation with source depth 60 m. The variable water depth (between 71 and 77 m) does not strongly affect the acoustic field. Figure 16(b) shows the magnitude of the complex field at 40 m depth, where the NNIWs are clearly seen to interrupt the modal interference pattern that prevails when sound does not encounter NNIWs. The depth-averaged acoustic energy [Fig. 16(c)] shows effects of horizontal refraction in the form of sound focusing by temporary ducting between internal waves, and of shadowing behind the internal waves at (x, y) near (−500, 6200).

Figure 17 shows effects of propagation disruption by the NNIW on array gain examined with array-like horizontal sampling of the field. In this calculation, fields at arrays of 40 elements at 0.75 m spacing (∼30 m length) aligned with the y axis at 40 m depth are examined for NNIW-induced coherence changes and beamformer array gain degradation. Figure 17(a) shows the coherence length scale, which is the point where the normalized y-lagged correlation function of the complex field falls to exp(−1) from unity at zero lag, or the array length, whichever is smaller. The correlation calculation is made as in Duda et al. (2012). Except for a few outliers, the correlation scale away from the NNIW exceeds the 30-m array length (20 wavelengths), but it can drop to as low 4 m (less than three wavelengths) where affected by the NNIW. The physics effects leading to the decorrelation are discussed in Lynch et al. (2010) and Duda et al. (2011). Figure 17(b) shows array-average single-element acoustic signal excess (SE) at the source in dB, simply SE = SL-TL-NL, which is a scaled form of transmission loss (TL) equal to sound pressure level at the source (SL) minus TL minus noise level (NL) at the receiver, with SL 130 dB and NL 85 dB chosen to illustrate array performance. 3D propagation effects are seen in SE, as they are in the coherence length. Figure 17(c) shows simulated array performance. Simulated field values are coherently summed with beam-steering in the direction of the source (conventional beamforming), with the plot showing SE = SL+AG-TL-NL, where AG is the array gain in dB from the summation. Figure 17(d) shows array gain degradation (see caption). The degradation of horizontal array gain within and behind the NNIW is apparent and is quantified. Comparable coherence changes and episodic array gain degradation within internal waves in oceanic experiments have been reported by Orr et al. (2004), Collis et al. (2008), and Duda et al. (2012). Those field measurements of array performance sensitivity to internal wave presence motivated this research into the factors that control internal wave energy and location. The simulations show that degradation can be correlated with NNIW presence.

The variability of AG, diagnosed from the 3D simulations, can be entered into the predictive probability of detection (PPD) framework (Robinson et al., 2002; Lermusiaux et al., 2002; Robinson and Lermusiaux, 2004; Emerson et al., 2015). In PPD, the quantities in the sonar equation are treated as random variables with probability distributions so that SE can also be expressed as a probability distribution. Emerson et al. (2015) examine TL and NL distributions in experimental data, with TL variability being near 2 dB rms and appearing lognormal. The AG in our example (Fig. 18) shows similar standard deviations near 2 dB in small areas inside the internal waves and in the area in the shadow of the waves, although the distributions are not lognormal. Including the AG variability in the PPD SE calculations can quantify the increased uncertainty of detection when using array sensors spanning tens of meters in or near NNIW.

The 3D ocean acoustic simulation system presented here builds on many well-developed component models that have been carefully investigated and validated when used individually for the purposes they were designed for. The challenges lie at the interfaces of the models, where physics inconsistencies force the simulation outputs from one modeling element to be adjusted to fit the framework of the interfaced modeling element. A few challenges have been pointed out in Secs. II–VII, and the chosen solutions (many ad hoc) were explained. However, much work remains to be done.

The most basic inconsistency at a model interface is the hydrostatic-pressure/nonhydrostatic pressure inconsistency. The regional models are restricted to hydrostatic pressure, by computational necessity, whereas the nonhydrostatic pressure is a key element of tidal bore and NNIW physics. To illustrate the regional model restriction, two simulations with nonhydrostatic pressure were performed using the MITgcm. One had tidally forced flow over an idealized canyon with Gaussian shape and was presented by Duda et al. (2016). The simulation required 720 processors operating for four weeks of clock time to simulate six tidal cycles of internal waves. Running a model like this in forecast mode with data assimilation is not plausible because the model pace is so much slower than real time. Thus, the insertion of a KdV-type nonhydrostatic wave model into a domain modeled with a hydrostatic model is likely to persist as a strategy.

This inconsistency means that the regional ocean model output is unreliable for some purposes where significant nonhydrostatic pressure exists and blending regional model output with KdV-type wave models will give unreliable NNIW output in these areas because the KdV initial conditions will be flawed and unphysical. The practice in our model is to use output only from areas with no steep internal waves (consistent with hydrostatic pressure throughout) as input to the internal-tide raytracing and eKdVf modeling. The outer shelf edge (shelf break) site satisfies this criterion for simulations with MSEAS regional model input examined here, and the outer shelf site satisfies this for Inner-Shelf program wave modeling (not shown here) with data-driven ROMS input.

The separation of internal-tide displacements from the background using isopycnal tracking has difficulties when advection brings new water of new density into the uppermost or lowermost depths. This is more of an operational problem than a physical inconsistency, and only occurs when isotherms are tracked independently by site. A 2 D or 3D method of isotherm tracking and time-averaging (to separate waves from evolving background) should be implemented. The methods explained here give reasonable internal-tide initial conditions (Figs. 7, 8, and 15), although the apparent outliers in ray initial direction (Fig. 15) suggest a possible inadequacy under some conditions. The initial conditions provided to the mode-one eKdVf model for NNIW calculation are more challenging, as along-ray radiating energy that does not pass through the tidal harmonic analysis may be essential to include but is difficult to isolate from other time-dependent motions. Effects of coupled-mode internal tide propagation in the ray starting area are not considered. Only near-seabed internal tide (multi-mode) displacements are used, while thermocline-trapped mode-one incident wave displacements (evident in data and models) have not yet been included in the ray and eKdVf modeling.

The results shown here do not include full effects from barotropic tidal currents on the solution, as explained in Sec. III. Instead, these oscillating currents were averaged away in the background state estimate, with a small residual remaining. These currents would serve to advect the entire set of solution fields from Sec. V onward (internal-tide rays, positions of NNIW along rays, 3D sound-speed field). These currents are knowable, and this correction can be made. Note that effects of barotropic tides on internal tide phase, pointed out by Zhang and Duda (2013) are included in the regional model physics, and thus in eKdVf initial condition displacements at ray origins, η₁(t, s₀).

Another inconsistency appears when internal-tide rays cross, and eKdVf wave simulations cannot be justifiably merged into a single result for the crossing site (Figs. 9, 11, and 15 show crossing rays). The two-dimensional evolution can be studied with the Kadomtsev-Petviashvili (KP) equation (Apel et al., 2007), which has only the first-order nonlinear term and has no rotation, consistent with A₂ and f equal to zero in Eq. (8). This equation is consistent only with a weakly-refracted beam of waves, as seen in Fig. 9(c). It is notable that internal-tide refraction on the shelf as in Fig. 9(c) results in ray deflections that are generally less than a few degrees after tens of km of propagation. This refraction is much weaker than can occur in extreme cases of strong currents in deep water such as the Gulf Stream (Duda et al., 2018), so this restriction is tolerable. On the other hand, the importance of the cubic nonlinear term in locations where A₁ (the quadratic nonlinear term coefficient) is zero is well known, and this situation is common on shelves.

Because one of the main objectives of this modeling is to locate NNIW packets for acoustic prediction purposes, and to provide estimates of packet speed, direction, wave height, and quantity of waves, the accuracy of such quantities should be determined using data. The effects of the packets on sound propagation can be strong, as explained in the work cited in the introduction, and the effects can be anisotropic, so knowing NNIW geometries and performing acoustic field predictions may be valuable in some situations. Nonhydrostatic modeling has been shown to give realistic nonlinear internal tides that compare well with data in a setting where bathymetry and tidal flow dominate the response (Rayson et al., 2018), but that effort had horizontal resolution of 100-m or coarser on an unstructured grid, and was not capable of resolving all NNIW. A shortcut approach to the NNIW prediction problem was explored by Min et al. (2019) who used nonhydrostatic MITgcm in an idealized way to study their formation in the South China Sea. Despite some uncertainty that must be assigned to our simulation results, our approach can provide objective estimates of NNIW properties in areas where the required supporting information is available.

We have shown that the modeling and prediction of NNIW can be pursued by a hierarchical composite modeling system. One common feature of the several wave models that make up the system is their restriction to one-way simulated waves. These are, namely, the extended Korteweg-de Vries model, the internal tide ray model, and the acoustic parabolic equation model. The observation of features at sea and analysis of more general (hyperbolic) forms of the governing equations validate the usefulness of the more easily obtained one-way solutions. The fact that the one-way wave methods give usable solutions for actual conditions makes possible the method of multiscale multiphysics ocean acoustic simulation explained here.

Despite the successful modeling to date, it is recognized that each dynamical model in the system has limitations. Together, however, they can reproduce with good reliability:

1. Data-informed tidally-forced (sub)-mesoscale regional conditions.

2. Tidally driven internal waves (internal tides) in their source region.

3. Trajectories of internal tide energy flux.

4. Nonlinear nonhydrostatic internal wave formation and propagation.

The volumes within which these processes are simulated using the methods outlined in the paper contain predicted features that are important for 3D acoustic propagation. An ability to simulate the acoustically relevant features like NNIW has been demonstrated. The quality and uncertainty of the NNIW predictions, and thus acoustic propagation predictions, remains to be examined in follow-up work.

This work was supported by Department of Defense Multidisciplinary University Initiative (MURI) Grant No. N00014-11-1-0701, managed by the Office of Naval Research Ocean Acoustics Program, and National Science Foundation Grant No. OCE-1060430. Final manuscript preparation was supported by ONR Ocean Acoustics Grant Nos. N00014-17-1-2624 and N00014-17-1-2692. P.F.J.L. also thanks ONR and NSF for research support under Grant Nos. N00014-13-1-0518 (Multi-DA) and OCE-1061160 (ShelfIT) to MIT, respectively. The MSEAS-based series of simulations for the New Jersey shelf region examined here was accelerated toward completion by the interest in realistic 3D acoustic fields expressed by Dr. Ivars Kirsteins at the Naval Undersea Warfare Center.