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.

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 system2 and develop convenient formulas for extracting wavenumber parameter dependence3 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 trains11 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 Weston15 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.

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 shown17 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 cW(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 hI and hII. The ocean surface is assumed to be pressure-release, and water has density ρW = 1000 kg/m3 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 > ρW, sound speed cB > cW, and attenuation αB.

FIG. 1.

Idealized model of two nonlinear internal wavefronts and duct between them. Upper shaded layer with constant sound speed c1 > c2 in lower unshaded layer. Duct region denoted by index I, wave region by II. Parameter values in ocean bottom are constant. Star indicates source location (rs, 0, zs). Adapted from Fig. 1 in Ref. 1.

FIG. 1.

Idealized model of two nonlinear internal wavefronts and duct between them. Upper shaded layer with constant sound speed c1 > c2 in lower unshaded layer. Duct region denoted by index I, wave region by II. Parameter values in ocean bottom are constant. Star indicates source location (rs, 0, zs). Adapted from Fig. 1 in Ref. 1.

Close modal

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

1rr(rPr)+1r22Pθ2+ρ(z)z(1ρ(z)Pz)+ω2c2(r,z)P=4πδ(rrs)rδ(θ)δ(zzs).
(1)

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,θ,z,t)=P(r,θ,z)eiωt, with time dependence from a single-frequency point source factored out of Eq. (1). The source is located in the duct at (r,θ,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

ρ(z)ddz(1ρ(z)dψmdz)+(k2(r,z)ζm2(r))ψm=δ(zzs),
(2a)
1rr(rAmr)+1r22Amθ2+ζm2(r)Am=4πδ(rrs)rδ(θ),
(2b)

where ψm=ψm(z),Am=Am(r,θ), and the medium wavenumber is given by k(r,z)=ω/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 ζ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

ζm(r)={ζm,I,r1<r<r2,ζm,II,r<r1,r>r2.
(3)

Subscripts k = 1, 2 identify radial interface locations r1 = D and r2 = 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 rk notation is convenient.

FIG. 2.

(Color online.) (a) Cross-section θ = constant of Fig. 1, showing an upper water layer (light blue, higher sound speed) over a lower (dark blue, lower sound speed). Ocean layers overlay attenuating bottom sediment (brown). Reference values indicated for sound speed, ocean depth, water layer depths, and bottom parameters. Well approximation (Refs. 7, 18) indicated by solid lines, wave profile by dashed; interior and exterior j = I regions (far left and right) replaced with sound speed of j = II region. (b) First unnormalized vertical mode m = 1 in the duct j = I (dashed blue) and wave j = II (solid red) regions, for source frequency f = 75 Hz. Thin horizontal lines indicate top water layer depth in duct (dashed blue) and wave (solid red). (c) Same as (b), but for m = 2. Adapted from Fig. 7 in Ref. 1.

FIG. 2.

(Color online.) (a) Cross-section θ = constant of Fig. 1, showing an upper water layer (light blue, higher sound speed) over a lower (dark blue, lower sound speed). Ocean layers overlay attenuating bottom sediment (brown). Reference values indicated for sound speed, ocean depth, water layer depths, and bottom parameters. Well approximation (Refs. 7, 18) indicated by solid lines, wave profile by dashed; interior and exterior j = I regions (far left and right) replaced with sound speed of j = II region. (b) First unnormalized vertical mode m = 1 in the duct j = I (dashed blue) and wave j = II (solid red) regions, for source frequency f = 75 Hz. Thin horizontal lines indicate top water layer depth in duct (dashed blue) and wave (solid red). (c) Same as (b), but for m = 2. Adapted from Fig. 7 in Ref. 1.

Close modal

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 ζ2,I=0.3092+5.3803×105i and ζ2,II=0.3056+6.0525×105i. 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).

Equation (2b) can be solved using a wavenumber integration technique.21 The solutions are

Am(r,θ)=12πG(r,η)eiηθdη,
(4)

where G solves

1rddr(rdGdr)+(ζm2(r)η2r2)G=4πδ(rrs)r.
(5)

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 ηmn is an angular wavenumber associated with radial mode 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 field22 and are neglected here.

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

k2=kzm2+ζm2=kzm2+krmn2+ηmn2r2,
(6)

where kzm is the vertical wavenumber of vertical mode m, and krmn is the radial wavenumber of horizontal mode n associated with vertical mode m. With approximate expressions for solutions of Eqs. (2a) and (2b), the modal pressure solution is constructed as

P(r,θ,z)=mnγmn(rs,zs)Gmn(r)eiηmnθψm(z),
(7)

where Gmn(r)=G(r,ηmn). Mode shape functions Gmn(r) and ψm(z) are normalized using procedures described in Refs. 1 and 23, and γmn(rs,zs)=iGmn(rs)ψ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.

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

G(r)={A0Hη(1)(ζIIr)+A0Hη(2)(ζIIr),r<r1,A1Hη(1)(ζIr)+B1Hη(2)(ζIr),r1<r<r2,A2Hη(1)(ζIIr),r>r2.
(8)

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 > r2 by excluding an incoming Hankel function. Boundedness for small r is ensured by combining outgoing and incoming waves for r < r1 to form a Bessel function of the first kind using the formula 2Jη=Hη(1)+Hη(2). The four remaining coefficients are found using continuity conditions of G and its first derivative at r = r1 and r = r2, 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=(ζ2η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

νjk=Hη(1)(ζjrk)Hη(2)(ζjrk),νjk()=Hη()(ζjrk)Hη()(ζjrk),
(9)

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

νII1νI1(ζIIνII1(1)ζIνI1(1))(ζIIνII2(1)ζIνI2(2))νI2(ζIIνII1(1)ζIνI1(2))(ζIIνII2(1)ζIνI2(1))νI1(ζIIνII1(2)ζIνI1(1))(ζIIνII2(1)ζIνI2(2))νI2(ζIIνII1(2)ζIνI1(2))(ζIIνII2(1)ζIνI2(1))=1.
(10)

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, A0) 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 A0 = 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

f(rk,η)g(rk,η)U(rk,η)+dU(rk,η)dr=0,
(11)

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

f(r1,ηn)g(r1,ηn)=ζIIJηn+1(ζIIr1)Jηn(ζIIr1)ηnr1,
(12a)
f(r2,ηn)g(r2,ηn)=ζIIHηn+1(1)(ζIIr2)Hηn(1)(ζIIr2)ηnr2,
(12b)

which are found using the appropriate region II solutions. As shown in the Appendix of Ref. 1, any of the mode coefficients A0 is given by

A02=2ηnr1r2Un2(r,ηn)rdrr1ηf(r1,η)g(r1,η)|ηnUn2(r1,ηn)+r2ηf(r2,η)g(r2,η)|ηnUn2(r2,ηn).
(13)

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 = A0U with A0 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 Gmn in Eq. (5) is readily shown from Eq. (6) to be krj2. The standard substitution Gmn(r)=r1/2Qmn(r) reduces Eq. (5) with source term omitted approximately to

d2Qdr2+krj2Q=0,
(14)

because 1/4r2 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 model1 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.

FIG. 3.

(Color online.) (a) Example radial mode types: whispering gallery (WG), fully bouncing (FB), and leaky (L-WG, L-B, L-T). Turning points (TPs) indicated by purple squares. (b) Schematic showing principal geometric features (oscillatory or exponential) of radial mode types, with vertical and horizontal axes as scaled azimuthal wavenumber η/D and radius r/D. Superscript denotes real part of complex wavenumbers; mode numbers m, n omitted on all quantities. On horizontal axis 1 and X = 1 + W/D are wave-duct interfaces; on right vertical axis, existence regions of four mode types shown for case XζII<ζI. The L-T mode region continues below horizontal axis. Sloped lines of radial wavenumbers krj=(ζj2η2/r2)1/2=0 show TP locations. (c) Same as (b), except for case ζI<XζII. Here L-WG modes occur in same relative position as FB modes in (b).

FIG. 3.

(Color online.) (a) Example radial mode types: whispering gallery (WG), fully bouncing (FB), and leaky (L-WG, L-B, L-T). Turning points (TPs) indicated by purple squares. (b) Schematic showing principal geometric features (oscillatory or exponential) of radial mode types, with vertical and horizontal axes as scaled azimuthal wavenumber η/D and radius r/D. Superscript denotes real part of complex wavenumbers; mode numbers m, n omitted on all quantities. On horizontal axis 1 and X = 1 + W/D are wave-duct interfaces; on right vertical axis, existence regions of four mode types shown for case XζII<ζI. The L-T mode region continues below horizontal axis. Sloped lines of radial wavenumbers krj=(ζj2η2/r2)1/2=0 show TP locations. (c) Same as (b), except for case ζI<XζII. Here L-WG modes occur in same relative position as FB modes in (b).

Close modal

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=ζ2η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 ζI<η/D<Xζ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 η<XζII oscillations leak across the interface at D + W, and thus the value XζII is a cutoff for leaky modes. Because XζII<ζ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 ζI<Xζ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 ζI and Xζ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 r2/r1), 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 104 for parameter values of interest.

As noted, this mode type classification extends previous work1 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.

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 integrals14 and were previously applied to vertical waveguides.15,16,25,26 This concept was later extended3 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 rr2. 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 kr, using Eq. (6):

τnD+Wζ2(r)ηn2r2dr=(n1Φn)π,
(15)

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),

