Nonlinear internal waves in shallow water have significant acoustic impacts and cause three-dimensional ducting effects, for example, energy trapping in a duct between curved wavefronts that propagates over long distances. A normal mode approach applied to a three-dimensional idealized parametric model [Lin, McMahon, Lynch, and Siegmann, J. Acoust. Soc. Am. **133**(1), 37–49 (2013)] determines the dependence of such effects on parameters of the features. Specifically, an extension of mode number conservation leads to convenient analytical formulas for along-duct (angular) acoustic wavenumbers. The radial modes are classified into five types depending on geometric characteristics, resulting in five distinct formulas to obtain wavenumber approximations. Examples of their dependence on wavefront curvature and duct width, along with benchmark comparisons, demonstrate approximation accuracy over a broad range of physical values, even including situations where transitions in mode types occur with parameter changes. Horizontal-mode transmission loss contours found from approximate and numerically exact wavenumbers agree well in structure and location of intensity features. Cross-sectional plots show only small differences between pattern phases and amplitudes of the two calculations. The efficiency and accuracy of acoustic wavenumber and field approximations, in combination with the mode-type classifications, suggest their application to determining parameter sensitivity and also to other feature models.

## I. INTRODUCTION

Curved nonlinear internal waves (NIWs), which occur widely and frequently in shallow-water oceans, can produce dramatic acoustic propagation effects. As only one example, a curved NIW train allows acoustic energy to become trapped in ducts between two waves and propagated over long distances. Many such phenomena are observed in the ocean and are investigated by computational and analytical models. The objective of this paper is to determine, using a modal approach, how acoustic quantities vary with changes in the environmental parameters of an idealized three-dimensional (3-D) curved duct feature model (Lin *et al.*^{1}). This work builds upon papers by DeCourcy *et al.* that define a mode-type classification system^{2} and develop convenient formulas for extracting wavenumber parameter dependence^{3} in a curved shelf-slope front model. These techniques are adapted here for cases of shallow-water NIWs over flat ocean bottoms. Using formulas based on mode number conservation, it is possible to obtain wavenumber approximations efficiently that are accurate enough for sensitivity computations of acoustic quantities including transmission loss.

The effects of NIWs are important because of the prevalence of mechanisms by which they are formed.^{4} Propagation in curved waveguides has a long history of modeling interest, and has been explored using several techniques as in Ref. 5. The 3-D NIW acoustic effects gained considerable attention as a result of field data, such as the SWARM95 experiment.^{6} Observations and analysis from that experiment and others, along with numerical simulations, provide striking evidence of 3-D acoustic effects produced by NIWs and ducting.^{7–9} Additional data with NIWs and further analysis followed results from a 2007 experiment on the northern shelf of the South China Sea,^{10} indicating that near-surface ducting plays a significant role in the long-distance propagation observed. This led to feature models for NIW trains^{11} and specification of effects from individual wavefronts. Mirror and 3-D ducting behaviors from a single curved NIW are quantified in Ref. 12, which treats a simpler geometry than the idealized two-wave ducting model in Ref. 1 and this paper.

Long-distance propagation in the Ref. 1 model is primarily the result of ducting by two acoustic mode types that can form in NIW trains. The first are whispering gallery (WG) modes that carry acoustic energy along the concave interface of curved wavefronts. The second are fully bouncing (FB) modes, which are essentially trapped between two wavefronts with very little energy leaking out. Although leaky (L) modes contribute less to in-duct propagation, retaining some of them may be necessary for representing the fields near a source or wavefronts. This NIW model has three types of L modes that are classified according to locations of a turning point, which determines a boundary of partial energy trapping. Conditions for acoustic ducting were previously derived by Lynch *et al.*^{13} and this paper extends their results by predicting wavenumber changes as feature parameters vary. These predictions in turn can be used to estimate the relevant modes and field quantities of interest for a particular set of feature parameters.

Computational approaches to this problem, such as PE methods, are extremely valuable. However, they are usually limited in providing broad quantitative insight into the effects from environmental parameter variations without a large computational burden. Normal mode approaches are helpful for insight, but there is a significant numerical challenge for the present model. The angular wavenumbers are complex-valued and determined from a complicated dispersion relation, in which they appear as (large) orders of Hankel functions. For this and related problems, numerical root-finding methods are known to be slow and error-prone, which inhibits examination of ducting effects across wide ranges of parameter values. Thus, the results here are directed toward finding approximate and still accurate wavenumbers, from formulas that provide physical insight into effects of parameter variations, and are significantly more computationally tractable.

The primary result of this paper is that useful wavenumber approximations can be obtained for the NIW feature model through the principle of mode number conservation. It is well known (see Ref. 14, for example) that acoustic normal modes are analogous to quantum mechanical waves with discrete energy levels. The energy levels are expressed in terms of mode number through Bohr quantization, while the acoustic equivalent is a conservation law previously applied to vertical modes in a Pekeris waveguide by Weston^{15} and Pierce.^{16} Because regions of modal energy trapping are often bounded by turning point locations, they appear as the limit(s) of integration over a modal wavenumber. The invariance of energy levels and mode numbers allows estimating wavenumber responses to parameter changes, provided phase effects from any interface(s) and turning point(s) are incorporated. The approach here extends work by DeCourcy *et al.*^{2,3} for the parameter dependence in an idealized single-interface coastal front over a sloping bottom. The NIW case requires handling pairs of interfaces and subsequent radial mode phase changes, and also new approximation formulas. Treating more than two interfaces would be a natural extension, suggesting that the approach could provide useful approximations for other feature models of interest. Further, the techniques in Ref. 3 can also be applied to the NIW problem in order to identify the environmental parameters to which propagation is most sensitive, and to quantify the effects of the variations. It is possible that an analogous approach could be used to examine consequences of uncertainty in parameter values on predictions obtained from data modeling or numerical computations. Direct applications of these methods require negligible mode coupling effects; possible extensions to handle weak coupling are not considered. Asymptotic methods are also needed, and parameter-magnitude conditions for the present case are provided.

