Horizontal Refraction of Acoustic Waves in Shallow-Water Waveguides Due to an Inhomogeneous Bottom Structure

: Three-Dimensional (3-D) sound propagation in a shallow-water waveguide with a constant depth and inhomogeneous bottom is studied through numerical simulations. As a model of inhomogeneity, a transitional region between an acoustically soft and hard bottom is considered. Depth-averaged transmission loss simulations using the “horizontal rays and vertical modes” approach and mode parabolic equations demonstrate the horizontal refraction of sound in this region, even if the water column is considered homogeneous. The observed wave effect is prominent at low frequencies, at which the water depth does not exceed a few acoustic wavelengths. The obtained results within the simpliﬁed model are veriﬁed by the simulations for a real seabed structure in the Kara Sea.


Introduction
Studying the sound propagation in shallow water with different kinds of inhomogeneities is still an actual problem in underwater acoustics. The most important part of it today is the investigation of 3-D effects [1], which restrict the applicability of the N × 2-D approximation when calculating the sound pressure field. One of these effects is the horizontal refraction.
The effect of horizontal refraction is observed in shallow water waveguides with inhomogeneous seafloor or sloping bottom. It also occurs in the presence of intense internal waves or temperature fronts. This wave effect is strongly dependent on the parameters of inhomogeneities. Its manifestation is not limited by the bending of the sound propagation path, but also includes the spatial redistribution of acoustic intensity in the horizontal plane, e.g., sound focusing and defocusing. More detailed description of the horizontal refraction effect due to the bottom or water column inhomogeneities can be found in [2][3][4][5][6][7][8][9][10].
Although the horizontal refraction has been studied extensively in different environments, there is one more possible scenario that needs to be examined. The scenario suggests the presence of volume inhomogeneities in the upper layer of bottom sediments [11]. Probably, the first mention of this factor is given by Ballard et al. in [12]. However, the effect of horizontal refraction due to inhomogeneous bottom sediments appeared to be negligible in comparison with the bottom roughness factor. The conclusion was supported by the results of numerical modelling for the New Jersey shelf and cannot be extended to other regions. In particular, the areas with a rather smooth seafloor and a highly inhomogeneous seabed can be encountered on the Arctic shelf, where the bottom sediments are usually unconsolidated and have the sound speed close to that in water. In this region, new 3-D effects of sound propagation, including horizontal refraction associated solely with an inhomogeneous bottom structure, can be expected.
This paper studies the effect of horizontal refraction due to an inhomogeneous seabed structure in shallow water within numerical simulations employing the mode 1.
A simplified model, where the sound speed is only dependent on one spatial coordinate (in the horizontal plane); 2.
A realistic model, based on the velocity analysis of 3-D seismic survey data from the Kara Sea [11].
Because of the lack of knowledge on spatial distribution of bottom density and sound attenuation coefficient, these parameters are assumed to be constant in both models. Moreover, shear waves are neglected owing to the fluid nature of the sediments in the region.

