Ray Trace Modeling of Underwater Sound Propagation

1. Introduction

Modeling acoustic propagation conditions is an important issue in underwater acoustics and there exist several mathematical/numerical models based on different approaches. Some of the most used approaches are based on ray theory, modal expansion and wave number integration techniques. Ray acoustics and ray tracing techniques are the most intuitive and often the simplest means for modeling sound propagation in the sea. Ray acoustics is based on the assumption that sound propagates along rays that are normal to wave fronts, the surfaces of constant phase of the acoustic waves. When generated from a point source in a medium with constant sound speed, the wave fronts form surfaces that are concentric circles, and the sound follows straight line paths that radiate out from the sound source. If the speed of sound is not constant, the rays follow curved paths rather than straight ones. The computational technique known as ray tracing is a method used to calculate the trajectories of the ray paths of sound from the source.

Ray theory is derived from the wave equation when some simplifying assumptions are introduced and the method is essentially a high-frequency approximation. The method is sufficiently accurate for applications involving echo sounders, sonar, and communications systems for short and medium short distances. These devices normally use frequencies that satisfy the high frequency conditions. This article demonstrates that ray theory also can be successfully applied for much lower frequencies approaching the regime of seismic frequencies.

This article presents classical ray theory and demonstrates that ray theory gives a valuable insight and physical picture of how sound propagates in inhomogeneous media. However, ray theory has limitations and may not be valid for precise predictions of sound levels, especially in situations where refraction effects and focusing of sound are important. There exist corrective measures that can be used to improve classical ray theory, but these are not discussed in detail here. Recommended alternative readings include the books. [1-4] and the articles [5-6].

A number of realistic examples and cases are presented with the objective to describe some of the most important aspects of sound propagation in the oceans. This includes the effects of geographical and oceanographic seasonal changes and how the geoacoustic properties of the sea bottom may limit the propagation ranges, especially at low frequencies. The examples are based on experience from modeling sonar systems, underwater acoustic communication links and propagation of low frequency noise in the oceans. There exist a number of ray trace models, some are tuned to specific applications, and others are more general. In this chapter the applications and use of ray theory are illustrated by using Plane Ray, a ray tracing program developed by the author, for modeling underwater acoustic propagation with moderately range-varying bathymetry over layered bottom with a thin fluid sedimentary layer over a solid half with arbitrary geo-acoustic properties. However, the discussion is quite general and does not depend on the actual implementation of the theory.

2. Theory of ray acoustics

The theory of ray acoustics can be found in most books and [1-4] will not be repeated here, but instead we follow a heuristic approach based on Snell’s law, which is expressed by.

$$\xi =\frac{cos\theta \left(z\right)}{c\left(z\right)}=\frac{cos{\theta }_0}{c_0}\text{}.$$E1

Figure 1 shows a small segment of a ray path and the coordinate system. The segment has horizontal and vertical components dz and dr, respectively, and has the angleθ with the horizontal plane. When the speed of sound varies with depth the ray paths will bend and the rays propagate along curved paths. The radius of curvature R is defined as the ratio between an increment in the arc length and an increment in the angle

$$R=\frac{ds}{d\theta }.$$E2

Figure 1 shows that the radius of curvature is

$$R=\frac{1}{sin\theta }\frac{dz}{d\theta }.$$E3

When the sound speed varies with depth the ray angle θ is a function of depth according to Snell’s law. Taking the derivative of Eq. (3) with respect to gives the ray’s radius of curvature at depth z expressed as

$$R\left(z\right)=-\frac{c\left(z\right)}{cos\theta \left(z\right)}\frac{1}{g\left(z\right)}\text{=}-\frac{\text{1}}{\xi \text{g(z)}}.$$E4

The ray parameter ξ is defined in Eq. (1) and g(z) is the sound speed gradient.

$$g\left(z\right)=\frac{dc\left(z\right)}{dz}.$$E5

At any point in space, the ray curvature is therefore given by the ray parameter ξ and the local value of the sound speed gradient g(z). The positive or negative sign of the gradient determines whether the sign of R is negative or positive, and thereby determines if the ray path curves downward or upward.

A ray with horizontal angle θ_in strikes a plane with inclination α, the reflected ray is changed toθ _out.

Figure 1.

A small segment of a ray path in a isotropic medium with arc length ds.

Figure 2.

A ray with horizontal angle θ_in strikes a plane with inclination α, the reflected ray is changed toθ _out.

The ray parameter is not constant when the bathymetry varies with range. The change in ray direction is illustrated in Figure 2 showing that after reflection the angle θ_in of an incoming ray is increased by twice the bottom inclination angle α.

$${\theta }_{out}={\theta }_{in}+2\alpha .$$E6

Consequently, after the ray is reflected, its ray parameter must change from ξ_in to ξ_out, which is expressed as

$$\begin{array}{c}{\xi }_{out}=\frac{cos({\theta }_{out})}{c}=\frac{cos({\theta }_{in}+2\alpha )}{c}\\ ={\xi }_{in}cos(2\alpha )\pm \frac{\sqrt{1-{\xi }_{in}^2{c}^2}}{c}sin(2\alpha ).\end{array}$$E7