In Sec. II the feature model formulation and its parameters are reviewed, and its normal mode components and wavenumbers defined. Section III introduces a classification scheme for radial duct modes and illustrates their spatial behavior. Section IV describes the wavenumber conservation approach that incorporates both radial-direction and endpoint modal phase changes to specify the responses of angular wavenumbers to changes in model parameters. Angular wavenumber variations from two different parameters are illustrated, and accuracy of the approximations are benchmarked with results from the modal dispersion relation. In Sec. V, modeling assumptions are briefly explored and wavenumber prediction quality is demonstrated. Approximate and benchmark-accurate wavenumbers, along with PE results, are used to generate horizontal transmission loss contour plots, which are compared qualitatively and quantitatively. Finally, conclusions are summarized in Sec. VI. The accuracy and efficiency of these approximations suggest possible extensions and other applications for specifying acoustic dependence and sensitivity.

## II. MODEL FORMULATION AND MODAL SOLUTION

In shallow-water ocean environments, surface warming leads to higher temperatures and sound speeds in the upper ocean. This produces variations that can often be modeled by one or more constant-property layers, and results in the formation of nonlinear internal waves of depression. An idealized 3-D feature model (Lin *et al.*^{1}) of this type is used to quantify effects of wave parameters on acoustic energy propagation in NIW ducts. The feature model differs from some others, for example, one proposed by Katsnelson and Pereselkov,^{11} which has a train of internal waves and a thermocline with continuous depth variation. The environmental complexity and consequent mode coupling requires parabolic equation (PE) methods to efficiently solve the horizontal acoustic problem. However, it was shown^{17} that a continuous thermocline for a single NIW can often be approximated by a square wave with two sharp interfaces. This treatment was extended to a pair of waves and the intermediate duct,^{18} which is used for the feature model here. With this formulation the acoustic propagation is adiabatic, and individual mode behaviors provide useful wavenumber approximations for any frequency, unlike ray theory employed by Ref. 11.

A schematic of the 3-D NIW duct geometry and its cylindrical coordinates is shown in Fig. 1. The ocean is modeled by two layers, where sound speed *c _{W}*(

*r*,

*z*) is constant within the unshaded region (labeled region I) containing a duct, and another (larger) constant within the shaded region (region II) containing the waves and the warm surface layer. The corresponding constant depths of the top layer in these two regions are

*h*

_{I}and

*h*

_{II}. The ocean surface is assumed to be pressure-release, and water has density

*ρ*= 1000 kg/m

_{W}^{3}in both layers. The constants

*H*,

*D*,

*W*, and

*w*represent total water depth, wavefront radius of curvature, duct width, and wave width. The geometry and the constant parameter values contribute to keeping the acoustic propagation problem not only adiabatic but also mathematically separable. No termination conditions in

*θ*are specified for ducts or wavefronts in this model, allowing free propagation in the angular direction. The ocean bottom is modeled as a semi-infinite fluid layer with constant density

*ρ*>

_{B}*ρ*, sound speed

_{W}*c*>

_{B}*c*, and attenuation

_{W}*α*.

_{B}Acoustic pressure for a point source of frequency *f* = *ω*/2*π* in a NIW duct is specified by the 3-D Helmholtz equation

The cylindrical coordinates *r* and *θ* are horizontally oriented polar coordinates, with depth *z* measured from the ocean surface as in Fig. 1. Acoustic pressure is $p(r,\theta ,z,t)=P(r,\theta ,z)e\u2212i\omega t$, with time dependence from a single-frequency point source factored out of Eq. (1). The source is located in the duct at $(r,\theta ,z)=(rs,0,zs)$, where $D<rs<D+W$ and $hI<zs<hII$.

The model coordinate geometry and environmental assumptions lead to a separable boundary value problem, because horizontal and vertical variability only occur at interfaces coinciding with lines of constant *r* and *z*. Thus, mode coupling need not be considered here. The separation parameter is *ζ _{m}*(

*r*), which is the horizontal wave number corresponding to vertical mode number

*m*. The vertical and horizontal components of the pressure solution satisfy

where $\psi m=\psi m(z),\u2009Am=Am(r,\theta )$, and the medium wavenumber is given by $k(r,z)=\omega /c(r,z)$.

In the duct and wave regions, vertical modes satisfy a three-layer Pekeris waveguide. The usual continuity conditions of pressure and normal velocity component are applied at each horizontal interface. Details of solving this familiar problem can be found in Refs. 19 and 20. Because the upper-layer thickness differs in the duct and wave regions, different vertical mode problems are solved in each region. The eigenvalues of Eq. (2a) for the two problems are complex (because of bottom attenuation) and are denoted by *ζ _{m}*(

*r*), which means this quantity is constant for

*r*within each region I and II. The two constant values are denoted for convenience by $\zeta m,j$ where the subscript

*j*is

*j*= I or

*j*= II. These values are needed for solutions of the horizontal problem given by Eq. (2b).

