Sound-pressure level calculations using the RRA algorithm for depth-dependent speeds of sound valid at turning points and focal points

Sound-pressure level (SPL) calculations are made along individual ray paths for arbitrary, one-dimensional, depth-dependent speeds of sound using an enhanced version of the RRA (recursive ray acoustics) algorithm. The SPL calculations are valid (i.e., finite) at turning points and focal points and do not require the use of Airy functions. The SPL calculations include the effects of frequency-dependent volume attenuation and frequency-dependent attenuation due to surface and bottom reflections. The ocean surface and bottom are treated as boundaries between viscous fluid media. Although the ocean surface is modeled as a planar boundary, the bathymetry is an arbitrary function of horizontal range. Sound speed versus depth and bathymetric data are represented by orthogonal function expansions. Computer simulation results from preliminary test cases are presented. >


I. INTRODUCTION
N a recent paper, Ziomek and Polnicky [l] introduced the I RRA (recursive ray acoustics) algorithm. The RRA algorithm is simple, fast, and accurate (based on the preliminary test case results presented in [l]) and uses arc length as the independent variable. It can be used to compute the position, angles of propagation, travel time, and path length along a ray path for three-dimensional speeds of sound. Ziomek and Polnicky [l] did not discuss sound-pressure level (SPL) calculations.
The main purpose of this paper is to demonstrate how to use the output from the RFL4 algorithm to calculate the SPL along individual ray paths for arbitrary, one-dimensional, depthdependent speeds of sound. Standard ray acoustics "correction formulas" for evaluating the acoustic field at a turning point or focal point for arbitrary, one-dimensional, depth-dependent speeds of sound require the use of Airy functions (e.g., see [2] and [3]). In contrast, the method used to perform SPL calculations in this paper is based on the approach presented by Krol [4]. Krol's method [4] does not require the use of Airy functions. Krol [4] demonstrated, via computer simulation results, intensity calculations along individual ray paths that were finite in value at both turning points and focal points. The computer simulation results reported by Krol [4] were for free-space propagation only, that is, the ray paths in his example encountered no surface and bottom reflections. Although the RRA algorithm can handle three-dimensional speeds of sound, the method described by Krol [4] to calculate the intensity along a ray path is only applicable to onedimensional, depth-dependent speeds of sound. However, the method of Gaussian beam tracing (e.g., see [5] and [6]) can be used to calculate the SPL not only for one-dimensional speeds of sound [5], but also for two-and three-dimensional speeds of sound [6]. The Gaussian beam method associates with each ray a beam with a Gaussian intensity profile normal to the ray. Auxiliary differential equations that govern the evolution of both the beam curvature and width are integrated along with the usual ray equations. In contrast, Krol's method [4] is based solely on one-dimensional ray acoustics-Gaussian intensity profiles are not artificially associated with individual ray paths.
The Gaussian beam method has the advantages of not requiring eigenray computations and not producing infinitely high intensity levels at caustics and zero intensity levels in shadow zones. However, the main disadvantage is that, at present, there is no agreed-upon method to compute the initial values for the beam width and curvature [5], [6]. Adjusting the initial conditions for the Gaussian beams can significantly alter the SPL at a receiver [5].
The SPL calculations discussed in this paper (see Section 11) include the effects of frequency-dependent volume attenuation and frequency-dependent attenuation due to surface and bottom reflections. The ocean surface and bottom are treated as boundaries between viscous fluid media. Although the ocean surface is modeled as a planar boundary, the bathymetry is an arbitrary function of horizontal range. Sound speed versus depth and bathymetric data are represented by orthogonal function expansions (see Section 111). Computer simulation results from preliminary test cases are presented in Section IV.

SPL CALCULATIONS
For an omnidirectional point source located at depth yo (see Fig. 10 in the Appendix), and for a speed of sound c(y) an arbitrary function of depth y, the magnitude squared of the acoustic pressure along a ray path is given by (see Appendix)

