A depth-dependent formula for shallow water propagation.

A depth-dependent formula for shallow water propagation € €seyin Ozkan Hu Sertleka) Behavioural Biology, Institute of Biology, Leiden University, Sylviusweg 72, 2333 BE, The Netherlands

Michael A. Ainslieb) TNO, Acoustic and Sonar Group, P.O. Box 96864, 2509 JG, The Hague, The Netherlands

(Received 15 October 2013; revised 9 May 2014; accepted 5 June 2014) In shallow water propagation, the sound field depends on the proximity of the receiver to the sea surface, the seabed, the source depth, and the complementary source depth. While normal mode theory can predict this depth dependence, it can be computationally intensive. In this work, an analytical solution is derived in terms of the Faddeeva function by converting a normal mode sum into an integral based on a hypothetical continuum of modes. For a Pekeris waveguide, this approach provides accurate depth dependent propagation results (especially for the surface decoupling) without requiring complex calculation methods for eigenvalues and corresponding eigenfunctions. C 2014 Acoustical Society of America. [http://dx.doi.org/10.1121/1.4884762] V PACS number(s): 43.30.Bp, 43.30.Dr, 43.20.Bi [ADP]

I. INTRODUCTION

The accurate calculation of propagation loss (PL) is needed for many acoustic problems such as sonar performance modeling, underwater acoustic communication, environmental risk assessment, and oceanography. These applications can require the simulation of propagation characteristics at different ranges and depths. PL can be defined as PL ¼ 10 log10

F1 dB re 1 m2 ; 2 rref

(1)

where F is known as the “transmission factor”1 or “propagation factor”2 and rref is the reference distance, equal to 1 m. PL can be estimated with different methods such as normal modes, ray theory, parabolic equation, and flux theory. Variation of PL with range, water depth, frequency, source depth, and receiver depth can be critical for determining detection ranges, isoclines for environmental risk assessment, etc. Weston introduced an energy flux approach for the calculation of depth-averaged PL without tracing rays or summing normal modes, and thus applicable to large scale problems.3 Weston’s flux equation can be derived from ray and mode theories.4 In Weston’s approach, four propagation regions are described as the spherical, cylindrical spreading, mode-stripping, and single mode regions. While permitting an analytical solution even for a range dependent environment,5 Weston’s flux formulation may not provide satisfactory simulation results when the source or receiver depths are located close to each other, to the seabed, or (especially) to the sea surface. Weston6 addressed this weakness by deriving depth dependence using wave theory concepts,

Pages: 573–582

applying these to a selection of predetermined functional forms for the distribution of energy with angle, making an implicit assumption that the frequency is sufficiently high for the various depth-dependent corrections not to overlap with one another. Harrison7 derived a range- and depthcorrection to Weston’s flux theory that includes ray convergence terms, while neglecting the depth dependence close to source depth, complementary depth, and boundaries. In the present paper, a method is developed to calculate the depth dependence by deriving a generalization of Weston’s approach that does not rely on prior knowledge of the angular energy distribution, and is valid even when the correction terms overlap. The solution, in terms of the Faddeeva function, is derived from an incoherent normal mode sum by using some trigonometric transformations, and is valid for any angular distribution compatible with propagation from a point source in a Pekeris waveguide, and for any sourcereceiver depth combination. Results computed in this way for different frequencies and receiver depths are compared with results using the normal mode program KrakenC.8 For the comparisons, a range independent test case (Scenario A2.I) from the Weston Memorial Workshop9,10 is used. The detailed description of this test case is given in Sec. III. II. DEPTH DEPENDENCE OF PL

The depth dependence of PL can have different characteristics due to the depth of the receiver relative to the sea surface, seabed, source depth, and complementary source depth.6 The complementary source ðzcs Þ or receiver depth ðzcr Þ can be defined for source depth zs (or receiver depth zr ) as zcs;cr ¼ D zs;r ;

(2)

where a)

Author to whom correspondence should be addressed. Also at: Gebze Institute of Technology, Electronics Engineering Department, P.O. Box 141, 41400, Gebze, Turkey. Electronic mail: [email protected] b) Also at: Institute of Sound and Vibration Research, University of Southampton, Highfield, Southampton SO17 1BJ, United Kingdom. J. Acoust. Soc. Am. 136 (2), August 2014

D ¼ h þ DW

(3)

and DW ¼ ðq2 =q1 Þ=k1 sin hc is the vertical wave shift and described by Weston,11,12 h is water depth, q2 is the density

0001-4966/2014/136(2)/573/10/$30.00

C 2014 Acoustical Society of America V

573

Redistribution subject to ASA license or copyright; see http://acousticalsociety.org/content/terms. Download to IP: 140.117.111.1 On: Mon, 18 Aug 2014 09:28:19

of sediment, q1 is the density of sea water, k1 is the wavenumber in the water layer, and hc is the sediment critical angle. When PL is plotted vs zs for a constant receiver depth, similar depth dependent properties are observed around the true receiver depth zr as the complementary receiver depth D zr .

ln jV ðhn Þj ð Þ dn ¼ Im k1 cos hn ¼ : rcM

(11)

The following relations can be obtained with the relation of c1n ¼ np=D in Eqs. (A4) and (A5) as shown in Appendix A: Fð D zr ; zs Þ ¼ Fðzr ; zs Þ ¼ Fðzr ; D zs Þ:

(12)

A. Derivation