Figure 2(a) is the schematic side view of Fig. 1, and also contains the reference parameter values used for example computations in this paper. The focus here is on propagation in a single duct region, so the width *w* of the waves (reference value *w* = 300 m) is treated as infinite, which is the so-called “well approximation.”^{1,13} As in work by Finette and Oba,^{7} the extended wave regions are assigned the wave sound speed value. The in-duct modal behavior is preserved and the model analysis becomes more tractable (half as many interfaces), at the cost of neglecting effects of acoustic tunneling (Refs. 1 and 18, and Sec. V). The well approximation is suggested in Fig. 2(a) by the solid lines at depth 40 m, and corresponding wavenumber results are expressed as

Subscripts *k* = 1, 2 identify radial interface locations *r*_{1} = *D* and *r*_{2} = *D* +* W*, which are the first (or inshore) and second (or offshore) duct interfaces. Subscript commas are added after vertical and horizontal mode number subscripts *m* and *n* and before any *j* and *k* subscripts. In this paper *D* and *W* are used to specify wave-duct interfaces, except for derivations where the *r _{k}* notation is convenient.

Numerical solutions at *f* = 75 Hz for the first two vertical modes in duct (dashed blue) and wave (solid red) regions appear in Figs. 2(b) and 2(c). Mode *m* = 1 [Fig. 2(b)] has exponential behavior in the upper water layer and oscillatory behavior in the lower water layer for both duct and wave regions. Mode *m* = 2 [Fig. 2(c)] oscillates in both water layers for both regions. All modes shown have exponential decay in the ocean bottom. For *m* > 2 vertical mode solutions differ little between duct and wave regions, and ducting effects vanish. The real parts of *ζ _{m}*(

*r*) are several orders larger than their corresponding imaginary parts; for example, mode

*m*= 2 has complex horizontal wavenumbers $\zeta 2,I=0.3092+5.3803\xd710\u22125i$ and $\zeta 2,II=0.3056+6.0525\xd710\u22125i$. This mode has the greatest difference (both absolute and relative) in phase speeds between the duct and wave regions, making it the lowest to show significant ducting effects;

^{1}it is used for numerical examples for the remainder of this paper. Note that the wavenumbers in both regions specify

*ζ*(

_{m}*r*) from Eq. (3), which appears in the horizontal mode Eq. (2b).

where *G* solves

Continuity conditions on *G* and its first derivative are enforced at each interface between regions I and II. Finiteness at *r* = 0 and the Sommerfeld radiation condition for large *r* are also applied. The integral in Eq. (4) can be evaluated by contour integration in the *η*-plane. Application of all the conditions leads to a formula for *G*, with the result that isolated singularities occur at *η* = *η _{mn}*, where

*η*is an angular wavenumber associated with radial mode

_{mn}*n*corresponding to vertical mode

*m*. Contributions at these points include all the proper normal mode solutions of Eq. (2b), and may include improper modes depending on branch cut choices. For the often-used Pekeris cuts, contributions from the branch line integrals are relatively small except possibly in the source near field

^{22}and are neglected here.

Local wavenumbers corresponding to the modal solutions of Eq. (1) satisfy the relation

where *k _{zm}* is the vertical wavenumber of vertical mode

*m*, and

*k*is the radial wavenumber of horizontal mode

_{rmn}*n*associated with vertical mode

*m*. With approximate expressions for solutions of Eqs. (2a) and (2b), the modal pressure solution is constructed as

where $Gmn(r)=G(r,\eta mn)$. Mode shape functions *G _{mn}*(

*r*) and

*ψ*(

_{m}*z*) are normalized using procedures described in Refs. 1 and 23, and $\gamma mn(rs,zs)=iGmn(rs)\psi m(zs)$ is an amplitude coefficient depending on source location. For the remainder of the paper, mode numbers

*m*and

*n*are suppressed except where defining new quantities or where their omission may be confusing.

## III. RADIAL MODE SOLUTIONS AND CLASSIFICATION

The primary objective of this paper is to develop useful approximations for acoustic quantities that characterize normal-mode sound propagation through models of NIW ducts. These rely on finding efficient and accurate approximations to normal mode phases in terms of environmental parameters. That is accomplished in Sec. IV, while this section provides necessary background for classifying radial mode properties.

The radial modes $Gmn(r)$ are expressed in terms of solutions to Eq. (5) as

The constants in Eq. (8) use the notation of Fig. (5) in Ref. 1, and are determined by applying subsidiary conditions.^{21} The Sommerfeld radiation condition is already incorporated in Eq. (8) for *r* > *r*_{2} by excluding an incoming Hankel function. Boundedness for small *r* is ensured by combining outgoing and incoming waves for *r* < *r*_{1} to form a Bessel function of the first kind using the formula $2J\eta =H\eta (1)+H\eta (2)$. The four remaining coefficients are found using continuity conditions of *G* and its first derivative at *r* = *r*_{1} and *r* = *r*_{2}, along with a normalization condition. The dispersion relation provides the essential fifth condition to determine the constant angular wavenumbers *η* = *η _{mn}*, which were introduced below Eq. (5). They specify the radial wavenumbers $krmn=(\zeta 2\u2212\eta 2/r2)1/2$ and the geometric properties of the radial modes.