Acoustic Pressure Field in an Inhomogeneous Waveguide
Complex amplitude of the 3-D pressure field P(x, y, z) in a waveguide with a rangeand depth-dependent bottom sound speed structure c b (x, y, z) is calculated by solving the mode parabolic equations. In simulations, the scenarios of a negligible mode coupling are selected. It also means that horizontal gradients of bottom parameters along the acoustic track is much smaller than those across the track. The water column is assumed to be a homogeneous layer of a constant depth H, with a sound speed c w and density ρ w .
To analyze possible 3-D effects, adiabatic normal modes are considered either separately or as a superposition. The total pressure field at a given point (x, y, z) can be written as a sum of M local normal modes Ψ l (z; x, y) [14] P(x, y, z) = The l-th mode amplitude A l (x, y) satisfies the equation where ξ l (x, y) = q l (x, y) + iα l (x, y) is the horizontal wavenumber of the l-th mode. Note that both the real and imaginary parts of eigenvalues are range-dependent.
To obtain the local horizontal wavenumbers ξ l (x, y) and mode profiles Ψ l (z; x, y), the following Sturm-Liouville problem is solved numerically Here, k w = ω/c w is the wavenumber in water, ω = 2π f is the sound source frequency, g(ξ l ; x, y) = iZ 1 (ξ l ; x, y)/ωρ w , Z 1 (ξ l ; x, y) is the total effective impedance of the bottom, which is defined by the bottom sound speed distribution c b (x, y, z), density ρ b , and attenuation coefficient β λ . To calculate the total effective impedance, one can use Equation (3.3.12) in [15] for the layered bottom.
Following the derivation by Collins et al. [16], by assuming outgoing waves only (no backscattering) and using Padé approximation of the operator square root, Equation (2) can be reduced to the parabolic equation in the form The coefficients a j,n and b j,n are selected to provide accuracy and stability constraints; where q 0 l is the real part of the l-th mode horizontal wavenumber at the sound source location. Numerical solution of Equation (4) is found by employing the split-step Padé algorithm implemented in the RAM code [17], which can be carefully configured and adapted for calculations in the horizontal plane.
The acoustic propagation direction is primarily along the x-axis. The gradients of bottom parameters are oriented mainly along the y-axis. At an initial range, the modal excitation is proportional to the amplitude of the mode at the sound source location (0, 0, z s ), A l (0, y) ∼ Ψ l (z s ; 0, 0), where z s is the source depth.
Sound attenuation in the horizontal plane is characterized by two-dimensional (2-D) transmission loss of individual modes and by 2-D depth-averaged transmission loss of the total acoustic field P(x, y, z) where x 1 = 1 m.

Numerical Simulations within a Simplified Model
Consider the idealized shallow-water waveguide shown in Figure 1, which represents a homogeneous water layer of a constant depth H over an inhomogeneous bottom halfspace. The waveguide parameters are given in Table 1.
The coefficients , and , are selected to provide accuracy and stability constraints; = 0 −2 ( 2 2 ⁄ + 2 − 0 2 ), where 0 is the real part of the l-th mode horizontal wavenumber at the sound source location. Numerical solution of Equation (4) is found by employing the split-step Padé algorithm implemented in the RAM code [17], which can be carefully configured and adapted for calculations in the horizontal plane.
The acoustic propagation direction is primarily along the x-axis. The gradients of bottom parameters are oriented mainly along the y-axis. At an initial range, the modal excitation is proportional to the amplitude of the mode at the sound source location (0,0, ), where 1 = 1 m.