ζ2ηn2r2dr=ηn[ζ2r2ηn21cos1(ηnζr)].
(16)

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

μ(ξ)=ξ21cos1(1ξ),
(17a)
μn,jk=μ(ζjrkηn).
(17b)

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 μ(ζjτn/η)=μ(1)=0. Thus, the only lower-limit contributions are from TPs at the inshore interface for FB and L-B modes, which have μ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 μI1 and μII1. All five mode types contain a term μ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 μn contributions from column three for any mode type of number n as Σμn. Then including the phase correction Φn, dividing by π, and rearranging terms, Eq. (15) is simply expressed (terms in cycles) as

n1=ηnπμn+Φn.
(18)

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,

R=krIkrIIkrI+krII.
(19)

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 ϕ = –ϕk value in cycles, where

ϕk=1πtan1(η2/rk2ζII2ζI2η2/rk2).
(20)

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 = r2 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 = r1 interface, but the Rayleigh reflection coefficient evaluates to zero there. The resulting total phase corrections Φn are summarized in column four of Table I.

TABLE I.

Definitions of quantities in Eqs. (15) and (18) for each mode type. Second column τn is turning point location for lower integration limit. Third column is an analytic expression obtained by integrating left-hand side of Eq. (15) and using Eqs. (16), (17a), and (17b). Fourth column is sum Φn of phase changes occurring at turning points that occur with mode oscillations in-duct, in-wave, or at in-shore and off-shore interfaces, with minus sign common to all terms factored out. Mode number subscripts suppressed in table entries.