Because one objective of this paper is accurate approximations for angular wavenumbers, results must be tested against benchmarks for validity (see Sec. IV). For any set of environmental parameters, the most accurate wavenumber values can be obtained from numerical solutions (referred to as “numerically exact”) of the radial-mode dispersion relation, which is given by Eq. (11) of Ref. 1. The bisection method is used to obtain most of the numerically exact *η* values in this paper. Because of challenges in the numerical evaluation of the complex-order Bessel functions, the algorithm did not converge for some higher radial mode numbers. For these an equivalent dispersion relation formula, expressed in terms of ratios of Hankel functions and their derivatives, is used as in Refs. 2 and 3. For convenience simpler quantities are defined

in which ℓ = 1, 2 indicates Hankel function type. After applying all the subsidiary conditions as before, the dispersion relation formula can be written as

Comparing curves of phase and amplitude for Eq. (10) in the complex plane locates additional roots. This approach is detailed in Ref. 24 and used for a shelf-slope front model. It is very helpful for leaky modes and complements the bisection method used for low *n* values.

Mode normalization is performed after *η* and all coefficients except one (for specificity, *A*_{0}) are determined. The procedure is described in the Appendix of Ref. 1 and is analogous to many other mode treatments, such as for vertical modes in the kraken algorithm.^{23} An initial value of *A*_{0} = 1 could be selected, and corresponding unnormalized radial modes *U* defined for the three regions of *G* in Eq. (8). The conditions at the two vertical wave-duct interfaces can be expressed as

for *k* = 1, 2. Exact expressions for ratios *f* and *g* of radial mode *n* at the interfaces are

which are found using the appropriate region II solutions. As shown in the Appendix of Ref. 1, any of the mode coefficients *A*_{0} is given by

Equations (11) and (13) are equivalent to Eqs. (A4) and (A13) in Ref. 1, using the interface location notation from this paper. Finally, normalized modes can be written as *G* = *A*_{0}*U* with *A*_{0} from Eq. (13).

Angular wavenumbers *η* are important in this paper for specifying quantities of acoustic interest, such as the transmission loss fields illustrated in Sec. V. They are also useful in determining turning point locations, and consequently the geometric properties of radial modes.^{2} The coefficient of *G _{mn}* in Eq. (5) is readily shown from Eq. (6) to be $krj2$. The standard substitution $Gmn(r)=r\u22121/2Qmn(r)$ reduces Eq. (5) with source term omitted approximately to

because 1/4*r*^{2} is small compared with $krj2$ for *r* values of interest here. Radial mode behavior can be ascertained from Eq. (14) in analogy to the constant coefficient case. Positive values of $kr2$ produce oscillatory solutions, while negative coefficients lead to exponential behavior. The case $kr2=0$ (with a simple zero) is a turning point (TP) across which *Q*(*r*) solutions change between oscillatory and exponential. The behavior in range of radial modes can thus be classified according to turning point locations, which along with boundaries or interfaces signify where acoustic energy is trapped.

Figure 3(a) shows the five different radial mode types which emerge in the idealized NIW model. The lowest radial mode numbers *n* are associated with whispering gallery (WG) modes, which have one turning point in the duct (*D* <* r* < *D* +* W*) and a second at the *r* =* D* +* W* interface. The familiar whispering gallery effect occurs in this feature model for shallow grazing angles, which cause acoustic energy to reflect repeatedly along the outer interface of the duct, trapping it over long distances along the front. Increasing grazing angles slightly may produce fully bouncing (FB) modes, which also have two TPs, one at each interface. These modes oscillate across the entire width of the duct, and again play a critical role for long-distance propagation in the along-front direction. Radial modes with steeper grazing angles leak energy across the interface at *r* =* D* +* W*, and are characterized by the occurrence of a single turning point. Although previous work on this model^{1} grouped all three types of leaky modes in a single category, in contrast here the wavenumber approximations require consideration of these cases separately. The first two share properties with WG and FB modes. Leaky whispering gallery (L-WG) modes, which may arise instead of FB modes as model parameters change, have an in-duct turning point (*D* <* r* < *D* +* W*), with acoustic energy again propagated along the outer edge of the duct. Trapping effects are naturally less than those of true WG modes, because their energy penetrates through the *r* =* D* +* W* interface. Similarly, leaky bouncing (L-B) modes have a turning point at the *r* =* D* interface and contribute to the modal solution in the duct and beyond the right interface, while experiencing more along-front energy loss than L-WG and FB modes. Finally, this feature model includes totally leaky (L-T) modes, which have the steepest grazing angles and a turning point for *r* <* D*. Because energy is lost through both wave-duct interfaces, their along-front propagation attenuates more rapidly than for other modes. Nonetheless, it may still be necessary to include some L-T modes in pressure-field computations that approximate the near field.

Figure 3(b) is a schematic that is primarily designed to illustrate the oscillatory and exponential range behavior of the radial modes, and the corresponding regions of ducted and transient energy. The figure has its vertical axis as a scaled azimuthal wavenumber *η*/*D* and its horizontal axis as a scaled range *r*/*D*. The middle region (green online) represents the duct, while the outer (red and blue online) are the two wave regions. Turning points occur where $kr2=\zeta 2\u2212\eta 2/r2=0$, which is equivalent to *η* = *ζr*. Dividing the latter equation by *D* specifies the sloped lines on the figure; there are two because *ζ* differs in the wave and duct regions. These lines are significant because for azimuthal wavenumbers below (or above) them, oscillatory (or exponential) behavior arises from $kr2$ being positive (or negative). For example, square symbols in Fig. 3(a) that correspond to in-duct and in-wave TPs of WG, L-WG, and L-T modes lie on these lines. Also, the two different *ζ* values in the wave and duct may cause one *η* value to produce different mode behavior on each side of a wave-duct interface, as for WG, FB, and L-B modes. Only the real parts of *η* and *ζ* are needed for this figure, because the imaginary parts for both parameters are much smaller than the real. To illustrate use of the figure, select a mode type and a value for *η*/*D* on the vertical axis, and track the behavior as *r* changes. For instance, considering a WG mode (top row of boxes), choose *η*/*D* that lies in the interval $\zeta I<\eta /D<X\zeta I$. Moving from left to right, there are exponentially decaying solutions in the left wave region and left side of the duct, oscillations in the right side of the duct up to the *r* =* D* +* W* interface, and exponential behavior in the right wave region before becoming oscillatory at large enough ranges. An analogous schematic appears as Fig. 1(b) in Ref. 2, to assist in classifying mode types and behaviors in a coastal shelf-slope front model.