The coordinates of a ray, starting with the angle θ₁ at the point (r₁,z₁), where the sound speed is c₁, as shown in Figure 3. For the coordinates of the running point at (r₂, z₂) along the ray path, the horizontal distance is

$$r_2-r_1={\displaystyle \underset{z_1}{\overset{z_2}{\int }}\frac{dz}{tan\theta \left(z\right)}}={\displaystyle \underset{z_1}{\overset{z_2}{\int }}\frac{\cos \theta \left(z\right)dz}{\sqrt{1-{\cos }^2\theta \left(z\right)}}}\text{ = }{\displaystyle \underset{z_1}{\overset{z_2}{\int }}\frac{\xi c\left(z\right)dz}{\sqrt{1-{\xi }^2{c}^2\left(z\right)}}}.$$E8

The travel time between the two points is obtained by integrating the quantity 1/c, the slowness, along the ray path:

$${\tau }_2-{\tau }_1={\displaystyle \underset{z_1}{\overset{z_2}{\int }}\frac{ds}{c\left(s\right)}}={\displaystyle \underset{z_1}{\overset{z_2}{\int }}\frac{dz}{c\left(z\right)sin\theta \left(z\right)}}={\displaystyle \underset{z_1}{\overset{z_2}{\int }}\frac{dz}{c\left(z\right)\sqrt{1-{\xi }^2{c}^2(z)}}}.$$E9

Figure 3.

Left: The sound speed profile. Right: A portion of a ray traveling from point (r₁, z₁) to (r₂, z₂).

The acoustic intensity of a ray can, according to ray theory, be calculated using the principle that the power within a ray tube remains constant within that ray tube. This is illustrated in Figure 4 showing two rays with a vertical angle separation of dθ₀ that define a ray tube centered on the initial angle θ₀. At a reference distance r₀ from the source, the intensity is I₀. Taking into consideration the cylindrical symmetry about the z axis, the power ΔP₀ within the narrow angle dθ₀ is

$$\Delta P_0=I_{0}2\pi r_0^{2}cos{\theta }_{0}d{\theta }_0.$$E10

At horizontal distance r, the intensity is I. In terms of the perpendicular cross section dL of the ray tube, the power is

$$\Delta P=I2\pi rdL.$$E11

Since the power in the ray tube does not change, we may equate Eq.(9) and eq. (10), and solve for the ratio of the intensities:

$$\frac{I}{I_0}=\frac{r_0^2}{r}cos{\theta }_0\left|\frac{d{\theta }_0}{dL}\right|.$$E12

Figure 4.

The principle of intensity calculations: energy radiated in a narrow tube remains inside the tube; r₀ represents a reference distance and θ₀ is the initial ray angle at the source; dθ₀ is the initial angular separation between two rays; dr is the incremental range increase; θ is the angle at the field point; dz is the depth differential; and dL is the width of the ray tube.

Instead of using Eq.(11), it may be more convenient to use the vertical horizontal ray dr, which is

$$dr=\left|\frac{dL}{sin\theta }\right|,$$E13

resulting in

$$\frac{I}{I_0}=\frac{r_0^2}{r}\left|\frac{cos{\theta }_0}{sin\theta }\frac{d{\theta }_0}{dr}\right|=\left(\frac{r_0^2}{r}\right)\left(\frac{c_0}{c}\right)\left|\frac{cos\theta }{sin\theta }\frac{d{\theta }_0}{dr}\right|.$$E14

The last expression in Eq.(13) is obtained by assuming that the ray parameter is constant and by using Snell’s law. The absolute values are introduced to avoid problems with regard to the signs of the derivatives and of sinθ.

With respect to the reference distance r₀, the transmission loss TL is defined as

$$TL=-10log\left(I/I_0\right).$$E15

By inserting Eq.(13) into Eq.(14) the transmission loss becomes

$$TL=10log\left(r/r_0^2\right)+10log\left|\frac{dz}{d{\theta }_0}\right|+10log\left(\frac{c}{c_0}\right).$$E16

The term c₀/c is close to unity in water and can be ignored in most cases.

In this treatment the transmission loss includes only the geometric spreading loss. Therefore bottom and surface reflection losses and sea water absorption loss must be included separately.

The geometric transmission loss in Eq.(15) consists of two parts. The first term represents the horizontal spreading of the ray tube and results in a cylindrical spreading loss. The second and third terms represent the vertical spreading of the ray tube and are influenced by the depth gradient of the sound speed.

Eq.(13) predicts infinite intensity under either of two conditions: when θ = 0 or when dr / dθ₀ = 0. The first condition signifies a turning point where the ray path becomes horizontal; the second condition occurs at points where an infinitesimal increase in the initial angle of the ray produces no change in the horizontal range traversed by the ray. The locus of all such points in space is called a caustic. In both cases there is focusing of energy by refraction and where classical ray theory incorrectly predicts infinite intensity. Caustics and turning points will be discussed further in section 8.2.

3. A recipe for tracing of rays

A simple receipt for a ray tracing algorithm is to divide the whole water column into a large number of layers, each with the same thickness Δz. Within each layer, the sound speed profile is approximated as linear so that, in the layer z_i < z < z_i+1, the sound speed is taken to be

$$c\left(z\right)=c_i+g_i\left(z-z_i\right).$$E17