Mode typeτnμn−Φn
WG η/ζI μI 2 1/4 + ϕ2 
FB D μI2μI1 ϕ1 + ϕ2 
L-WG η/ζI μI 2 1/4 
L-B D μI2μI1 ϕ1 
L-T η/ζII μI2μI1+μII1 1/4 
Mode typeτnμn−Φn
WG η/ζI μI 2 1/4 + ϕ2 
FB D μI2μI1 ϕ1 + ϕ2 
L-WG η/ζI μI 2 1/4 
L-B D μI2μI1 ϕ1 
L-T η/ζII μI2μI1+μ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 ηn 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.3 

FIG. 4.

(Color online.) Parameter dependence of scaled azimuthal wavenumber ηmn/D for vertical mode m = 2. Numerical solutions of dispersion relation shown as solid shapes, with each mode type identified by shape and background shading. Approximation curves generated using solutions of Eq. (18). See text discussion for definitions of horizontal lines (heavy solid) and sloped curves (light solid). Mode numbers suppressed. (a) Plot versus radius of curvature D from 25 to 75 km, with reference value D = 50 km, for radial modes n = 1–10 (top–bottom). (b) Plot versus duct width W from 300 to 700 m, with reference value W = 500 m, for radial modes n = 1–8 (top–bottom). Approximation accuracy excellent for nearly all cases and mode types, even as they change type.

FIG. 4.