For $\eta <X\zeta II$ oscillations leak across the interface at *D* +* W*, and thus the value $X\zeta II$ is a cutoff for leaky modes. Because $X\zeta II<\zeta I$ as shown on the vertical axis of Fig. 3(b), it is not possible for an in-duct turning point to occur in the leaky region, so L-WG modes cannot arise. Similarly, if instead $\zeta I<X\zeta II$ as shown on the vertical axis of Fig. 3(c), it can be argued that FB modes do not arise. The relationship of these two mode types to each other and to WG modes is reinforced by the fact that the FB and L-WG regions are in the same relative positions on the two sub-figures. The applicable sub-figure depends on the relative positions of $\zeta I$ and $X\zeta II$, which in turn is determined by model parameters. Thus for fixed parameter values, only four of the five possible mode types occur.

The quantity *X* = 1 + *W*/*D* (also equal to *r*_{2}/*r*_{1}), appearing in cutoff values that separate mode types, provides physical intuition into the formation of FB and L-WG modes. The NIWs with smaller *X* values (duct width small relative to radius of curvature) support FB modes. This is because larger amounts of acoustic energy is more easily trapped in a narrow region, where interaction occurs between two relatively close interfaces. On the other hand, larger *X* values (relatively wider ducts) have less trapping, and FB modes are supplanted by L-WG modes. In the present problem *X* remains close to 1 because *W*/*D* is typically of order 10^{−2}; nonetheless, small *X* changes are sufficient to produce different modal types. Note that a similar scaling of wave, duct, or front widths by a large parameter such as radius of curvature would be necessary to apply the methodology to other cylindrical feature models. Another useful condition in this problem is that *η* and *ζD* are of similar size, both order 10^{4} for parameter values of interest.

As noted, this mode type classification extends previous work^{1} by distinguishing three leaky mode types. While all three have non-negligible energy passing through the offshore wave-duct interface, their oscillatory behavior begins at different locations: in the duct, at the inshore wave-duct interface, or in the inshore wave region. Consequently, the three types contribute differently to full acoustic field results. In addition, specification of modal TP locations is needed for developing wavenumber approximations in Sec. IV.

## IV. PARAMETER DEPENDENCE OF ANGULAR WAVENUMBERS

Accurate values of *η* that are useful for benchmarking subsequent approximations can be found numerically from the radial-mode dispersion relation. Determining the needed complex roots can be challenging and inefficient, because of complicated dispersion relation formulas such as Eqs. (9) and (10). Further, the results do not provide physical insight into the influence of feature-parameter variations. In contrast, finding the response of *η* values to parameter changes is feasible without solving the dispersion relation, by using mode number conservation that leads to convenient formulas for calculations and for understanding parameter dependence.

Mode number conservation integrals are an acoustic equivalent to Bohr quantization integrals^{14} and were previously applied to vertical waveguides.^{15,16,25,26} This concept was later extended^{3} to radial modes in a shelf-slope front model and used to predict accurate wavenumber variations from its feature-parameter changes. In adapting this approach to the NIW model, it should be mentioned that two different radial mode indexing conventions are in the literature. Reference 1 and this paper use *n* to count the number of local amplitude extrema, that is peaks or troughs, for *r* ≤ *r*_{2}. The formulations in Refs. 3 and 15 instead count the number of modal nodes in that region. Thus, the mode number in corresponding formulas [such as Eq. (8) in Ref. 3] are replaced here by *n* – 1. The left side of the radial mode number conservation Eq. (15) for the NIW model has an integral over a formula for radial wavenumber *k _{r}*, using Eq. (6):

and (fixed) vertical mode numbers *m* are suppressed in the remainder of this section. In Eq. (15)*τ _{n}* is the TP located at the smallest

*r*value for radial mode

*n*. The total phase correction Φ

_{n}is the sum of partial-oscillation phase changes in cycles at any TPs for a mode that occur at endpoints of the integration integral. Thus, for any selected radial mode, the integral accounts for the number of complete and partial oscillations between the most-inshore TP and the offshore duct interface. That value plus Φ

_{n}

*π*produces an excellent approximation to

*π*times the integer

*n*– 1.

It is valuable to make use of the antiderivative of the wavenumber integral in Eq. (15),

As in Ref. 3, define the function *μ*(*ξ*)

The contributions from evaluating the wavenumber integral are summarized next. The most-inshore TP at *r* = *τ*_{n} in Eq. (15) for each mode type is specified in column two of Table I. Contributions from the lower limit of integration should be considered for all five mode types. However, for WG, L-WG, and L-T modes [see Fig. 3(a)], it follows from Eqs. (17a) and (17b) and *τ _{n}* that $\mu (\zeta j\tau n/\eta )=\mu (1)=0$. Thus, the only lower-limit contributions are from TPs at the inshore interface for FB and L-B modes, which have $\u2212\mu I1$ terms in their rows of the third column. The L-T mode type has two terms from piecewise evaluation of the integral on both sides of the interface