Numerical Simulations within a Simplified Model
Consider the idealized shallow-water waveguide shown in Figure 1, which represents a homogeneous water layer of a constant depth over an inhomogeneous bottom halfspace. The waveguide parameters are given in Table 1.   Let us examine sound propagation in the transitional region between an acoustically soft (c b < c w ) and hard (c b > c w ) bottom, as depicted in Figure 1. The bottom density ρ b and attenuation coefficient β λ are constant. Sound speed in the bottom varies as a piecewise-linear function of y, To prevent undesirable reflections, artificial absorbing layers are placed on the sides of the model domain at |y| > 2000 m. Figure 2 shows the phase speed V ph l = ω/q l and attenuation coefficient α l of the first mode (l = 1) plotted vs. sound speed in the bottom c b for frequencies of 29, 55, 81, 108, and 134 Hz. For the purpose of illustration, we choose the bottom sound speed variation from 100 to 2500 m/s, which exceeds the interval of c b in the transitional region. Vertical dashed lines in Figure 2 demonstrate the boundaries of the region. Note that the phase speed V ph 1 has a distinct minimum between these lines. The relative "depth" of the minimum is lower at a higher frequency. By placing a low-frequency sound source in the transitional region in the vicinity of the modal phase speed minimum, one can expect the effect of horizontal refraction to be pronounced. It is also important that the modal attenuation coefficient for the soft bottom is two orders of magnitude greater than for the hard bottom. The imaginary part of ξ l can have a large value, comparable to the real part. This feature should be carefully treated when solving Equation (4). All these remarks are also true for other normal modes (l > 1).  Figure 2 shows the phase speed ℎ = ⁄ and attenuation coefficient of the first mode ( = 1) plotted vs. sound speed in the bottom for frequencies of 29, 55, 81, 108, and 134 Hz. For the purpose of illustration, we choose the bottom sound speed variation from 100 to 2500 m/s, which exceeds the interval of in the transitional region. Vertical dashed lines in Figure 2 demonstrate the boundaries of the region. Note that the phase speed 1 ℎ has a distinct minimum between these lines. The relative "depth" of the minimum is lower at a higher frequency. By placing a low-frequency sound source in the transitional region in the vicinity of the modal phase speed minimum, one can expect the effect of horizontal refraction to be pronounced. It is also important that the modal attenuation coefficient for the soft bottom is two orders of magnitude greater than for the hard bottom. The imaginary part of can have a large value, comparable to the real part. This feature should be carefully treated when solving Equation (4). All these remarks are also true for other normal modes ( > 1).  Acoustic field is calculated for a single sound source placed at the depth of z s = 15 m in the center of the transitional region, i.e., near the phase speed minimum. To demonstrate possible 3-D effects, sound propagation is simulated for a frequency of 55 Hz, at which the first mode dominates the total field. Figure 3 shows the modal phase speed V ph 1 along the y-axis and the corresponding trajectories of modal rays calculated by a ray tracing code, BELLHOP [18], in the horizontal plane. The bending of the rays is clearly seen, which reveals the effect of horizontal refraction. The maximum angle of refraction is 7.5 • . (The angle of refraction is defined as the angle between the source-receiver direction and the actual ray trajectory at the receiver.) Acoustic field is calculated for a single sound source placed at the depth of = 15 m in the center of the transitional region, i.e., near the phase speed minimum. To demonstrate possible 3-D effects, sound propagation is simulated for a frequency of 55 Hz, at which the first mode dominates the total field. Figure 3 shows the modal phase speed 1 ℎ along the y-axis and the corresponding trajectories of modal rays calculated by a ray tracing code, BELLHOP [18], in the horizontal plane. The bending of the rays is clearly seen, which reveals the effect of horizontal refraction. The maximum angle of refraction is 7.5. (The angle of refraction is defined as the angle between the source-receiver direction and the actual ray trajectory at the receiver.) Now we estimate the effect of the refraction on the acoustic intensity distribution in the horizontal plane. Figure 4a shows the calculated 2-D transmission loss at 55 Hz. (In this case, transmission loss of the first mode approximates the total transmission loss, ≈ 1 .). We notice a strong azimuthal anisotropy of , that is explained by a significant spatial variability of modal attenuation, but not by the horizontal refraction itself. To highlight the effect of horizontal refraction, the transmission loss is compared to the transmission loss 0 calculated using the N × 2-D approximation. As , 0 can also be defined by Equation (6) and Equation (1), but the modal amplitude in Equation (1) is found from Equation (4) with the complex wavenumber ( , ) = ( , ) + ( , ) replaced by ̃( , ) = 0 + ( , ), i.e., by keeping the spatial variability only in the imaginary part. In fact, it means that we use the adiabatic N × 2-D approximation, where mode coupling is neglected. (The proposed approach is valid only for the depth-averaged intensity estimates.) According to [12], the effects of mode coupling are secondary to those of horizontal refraction even at higher frequency, 300 Hz, in a ~80 m deep waveguide with a rough seafloor. In [19], the applicability of adiabatic approximation to depth-averaged transmission loss calculations is demonstrated for a single source in shallow water with a range-dependent bottom structure.
The spatial distribution of the transmission loss difference Δ = 0 − at a frequency of 55 Hz is shown in Figure 4b, which exhibits an additional insonification of the transitional region due to the horizontal refraction effect. This is explained by the fact that a curved ray, despite its longer trajectory, travels a shorter distance in the region of strong attenuation compared to a straight ray. The maximum increase in amplitude under given conditions reaches Δ = 22 dB at 10-km range from the sound source. In the multimodal regime at higher frequencies, the effect of horizontal refraction weakens for the first mode, but can still be exhibited in the acoustic field of other modes. However, the amplitudes of those modes are attenuated with range more rapidly than the amplitude of the first mode. Figure 4c shows the 2-D transmission loss at 134 Hz, at Now we estimate the effect of the refraction on the acoustic intensity distribution in the horizontal plane. Figure 4a shows the calculated 2-D transmission loss TL at 55 Hz. (In this case, transmission loss of the first mode approximates the total transmission loss, TL ≈ TL 1 ). We notice a strong azimuthal anisotropy of TL, that is explained by a significant spatial variability of modal attenuation, but not by the horizontal refraction itself. To highlight the effect of horizontal refraction, the transmission loss TL is compared to the transmission loss TL 0 calculated using the N × 2-D approximation. As TL, TL 0 can also be defined by Equation (6) and Equation (1), but the modal amplitude A l in Equation (1) is found from Equation (4) with the complex wavenumber ξ l (x, y) = q l (x, y) + iα l (x, y) replaced by ξ l (x, y) = q 0 l + iα l (x, y), i.e., by keeping the spatial variability only in the imaginary part. In fact, it means that we use the adiabatic N × 2-D approximation, where mode coupling is neglected. (The proposed approach is valid only for the depth-averaged intensity estimates.) According to [12], the effects of mode coupling are secondary to those of horizontal refraction even at higher frequency, 300 Hz, in a~80 m deep waveguide with a rough seafloor. In [19], the applicability of adiabatic approximation to depth-averaged transmission loss calculations is demonstrated for a single source in shallow water with a range-dependent bottom structure.
The spatial distribution of the transmission loss difference ∆TL = TL 0 − TL at a frequency of 55 Hz is shown in Figure 4b, which exhibits an additional insonification of the transitional region due to the horizontal refraction effect. This is explained by the fact that a curved ray, despite its longer trajectory, travels a shorter distance in the region of strong attenuation compared to a straight ray. The maximum increase in amplitude under given conditions reaches ∆TL = 22 dB at 10-km range from the sound source.
In the multimodal regime at higher frequencies, the effect of horizontal refraction weakens for the first mode, but can still be exhibited in the acoustic field of other modes. However, the amplitudes of those modes are attenuated with range more rapidly than the amplitude of the first mode. Figure 4c shows the 2-D transmission loss TL at 134 Hz, at which the acoustic field is effectively formed by the first three normal modes. Owing to the dominance of the first mode in the full field, the maximum difference in propagation loss between the three-dimensional case and the N × 2-D approximation does not exceed ∆TL = 12 dB at 10-km range (see Figure 4d). Additional calculations for other positions of the sound source along the y-axis show that the transmission loss difference ∆TL decreases to 6 dB, if the source is shifted by ±500 m relative to the initial location.
which the acoustic field is effectively formed by the first three normal modes. Owing to the dominance of the first mode in the full field, the maximum difference in propagation loss between the three-dimensional case and the N × 2-D approximation does not exceed Δ = 12 dB at 10-km range (see Figure 4d). Additional calculations for other positions of the sound source along the y-axis show that the transmission loss difference Δ decreases to 6 dB, if the source is shifted by ±500 m relative to the initial location. Note that the indicated effect of horizontal refraction, in principle, can be observed not only in the transitional region between an acoustically soft and hard bottom, but also, for example, in a region of hard bottom with a variable sound speed, which is close to the sound speed in water.