where c_i is the speed at depth z_i, and g_i is the sound speed gradient in the layer. From Eq. (7) and Eq. (8) the range and travel time increments in the layer are given by

$$r_{i+1}-r_i=\frac{1}{\xi g_i}\left[\sqrt{1-{\xi }^2{c}^2(z_i)}-\sqrt{1-{\xi }^2{c}^2(z_{i+1})}\right],$$E18

and

$${\tau }_{i+1}-{\tau }_i=\frac{1}{\left|g_i\right|}ln\left(\frac{c(z_{i+1})}{c(z_i)}\frac{1+\sqrt{1-{\xi }^2{c}^2(z_i)}}{1+\sqrt{1-{\xi }^2{c}^2(z_{i+1})}}\right).$$E19

When$${\xi }^2{c}^2(z_{i+1})\ge \text{1}$$, the ray path turns at a depth between z_i and z_i+1, and Eq.s (17) and (18) must be replaced by the following expressions:

$$r_{i+1}-r_i=\frac{2}{\xi g_i}\sqrt{1-{\xi }^2{c}^2(z_i)},$$E20

and

$${\tau }_{i+1}-{\tau }_i=\frac{2}{\left|g_i\right|}ln\left(\frac{1+\sqrt{1-{\xi }^2{c}^2(z_i)}}{\xi c(z_i)}\right).$$E21

These equations give the trajectories and the travel times for any ray’s path to the desired range. By applying Eqs. (13) and (14), the geometrical transmission loss is also determined. The simplicity of this method lies in the approximation of the sound speed profiles with straight-line segments and the ray path’s subsequent decomposition into circular segments. The method’s accuracy is determined by how well the linear fit matches the actual profile. In practice, the sound speed profile is often given as measured sound speeds at relatively few depth points. It is therefore advisable to use an interpolation scheme that is consistent with the usual behavior of the sound speed profile to increase the number of depth points to an acceptable high density.

The examples in this article are generated using the ray trace program PlaneRay that has been developed by the author [7-8]. However, any other ray programs with similar capabilities could have been use and the discussion is therefore valid for ray modeling in general. Other models frequently used and are the Bellhop model [9], and the models [10-11].

Figure 5 shows an example of ray modeling. The sound speed profile is shown at the left panel and the rays from a source at 50 m depth is shown in the right panel, which also shows the bathymetry and the thickness of the sediment layer over the solid half space.

Figure 5.

Sound speed profile and ray traces for a typical case. The source depth is 150 m and the red dotted line indicates a receiver line at a depth of 50 m. The initial angles of the rays at the source are from –30º to 30º.

4. Eigenray determination

To calculate the acoustic field it is necessary to have an efficient and accurate algorithm for determination of eigenrays. An eigenray is defined as a ray that connects a source position with a receiver position. In most case with multipath propagation there are many eigenrays for a given source/receiver configuration, which means that finding all eigenrays is not a trivial task.

The PlaneRay model uses a unique sorting and interpolation routine for efficient determination of a large number of eigenrays in range dependent environments. This approach is described by the two plots in Figure 6, which displays the ray history as function of initial angle at the source. All facts and features of the acoustic fields such as the transmission loss, transfer function and time responses are derived from the ray traces and their history The two plots show the ranges and travel times to where the rays cross the receiver depth line (marked by the red dashed line in Figure 5). A particular ray may intersect the receiver depth line, at several ranges. For instance at the range of 2 km, there are 11 eigenrays and from Figure 6 the initial angles of these rays are approximately found to be 5.9°, 9.6°, 22°, 24° for the positive (down going) rays and−2.0°,−3.6, °− 7-4°− 15.0°− 17.0°− 25.0°,−27.0°, for the negative (up-going waves). However, the values found in this way are often not sufficiently accurate for the determination of the sound field. Further processing may therefore be required to obtain accurate results.

The graphs of Figure 6 are composed of independent points, but it is evident that the points are clustered in independent clusters or groups. This property is used for sorting the points into branches of curves that represents different ray history. These branches are in most case relatively continuous and therefore amenable to interpolation. An additional advantage of this method is that the contribution of the various multipath arrivals can be evaluated separately, thereby enabling the user to study the structure of the field in detail.

Figure 6.

Ray history of the initial ray tracing in Figure 5 showing range (left) and travel time (right) to the receiver depth as function of initial angel at the source.

In most cases the eigenrays are determined by one simple interpolation yields values that are sufficiently accurate for most application, but the accuracy increases with increasing density of the initial angles at the cost of longer computation times.

Figure 7 shows examples of eigenrays traces with rays a receiver located at 2.5 km from the source for the scenario shown in Figure 5. To this receiver there are a total of 12 eigenrays, spanning the range of initial angles from -30° to 29°.

Figure 7.

Eigenrays from a source at 150 m depth to a receiver at 50 m depth and distance of 2.5 km from the source.

5. Acoustic absorption in sea water

Sound absorption is important for long range propagation especially at higher frequencies. The absorption increases with frequencies and is dependent on temperature, salinity, depth and the pH value of the water. There exists several expressions for acoustic absorption in sea water; one of the preferred options is the semi-empirical formulae by Francoise and Garrison [12]. Figure 8 shows sound absorption as function of frequency in sea water using this expression for the values given in the figure caption.