*r*=

*D*, which produces $\u2212\mu I1$ and $\mu II1$. All five mode types contain a term $\mu I2$ in column three that arises from the upper limit of integration. At this point it is easy to express Eq. (15) in a concise form to make its terms clear. Denote the sum of the

*μ*contributions from column three for any mode type of number

_{n}*n*as $\Sigma \mu n$. Then including the phase correction Φ

_{n}, dividing by

*π*, and rearranging terms, Eq. (15) is simply expressed (terms in cycles) as

It remains to determine the phase contributions (column four of Table I) from the TPs. These serve as endpoints of the mode oscillation regions that are inshore of the interface at *D* +* W*. It follows directly from Fig. 3(a) that there are seven such TPs: three not at interfaces (two in duct for WG, L-WG; one in wave for L-T), and four at interfaces (one for WG, two for FB, one for L-B). The former three produce terms from the WKBJ connection formulas for Airy functions.^{27,28} These generate a –*π*/2 shift that contributes a cycle shift of *ϕ* = –1/4 to Φ_{n}, as seen in column four. The latter four TPs are modeled using Rayleigh reflection theory, and the phase follows from using radial wavenumbers in the reflection coefficient *R*,

The ratio of region I and II densities is of course one, and only the real parts of *η _{n}* and

*ζ*are used. As modes transition from oscillatory to exponential behavior at an interface, $krI2>0$ and $krII2<0$, so total internal reflection occurs from Eq. (19). Its phase can be represented in polar coordinates using an arctangent, and produces a

*ϕ*= –

*ϕ*value in cycles, where

_{k}At an interface for a leaky mode, oscillations occur on both sides because $krj2>0$. As expected *R* is real and positive, so there is no phase contribution. In summary of results in column four, it is easily seen from Fig. 3(a) that the in-duct TP of WG modes accounts for the phase correction *ϕ* = –1/4 and the interface TP accounts for *ϕ* = –*ϕ*_{2} from Eq. (20). Similar reasoning applies to FB modes with two interface turning points, giving phase corrections *ϕ* = –*ϕ*_{1} and *ϕ* = –*ϕ*_{2}. All three leaky modes have no phase change at the *r* = *r*_{2} interface, and the single TPs contribute *ϕ* = –1/4 for L-WG and L-T modes, and *ϕ* = –*ϕ*_{1} for L-B modes. The L-T modes cross the *r* = *r*_{1} interface, but the Rayleigh reflection coefficient evaluates to zero there. The resulting total phase corrections Φ_{n} are summarized in column four of Table I.

Mode type . | τ
. _{n} | $\u2211\mu n$ . | −Φ_{n}
. |
---|---|---|---|

WG | η/ζ_{I} | μ_{I 2} | 1/4 + ϕ_{2} |

FB | D | $\mu I2\u2212\mu I1$ | ϕ_{1} + ϕ_{2} |

L-WG | $\eta /\zeta I$ | μ_{I 2} | 1/4 |

L-B | D | $\mu I2\u2212\mu I1$ | ϕ_{1} |

L-T | $\eta /\zeta II$ | $\mu I2\u2212\mu I1+\mu II1$ | 1/4 |

Mode type . | τ
. _{n} | $\u2211\mu n$ . | −Φ_{n}
. |
---|---|---|---|

WG | η/ζ_{I} | μ_{I 2} | 1/4 + ϕ_{2} |

FB | D | $\mu I2\u2212\mu I1$ | ϕ_{1} + ϕ_{2} |

L-WG | $\eta /\zeta I$ | μ_{I 2} | 1/4 |

L-B | D | $\mu I2\u2212\mu I1$ | ϕ_{1} |

L-T | $\eta /\zeta II$ | $\mu I2\u2212\mu I1+\mu II1$ | 1/4 |

For given mode numbers *m* and *n*, the expressions in the third and fourth columns of Table I specify all quantities in Eq. (18) except for the value of *η*. The resulting formula provides an implicit equation for *η* as a function of model parameters that is far easier to solve numerically than the dispersion relation. For example, Fig. 4(a) [or Fig. 4(b)] is generated by allowing the parameter of interest *D* (or *W*) to vary while all others are fixed at the reference values. The solutions produce curves (shown dashed) of azimuthal wavenumbers for 10 (or 8) modes versus the varying parameter. The algorithm for obtaining these solutions typically converges to complex *η _{n}* values, of which only the real parts are plotted. Numerical solutions of the dispersion relation for

*η*are shown as symbols, different for each mode type, and are benchmarks for the wavenumber approximations. Overall the approximations show excellent accuracy (lines pass near centers of symbols), with the lowest accuracy for L-T modes. It is striking that even as modes transition from one type to another, approximation accuracy is maintained, as noted for a shelf-slope model.

_{n}^{3}