By using the incoherent mode sum, the propagation factor F for an isovelocity and isodensity waveguide can be written as 1 w n ðzs Þw n ðzr Þ 2p X expð2dn rÞ: jn r n¼1 2

Fðzr ; zs Þ ¼

2

(4)

This property is useful because if Fðzr ; zs Þ is known for 0 < zr ; zs < D=2, Eq. (30) permits straightforward evaluation of Fðzr ; zs Þ for D=2 < zr ; zs < h. The discrete mode sum ðFÞ can be approximated by the following integral which is denoted by F0 : 2p F0 ðzr ; zs Þ ¼ r

ð1 0

w2n ðzs Þw2n ðzr Þ ð2r=rcM Þ jVj dn : jn

(13)

Here wn ðzÞ is the eigenfunction wn ðzÞ ¼ An sinðc1n zÞ;

(5)

where An is the amplitude of the nth eigenfunction, cn is the vertical wavenumber, and jn is the horizontal wavenumber. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ð Defining the quantity of Ln ¼ ð1= k12 ðzÞ j2n Þdz ¼ rcM =2jn , the normalization constant A2n can be written13 1=2 4 cot hn 2 2 A2n ¼ k1 j2n ¼ ; rcM Ln

(6)

where rcM is the modal cycle distance, and hn is the nth mode grazing angle given by nk ; n ¼ 1; 2; 3; …: (7) hn ¼ arcsin 2D The relation between the modal cycle distance and geometric ray cycle distance (rcG ) is13 rcM ¼ rcG þ 2DB cot h;

(8)

where DB ¼ c21n þ c22n =c1n c22n ½ðq1 c2n =q2 c1n Þþ ðq2 c1n =q1 c2n Þ is the vertical beam shift12 (converted from the corresponding horizontal shift stated by Weston and Tindle14). This relation can be written for the isovelocity case by substituting for rcG ¼ 2hcotðhÞ

By this integral representation, the integer n is replaced by a continuous function of the grazing angle h as n nðhÞ. The variation of nðhÞ versus h can be written as dnðhÞ rcM k1 sinðhÞ: ¼ 2p dh

(14)

However, the integrand is defined only at the discrete eigenvalues. In the following steps, this integrand will be interpolated between the eigenvalues as a continuous function of h, thus permitting evaluation of the integral, F 0 ðzr ; zs Þ ¼

1 r

ð p=2 0

A4n sin2 ðc1n zs Þsin2 ðc1n zr ÞrcM

tanðhÞjVjð2r=rcM Þ dh;

(15)

where An and rcM are given by Eqs. (6) and (9). The integral can be written in different forms by using trigonometric identities. First, using the trigonometric identity sin2 ðc1n zÞ 12 12 cosð2c1n zÞ for eigenfunctions, this integral becomes 4 F0 ðzr ; zs Þ ¼ r

ð p=2 0

ð1 cosð2c1n zs ÞÞð1 cosð2c1n zr ÞÞ rcM tanðhÞ

jVjð2r=rcM Þ dh:

(16)

(9)

Then, changing the variable of integration to continuous grazing angle and using the trigonometric identity

If one uses the relations of c21n ¼ k12 sin2 h and c22n ¼ k12 c21n k22 , the beam shift takes the form q1 q2 k12 k22 ; DB ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi k12 cos2 h k22 q22 k12 sin2 h þ q21 k12 cos2 h q21 k22 (10)

cosð2k1 zs sin hÞcosð2k1 zr sin hÞ cos 2k1 ðzs zr Þsin h þ cosð2k1 ðzs þ zr Þsin hÞ ¼ ; 2 (17)

rcM ¼ 2ðh þ DB Þcot h:

results in the following integral form for the propagation factor: where k2 ¼ x=c2 . The attenuation term in the horizontal direction is added by a perturbation approach with the contribution of a small imaginary part dn to jn .15 This small imaginary part leads to an exponential attenuation with range as expð2dn r Þ; where dn can be related to the reflection coefficient of sea bottom V ðhÞ according to 574

J. Acoust. Soc. Am., Vol. 136, No. 2, August 2014

F 0 ðzr ; zs Þ ¼

4 r

ð p=2 0

ð1 W ðzr ; zs ; hÞÞ ð2r=rcM Þ jVj dh; rcM tanh

(18)

where the term involving source and receiver depths can be written

€ Sertlek and M. A. Ainslie: A depth-dependent formula for shallow water H. O.


W ðzr ; zs ; h Þ ¼ cosð2k1 zs sin hÞ þ cosð2k1 zr sin hÞ cos 2k1 sin hðzs zr Þ þ cosð2k1 sin hðzs þ zr ÞÞ : 2 (19) If W ðzr ; zs ; hÞ were sufficiently close to zero on average, Eq. (18) would have the same form as the flux integral, which has an analytical solution in the form of an error function if the reflection loss increases linearly with angle.2,16 This gives the depth-averaged propagation factor. B. Reflections from the seabed

Equation (18) can be evaluated numerically for any given seabed reflection coefficients. The possibility of analytical evaluation for special forms of the reflection coefficient is explored next.

The vertical beam shift is important for low frequencies and is itself a function of angle.12 Equation (22) can be used for low frequencies by evaluating the integral numerically. The beam shift can be ignored for higher frequencies, in which case the following equation results: 2 F0 ðzr ; zs Þ ¼ rh

ð hc 0