Figure 8.

Acoustic absorption (dB/km) for fresh water and saltwater, plotted as a function of frequency (kHz) for water temperature of 10ºC, atmospheric pressure of one atmosphere (surface), salinity of 35 pro mille, and pH value of 7.8. The various contributions to the absorption are also indicated.

6. Boundary conditions at the surface and bottom interfaces

Ray tracing is greatly simplified when no rays are traced into the bottom, but stops at the water-bottom interface. This avoids tracing of multiple reflections in layered bottoms. Instead the boundary conditions at the sea surface and the bottom can be approximately satisfied by the use of plane wave reflection coefficient.

A simple and useful bottom model is assuming a fluid sedimentary layer over a homogeneous solid half space. The reflection coefficient of a bottom with this structure is

$$R_b=\frac{r_{01}+r_{12}\exp \left(-2i{\gamma }_{p1}D\right)}{1+r_{01}{r}_{12}\exp \left(-2i{\gamma }_{p1}D\right)},$$E22

where γ_p1 is the vertical wave number for sediment layer and D is the thickness of the sediment layer. The reflection coefficient between the water and the sediment layer, r₀₁, is given as

$$r_{01}=\frac{Z_p{}_1-Z_p{}_0}{Z_p{}_1+Z_p{}_0},$$E23

and r₁₂ is the reflection coefficient between the sediment layer and the solid half space,

$$r_{12}=\frac{Z_{p2}{\cos }^{2}2{\theta }_{s2}+Z_{s2}{\sin }^{2}2{\theta }_{s2}-Z_{p1}}{Z_{p2}{\cos }^{2}2{\theta }_{s2}+Z_{s2}{\sin }^{2}2{\theta }_{s2}+Z_{p1}}\text{}.$$E24

In Eq (15) and (16) Z_ki is the acoustic impedance for the compressional (k = p) and shear (k = s) waves in water column (i = 0), sediment layer (i = 1) and solid half-space (i = 2), respectively. The grazing angle of the transmitted shear wave in the solid half-space is denoted θ_s2.

Figure 9 shows an example of the bottom reflection loss as function of angle and frequency for a bottom with a sediment layer with the thickness D = 2 m with sound speed of 1700 m/s and density 1800 kg/m³ over a homogenous solid half space with compressional speed 3000 m/s, shear speed 500 m/s and density 2500 kg/m³. The wave attenuations are 0.5 dB/ wavelength. The critical angle changes from 60° at very low frequencies to about 28° at high frequencies, the two angles are given by the sound speed in the water and the two bottom sound speed of 3000 m/s and 1700 m/s. The small, but significant, reflection loss at lower angles is caused by shear wave conversion and bottom absorption In this case the attenuation is about 1 dB in the frequency band around 50 Hz to 100 Hz.

The reflection coefficient of a flat even sea surface is −1 for. For a sea surface with ocean waves there will be diffuse scattering to all other direction than the specular direction, which result in a reflection loss that in the first approximation can be modeled by the coherent rough surface reflection coefficient

$$R_{coh}=exp\left[-2{\left(\frac{2\pi }{\lambda }{\sigma }_{h}sin\theta \right)}^2\right].$$E25

In this expressionθ is the grazing angle and σ_h is the rms. wave height and λ, is the acoustic wavelength, both in meters.

Figure 9.

Bottom reflection loss (dB) as function of frequency and incident angle for a 2 m sediment layer over solid rock. The parameters are given in the text.

The reflection loss associated with reflection from a rough sea surface is

$$RL=-20\log 10\left(R_{coh}\right).$$E26

The same rough surface reflection coefficient may also be applied to a rough bottom.

Figure 10 shows the rough surface reflection loss as function of grazing angle, calculated for a wave height of 0.5 m and the frequencies of 50 Hz, 100 Hz, 200 Hz and 400 Hz.

Figure 10.

Reflection loss (dB) of rough surface with rms. wave height of 0.5 m as function of grazing angle, for the frequencies in the legend

7. Synthesizing the frequency domain transfer function and the time responses

The total wave field at any receiving point is calculated in the frequency domain by coherent summation of all the eigenray contributions. The first step in the calculation is to determine the geometrical transmission loss of each of the multipath contributions by applying Eq. (13) and Eq.(14) to the sorted and interpolated range-angle values. The frequency domain transfer function and the transmission loss are obtained by adding the multipath contributions coherently in frequency domain taken into account the phase shifts associated the travel times from the interpolated history of the travel times. The frequency dependent acoustic absorption of sound in water is included at this point in the process. The transfer function H(ω, r) can be expressed as

$$H\left(\omega ,r\right)={\displaystyle \sum _n{A}_n{B}_n}\left(\omega \right)S_n\left(\omega \right)T_n\exp \left(i\omega {\tau }_n\right).$$E27

Eq. (26) expresses the transfer function H(ω, r) to a distance r from the source at the at angular frequency ω as a sum over the n eigenrays that are included in the synthesis. A_n is the geometrical spreading loss factor, defined as the square root of the expression in Eq. (13). B_n, and S_n, are the combined effects of all bottom reflections and surface reflections, respectively, T_n, is −90º phase shift associated with caustics and turning points, and τ_n is the travel time.