Figure 4(a) shows responses of scaled azimuthal wavenumber *η*/*D* to changes in radius of curvature *D*, with a fixed duct width *W* = 500 m and all other parameters at their reference values. Horizontal lines (heavy solid) *η*/*D* = *ζ*_{I} and *η*/*D* = *ζ*_{II} indicate two boundaries where changes in mode type and TP (*τ*) locations occur. Although approximation accuracy may slightly decrease close to these curves, improvements are possible by including additional terms. However, this approach sacrifices the convenience of the current versions, especially because accuracy improves to its typical high level by moving away from the curves. One light solid curve *η*/*D* =* Xζ*_{I} is an upper bound for azimuthal wavenumbers, and the other, *η*/*D* =* Xζ*_{II}, separates leaky modes from WG and FB modes. In summary, approximate results show high accuracy across a wide range of parameter values; specifically, in this figure, increasing or decreasing *D* up to 50% from its reference value shows that scaled wavenumber curves have non-linear variations with *D*.

An analogous example is Fig. 4(b) for duct width *W*, with radius of curvature fixed at its reference value of *D* = 50 km. The same overall conclusions apply to this case as to Fig. 4(a), with accurate approximate results when *W* varies by up to ± 40% from its reference value. One feature of Fig. 4(b) is the approximation-benchmark comparisons as duct width decreases and L-B modes move into the L-T region. A few such modes appear in Fig. 4(b), where the solid symbols oscillate slightly around the approximation curves. Because in-duct acoustic propagation is of primary interest here, including additional term(s) to improve that accuracy is not pursued. Overall the wavenumber approximations in this figure appear sufficiently accurate for determining acoustic field properties, as, for example, transmission loss in Sec. V.

## V. TRANSMISSION LOSS APPROXIMATIONS

The objective of this section is to use approximate and numerical-benchmark angular wavenumbers *η* to compute transmission loss (TL) plots for qualitative (contours) and quantitative (level curves) comparisons. Because horizontal field properties are of particular interest in NIW environments, the TL is examined for a single vertical mode *m*. Its contribution to the time-independent pressure solution is

and horizontal-mode TL normalized with respect to a source value (1 m away) is

For a horizontal cross-section at depth *z* = *z _{s}*, Eq. (22) can be expressed as

where

is the normal mode approximation to the integral in Eq. (4).

To demonstrate the capabilities of wavenumber approximations from Sec. IV, TL is computed twice from Eqs. (23) and (24) and twice using a documented parabolic equation (PE) algorithm.^{1} The same model parameter values are used: *D* = 60 km, *W* = 500 m, reference values for others, and frequency *f* = 75 Hz. The second vertical mode *m* = 2 is selected for all examples, because for the reference parameter values (at least), this mode is strongly trapped in the NIW duct. The source is placed inside the duct at $(r,\theta ,z)=(rs,0,zs)$, where $D<rs<D+W$ and $hI<zs<hII$, to further highlight trapping. Figure 5(a) is obtained from Eq. (23) using numerically exact *η* values for *n* = 1 to 10 obtained by solving the full dispersion relation. Radial source coordinate *r _{s}* = 60.425 km is selected, and a specific value for

*z*does not appear in horizontal TL Eqs. (23) and (24). Next, approximate values for Re(

_{s}*η*) at

*D*= 60 km are taken from the curves shown in Fig. 4(a), which are then used in the same TL computation to generate Fig. 5(b). The small imaginary parts of

*η*found from the dispersion relation are used in both calculations, to ensure that differences in Figs. 5(a) and 5(b) arise only from the real parts of

*η*. Figure 5(c) shows the TL in the same waveguide calculated using the PE algorithm, again with vertical mode

*m*= 2. Finally, Fig. 5(d) is a second PE computation, this time with the well approximation removed and the NIW fronts as sketched by the dashed vertical lines in Fig. 2(a), using wave width

*w*= 300 m.

Comparing Figs. 5(a) and 5(b), it is clear that the wavenumber approximations preserve the pattern features of the transmission loss field, including effects from different mode types. For example, mechanisms responsible for long-distance propagation are clearly visible in both images; very similar whispering gallery patterns are observable along the inshore side of the outer duct interface from the four WG modes included. Features from the four leaky modes also remain intact using the approximate wavenumbers; beams leaking across the offshore interface show very minor differences in positioning and intensity. The in-duct results from Fig. 5(b), near the inshore interface region shown in the inset plot, slightly misestimate TL compared with Fig. 5(a). The relatively intense approximate and benchmark fields near the offshore interface are in strong agreement along the entire duct length shown in Figs. 5(a) and 5(b). The magnified insets show the same number and positions of interference beams, although the benchmark field is slightly better defined than the approximation. A key result is the striking pattern agreement, particularly in the duct, between Figs. 5(a), 5(b), and the PE computation in 5(c) that confirms the approximation accuracy compared with both the benchmark modal solution and the high-accuracy PE computation using the well approximation. Finally, Figs. 5(c) and 5(d) illustrate the effectiveness of the well approximation. Using the notation *w* for wave width, the propagation patterns in the wave regions $D\u2212w<r<D$ and $D+W<r<D+W+w$ cannot be accurately captured by the well approximation, but agreement is very good with Fig. 5(d) elsewhere. For example, the inset box in Fig. 5(d) shows partially trapped energy in the inner wave region that is transmitted back into the duct; although the well approximation does not handle this behavior, the in-duct TL patterns show only minor differences. In summary, this example supports the conclusion that the TL field of Fig. 5(a), constructed from approximations using the approach of Sec. IV, preserves almost every TL feature inside and outside the duct, other than in the adjacent wave regions without the well approximation.

Figure 6(a) plots cross sections, along the line *y* = 59.4 km shown (blue) in Fig. 5 subfigures. The solid (blue) curve (labeled “Numerically Exact”) is from Fig. 5(a), and the dashed (red) curve (“Approximate”) is from Fig. 5(b). There are six distinct peaks in the duct interference pattern. Counting from inshore to offshore, peak 1 shows the largest differences in amplitude and pattern phase; peaks 4, 5, and 6 (strongest) have excellent amplitude and phase agreement; and peaks 2 and 3 have very good maximum amplitude agreement (differences of order 1 dB) with small discrepancies in pattern phase. The most observable differences, other than peak 1, are the pattern phase shifts in peaks 2 and 3, and the two smaller fades between peaks 3 and 5. In summary, the TL approximation provides excellent predictions for the three peaks in the low-loss region near the offshore interface, and good predictions elsewhere in the duct except for the lowest peak near the inshore interface. Within both wave regions, the error of the approximation is low.

Figure 6(b) shows incoherent transmission loss along the same cross-sectional line *y* = 59.4 km. The same expression Eq. (23) is used, but now with $Am(r,\theta )$ given instead by

Incoherent loss predictions are expected to agree well, because the influence of benchmark and approximate angular wavenumbers is confined to radial mode functions and coefficients. Figure 6(b) shows excellent agreement in the whispering gallery region (including the highest peak) and everywhere in the wave regions, along with less than 1 dB difference elsewhere in the duct. It should be mentioned that wavenumber approximations which are even a little less accurate than those from Fig. 4(b) can significantly change incoherent TL comparisons. Overall, coherent loss predictions from this paper's approach are very good.

Figure 6(c) represents a benchmark test for the well approximation, and compares results from Figs. 5(c) and 5(d) along the line *y* = 59.4 km. Transmission loss from a PE computation using finite wave width of *w* = 300 m (solid, green online, labeled “PE Exact”) is compared with TL from another treating wave width as infinite (dashed, purple, “PE Well”). The two curves in Fig. 6(c) have excellent agreement in the duct, except near the *r* =* D* interface. This discrepancy arises from energy that is reflected into the duct from the inner wave, as is visible in Fig. 5(d), and also affects TL in the inner wave region. The outer wave regions have pattern oscillations without the well approximation, whereas the solution using the well approximation tracks closely the mean of those oscillations. An important observation is that both PE results in the duct are very similar to Fig. 6(a); the solid and dashed curves in Figs. 6(a) and 6(c) almost share ink. This observation has several implications: first, modes indeed represent a valid high-accuracy solution approach for this NIW model; second, the well approximation has little effect on in-duct propagation, and often provides good estimates outside of the duct region; and third, the wavenumber prediction methods outlined in Sec. IV are accurate enough to be useful in comparison with two high-accuracy modeling approximations.

## VI. SUMMARY AND CONCLUSIONS

The 3-D shallow-water nonlinear internal wave duct model of Lin *et al.*^{1} provides a computational approach for examining how acoustic quantities in modal propagation depend on the environmental feature parameters. A theoretical and analytical approach centered on mode number conservation, as applied to a 3-D shelf-slope model,^{2,3} is shown to be useful for formulating the dependence of sound propagation on nonlinear internal wave parameters. A previous normal mode classification, which specified whispering gallery, fully bouncing, and leaky modes, is extended to distinguish via turning point locations three categories of leaky modes. The complete classification scheme shows the regions of oscillation and phase changes for each possible mode type, leading to five distinct formulas for mode number conservation.

From these five formulas, angular wavenumber responses to changes in NIW parameters can be predicted efficiently and without using any information from the complicated modal dispersion relation. Example plots are shown for wavenumber dependence on parameters specifying radius of wavefront curvature and duct width. Wavenumber approximations from mode conservation are compared to those determined numerically from the dispersion relation, and the accuracy agreement is excellent. The agreement is particularly important for whispering gallery and fully bouncing modes that are responsible for long-distance acoustic propagation in the duct. Wavenumber approximations remain accurate even when parameter changes result in transitions between different mode types.

Because of the wavenumber approximation accuracy, it is feasible to perform sensitive computations including transmission loss. Furthermore, benchmark tests against PE solutions show this method can be used in conjunction with model simplifications such as the well approximation. An example of a horizontal-plane coherent loss using the approximate wavenumbers preserves all the main features of the numerically exact modal pressure field, particularly for long-distance propagation patterns. A cross-sectional plot shows only minor differences in pattern phase and amplitude. As expected, incoherent loss pattern agreement is excellent. The efficiency and accuracy of these approximations suggest usefulness in other applications. They can provide computational efficiency for insight into environmental parameter influences on acoustic field quantities. In addition, mode number conservation formulas can quantify the sensitivity of propagation predictions to model parameter changes, which are proxies for uncertainty in field data.

Future work can further explore influences of modeling assumptions in this paper, by avoiding the well approximation and by including more NIW wavefronts. Taking advantage of mode number conservation permits handling of additional wave-duct interfaces for such purposes. The approach can also be extended to models for other features, including eddies and variable bathymetry. The approach and results of this paper provide tools for investigating the dependence and sensitivity of modal field properties on 3-D physical oceanographic processes.

## ACKNOWLEDGMENTS

This work was supported by the Office of Naval Research under grants to Rensselaer Polytechnic Institute (Grant Nos. N00014-14-1-0372 and N00014-17-1-2370) and to Woods Hole Oceanographic Institution (Grant Nos. N00014-11-1-0701 and N00014-17-1-2692). Additional funding was provided by Naval Undersea Warfare Center Division Newport through the SMART Scholarship for the first author's doctoral degree program. The authors also thank Dr. Timothy F. Duda of WHOI and Dr. David Wells of University of North Carolina at Chapel Hill for their assistance with this paper.

## References

*GSA Memoirs*(