(Color online.) Parameter dependence of scaled azimuthal wavenumber ηmn/D for vertical mode m = 2. Numerical solutions of dispersion relation shown as solid shapes, with each mode type identified by shape and background shading. Approximation curves generated using solutions of Eq. (18). See text discussion for definitions of horizontal lines (heavy solid) and sloped curves (light solid). Mode numbers suppressed. (a) Plot versus radius of curvature D from 25 to 75 km, with reference value D = 50 km, for radial modes n = 1–10 (top–bottom). (b) Plot versus duct width W from 300 to 700 m, with reference value W = 500 m, for radial modes n = 1–8 (top–bottom). Approximation accuracy excellent for nearly all cases and mode types, even as they change type.

Close modal

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.

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

Pm(r,θ,z)=iψm(zs)ψm(z)nGmn(rs)Gmn(r)eiηmnθ,
(21)

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

TLm=20log10|Pm(r,θ,z)Pm(rs,0,zs)|.
(22)

For a horizontal cross-section at depth z = zs, Eq. (22) can be expressed as

TLm=20log10|Am(r,θ)Am(rs,0)|,
(23)

where

Am(r,θ)=inGmn(rs)Gmn(r)eiηmnθ
(24)

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,θ,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 rs = 60.425 km is selected, and a specific value for zs does not appear in horizontal TL Eqs. (23) and (24). Next, approximate values for Re(η) 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.

FIG. 5.

(Color online.) Horizontal TL comparisons in (x, y) for D = 60 km, W = 500 m, and reference values of other parameters. Point source (blue circle) at (0, 60.425), f = 75 Hz, and m = 2. (a) Numerically exact (benchmark) computations of η from dispersion relation Eq. (10) using ten radial modes n = 1–10 (four WG, two FB, three L-B, one L-T). (b) Same as (a), with approximate values of Re(η) as in Fig. 4(a). (c) PE computation using same environment as (a) and (b), including well approximation with a duct between two NIW wavefronts [solid line, Fig. 2(a)]. (d) Same as (c) but without well approximation [dashed line, Fig. 2(a)] using wave width w = 300 m. In all subfigures, larger rectangular box magnifies region in smaller rectangular box to highlight details of TL field. Key pattern features of all four subfigures correspond very well in the duct region, with disagreements mainly in wave regions of Fig. 5(d) because well approximation is not used.

FIG. 5.

(Color online.) Horizontal TL comparisons in (x, y) for D = 60 km, W = 500 m, and reference values of other parameters. Point source (blue circle) at (0, 60.425), f = 75 Hz, and m = 2. (a) Numerically exact (benchmark) computations of η from dispersion relation Eq. (10) using ten radial modes n = 1–10 (four WG, two FB, three L-B, one L-T). (b) Same as (a), with approximate values of Re(η) as in Fig. 4(a). (c) PE computation using same environment as (a) and (b), including well approximation with a duct between two NIW wavefronts [solid line, Fig. 2(a)]. (d) Same as (c) but without well approximation [dashed line, Fig. 2(a)] using wave width w = 300 m. In all subfigures, larger rectangular box magnifies region in smaller rectangular box to highlight details of TL field. Key pattern features of all four subfigures correspond very well in the duct region, with disagreements mainly in wave regions of Fig. 5(d) because well approximation is not used.

Close modal

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 Dw<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.

FIG. 6.

(Color online.) Transmission loss curves, from cross sections along horizontal line y = 59.4 km shown (solid, blue) in Fig. 5 subfigures. Vertical lines are wave-duct interfaces. (a) Coherent TL curves: solid (blue, labeled “Numerically Exact”) is from Fig. 5(a), and dashed (red, labeled “Approximate”) is from Fig. 5(b). Both coherent loss curves share close similarities in phase and amplitude of interference patterns, such as in the whispering gallery region near the offshore interface. Some differences occur, as from the inshore interface to the middle of the duct. (b) Incoherent TL curves: same as (a) but calculated using Eq. (25). Incoherent loss results from both curves agree closely in the duct and wave regions shown. (c) High-accuracy TL curves from PE calculations: dashed (purple, labeled “PE Well”) is from Fig. 5(c), and solid (green, labeled “PE Exact”) is from Fig. 5(d). PE calculations without well approximation show excellent agreement with normal mode solutions, except in two wave regions adjacent to duct interfaces.

FIG. 6.