The synthesis of the received signals is performed in the frequency domain by multiplying the frequency spectrum of the source signal with the transfer function of each of the eigenrays and summing the contributions. The time domain response is obtained after multiplication with the frequency function of a source signal followed an inverse Fourier transform of the product. This requires the choice of a source signal, a sampling frequency (f_s) and a block length (N_fft) of the Fourier transform.

The total duration of the time window (T_max) after Fourier transform is

$$T_{\max }=\frac{N_{fft}}{f_s}.$$E28

It is important to select the values of N_fft and f_s such that Fourier time window, T_max, is larger than the actual length or duration of the signal. In reality the real time duration of the received signal is often not known in advanced and therefore the user may have to experiment with different values to find appropriate values for of N_fft and f_s.

Figure 11 shows an example where the transmission loss (in dB) as function of range has been calculated for the frequencies of 100 Hz and 200 Hz. The dashed black line indicates the geometrical spreading loss, which is added for comparison and given by,

$$T{L}_{geo}\left(r\right)=10{\log }_{10}\left[r^2{\left(1+\frac{r^2}{r_t^2}\right)}^{-1/2}\right].$$E29

This expression yields a transmission loss proportional to 20log(r) when r< r_t and proportional to 10 log(r) for r> r_t. This approximation to the geometrical transmission loss may be used for approximate calculations of transmission loss for flat bottom and simple sound speed profiles. In the case shown in Figure 11 r_t is set equal to the water depth at source location, which in this case is 200 m.

Figure 11.

Transmission loss as function of range calculated for 100 Hz and 200 Hz The dashed black line is values of Eq.(28)

Figure 12 shows the synthesized time response at receivers spaced at 200 m separations in range up to 6 km. The sound speed and bathymetry is the same as in Figure 5 with the source at 150 m and all receivers at 50 m depth. The time scale is in reduced time to remove the gross transmission delay between the source and receiver. The reduced time is defined as

$$t_{red}=t_{real}-\frac{r}{c_{red}}.$$E30

In Eq. (29), t_real and t_red are the real and reduced times, respectively, r is range and c_red is the reduction speed. The actual value of c_red is not important as long as the chosen value results in a good display of the time responses.

Figure 12.

Received time signals as function of range and reduced time.

In the example shown above, the time signal and calculated assuming a narrow band-limited source signal in the form of a Ricker pulse. An example of a Ricker pulse and its frequency spectrum are shown in Figure 13.

Figure 13.

Ricker time pulse and frequency function the center frequency of 100 Hz

The time responses in ‎Figure 12 are sorted according to the history of their eigenrays and color coded to allow for studying the various multipath contributions. This is particularly useful when dealing with transient signal and broad band signal, especially when knowledge of the multipath structure is important. In many such situations only the direct arrival or the refracted arrivals in the water column may carry the useful signals and all the other arrivals represent interference. In this case there are direct arrivals, followed by surface reflected and refracted arrivals at the turning points. Notice the high sound pressure values caused by the caustics at 3 km, 6, km and 7 km, which are apparent in both plots, this issue is discussed in section 8.2.

The red dotted line in Figure 12 represents an estimate of the duration of the cannel impulse response. This time duration is mainly given by the bottom reflection coefficient and the critical angle. Rays that propagate at angles closer to the horizontal plane than the critical angle experience almost no bottom reflection loss and may therefore propagate to long distances. Rays with steeper angles will experience higher reflection losses and die out more rapidly with range. Thus the time duration of the impulse is directly determined by the ratio of sound speeds in the water and the bottom as

$$t_{red}=\frac{r}{c_0}\left(\frac{1}{\cos {\theta }_{crit}}-1\right)=r\left(\frac{c_b-c_0}{c_0^2}\right)$$E31

This estimate of the time duration of the channel impulse response assumes that the bottom is fluid, homogenous and flat, but the estimate may also be useful in other cases with moderately range dependent depth and with solid or layered bottom.

8. Special considerations

8.1. Frequency of applications

Ray tracing is a high frequency approximation to the solution of the wave equation and in principle more valid for high than for low frequency applications. However, high resolution prediction of higher frequency acoustic fields is difficult both for numerical and physical reasons. Principally most important is the physical limitation caused by the fact that the sound speed and the environment are generally not known in sufficient detail. This can be illustrated by a simple example. Consider coherent communication using a frequency of 10 kHz with wavelength of 10 mm. The required accuracy in order to be correct at a distance of 1 km is that the sound speed is known and stable with a relative error less 10^-5, an impossible requirement to satisfy in practice regardless of the numerical accuracy of the computer model.

8.2. Caustics and turning points

As mentioned before, the locations where dr/dθ₀=0 are called caustics where the ray phase is decreased by 90º and the where the intensity, according to ray theory, goes to infinity. In reality the intensity is high, but finite, and the basic ray theory breaks down at these points. There exists theories to amend and repair the defects of ray theory at these points [1, 2, 13], but that is not discussed here.