where wð x þ iyÞ is the Faddeeva function17,18 which is defined as ! ð xþiy 2i expðt2 Þdt wðx þ iyÞ ¼ expðð x þ iyÞ2 Þ 1 þ pffiffiffi p 0 (25) and is related to the error function18 according to

wðx þ iyÞ ¼ exp ð x þ iyÞ2 ð1 þ erf ðiðx þ iyÞÞÞ:

ð hc

ð1 W ðzr ; zs ; hÞÞ h þ DB 0 gr h2 dh: exp ðh þ DB Þ

2 F 0 ðzr ; zs Þ ¼ r

gr 2 ð1 W ðzr ; zs ; hÞÞexp h dh; h (23)

cosðZhÞexpðRh2 Þdh 2 Z rffiffiffi exp pffiffiffi p Z 4R w pffiffiffi i Rhc ¼ 4 R 2 R ! pffiffiffi 2 Z pffiffiffi i Rhc exp 2 R pffiffiffi Z w pffiffiffi þ i Rhc 2 R !! pffiffiffi 2 Z pffiffiffi þ i Rhc ; (24) exp 2 R

U½Z; R; hc ¼

One reflection coefficient assumption that permits an analytical evaluation of the integral is exponential reflection coefficient in the form (20) jVj ¼ exp gh2 =tan h ;

This integral can be solved analytically for the isovelocity and isodensity case. For the isovelocity case, ray cycle distance is rcM ¼ 2ðh þ DB Þcot h. Then, the propagation factor becomes

0

where W ðzr ; zs ; hÞ is the sum of cosine functions from Eq. (19). Thus, the analytical solution for the following type of integral leads to an analytical solution for the propagation factor:

1. Exponential reflection coefficient

where g is the rate of increase with grazing angle of the sediment reflection loss,2 in nepers per radian (Np/rad). This form of reflection coefficient simplifies the derivation without requiring small angle approximation. Specifically, the propagation factor becomes ð 4 hc ð1 W ðzr ; zs ; hÞÞ F 0 ðzr ; zs Þ ¼ r 0 r tan h cM 2gr h2 dh: (21) exp rcM tan h

ð hc

(26)

(22)

The analytic expression for the propagation factor for the isovelocity case can be written in terms of the function U½Z; R; hc as

8 2 " # " # rffiffiffiffiffi ! > rffiffiffiffiffi :rh 39 gr gr = U 2k1 ðzs zr Þ; ; hc þ U 2k1 ðzs þ zr Þ; ; hc 5 h h ;; 2



(27)

575


where the first term is the error function based solution of the classic flux integral.2,16 The term in curly parentheses describes the depth dependence for an isovelocity waveguide. In Fig. 1, U½Z; R; hc is plotted as a function of R, a dimensionless range variable equal to gr=h, and a dimensionless depth variable Z that can be any one of 2k1 zr ; 2k1 zs ; 2k1 ðzr þ zs Þ and 2k1 ðzr zs Þ: Specifically, the graph on the left shows the depth dependence of U½Z; R; hc for selected constant values of R, as stated in the figure legend. Similarly, the lower graph shows the range dependence of U½Z; R; hc for selected values of Z. The two-dimensional (2D) contour plot (upper right panel) shows the combined R and Z dependence in a single graph. In all three graphs, the value of hc ¼ arccosð1500=1700Þ ¼ 0:489957 rad is used for consistency with the parameters of the selected scenario (A2.I) from Weston Memorial Workshop 2010. A more detailed description of the test cases are given in Sec. III. According to these figures, we expect more significant variations at low values of Z (low frequencies) and R (at small ranges or deeper water). The function U½Z; R; hc pffiffiffiffiffiffiffiffiffiffiffiffi asymptotically approaches 12 ðp=RÞexp Z 2 =4R for large R [by using Eq. (36)]. For small R, limR!0 U½Z; R; hc ¼ sinðZhc Þ=Z: For large Z, the asymptotic form is R; hc ¼p0. limp Z!1 pffiffiffi ffiffiffi For small Z, limZ!0 U½Z; R; hc ffiffiffi U½Z; ¼ pðerfð Rhc Þ=2 RÞ. When they are both small, limZ;R!0 U½Z; R; hc hc . 2. Harrison’s approximation to Rayleigh reflection coefficient

Use of an exponential reflection coefficient, as shown in Sec. II B 1, provides an analytical solution. This analytical solution [Eq. (27)] is based on a linear approximation to the seabed

reflection loss. As it will be shown in Fig. 2, the Rayleigh reflection coefficient can provide a better approximation for the decay rate.19 Expressing the Rayleigh reflection coefficient in an exponential form can be computationally more efficient because of having less complicated equations in the final form. This exponential form also enables to facilitate comparisons with the previous derivations. A useful approximation to the Rayleigh reflection coefficient that also leads to an analytical expression in a similar form to Eq. (18) is19 0 1 g sin h ! !C; (28) jVj ¼ expB 2 B pffiffiffiffiffiffiffiffiffiffiffi C q2 @ 1 v A 1v 1þ q1 where v ¼ ðsin h=sin hc Þ2 , q1 and q2 are the densities of water and sediment. Thus, the propagation factor can be written as F 0 ðzr ; zs Þ ¼

ð hc

ð1 W ðz r ; zs ; h ÞÞ ðh þ DB Þ 0 ggr sin h tan h dh; exp ðh þ DB Þ 2 r

(29)