(Color online.) Transmission loss curves, from cross sections along horizontal line y = 59.4 km shown (solid, blue) in Fig. 5 subfigures. Vertical lines are wave-duct interfaces. (a) Coherent TL curves: solid (blue, labeled “Numerically Exact”) is from Fig. 5(a), and dashed (red, labeled “Approximate”) is from Fig. 5(b). Both coherent loss curves share close similarities in phase and amplitude of interference patterns, such as in the whispering gallery region near the offshore interface. Some differences occur, as from the inshore interface to the middle of the duct. (b) Incoherent TL curves: same as (a) but calculated using Eq. (25). Incoherent loss results from both curves agree closely in the duct and wave regions shown. (c) High-accuracy TL curves from PE calculations: dashed (purple, labeled “PE Well”) is from Fig. 5(c), and solid (green, labeled “PE Exact”) is from Fig. 5(d). PE calculations without well approximation show excellent agreement with normal mode solutions, except in two wave regions adjacent to duct interfaces.

Close modal

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,θ) given instead by

Am(r,θ)=n|Gmn(rs)Gmn(r)eiηmnθ|2.
(25)

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.

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.

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.

1.
Y.-T.
Lin
,
K. G.
McMahon
,
J. F.
Lynch
, and
W. L.
Siegmann
, “
Horizontal ducting of sound by curved nonlinear internal gravity waves in the continental shelf areas
,”
J. Acoust. Soc. Am.
133
(
1
),
37
49
(
2013
).
2.
B. J.
DeCourcy
,
Y.-T.
Lin
, and
W. L.
Siegmann
, “
Approximate formulas and physical interpretations for horizontal acoustic modes in a shelf-slope front model
,”
J. Acoust. Soc. Am.
140
(
1
),
EL20
EL25
(
2016
).
3.
B. J.
DeCourcy
,
Y.-T.
Lin
, and
W. L.
Siegmann
, “
Estimating the parameter sensitivity of acoustic mode quantities for an idealized shelf-slope front
,”
J. Acoust. Soc. Am.
143
(
2
),
706
715
(
2018
).
4.
J. R.
Apel
,
L. A.
Ostrovsky
,
Y. A.
Stepanyants
, and
J. F.
Lynch
, “
Internal solitons in the ocean and their effect on underwater sound
,”
J. Acoust. Soc. Am.
121
(
2
),
695
722
(
2007
).
5.
D. S.
Ahluwalia
,
J. B.
Keller
, and
B. J.
Matkowsky
, “
Asymptotic theory of propagation in curved and nonuniform waveguides
,”
J. Acoust. Soc. Am.
55
(
1
),
7
12
(
1974
).
6.
R. H.
Headrick
,
J. F.
Lynch
,
J. N.
Kemp
,
A. E.
Newhall
,
K.
von der Heydt
,
J.
Apel
,
M.
Badiey
,
C.-S.
Chiu
,
S.
Finette
,
M.
Orr
,
B.
Pasewark
,
A.
Turgot
,
S.
Wolf
, and
D.
Tielbuerger
, “
Acoustic normal mode fluctuation statistics in the 1995 SWARM internal wave scattering experiment
,”
J. Acoust. Soc. Am.
107
(
1
),
201
220
(
2000
).
7.
S.
Finette
and
R.
Oba
, “
Horizontal array beamforming in an azimuthally anisotropic internal wave field
,”
J. Acoust. Soc. Am.
114
(
1
),
131
144
(
2003
).
8.
S. D.
Frank
,
M.
Badiey
,
J. F.
Lynch
, and
W. L.
Siegmann
, “
Experimental evidence of three-dimensional acoustic propagation caused by nonlinear internal waves
,”
J. Acoust. Soc. Am.
118
(
2
),
723
734
(
2005
).
9.
R.
Oba
and
S.
Finette
, “
Acoustic propagation through anisotropic internal wave fields: Transmission loss, cross-range coherence, and horizontal refraction
,”
J. Acoust. Soc. Am.
111
(
2
),
769
784
(
2002
).
10.
T. F.
Duda
,
Y.-T.
Lin
, and
D. B.
Reeder
, “
Observationally constrained modeling of sound in curved ocean internal waves: Examination of deep ducting and surface ducting at short range
,”
J. Acoust. Soc. Am.
130
(
3
),
1173
1187
(
2011
).
11.
B. G.
Katsnel'son
and
S. A.
Pereselkov
, “
Low-frequency horizontal acoustic refraction caused by internal wave solitons in a shallow sea
,”
Acoust. Phys.
46
(
6
),
684
691
(
2000
).
12.
K. G.
McMahon
,
L. K.
Reilly-Raska
,
W. L.
Siegmann
,
J. F.
Lynch
, and
T. F.
Duda
, “
Horizontal Lloyd mirror patterns from straight and curved nonlinear internal waves
,”
J. Acoust. Soc. Am.
131
(
2
),
1689
1700
(
2012
).
13.
J. F.
Lynch
,
Y.-T.
Lin
,
T. F.
Duda
, and
A.
Newhall
, “
Acoustic ducting, reflection, refraction, and dispersion by curved nonlinear internal waves in shallow water
,”
IEEE J. Ocean Eng.
35
(
1
),
12
27
(
2010
).
14.
D. M.
Milder
, “
Ray and wave invariants for SOFAR channel propagation
,”
J. Acoust. Soc. Am.
46
(
5B
),
1259
1263
(
1969
).
15.
D. E.
Weston
, “
Guided propagation in a slowly varying medium
,”
Proc. Phys. Soc.
73
(
3
),
365
384
(
1959
).
16.
A. D.
Pierce
, “
Extension of the method of normal modes to sound propagation in an almost stratified medium
,”
J. Acoust. Soc. Am.
37
(
1
),
19
27
(
1965
).
17.
J. C.
Preisig
and
T. F.
Duda
, “
Coupled acoustic mode propagation through continental-shelf internal solitary waves
,”
IEEE J. Ocean Eng.
22
(
2
),
256
269
(
1997
).
18.
Y.-T.
Lin
,
T. F.
Duda
, and
J. F.
Lynch
, “
Acoustic mode radiation from the termination of a truncated nonlinear internal gravity wave duct in a shallow ocean area
,”
J. Acoust. Soc. Am.
126
(
4
),
1752
1765
(
2009
).
19.
C. L.
Pekeris
, “
Theory of propagation of explosive sound in shallow water
,” in
Propagation of Sound in the Ocean
, edited by
M.
Ewing
,
J. L.
Worzel
, and
C. L.
Pekeris
, Vol.
2
of GSA Memoirs (
Geological Society of America
,
New York
,
1948
).
20.
F. B.
Jensen
,
W. A.
Kuperman
,
M. B.
Porter
, and
H.
Schmidt
,
Computational Ocean Acoustics
, 2nd ed. (
Springer
,
New York
,
2011
), pp.
112
126
, 269–276.
21.
G. V.
Frisk
,
Ocean and Seabed Acoustics
(
Prentice-Hall
,
Englewood Cliffs, NJ
,
1994
), pp.
33
, 38.
22.
D. C.
Stickler
, “
Normal-mode program with both the discrete and branch line contributions
,”
J. Acoust. Soc. Am.
57
(
4
),
856
861
(
1975
).
23.
M. B.
Porter
,
The KRAKEN Normal Mode Program
,
SACLANT Undersea Research Centre, currently Centre for Maritime Research and Experimentation
,
La Spezia, Italy
(
1991
).
24.
B. J.
DeCourcy
,
Y.-T.
Lin
, and
W. L.
Siegmann
, “
Effects of front width on acoustic ducting by a continuous curved front over a sloping bottom
,”
J. Acoust. Soc. Am.
146
(
3
),
1923
1933
(
2019
).
25.
C. H.
Harrison
, “
A relation between multipath group velocity, mode number, and ray cycle distance
,”
J. Acoust. Soc. Am.
132
(
1
),
48
55
(
2012
).
26.
I.
Tolstoy
and
C. S.
Clay
,
Ocean Acoustics: Theory and Experiment in Underwater Sound
(
AIP
,
New York
,
1987
), pp.
16
20
, 77–78.
27.
I.
Tolstoy
, “
Phase changes and pulse deformation in acoustics
,”
J. Acoust. Soc. Am.
44
(
3
),
675
683
(
1968
).
28.
M. H.
Holmes
,
Introduction to Perturbation Methods, Texts in Applied Mathematics
(
Springer
,
New York
,
2013
), pp.
236
243
.