U.S. Government Work Not Protected by U.S. Copyright
is the horizontal-range, z and z are the cross-range and downrange coordinates, respectively, Po is the initial angle of propagation at the source, po(y) is the depth-dependent ambient (equilibrium) density, and P( y ) is the angle of propagation (arrival) at depth y (see Fig. 10 in the Appendix). In order to evaluate ( l ) , an equation for the denominator term cosP(y) dT/d,& must be obtained. Following the approach used by Krol [4], it can be shown that  when where C"(Y) = ,C(Y)l d2 dY (9) CO = c(yo), cb = c'(yo), and where it is understood that y = IJ(T), that is, the horizontal range T along a ray path is the independent variable. The denominator term D given by (5) and (7) are exact-no additional simplifying assumptions or approximations were made. According to (€9, (7) would be used, for example, whenever a ray path approached a turning point, that is, whenever P(y) -+ go", resulting in cosP(y) + 0. Note that if we let RO = 1 m, then ~~( T O I , yO1)l2 appearing in (1) is related to the source level (SL) of the omnidirectional point source.
In order to make the calculation of I~(T, y)I2 given by An infinitesimal element of arc length ds at an arbitrary point P integrations over the arc length s along a ray path (see Fig.  1). Recall that Snell's law is given by where is the ray parameter. Therefore, by referring to Fig. 1 and using ( l l ) , Also note that is the direction cosine with respect to the Y axis. With the use of (1 1)-(14), (5)-(8) can be rewritten as follows: where it is understood that y = y(s), that is, arc length s along a ray path is the independent variable.
Note that the RRA algorithm [ l ] already provides values for the cross-range z, depth y, down-range z , and the direction cosine v as functions of arc length s along a ray path, that is, z = ~( s ) , y = y(s), z = z ( s ) , and v(y) = v[y(s)]. Therefore, I the horizontal range T along a ray path can also be expressed as a function of arc length, that is,