pffiffiffiffiffiffiffiffiffiffiffi where g ¼ ð 1 vð1 þ ððq2 =q1 Þ2 1ÞvÞÞ1 . In Sec. III, propagation factors based on exponential and Rayleigh reflection coefficients are compared in Fig. 2. C. Symmetry of depth dependent solution

The approximation F0 is based on a continuous angle assumption for the eigenvalues as c1n ¼ kw sinðhÞ. For the

FIG. 1. (Color online) U½Z; R; hc function is plotted versus v (on the upper left) and versus R (on the bottom right) for critical angle hc ¼ 0.489957 rad, corresponding to the scenario described in Sec. III. 2D representation of U½Z; R; hc versus Z and R (on the upper right).

576




FIG. 2. (Color online) PL vs depth calculated by KrakenC, Weston’s approach [Eq. (33)] for the various energy distributions with angle (low pass and Gaussian) and Flux integral without depth dependent properties [Eq. (34)] for 250 Hz at 1, 5, 25, and 125 km (upper graphs). PL vs depth calculated by KrakenC, Eqs. (23) and (29) for 250 Hz at 1, 5, 25, and 125 km. Beam shift is neglected in these comparisons. For the dashed and dash-dotted lines, PL is approximated using F F0 (lower graphs). The source depth is 30 m.

continuous angle approximation, the symmetry properties of PL are not satisfied [see Eq. (19) and Appendix B]. Therefore, the derived solution can be used to model the features at the sea surface and source depth. This derivation gives poor results for the features at sea bottom and complementary source depth as it could be seen in Fig. 2. By using these symmetry properties from Eqs. (A4) and (A5), the features at sea bottom and complementary source depth can be obtained from the features at sea surface and source depth as 8 < F ðzr ; zs Þ; zs < D=2 (30) F ð zr ; z s Þ ¼ : Fþ ðzr ; zs Þ; zs > D=2; J. Acoust. Soc. Am., Vol. 136, No. 2, August 2014