Figure 14 shows details of the field at a showing the rays with initial angles in the range of −6° to −1°. The scenario is the same as in of Figure 5, but for clarity the tracing of rays have been stopped after the first bottom reflection and the figure concentrates on the details the field at the caustic at 1760 m range for a ray with initial angle of −5.6°. Figure 15 shows the time responses for ranges in the interval from 1.6 km to 1.9 km. In this case, the source signal is a Ricker pulse with a peak frequency of 200 Hz. There is a first direct arrival (black color) at all ranges. From the range 1760 there is also a refracted arrival a little later than the direct, but with higher amplitude, in particular near the range of 1760 m. Notice the effect of the 90° phase shift for ranges beyond the caustic at 1760 m and that the amplitude at this range is considerable higher than at the other ranges.

Figure 14.

Rays through a caustic

Figure 15.

Time responses around the caustic at 1.76 km. The transmitted signal is a Ricker pulse with peak frequency of 200 Hz.

8.3. The principle of reciprocity and its validity in ray modeling

The principle of reciprocity is an important and useful property of linear acoustics and systems theory. The principle is very general and valid also in cases where the wave undergoes reflections at boundaries on its path from source to receiver [14]. The reciprocity principle is correctly represented in ray modeling, as can easily be understood from the eigenray plots of Figure 7. The eigenrays from a source position to the receiver position are the same as when source and receiver changes positions. The reflections at the bottom and at the sea surface are also symmetric in angles and consequently the acoustic fields are the same. However, it should be noted that the reciprocity principle applies to a point-to-point situation. This means that, for instance, that the development of the transmission loss as function source-receiver separation is generally not the same for the two directions.

8.4. The validity of using plane wave reflection coefficients

The accuracy of any ray model depends on the validity and limitation of ray theory and the implementation. A fundamental assumption of model is that the interactions with the boundaries are adequately described by plane wave reflection coefficient. In this section the validity of this assumption is investigated.

The general expression for the reflected field is given in text books, for instance in [13], over horizontal wave numbers k, as

$${\Phi }_{ref}\left(r,z_r,\omega \right)=\frac{S(\omega )}{8\pi i}{\displaystyle \underset{0}{\overset{\infty }{\int }}\Re \left(k\right)\frac{\exp \left(i\gamma \left|z_r+z_s\right|\right)}{\gamma }k{H}_0^{(1)}\left(kr\right)dk\text{.}}$$E32

Φ_ref(r, z_r ω) is the reflected field due to point source with frequency ω and source strength S(ω).$\Re \left(k\right)$ is the reflection coefficient,$H_0^1\left(kr\right)$ is the Hankel function of first kind, which represents a wave progressing in the positive r-direction. The horizontal wave number k and the vertical wave numberγ are related to the sound speed, frequency and the angle θ by

$$\begin{array}{l}k=\frac{\omega }{c}\cos \theta ,\\ \gamma =\frac{\omega }{c}\sin \theta .\end{array}$$E33

Eq.(31) states that the field is given as an integral over all horizontal wave numbers, or as consequence of Eq.(32), integration over all angles both real and the imaginary.

Consider now the situation where $\Re \left(k\right)=\Re$ is constant and independent of k or the angle. The integral in Eq.(31) becomes a standard integral and

$${\Phi }_{ref}\left(r,z_r,\omega \right)=\frac{S(\omega )}{4\pi R}\Re \exp (ikR).$$E34

$$R=\left[r^2+{\left(z_s+z_r\right)}^2\right].$$E35

According to Eq. (33) the reflected wave is the same as the outgoing spherical wave from the image of the source in the mirror position of the real source and modified by the constant reflection coefficient$\Re$. The situation with a constant reflection coefficient is valid for perfectly flat sea surface where the reflection coefficient is equal to -1 for all angles of incidence. Thus the reflection from a smooth sea surface is accurately described plane wave reflection coefficients.

In the general case, and for reflections from the bottom, the reflection coefficient $\Re \left(k\right)$ is not constant and the integral can only be solved approximately or numerically. In order to obtain an approximation of the integral in Eq.(31), the Hankel function is expanded in a power series with the first terms being

$$H_0^1\left(kr\right)\approx \sqrt{\frac{2}{\pi kr}}\exp \left[i\left(kr-\raisebox{1ex}{$\pi $}\!\left/ \!\raisebox{-1ex}{$4$}\right.\right)\right]\left[1+\frac{1}{8ikr}+\cdot \cdot \cdot \cdot \right].$$E36

Restricting the integral of Eq.(31) to the first term yields

$$\begin{array}{l}{\Phi }_{ref}\left(r,z_r,\omega \right)=\\ \text{= }\frac{S\left(\omega \right)}{4\pi }\frac{1}{\sqrt{2\pi r}}{\displaystyle \underset{-\infty }{\overset{\infty }{\int }}\Re \left(k\right)\frac{\sqrt{k}}{\gamma }\exp \left[ikr+i\gamma \left(z_r+z_s\right)\right]dk}.\end{array}$$E37

The exponential in the integrand will normally be a rapid varying function and therefore the value of the integral will be small except when the phase term of Eq.(36) is nearly constant. The phase term of Eq.(36) is

$$\alpha =i\gamma \left(z_r+z_s\right)+ikr.$$E38

The stationary points are defined to the values of the horizontal wave number k where the derivative of the phase with respect to k is equal to zero, that is where dα/dk=0, giving the stationary point as