(19)
However, to be more specific, the RRA algorithm [l] actually provides values for the aforementioned quantities as functions of step number i, that is, x; = ~(s;), y; = y(s;), z; = z(s;), and vi = v(y;) = w [ y ( s ; ) ] . Therefore, the magnitude squared of the acoustic pressure along a ray path can finally be expressed as follows: where Equation (15) is used when and (17) is used when It is often convenient to express the values of acoustic pressure in terms of sound-pressure levels. The sound-pressure level (SPL) is defined as follows: where (pl is the magnitude of the acoustic pressure, Pref is the root-mean-square (rms) reference pressure, and "re" means "relative to." For problems in underwater acoustics, Pr,f = 1 pPa (rms) is common. Therefore, the SPL along a ray path as a function of step number i can be expressed as The effects of frequency-dependent attenuation due to surface and bottom reflections can be modeled by using appropriate equations for the reflection coefficients at the ocean surface and bottom. If we designate the air as medium I, the ocean medium as medium 11, and the ocean bottom as medium 111, and since the sound source is in medium 11, then the reflection coefficient at the ocean surface Rzl and the reflection coefficient at the ocean bottom R23 are given by p2B2m, where pm, m = 1, 2, 3 is the ambient (equilibrium) density in kilograms per cubic meter of media I, 11, and 111, respectively, Gnc is the angle of incidence measured with respect to the local normal to the boundary at the point of incidence, is the complex index of refraction between media I1 and I, and I1 and 111, respectively, is the complex wavenumber, is the real wavenumber in radians per meter, cm is the speed of sound in meters per second, and a,(f) is the real, frequencydependent attenuation coefficient in nepers per meter in media I, 11, and 111, respectively. Note that for ideal, nonviscous fluids, a m ( f ) = 0, m = 1, 2 , 3. The right-hand side of (20) is multiplied by the magnitude squared of either Rzl or R23 in order to take into account the effects of a surface or bottom reflection, respectively. In addition, the effect of frequency-dependent vdume attenuation can be modeled by multiplying the right-hand side of (20)

by exp [ -2 a~( f ) s ; ] ,
where a 2 ( f ) is the frequencydependent attenuation coefficient of the ocean medium in Np/m and si is the value of the arc length along a ray path at step number i (an output of the RRA algorithm). The factor of 2 is present in the exponent because (20) is an expression for the magnitude squared of the acoustic pressure.

ORTHOGONAL FUNCTION EXPANSIONS OF SOUND SPEED AND BATHYMETRIC DATA
In this section, we shall discuss why and how sound speed versus depth data and bathymetry versus horizontal range data are represented by orthogonal function expansions. Measured data are usually noise corrupted. Even in the absence of noise, data are uncertain due to the lack of precision in the measurement equipment. As a result, if an interpolation scheme such as cubic splines is used to fit noise corrupted and/or uncertain measured data, then the noise and/or the uncertainty in the data are being fit rather than the trend in the data since an interpolation scheme will pass a curve through each data point.
When dealing with noise-corrupted and/or uncertain data, it is better to fit the data using an orthogonal function expansion. An orthogonal function expansion will fit the data in a minimum mean-squared error (MMSE) sense using orthogonal polynomials. It is analogous to, but not the same as, the method of nonlinear least squares estimation. The order of the orthogonal function expansion of the data can be increased by simply computing the next higher order coefficient-all previous coefficients do not have to be recomputed. In addition, the data do not have to be evenly spaced. In order to determine the order N, and Nb to be used for the orthogonal function expansion of both the sound speed and bathymetric data, respectively, the following procedure is used. Since the procedure is the same for both the sound speed and bathymetric data, we shall explain it for the soundspeed data only. Using the measured sound speed versus depth data, an estimate of the speed of sound t(y) given by (32) From Fig. 4, it can be seen that both rays "compress" as they propagate up-slope, as predicted by theory. However, the first ray (PO = 55') turns back on itself at T M 2.3 km and propagates back toward the source, arriving at T = 0 km with a depth of y M 50 m. This phenomenon of rays turning back on themselves and propagating back toward the sound source is also predicted by theory for up-slope propagation. Although both rays shown in Fig. 4 underwent many surface and bottom reflections, since the ocean surface and bottom were modeled as ideal pressure-release and rigid boundaries, the corresponding SPL plots shown in Fig. 5 are smooth with no discontinuities. In addition, since the ocean medium was modeled as being nonviscous with a constant speed of sound, the SPL decreased at the expected rate of 6 dB for each doubling of the range.
Test case 2 is identical to test case 1, with the exception that medium I1 (ocean water) and medium I11 (ocean bottom) were modeled as viscous fluid media using the following parameter values: Gnc < cc = 60.1' (see Fig. 7). Note that SPL values below -25 dB were set equal to -25 dB for simulation purposes.
In contrast, the angle of incidence Gnc > cc for each bottom reflection that ray 2 (PO = 80') underwent, except for the lust two bottom reflections where Gnc = 59' and 53", respectively, which are less than cc = 60.1'. As a result, for those bottom one would find small, finite discontinuities since I R23 I < 1 for Gnc > cc because of frequency-dependent attenuation.
Test case 3 simulated ray propagation in an ocean waveguide with a focal point problem for the calculation of SPL. The ocean surface and bottom were once again modeled as ideal pressure-release and rigid boundaries, with reflection coefficients R21 = -1 and R23 = 1, respectively.
The ocean medium was modeled as being nonviscous with the sound-speed profile shown in Fig. 8. The sound source was placed at a depth of yo = 55 m with a SL = 180 dB re 1 pPa (rms). Three rays were transmitted from the source with initial angles of propagation PO = 94.2', 94.3', and 94.4'.
An arc length step size of 0.5 m was used. Figs. 8 and 9 illustrate the corresponding ray trace and SPL plots for test case 3. Fig. 9 correctly shows a dramatic increase in SPL for the three rays at the location of the focal points.   Fig. 10 illustrates a slice through a ray tube where T is the horizontal range. An omnidirectional point source is located at depth yo, producing an axisymmetric acoustic field. If we invoke conservation of energy, then where Iavg(.) is the time-average intensity, and dS1 and dS are the infinitesimal wavefront surface areas at the beginning and end of the ray tube, respectively, created by rotating the ray tube slice illustrated in Fig. 10, 360' around the Y axis. If both the speed of sound and the equilibrium (ambient) density are only functions of depth y, then In order to proceed further, we need expressions for dS1 and dS.
By refemng to Fig. 10, we can write that where 2m-l is the circumference of a circle with radius is an infinitesimal element of arc length along the wavefront that passes through the point ( r , y). Substituting (A9) into (A8) yields d S = 21rr cos p(g/) dr.

W O )
dl = cos P(y) d r Therefore, upon substituting (A2), (A6), and (A10) into (Al), we obtain the following expression for the magnitude squared of the acoustic pressure along a ray path: where the magnitude of cosP(y) dr/dPo is taken in order to ensure that the magnitude squared of the acoustic pressure is positive.