where F and Fþ are defined as ( F0 ðzr ; zs Þ; F ðzr ; zs Þ ¼ F0 ðD zr ; zs Þ; and

( F þ ðzr ; zs Þ ¼

zr < D=2 zr > D=2

F0 ðzr ; D zs Þ;

zr < D=2

F0 ðD zr ; D zs Þ;

zr > D=2:

(31)

(32)

These terms give the propagation factor ðFÞ by using the symmetry property of source, receiver, complementary source, and complementary receiver depths.


577


III. COMPARISONS

In this section, results obtained using the methods described in Sec. II are compared with KrakenC. In these comparisons, a range independent test case (Scenario A2.I) from the Weston Memorial Workshop10,11 is used, based on Problem XI from the 2006 Reverberation Modeling Workshop at University of Texas at Austin.20 For Scenario A2.I, the water depth (h) is 100 m, source depth (zs ) is 30 m, sound speed in water (c1 ) and sediment (c2 ) are, respectively, are 1500 and 1700 m/s, corresponding to a critical angle hc ¼ 0.489957 rad. The sediment-seawater density ratio (q2 =q1 ) is 2 and resulting reflection loss gradient (g) is 0.273777 Np/rad, evaluated using Eq. (8.86) of Ref. 2. Weston6 derived the depth dependence for a number of specified distributions of energy with angle. For each angular distribution, which must be known in advance, four correction factors to the depth-averaged propagation factor are given of the form i 1 f ðx0 Þ (sea surface and seabed) and h 1 þ 12 f ðx0 Þ (source depth and complementary source depth), where x0 ¼ k1 d/0 , d is the distance to the depth in question and /0 is an angle that characterizes the distribution, equal to the critical angle for Weston’s case 2 p (“low pass”) and to the ffiffiffiffiffiffiffiffiffiffiffiffi effective propagation angle /0 h=2gr for cases 5 (Weston’s “Gaussian”) and 6 (“dipole Gaussian”). In each case, the functional form f(x0) is given by Weston’s Table I. In contrast to Weston’s approach, use of Eq. (23) in terms of the Faddeeva function permits direct computation of the depth-dependent propagation factor without the need for prior knowledge of the angular distribution of energy. Instead, Eqs. (23) and (29) automatically take into account all possible angular distribution cases for a Pekeris waveguide, dealing with transitions between different cases and combinations of

multiple cases in a natural way, and permitting evaluation of all possible combinations of source or receiver depth. In Appendix C, relevant cases from Weston’s table (Gaussian and dipole Gaussian) are derived from Eq. (23). Figure 2 compares the result of evaluating expressions for cases 2 (low pass) and 5 (Gaussian). In the cylindrical spreading region (e.g., 0.5 km), the low pass solution is applicable, whereas for mode stripping (e.g., 25 km), one expects a Gaussian distribution with angle, and these expectations are confirmed by Fig. 2. Less obvious is what happens for intermediate ranges (e.g., 2 km) at which it is necessary to make a choice between the low pass and Gaussian distributions. By contrast, Eq. (23) from the present paper is applicable to both cylindrical spreading and mode stripping regions, including any intermediate range. Further, if greater accuracy is required, this can be achieved by means of Harrison’s approximation to the Rayleigh reflection coefficient, as implemented here in the form of our Eq. (29). For the purpose of Figs. 2 and 3, to combine the effects of multiple depth correction factors from Weston’s paper, we have multiplied these together in the form (denoted F2;5 for Weston’s cases 2 and 5, applicable when the source is not close to the sea surface) f ðx0;zs Þ f ðx0;zcs Þ F2;5 ¼ Fref 1 þ 1þ 2 2 (33) 1 f ðx0;0 Þ 1 f ðx0;D Þ ; where the second subscript of the variable x0 indicates the depth relative to which d is measured, such that x0;f ¼ k1 ðzr fÞ/0 . Here Fref is the analytical solution of Eq. (23) for W ðzs ; zr ; hÞ ¼ 0 as16

FIG. 3. (Color online) PL vs depth calculated by Eq. (23), Eq. (37) Weston’s approach [Eq. (33)] for the different energy distribution with angle (Gaussian dipole Gaussian), for different source depths (1, 3, 10, and 30 m). 578




Fref ¼ r

3=2

rffiffiffiffiffi ! rffiffiffiffiffi p gr erf hc : gh h

(34)

In Eq. (3) of Weston’s paper,6 the low pass and Gaussian cases correspond to an assumption that sin2 ðc1n zs Þ may be replaced by its average value of 1/2. Consequently, Eq. (33) is not valid for source depths close to the sea surface or seabed on a wavelength scale. For small source depths, the depth corrections can be combined as (denoted F6 , for Weston’s case 6, dipole Gaussian–as used in Fig. 3)

F6 ¼ 2ðk1 zs /0 Þ2 Fref 1 f ðx0;0 Þ 1 f ðx0;D Þ :

(35)

A similar form of Eq. (35) for the large receiver depth is also shown by Denham.21 In Fig. 2, the PL results obtained with Eqs. (23) and (29), Weston’s cases [Eq. (33)], and KrakenC are compared at different ranges. For these comparisons, Eqs. (23) and (29) are evaluated numerically without applying the symmetry relations of Eqs. (30)–(32). The beam shift is ignored in Fig. 2. The classic form of the flux integral without depth dependent properties2,16 is also plotted, as obtained by substituting W ðzr ; zs ; hs Þ ¼ 0 in Eq. (23) gives this classical form of flux equation [Eq. (34)]. The features at the complementary depth and seabed are not visible in the graphs of proposed solutions by Eqs. (23) and (29) because the symmetry property of Eq. (30), Sec. II C, is not applied. They can be obtained, if desired, by applying the symmetry property described in Sec. II C. In Fig. 3, Eqs. (33) and (35) are evaluated for the different angle distributions (Gaussian and dipole Gaussian) for different source depths at 250 Hz and 25 km. In these comparisons, the same parameters of Scenario A2.I (described above) are used. The benefit of the Faddeeva function is that it includes all relevant limiting cases in one formula. To consider the large R limit illustrate pffiffiffi this pffiffiffi benefit, Z=2 R Rhc in Eq. (26). In this limit, U½Z; R; hc simplifies to

2 rffiffiffi exp Z p 4R ½ð1 þ erfðpffiffiffi U½Z; R; hc Rhc ÞÞ 4 R pffiffiffi ð1 erfð Rhc ÞÞ rffiffiffi 1 p Z2 ; (36) exp 4R 2 R pffiffiffi where erf Rhc 1 for large R. When Eq. (36) is substituted in Eq. (23), we obtain the following general-purpose formula for the mode stripping region (generalizing Weston’s case 5 and 6), F0 ðzr ; zs Þ ¼ Fref f1 expð2ðk1 zr /0 Þ2 Þ expð2ðk1 zs /0 Þ2 Þ ½1 expð2ðk1 zr /0 Þ2 Þ coshð4zr zs ðk1 /0 Þ2 Þg;

(37)

pffiffiffiffiffiffiffiffiffiffiffiffi where /0 ¼ h=2gr. The transition between cases 5 and 6 is illustrated in Fig. 3 by showing transition from 1 m (case 6) to 30 m (case 5) at 250 Hz. At the intermediate source depth of 10 m, neither case 1 nor case 5 predicts the correct solution. Equations (33) and (35) each satisfy the reciprocity principle for their respective regime of applicability. For instance, Eq. (33) satisfies reciprocity for large receiver depths because of the use of sin2 ðc1n zs Þ 1=2 approximation. Equation (35) satisfies reciprocity for a sufficiently small source and receiver depth. We apply Eq. (23) next (Fig. 4) to investigate the depth dependence for Scenario A2.I from the Weston Memorial Workshop, for frequencies in the range 50 to 3500 Hz, at a range of 5 km. The symmetry approach described by Eq. (30) is used to obtain features close to the complementary depth and the seabed. The expressions for F and Fþ from Eqs. (31) and (32) are evaluated numerically by the Rayleigh reflection coefficient based solution from Eq. (29), with and without beam shifts. There is good agreement with KrakenC, especially when the source or receiver is close to the sea surface. Although the

FIG. 4. (Color online) PL vs depth by KrakenC and Eq. (30) [F0 is calculated by Eq. (29)] for 50 Hz, 250 Hz, 1 kHz, and 3.5 kHz at a range of 5 km from the source. The source depth is 30 m.



579


FIG. 5. (Color online) Relative PL (PLrelative ¼PL PLref ) versus range at different receiver depths (1, 30, and 50 m) for 250 Hz (upper graph). Relative PL versus depth and range for 250 Hz (lower graph). PL is calculated using Eq. (30) (without beam shift) and KrakenC (incoherent mode sum). PLref is the reference PL [Eq. (34)].

results with the effect of beam shift are closer to KrakenC predictions than without beam shift, the difference is small. The proposed solution [Eq. (29)] can provide an accurate approximation to the incoherent normal mode sum, without the need to compute discrete eigenvalues. The classic flux integral without depth dependent properties [Eq. (34)] is selected as a reference PL ðPLref 2 ¼ 10 log10 ðF1 ref =1 m Þ dBÞ, with h ¼ 100 m, g ffi 0.273777 Np/rad and hc ¼ 0:489957 rad. The relative PL ðPLrelative ¼ PL PLref Þ is shown in Fig. 5. PL is calculated using Eq. (30) (with the Rayleigh reflection coefficient, without beam shift) and KrakenC. In Fig. 5, the range variation of the relative PL versus range at receiver depths 1, 30, and 50 m is shown for 250 Hz. Good agreement between Kraken and Eq. (30) is observed. The largest differences between the proposed solution and Kraken at the different receiver depths considered are 0.19 dB (at 1 m), 0.07 dB at (30 m), and 0.15 dB (at 50 m).

complementary source depth is equal to that for a source at the true source depth, and similarly for the receiver (for the isovelocity and isodensity case). Substituting c1n ¼ np=D into 1 A4n sin ðc1n zr Þsin ðc1n zs Þ 2p X expð2dn rÞ j r n¼1 2

Fðzs ; zr Þ ¼

2

(A1) the propagation factor becomes zs zr 2 2 4 sin np np 1 An sin D D 2p X F ðzs ; zr Þ ¼ jn r n¼1 expð2dn r Þ:

(A2)

By using the periodicity property of trigonometric functions as sin2 ðnpðD zr =DÞÞ ¼ sin2 ðnpðzr =DÞÞ, the equation F ðzs ; D zr Þ

IV. CONCLUSIONS

In this paper, the depth dependence of PL is investigated by deriving a depth-dependent correction term [as given by Eq. (19)] to the error function associated with the classical flux solution [Eq. (34)]. This equation is derived from an incoherent normal mode sum, and provides a depth-dependent solution for a lossy Pekeris waveguide without the need for (computationally expensive) calculation of eigenvalues. Equation (23) provides an analytical solution which includes all relevant depth corrections (Weston’s low pass, Gaussian, and dipole Gaussian cases) in one formula. The depth dependence of PL for long ranges can be estimated by a simpler formula Eq. (37). If it is desired, a numerical solution based on Rayleigh reflection coefficient [Eq. (29)] including beam shifts can also be used in order to increase the accuracy. The explicit form of the result provides insight into the behavior of the sound field close to the source depth, complementary source depth, sea surface, and seabed.

zs zr 2 2 4 A sin np np sin 1 n D D 2p X expð2dn r Þ ¼ r n¼1 jn (A3)

is obtained for the complementary receiver depth. The righthand side of this equation is equal to the propagation factor for zr thus proving that the propagation factors at receiver and complementary receiver depths are equal. By repeating the same steps for the complementary source depth as D zs , a similar identity can be obtained. The resulting identities can be written Fðzs ; zr Þ ¼ Fðzs ; D zr Þ;

(A4)

Fðzs ; zr Þ ¼ FðD zs ; zr Þ;

(A5)

and therefore Fðzs ; zr Þ ¼ FðD zs ; D zr Þ:

APPENDIX A: THE SYMMETRY OF PROPAGATION FACTOR BY MODE SUM

The purpose of this appendix is to demonstrate that, if Weston’s vertical wave shift is assumed to be independent of angle, the propagation factor for a source at the 580


(A6)

These symmetry properties can be exploited to avoid the need to calculate PL for all receiver and source depth combinations separately. For computational efficiency, PL for all receiver depths can be calculated by Eq. (30) based on these properties.



APPENDIX B: THE ASYMMETRY OF PROPAGATION FACTOR BY INTEGRAL SOLUTION

In this section, the symmetry of the integral representation is investigated. Figure 2 demonstrates that the integral over a continuum of modes fails to exhibit the symmetry property of the discrete mode sum characterized by Eqs. (A4) and (A5). The purpose of this appendix is to describe

the conditions for which the integral solution breaks down and explain the reason for this failure. The differences between the discrete mode summand and continuous angle integrand are also shown. The W ðzr ; zs ; hÞ term in the integral solution is derived from eigenfunctions, which lead to the symmetrical propagation factor as described in Appendix A. First, the symmetry of W ðzr ; zs ; hÞ is investigated by replacing zr with zr ,

W ð D zr ; zs ; hÞ ¼ cosð2zs k1 sin hÞ þ cosð 2Dk1 sin h 2zr k1 sin hÞ cosð2ðzs þ zr Þk1 sin h 2Dk1 sin hÞ þ cos 2ðzs zr Þk1 sin h þ 2Dk1 sin h : 2

The symmetry property shown by Eqs. (A5) and (A6) come from the periodicity of trigonometric functions. Thus, in order to ensure the symmetry property, 2Dk1 sin h must be equal to 2np, which holds only for the discrete angles corresponding to the eigenvalues. However, this relation is not valid for the non-integer values of n (continuous angles). Thus, having continuous angles in the integration can cause an asymmetric solution. The following equality can only be written for the discrete angles ðhn Þ: cn ¼ k1 sinðhn Þ ¼

np ; D

n ¼ 1; 2; 3; …:

(B2)

Equations (B1) and (18) for continuous angles are different. This means the symmetry relations formulated by Eqs. (A4)–(A6) are not satisfied for the continuous angle approximation. For an illustration of symmetry properties for discrete and continuous angles, the argument of sin2 ðcn zr Þsin2 ðcn zs Þ is plotted for discrete eigenvalues ðcn ¼ np=DÞ and continuous angles ðcn ¼ k1 sin hÞ in Fig. 6. The same comparison is also repeated for the complementary receiver depth as sin2 ðcn ðD zr ÞÞsin2 ðcn zs Þ. The symmetry property for the receiver [Eq. (A4)] requires sin2 ðcn zr Þsin2 ðcn zs Þ ¼ sin2 ðcn ð D zr ÞÞsin2 ðcn zs Þ. The validity of the continuous angle approach is investigated by the comparisons in Fig. 6. For this comparison, f ¼ 250 Hz, zr ¼ 10 m, D ¼ 104:06 m, and hc ¼ 0:49 rad. This comparison shows that the discrete mode solutions at receiver depth and complementary receiver depth are identical when calculated using Weston’s wave shift

(B1)

concept, if the vertical wave shift is assumed independent of grazing angle. According to Eq. (B2), the discrete angles are calculated as hn ¼ arcsinðnc1 =2Df Þ. The continuous angles can be any value between zero and hc , irrespective of the eigenvalues. In order for the continuous angle integral to approximate the discrete mode sum, the area under the curves in Fig. 6 must be equal to the sum of contributions from discrete eigenvalues, which is the case if the curve contains no energy above the Nyquist frequency corresponding to the rate at which the curves are sampled by the series of eigenvalues. Writing sinððnp=DÞzs Þsinððnp=DÞzr Þ ¼ ð1= 2iÞðexpðiðnp=DÞzs Þ expð iðnp=DÞ zs ÞÞð1 = 2iÞðexpðiðnp= DÞzr Þ expðiðnp=DÞzr ÞÞ, the highest frequency component is the term proportional to expð6i2npFÞ, where the highest frequency, F, is equal to ðzs þ zr =2DÞ. According to the Nyquist-Shannon theorem, the continuous function is well sampled, in the sense that it can be reconstructed from these samples, by the discrete modes if F is less than the Nyquist frequency (half the sampling rate), which is equal to 0.5. While it does not follow directly from this theorem, it seems reasonable to assume that the integral of a continuous function will be a good approximation to a discrete sum if the maximum frequency is below the Nyquist frequency corresponding to the sampling rate, and not otherwise. If the continuous function is sufficiently well sampled in the same sense, one can therefore expect the integral to exhibit the same symmetry property as the discrete sum. Given that the maximum frequency is F ¼ zs þ zr =2D and the Nyquist frequency is 0.5, it follows that the condition for symmetry is

FIG. 6. (Color online) Variation of sin2 ðcn zr Þsin2 ðcn zs Þ for discrete angles at eigenvalues (vertical bars) and continuous angles (solid and dashed lines). The source depth is 30 m.



581


zs þ zr < D. This condition shows that the proposed solution is valid up to the complementary source depth ðzr < zcs Þ. Thus, the proposed solution by the continuum approach is valid for zr < zcs 2r, where r is the width of the peak at the complementary source depth. In the mode stripping, this width is given by the simple formula as r ¼ 1=2k1 /0 . In the cylindrical spreading region, the width is less well defined. Even though the symmetry properties are not satisfied for the integral representation at the complementary source and receiver depths, the correct propagation factor can be obtained by exploiting the symmetry of Eq. (30) for D=2 < zr < D. APPENDIX C: DERIVATIONS OF WESTON’S CASES 5 AND 6 FROM THE PROPOSED SOLUTION

In this appendix, Weston’s case 5 (Gaussian) and case 6 (dipole Gaussian) are derived from Eq. (37). 1. Derivation of case 5

When zs is large enough and far from the waveguide’s boundaries (k1 zs /0 1, such that expð2ðk1 zs /0 Þ2 Þ) may be neglected, the propagation factor can be written as F0 ðzr ; zs Þ Fref ð1 expð2ðk1 zr /0 Þ2 ÞÞ

(C1)

consistent with Weston’s case 5. Equation (C1) can also be derived from Ref. 21 by applying the reciprocity principle to Eq. (25) from that paper. 2. Derivation of case 6

When zs is small and close to sea surface, the propagation factor can be written as F0 ðzr ; zs Þ ¼ Fref f1 expð2e2 Þ expð2ðk1 zr /0 Þ2 Þ ½1 expð2e2 Þcoshð4k1 zr /0 eÞg;

(C2)

where e ¼ k1 zs /0 . For small e, the exponential and hyperbolic cosine function can be replaced with expð2e2 Þ 1 2e2 and coshð4ek1 zr /0 Þ 1 þ ð4ek1 zr /0 Þ2 =2 ¼ 1 þ8ðek1 zr /0 Þ2 , F0 ðzr ; zs Þ ¼ Fref ð2e2 expð2ðk1 zr /0 Þ2 Þ ½1 ð1 2e2 Þð1 þ 8ðek1 zr /0 Þ2 ÞÞ: (C3) Then, neglecting terms of order e4 or higher

582


F0 ðzr ; zs Þ ¼ 2e2 ð1 ð1 4ðk1 zr /0 Þ2 Þ expð2ðk1 zr /0 Þ2 ÞÞFref ;

(C4)

consistent with Weston’s case 6. For large receiver depths, this expression may be approximated using expð2ðk1 zr /0 Þ2 Þ 0; in which case Eq. (27) from Ref. 21 is recovered. 1

D. E. Weston, “Acoustic flux methods for oceanic waveguides,” J. Acoust. Soc. Am. 68(1), 287–296 (1980). 2 M. A. Ainslie, Principles of Sonar Performance Modeling (Praxis Publishing, Chichester, UK, 2010). 3 D. E. Weston, “Intensity-range relations in oceanographic acoustics,” J. Sound Vib. 18(2), 271–287 (1971). 4 C. H. Harrison and M. A. Ainslie, “Fixed time versus fixed range reverberation calculation: Analytical solution,” J. Acoust. Soc. Am. 128, 28–38 (2010). 5 D. E. Weston, “Propagation in water with uniform sound velocity but variable depth lossy bottom,” J. Sound Vib. 47, 473–483 (1976). 6 D. E. Weston, “Wave theory peaks in range averaged channels of uniform sound velocity,” J. Acoust. Soc. Am. 68(1), 282–286 (1980). 7 C. H. Harrison, “Ray convergence in a flux-like propagation formulation,” J. Acoust. Soc. Am. 133, 3777–3789 (2013). 8 M. B. Porter, The KRAKEN Normal Mode Program, SACLANT Undersea Research Centre Technical Report (1990). 9 M. A. Ainslie, Editorial: Validation of Sonar Performance Assessment Tools, in Validation of Sonar Performance Assessment Tools (workshop held in memory of David E. Weston, April 7–9, 2010), in Proceedings of the Institute of Acoustics, Vol. 32, Part 2 (2010), pp. 1–8. 10 M. Zampolli, M. A. Ainslie, and P. Schippers, “Scenarios for benchmarking range-dependent active sonar performance models,” in Proceedings of the Institute of Acoustics, Vol. 32, Part 2 – Research Symposium (2010), pp. 53–63. 11 D. E. Weston, “A Moire fringe analog of sound propagation in shallow water,” J. Acoust. Soc. Am. 32(6), 647–654 (1960). 12 D. E. Weston, “Wave shifts, beam shifts, and their role in modal and adiabatic propagation,” J. Acoust. Soc. Am. 96, 406–416 (1994). 13 C. T. Tindle and D. E. Weston, “Connection of acoustic beam displacement, cycle distances and attenuations for rays and normal modes,” J. Acoust. Soc. Am. 67, 1614–1622 (1980). 14 D. E. Weston and C. T. Tindle, “Reflection loss and mode attenuation in a Pekeris model,” J. Acoust. Soc. Am. 66(3), 872–879 (1979). 15 C. T. Tindle, “The equivalence of bottom loss and mode attenuation per cycle in underwater acoustic,” J. Acoust. Soc. Am. 66(1), 250–255 (1979). 16 J. D. Macpherson and M. J. Daintith, “Practical model of shallow-water acoustic propagation,” J. Acoust. Soc. Am. 41(4), 850–854 (1967). 17 G. P. M. Poppe and C. M. J. Wijers, “More efficient computation of the complex error function,” ACM Trans. Math. Software 16, 38–46 (1990). 18 M. Abramowitz and I. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables (Dover Publications, New York, 1972). 19 C. H. Harrison, “An approximate form of the Rayleigh reflection loss and its phase: Application to reverberation calculation,” J Acoust. Soc. Am. 128(1), 50–57 (2010). 20 E. I. Thorsos and J. S. Perkins, “Overview of the reverberation modeling workshops,” J. Acoust. Soc. Am. 122(5), 3074 (2007). 21 R. N. Denham, “Intensity decay laws for near-surface sound sources in the ocean,” J. Acoust. Soc. Am. 79, 60–63 (1986).



Characteristics and Propagation of Airgun Pulses in Shallow Water with Implications for Effects on Small Marine Mammals.

Seismic footprints of shallow dyke propagation at Etna, Italy.

Ring-type multisoliton dynamics in shallow water.

Autoregressive model for high-resolution wavenumber estimation in a shallow water environment using a broadband source.

Developing Benthic Class Specific, Chlorophyll-a Retrieving Algorithms for Optically-Shallow Water Using SeaWiFS.

The shallow water equation and the vorticity equation for a change in height of the topography.

A theoretical model of linearly filtered reverberation for pulsed active sonar in shallow water.

Soliton turbulence in shallow water ocean surface waves.

Shallow-water zoantharians (Cnidaria, Hexacorallia) from the Central Indo-Pacific.

High Resolution Marine Magnetic Survey of Shallow Water Littoral Area.

Shallow-water soliton dynamics beyond the Korteweg-de Vries equation.

Nitrate removal in shallow, open-water treatment wetlands.

Solving global shallow water equations on heterogeneous supercomputers.

Methane production in shallow-water, tropical marine sediments.

Expansion shock waves in regularized shallow-water theory.

Marine mammal audibility of selected shallow-water survey sources.

Pneumomediastinum and bilateral pneumothorax following near drowning in shallow water.

Sediment color tool for targeting arsenic-safe aquifers for the installation of shallow drinking water tubewells.

A formula for success.

[A formula for health].

Theoretical analysis of warping operators for non-ideal shallow water waveguides.

Vitellibacter nionensis sp. nov., isolated from a shallow water hydrothermal vent.

Shallow water marine gammaridean amphipods of Pulau Tioman, Malaysia, with the description of a new species.

On shallow water waves in a medium with time-dependent dispersion and nonlinearity coefficients.