$$r=\frac{\left(z+z_s\right)}{\tan \left({\theta }_0\right)}.$$E39

The interpretation of this result is quite simple; the reflected wave field is equal to that of the image source multiplied with the reflection coefficient at the specular angleθ ₀.

There are however situations where this approximation is not sufficient in practice. This is discussed in [13] and in the following their results are cited without proof. The accuracy of the approximation depends on the source or receiver distance from the bottom interface. The result of the analysis is that the distance z from the bottom must satisfy

$$z>>\frac{\lambda }{2\pi }\frac{\frac{{\rho }_b}{{\rho }_w}}{\sqrt{{\left(\frac{c_b}{c_w}\right)}^2-1}}.$$E40

With the water parameters of ρ_w = 1000 kg/m³ and c_w=1500 m/s, and the bottom parameters of ρ_b = 1500 kg/m³ and c_b=1700 m/s. Equation (39) requires than the distance from the bottom satisfy z>> 0.5 λ for the validity of using plane wave reflection coefficient at the bottom interface. A harder bottom with ρ_b = 1800 kg/m³ and c_b=3000 m/s, gives the requirement that z>> 1.0 λ. Hence the condition for validity is somewhat easier to satisfy for a soft bottom than for a hard bottom.

8.5. Bench marking ray modeling

The wave number integration model OASES [15] has been used to validate the accuracy and the limitation of the ray trace model using the simple case with constant water depth of 100 m and constant sound speed of 1500 m/s.

Figure 16 show the calculated transmission loss for the frequencies of 25 Hz, 50 Hz, 100 Hz and 200 Hz. The agreements between the results are very god for the higher frequency, but with some discrepancies for the lower frequencies, in particular for 25 Hz. The discrepancy is mainly a phase shift in the interference patterns of the two results, most pronounced for low frequencies and long ranges. This observation agrees with the theory outlined earlier. The seriousness of this discrepancy or errors may not very important in practice since the mean level is nearly the same as shown by the comparison with the OASES model.

Figure 16.

Comparison of the transmission loss as function of range for selected frequencies by PlaneRay (solid blue line) and OASES (dotted red line) for Pekeris’ waveguide with a homogenous solid bottom with compressional wave speed of 3000 m/s and shear wave speed 500 m/s. Both wave attenuations have the values of 0.5 dB/wavelength.

9. Case studies

In the following we present two case studies that are relevant application of the modeling techniques descried in this article. The first if these is in connection with acoustic underwater communication and the transmission of digital information. In this case the multipath communication may be a significant problem causing intersymbol interference and significant degradation of reliability and performance. The second case is related to studies on the propagation of low frequency sound and the effect such noise may affect marine life, sea mammals and fish.

9.1. Seasonal variations of communication links

In connection with a study of underwater acoustic communication the propagation over a 6 km track has been modeled for the various seasonal sound speed profiles.

The sound speed some months are shown in Figure 17. The sound speed profiles depend on the sea water temperature, the salinity and the depth. In the present case the sea water temperature variation with depth and the seasons is the main reason for changes in sound speed profile. During winters the surface water is cold and the sound speed is low, in the summer the surface water temperature and the sound speed is higher. The seasonal heating and cooling of the surface water propagates also to deeper depths, but with diminishing temperatures changes. At very large depths the water temperature is nearly the same at all seasons and the sound speed increases linearly and slowly with depth.

Figure 17.

Sound speed profiles measured at specific dates for the months given in the figures

Figure 18 shows ray tracing results are for the same profiles as displayed in Figure 17. The purpose of the study was to investigate the possibility of communication to positions beyond a sea mount and to study the multipath arrival structure as function of range and depth.

There is a seamount with a peak at about 3 km from the transmitting station. In order to simplify the interpretation ray tracings in these plots have been terminated after 6 bottom reflections, but all rays are included in the calculation of the acoustic field, but rays with so many bottom reflections, or more, will in most case not be useful for data communication because of the reflection loss and reduced coherence.

Figure 18.

Ray tracing plots assuming a source depth of 15 meter for four monthly conditions at the Roberg test site. The sound speed profiles are the same as shown in Figure 17.

Figure 19 shows examples of received time responses at 25 m depth using a Ricker pulse as source signal. The different multipath contributions are color coded for clarity. At distances from the source over 1.5 km the first arrivals is follow paths surface reflected and upward refracted paths

Figure 20 shows the channel responses at a fixed range as function of depth down to 50 m. This figure shows the total response after adding all the individual multi path contributions. The plots demonstrate that the surface channel consists of deep refracted path and a number of paths reflected from the surface and deeper upwards refractions. The stability of these paths may be uncertain and subject to rapid changes in the environmental conditions near the surface due to temperature wind and current.

Figure 19.

Time responses as function of range for receivers at depths of 25 m with a source at 15 m.

Figure 20.

Time responses as function of receiver depth at a fixed horizontal distance of 3 km from a source at 15 m depth.

9.2. Seismic noise propagation