Numerical Simulations for a Real Seabed Model in the Kara Sea
The results of the previous section show the effect of horizontal refraction in an idealized waveguide of a constant depth, having the "soft bottom-hard bottom" transitional region. The key remaining question is, do such environments exist in a real ocean? Indeed, in order to observe horizontal refraction, the transitional region must be sufficiently long (at least several kilometers in length) and have a fairly regular structure.
An example of a shallow-water waveguide with an inhomogeneous bottom and almost constant water depth is presented in Figure 5, which shows the sound speed Note that the indicated effect of horizontal refraction, in principle, can be observed not only in the transitional region between an acoustically soft and hard bottom, but also, for example, in a region of hard bottom with a variable sound speed, which is close to the sound speed in water.

Numerical Simulations for a Real Seabed Model in the Kara Sea
The results of the previous section show the effect of horizontal refraction in an idealized waveguide of a constant depth, having the "soft bottom-hard bottom" transitional region. The key remaining question is, do such environments exist in a real ocean? Indeed, in order to observe horizontal refraction, the transitional region must be sufficiently long (at least several kilometers in length) and have a fairly regular structure.
An example of a shallow-water waveguide with an inhomogeneous bottom and almost constant water depth is presented in Figure 5, which shows the sound speed structure in one of the regions of the Kara Sea. This complicated range-and depth-dependent structure is obtained by the velocity analysis of 3-D seismic survey data [11]. To study three-dimensional acoustic effects, we choose a rectangular area of the bottom, marked in Figure 5a by the dashed ellipse and shown separately in Figure 5b. The selected area is characterized by strong horizontal gradients of the sound speed in the seabed. The water depth there is equal to 28 m over the whole area. structure in one of the regions of the Kara Sea. This complicated range-and depth-dependent structure is obtained by the velocity analysis of 3-D seismic survey data [11]. To study three-dimensional acoustic effects, we choose a rectangular area of the bottom, marked in Figure 5a by the dashed ellipse and shown separately in Figure 5b. The selected area is characterized by strong horizontal gradients of the sound speed in the seabed. The water depth there is equal to 28 m over the whole area. As shown in Figure 5, the sound speed in the seabed varies from 1460 m/s (blue areas) to 1700 m/s (yellow areas). It is important that there exist some extended zones, where the sound speed in the upper sediments is approximately equal to the sound speed in water ≈ = 1460 m/s. The existence of such zones can be explained by the presence of gas-saturated or water-saturated silty sediments. The bottom density is assumed to be constant and equal to 1850 kg/m 3 , which is supported by the core sample analysis results [11]. Due to the lack of information on the sound attenuation coefficient in the seabed, its value is taken equal to 0.33 dB/λ. Figure 6 shows the distributions of the phase speed ℎ and attenuation coefficient of the first mode at a frequency of 55 Hz in the selected region. One can see a strong variability of the modal phase speed in the zone of a low sound speed in the seabed, −2000 m < y < 0 m. The maximum phase speed gradient is equal to 0.06 (m/s)/m. The same zone is also characterized by the largest values of modal attenuation (up to 1.5 × 10 -3 1/m or 13 dB/km). As shown in Figure 5, the sound speed in the seabed c b varies from 1460 m/s (blue areas) to 1700 m/s (yellow areas). It is important that there exist some extended zones, where the sound speed c b in the upper sediments is approximately equal to the sound speed in water c b ≈ c w = 1460 m/s. The existence of such zones can be explained by the presence of gas-saturated or water-saturated silty sediments. The bottom density ρ b is assumed to be constant and equal to 1850 kg/m 3 , which is supported by the core sample analysis results [11]. Due to the lack of information on the sound attenuation coefficient in the seabed, its value is taken equal to 0.33 dB/λ. Figure 6 shows the distributions of the phase speed V ph l and attenuation coefficient α l of the first mode at a frequency of 55 Hz in the selected region. One can see a strong variability of the modal phase speed in the zone of a low sound speed in the seabed, −2000 m < y < 0 m. The maximum phase speed gradient is equal to 0.06 (m/s)/m. The same zone is also characterized by the largest values of modal attenuation α l (up to 1.5 × 10 -3 1/m or 13 dB/km).
The presence of a canyon-like distribution of the phase speed shown in Figure 6a, leads to the effect of horizontal refraction, which is demonstrated in Figure 7 for a 55-Hz sound source located at (0, 0, 15 m). As in the case of the idealized waveguide model, we note the curvature of modal rays in the horizontal plane (the modal rays are depicted by white lines in Figure 7a). The maximum angle of refraction is 6 • . The maximum difference in the transmission loss ∆TL between the mode parabolic equation solution and N × 2-D solution, is observed in the zone of a low sound speed in the sediments, and equals to 7 dB. If the sound source is moved along the y-axis to the point (0, −1000 m, 15 m), which is in the middle of the "canyon", the transmission loss difference increases up to 10 dB. The corresponding results at 134 Hz show the maximum difference of 4 dB. The main factor leading to the decrease of ∆TL in the real simulation with respect to the simplified model is that the modal phase speed gradient is not constant along the x-axis, and its mean value is less than that in the idealized waveguide. The presence of a canyon-like distribution of the phase speed shown in Figure 6a, leads to the effect of horizontal refraction, which is demonstrated in Figure 7 for a 55-Hz sound source located at (0,0,15 m). As in the case of the idealized waveguide model, we note the curvature of modal rays in the horizontal plane (the modal rays are depicted by white lines in Figure 7a). The maximum angle of refraction is 6°. The maximum difference in the transmission loss Δ between the mode parabolic equation solution and N × 2-D solution, is observed in the zone of a low sound speed in the sediments, and equals to 7 dB. If the sound source is moved along the y-axis to the point (0, −1000 m,15 m), which is in the middle of the "canyon", the transmission loss difference increases up to 10 dB. The corresponding results at 134 Hz show the maximum difference of 4 dB. The main factor leading to the decrease of Δ in the real simulation with respect to the simplified model is that the modal phase speed gradient is not constant along the x-axis, and its mean value is less than that in the idealized waveguide.