In many areas of the world anthropogenic noise often dominates over the natural ambient noise, especially in the low frequency band from approximately from 10 Hz and upwards to 1000 Hz, or more. This frequency band coincides approximately with the frequencies of perception of sea mammals and fish and may therefore be harmful to their natural activities, or even cause physical damages. An example is the case of the seismic exploration for oil and gas in certain areas where there is important commercial fishing interest. The propagation and distribution of acoustic noise depends the environmental conditions, in particular the oceanographic parameters, the topography of the seafloor and the acoustic properties of the bottom. In this section some of examples are presented to illustrate how the environment may affect the distribution of sound and noise. This study and discussion is also relevant for passive sonar applications to detect and track submerged vehicles and objects base emitted acoustic noise

The effects of bathymetric are illustrated in Figure 21 showing ray traces of upslope and downslope conditions for typical summer conditions at the Halten Bank in the Norwegian Sea. With downslope propagation there is a thinning the ray density with distance and upslope propagation gives a concentration of rays as the water depth diminishes.

Figure 21.

The effect of up and down sloping bottoms on the acoustic field distribution calculated for the typical summer condition in the month of July.

Figure 22 and Figure 23 show the calculated sound pressure level as function of range for the downslope and upslope propagation. The sound pulse from an airgun array is modeled as s a Ricker pulse with a peak pressure of 260 dB rel. 1µPa, centered on the frequency of 50 Hz, The horizontal dashed line is the assumed threshold value for fish reaction to sound. The bottom is modeled with a 2 m thick sedimentary layer over solid rock. The sound speed in the sediment layer is 1700 m/s and the density is 1800 kg/m³. The compressional sound speed in the rock is 3000 m/s, and density is 2500 kg/m^3. The results in Figure 22 and Figure 23 are obtained under two conditions: (a) with a shear speed of 500 m/s, and (b) with no shear wave in the rock, i.e. the shear speed is zero. The absorptions are assumed to be 0.5 dB per wavelength for all the waves in the sediment layer and the rock. In the first case (a) the bottom reflection loss is as shown in Figure 9 with a significant low frequency reflection loss at angles lower than the critical angle caused by absorptions and conversion to shear wave in the bottom, which draws energy for the reflected wave. In the case of Figure 22 this results in a low-frequency and low-angle reflection loss of about 1 dB. For long ranges and many reflections this adds up to a significant total propagation loss. With no shear conversion the reflection loss is considerably reduced and the sound propagates easier to long ranges. The difference between the sound level at 50 Hz and 100 Hz is partly a result of increase attenuation at the higher frequency and partly that the source level in this case is higher for 50 Hz than for 100 Hz.

Figure 22.

Sound pressure level as function of range for downslope propagation and July conditions. Left: With shear wave conversion (500 m/s). Right: No shear wave conversion.

Figure 24 and Figure 23 show similar results for downslope and upslope propagation for typical winter conditions represented by a sound speed profile measured in the month of February. For downslope conditions the sound level decrease rapidly with increasing depth and much more rapidly with shear wave conversion (Figure 24a) than without shear (Figure 24b). With upslope propagation (Figure 23) the sound levels are near independent of shear conversion except at the very long rages where the water depth becomes constant. The examples demonstrate that sound propagation in the ocean is strongly influenced by both by the oceanographic conditions and the geophysical properties of the bottom. Reliable prediction of acoustic propagation condition requires modeling tool that can that can handle both bottom and water properties.

Figure 23.

Sound pressure level as function of range for upslope propagation and July conditions. Left: With shear wave conversion (500 m/s). Right: No shear wave conversion.

Figure 24.

Sound pressure level as function of range for downslope propagation and February conditions. Left: With shear wave conversion (500 m/s). Right: No shear wave conversion.

Figure 25.

Sound pressure level as function of range for upslope propagation and February conditions. Left: With shear wave conversion (500 m/s). Right: No shear wave conversion

10. Summary

The article has outlined the theory of ray modeling and described how the theory can be applied to study acoustic wave propagation in the ocean. The complete acoustic fields are calculated by coherent addition of the contributions of a large number of eigenrays. In this method no rays are traced into the bottom, but the bottom interaction is modeled by plane wave reflection coefficients. Ray tracing is, by definition, frequency independent and therefore the ray trajectories through the water column are valid for all frequencies. Frequency dependency is introduced by reflections from the sea surface and the bottom, including loss associated with absorption and diffuse scattering of a rough ocean and bottom interfaces. Ray tracing is therefore a computational effective method for modeling broad of frequency band wave fields and for calculation of time responses.

Ray tracing is high-frequency approximation to the solution of the wave equation and the accuracy and validity at lower frequencies may be questioned, in particular the use of plane ray reflection coefficient to represent the bottom effects. This problem has been considered both theoretically and by simulations and comparison with more accurate model. The results of this study shows that source and receiver should be at a height above the bottom of at least half a wavelength, but there is no similar requirement to the distance from the sea surface. Less fundamental is the limitation of the numerical accuracy of the determination of the eigenrays, which is most serious in the calculation of the ray amplitude and the transmission loss. These inaccuracies are of more practical nature and can be reduced by refinements in the calculations.

Examples relevant for application in acoustic underwater communication and active sonar have been presented. The propagation of low frequency sound to large distances has been presented showing the effect of the bathymetry and the acoustic properties of the bottom. An important conclusion is the effect of bathymetry and the sound speed structure interacts and that accurate modeling of sound propagation requires information about the oceanography, the bathymetry and the geology of the bottom.