Discussion and Conclusions
The simulated depth-averaged intensity and transmission loss demonstrate the effect of horizontal refraction due to an inhomogeneous bottom structure, which can be significant even if the water column has a constant sound speed and depth. This effect is most pronounced at low frequencies in the zones, where the bottom sound speed is close to that

Discussion and Conclusions
The simulated depth-averaged intensity and transmission loss demonstrate the effect of horizontal refraction due to an inhomogeneous bottom structure, which can be significant even if the water column has a constant sound speed and depth. This effect is most pronounced at low frequencies in the zones, where the bottom sound speed is close to that in water.
For the realistic model of a shallow-water waveguide in the Kara Sea, the horizontal refraction at 55 Hz and 10-km range leads to the decrease in the depth averaged transmission loss by 7 dB relative to the N × 2-D approximation. In this case, the maximum angle of refraction is 6 • . This value is of the same order of magnitude as in environments with a typical sloping bottom (the slope angle is~0.5 • ) or intense internal waves (the amplitude is 10 m).
The horizontal refraction associated with an inhomogeneous bottom is characterized by the following features: 1.
The redistribution of sound intensity in the horizontal plane is determined not only by the modal phase speed structure, but also by the modal attenuation; 2.
The interference of a refracted and direct modal ray is weakly exposed due to the strong difference in the modal attenuation coefficient along the rays; 3.
The degree of horizontal refraction strongly depends on the sound source location; the maximum effect is observed for a sound source placed above the bottom, where the sound speed is close to that in water.
The results of this work can be important, for example, for more accurate estimates of safety zone boundaries introduced for marine mammals on the Arctic shelf in the areas with a high level of anthropogenic noise. Another possible application is an arrayperformance prediction when estimating source bearing and range in shallow water with inhomogeneous bottom, since their estimates can be affected by horizontal refraction.