Beaufort Sea observations of 11 to 12.5 kHz surface pulse reflections near 50 degree grazing angle from summer 2016 to summer 2017^a)

The state of sea ice in the Arctic Ocean has been of great concern for the last two decades because it has become much younger and thinner with a dramatic decline in multiyear ice (Frey et al., 2015; Krishfield et al., 2014; Kwok, 2018; Maslanik et al., 2011; Mikhalevsky et al., 2015; Stroeve and Notz, 2018; Serreze and Meier, 2019; Wadhams, 2012; Worcester et al., 2020). These fluctuations are expected to lead to significant changes in acoustic scattering from ice through changes in the sea ice acoustic properties, height and roughness spectrum (Gavrilov and Mikhalevsky, 2006), and morphological differences between the various types of ice (Wadhams and Doble, 2008). Therefore, as a remote sensing tool, acoustic ice scattering statistics offer some potential for monitoring ice evolution (Bassett et al., 2020). This paper is a step in this direction as it examines a yearlong record of the direct-path ice reflection data from the long-baseline (LBL) navigation systems of seven acoustic source and/or receiver moorings deployed in the Beaufort Sea from summer 2016 to summer 2017. These observations were part of the Canada Basin Acoustic Propagation Experiment (CANAPE; Kucukosmanoglu et al., 2021; Worcester et al., 2018; Worcester et al., 2020) and represent a novel view into the space and time scales of the single-bounce ice scattering statistics in the 11–12.5 kHz frequency range with incident grazing angles near 50°. In particular, the data are interpreted in terms of the identified temporal epochs defined as (1) open water, (2) initial ice formation lasting several weeks, (3) ice solidification also lasting several weeks, (4) ice thickening (IT) lasting roughly 6 months, and (5) ice melting lasting roughly 1.5–2 months. The observables analyzed here are arrival angle, moments of reflected intensity, its probability density function (PDF), and pulse time spread. The most variable epochs are observed to be during the transitions, which are the epoch of ice formation, ice solidification, and ice melting.

The statistics of acoustic scattering from ice as a function of the frequency and grazing angle (Duckworth et al., 2001; LePage and Schmidt, 1996) are mainly governed by the ice thickness; elastic and acoustic properties (Kuperman and Schmidt, 1986); underside and topside roughness (Hope et al., 2017); internal structure, such as salt channels, dendritic, and frazil ice structures; and keels (Jezek et al., 1990), all of which are spatially and temporally highly variable. In previous modeling studies, the variation in the ice layer acoustic properties have been shown to cause significant changes in the reflected energy (Alexander et al., 2012; McCammon and McDaniel, 1985). In the experiments, on the other hand, it has been shown that the acoustic energy loss can be attributed to compressional and shear wave attenuation (Jin et al., 1994) and bubbles within the ice (Bassett et al., 2020). In another experimental study, the loss was ascribed to variations in the temperature, salinity, porosity, and density, using the sea ice model of Laible and Rajan (1996). In addition, some experiments have found important differences in the high-frequency reflection properties due to slushy ice with its skeletal layer giving a low reflection and the consolidated thick ice giving a higher reflection (Garrison et al., 1991; Jezek et al., 1990; Stanton et al., 1986; Wen et al., 1991; Williams et al., 1992). The ice keels are yet another consideration, creating large-scale surface roughness, which can have a significant impact on the acoustic scattering from the ice (Bishop, 1989; LePage and Schmidt, 1994). Moreover, the effects from the rough sea ice surface and ice keels (which may be up to 25 m deep; Strub-Klein and Sudom, 2012), can lead to lossy out-of plane scattering and high angle scattering, where less energy is scattered in the specular direction toward the receiver (Ballard, 2019; Simon et al., 2018). This study is novel in looking at several acoustic observables over the duration of a year to try to piece together the relative effects of these aforementioned processes during the various ice epochs.

A few general comments can be made about how the geometry and environment set the scattering behavior of the measured data. We consider two dimensionless parameters that characterize the scattering from the point sources, namely, $Φ$ and Λ. $Φ$ is defined as the root mean square (RMS) phase fluctuation experienced by a ray with a grazing angle θ, impinging on a rough surface with the RMS height $⟨ h 2 ⟩ 1 / 2$⁠. Because the acoustic wavelengths are small (12–14 cm) compared to $⟨ h 2 ⟩ 1 / 2$ (which can be as large as 2 m), this means that $Φ ≫ 1$⁠. As a result, the mean pressure field will be very close to zero because the complex pressure phasor can be pointing in nearly any direction. Consequentially, the observed reflected intensity is entirely the incoherent component (Thorsos, 1984). Next, consider Λ, which characterizes the role of diffraction in the scattering process. It is defined as the ratio of the Fresnel zone radius, R_f, to the horizontal characteristic length of the rough surface, L_h. R_f is defined as the surface region around the specular point in which the paths reflecting off this area have phase differences of less than π compared to the specular ray. The paths within a Fresnel zone are considered to be smoothed by diffraction. For $Λ ≪ 1$⁠, the scattering is very geometric, that is, the sound field can sense the entire structure of the roughness. However if $Λ ≫ 1$⁠, then a lot of the roughness structure is smoothed over by diffraction. For the CANAPE geometry, R_f is of the order of several meters, whereas L_h is estimated to be tens of meters (Gavrilov and Mikhalevsky, 2006), thus, $Λ < 1$⁠. As a consequence of also having $Φ ≫ 1$⁠, the interference from the different scattering points on the ice are expected to be significant and, therefore, the acoustic field statistics should show the properties of the partial and fully saturated wave propagation regimes (Colosi, 2016). In this experimental configuration, the Kirchhoff approximation is expected to be accurate (Thorsos, 1988, 1990), and there are expressions for the reflected intensity and the scintillation index (SI), which could be evaluated under certain additional conditions (Jackson and Richardson, 2007; Yang and McDaniel, 1991), but this is not done here because the focus of this paper is the observations.

The paper is organized is as follows. Section II explains the CANAPE ice reflection observations and describes the methods of analysis. The observed surface scattering statistics, such as the RMS arrival angle, mean reflected intensity, SI, reflected intensity PDF, and pulse time spread, are presented in Sec. III. In Sec. IV, a review of some of the scattering theory is presented, and a discussion and synthesis of the results are given. Here, the focus is to identify the key acoustic scattering processes associated with the seasonal epochs and point the way to future research and modeling. The paper concludes in Sec. V with a summary of the most important results.

The CANAPE experiment was performed to assess the impact of the surface and water column processes on the low-frequency sound propagation on the varying space and time scales over the annual cycle in the Beaufort Sea (Worcester et al., 2018; Worcester et al., 2020). The fieldwork occurred between 10 September 2016 and 31 August 2017, using six transceiver moorings (T moorings with 200–300 Hz sources) and one wide aperture distributed vertical line array (DVLA) receiver mooring (Kucukosmanoglu et al., 2021). The overall experimental geometry/bathymetry and the locations for all seven moorings, arranged in a pentagon pattern, are shown in Kucukosmanoglu et al. (2021). In this section, the mooring acoustic LBL navigation systems used to obtain the ice reflection data are described, and important beampattern and response characteristics for the hydrophones, which were designed with lower frequencies in mind, are illustrated. Auxiliary information is also discussed, for example, the mooring motion behavior and high-frequency ice profiling sonars (IPS), which provide second by second estimates of the ice thickness.

The mooring LBL navigation systems consist of a 9 kHz pinger, located near the receiving array, which interrogates the near omnidirectional Benthos XT-6001-13 expendable transponders (North Falmouth, MA) that are attached to a tether 2 m above the seafloor. Then the transponders respond to the broadcasting at frequencies between 11 and 12.5 kHz with a transmit pulse width of 10 ms (±1 ms), which is then recorded on the receiving arrays. For the DVLA, the responses are recorded on 60 Scripps hydrophone modules (HMs; San Diego, CA) at depths between 50 and 590 m with a nominal 9-m depth separation, whereas for the T moorings there are only 15 HMs located at depths between 50 and 180-m, again with a 9-m separation. The DVLA had four transponders (frequencies of 11, 11.5, 12, and 12.5 kHz), and the T moorings were similarly navigated using only three transponders (frequencies of 11, 11.5, and 12 kHz). The matched filter processing was used to extract the broadband signals at each frequency with a subsequent complex demodulation to obtain the signal envelope. In this study, the signal intensity is analyzed but not the phase. Figure 1 shows an example of the geometry. The direct path is used for the mooring navigation, usually down to submeter accuracy but importantly, a surface bounce is also recorded several hundreds of ms after the direct arrival (Fig. 2). Because the goal of this system is hydrophone navigation, the transponders were situated to give roughly 45° arrival angles at the HMs, although the actual range of the angles was 50°–55°. This means that the surface grazing angles are also in this range (assuming straight-line propagation).

The HMs for this experiment were designed for the low-frequency signals (Dzieciuch, 2020) and they have idiosyncrasies at the higher frequencies, which are used by the LBL navigation system. Figure 1 shows a picture of the HMs in which the hydrophone part of the module is oriented toward the sea surface (blue section) and the electronics portion is below it (silver section with clamps). The HM is about 61-cm tall and the barrel has a diameter of roughly 6.4 cm. At low frequencies and for long-range propagation, the sound comes in mostly horizontally or broadside to the HM, but because the wavelength is much larger than the HM, it acts as a point receiver. However, for the navigation pings, the sound is coming in at steeper angles and the wavelength is comparable to the HM size. For the direct path arriving from the seafloor, the electronics portion of the HM partially blocks the signal and, therefore, a weaker arrival is recorded (Fig. 2). On the other hand, the surface-reflected arrival comes down onto the hydrophone section of the HM and, therefore, is somewhat less affected by the HM shape. To understand this beam pattern effect and quantify the variability between the different HMs, a test was performed at the Transducer Evaluation Center (TRANSDEC) in San Diego, CA. Four HMs were evaluated at the transponder frequencies. Figure 3 shows the results for the 11 kHz case, which was very similar for the other frequencies. Here, the blocking of the direct path by the electronics package is evident, and there is a flatter response for the surface-reflected path. This means that with the mooring motion (Fig. 4), the direct path will vary more than the surface paths as the angle changes: this makes the direct path intensity much less reliable. For this reason, whereas the direct arrival is used in a few places for this analysis, the focus is primarily on the variability of the surface-reflected arrival. The bottom reflections might possibly interfere with the direct and reflected paths as an irreducible noise source, explaining at least part of the variation in the direct and reflected path intensities. However, in comparison to the electronics package blocking the direct path, its impact is believed to be minor. Figure 3 also shows some significant variability phone to phone, likely due to the small wavelength of the transponder pings. This means that the analysis will be restricted to the statistics, which are only derived from a single phone. That is, the vertical correlations of the signals will not be examined.

Returning to Fig. 2, a time series of the arrival patterns for a given hydrophone and bottom transponder from the T1 mooring are displayed. The direct and reflected arrivals are evident, and the variability in the time of arrival, which is due to the mooring motion and a large pulldown from an eddy, is clear (Fig. 4). Also, the differences in the signal's intensity during the ice epochs are visible. Figure 5 shows an example of an 11 kHz time front observed on the DVLA. Open water, thickening ice, and ice melting conditions are shown and, again, the differences in the received signals are apparent for the various ice epochs.

Given the different slices that can be made through the data set, some notation is required. At a given mooring, the received intensity can be written as $I ( n , m , τ , t , r → m ( t ) )$⁠, where $n = 1 , 2 , 3 , 4$ is the transponder number, m is the hydrophone number, τ is the time in the arrival pattern, t is the geophysical time of each transmission, and $r → m ( t )$ is the hydrophone location, which changes due to the mooring motion. From these data, the intensity (proportional to the magnitude squared of the complex pressure) and arrival time are determined for the direct and surface-reflected paths, which are labelled as $I d ( n , m , t , r → m ( t ) )$ and $I r ( n , m , t , r → m ( t ) )$⁠, and $τ d ( n , m , t , r → m ( t ) )$⁠, and $τ r ( n , m , t , r → m ( t ) )$⁠. The direct path intensity is easily identified as the peak of the early arrival, which does not have any interactions with the surface boundaries. Because the path length is short, it does not experience much scattering by the water column (Fig. 2). For the surface-reflected path, the arrival pattern is more complex and, thus, the maximum intensity is chosen in the pattern. Other intensity metrics, such as integrated energy and mean energy, were experimented with and gave essentially the same results, so it was decided to apply the simple approach of using peaks.

The direct path travel times from all of the transponders, $τ d ( n , m , t , r → m ( t ) )$⁠, are important because they are used to determine the hydrophone positions $r → m ( t )$ and, in turn, correct the long-range transmission travel times for the mooring motion, leaving only the ocean sound-speed effects. The hydrophone navigation is obtained using the standard methods (Gaillard et al., 2006). The hydrophone positions are used in this analysis to compute the direct and reflected path lengths, which are written as $L d ( n , m , t , r → m ( t ) )$ and $L r ( n , m , t , r → m ( t ) )$⁠, respectively, where straight-line propagation is assumed (which is an excellent approximation). This allows a spherical spreading corrected intensity to be written as

The original intention of the work was to form an intensity reflection coefficient using the direct and reflected intensities, $R ( n , m , t ) = I r c ( n , m , t ) / I d c ( n , m , t )$⁠. However, although the time series of $I d c ( n , m , t )$ were expected to be very stable, they show considerable variability that can be attributed to two factors: (1) a strong variation in the hydrophone response as the arrival angle changes (Figs. 3 and 8) and (2) a weak out of band interference from the different transponders (Fig. 2). The phone to phone differences in the direct intensity is found to be roughly 15 dB. Therefore, the direct path was deemed unusable for the present analysis, therefore, the focus is exclusively on the reflected path. The weak out of band interference also applies to the reflected arrival and must be considered to be an irreducible noise. The key factor contributing to the reflected arrival intensity's high quality is the weak variation of the hydrophone response to small changes in the arrival angle (Fig. 8).

The signal-to-noise ratios (SNRs) for the reflected arrivals have been analyzed. The noise intensity is computed within a time window before the direct ray arrives as a function of the depth and yearday. The peak intensities at each hydrophone are divided by the noise intensities to calculate the SNR with

where $⟨ ⟩ t$⁠, and $⟨ ⟩ τ$ are the averages over t and τ, respectively, and $⟨ I ( n , m , τ , t ) ⟩ τ$ is the mean noise intensity in the time window away from the signal. The SNR statistics for the DVLA show no significant depth dependence with a mean of roughly 50 dB and RMS of 12 dB. The other T moorings show similarly high SNR results.

All of the T moorings were equipped with 420 kHz ASL Environmental Sciences IPS-5 pencil beam IPS (Victoria, British Columbia, Canada), which provided the ice draft estimates, h(t), every second. Due to the limited acoustic range of these devices, the mooring pulldowns from the eddies resulted in data gaps. For mooring T1, Fig. 6 shows the daily mean ice thickness, daily max ice thickness, roughness (i.e., RMS thickness), and, last, the percentage of leads (the percentage of time with $h ≃ 0$⁠) and percentage of ice keels (the percentage of time with $h > h ¯ + 2 σ$⁠). The percentage of ice keels is an estimate subject to our chosen threshold. For example, $h ¯$ and σ are 0.98 and 0.59 m between yeardays 404 and 405, respectively, and the keel threshold is found to be equal to 2.16 m. This approach was applied to all of the observations in a one-day window when h(t) > 0 because the ice keels range from less than 1 m to roughly 25 m deep (Strub-Klein and Sudom, 2012), and a new keel threshold is computed each day. The T1 mooring is shown here because it had suffered relatively few pull downs. The other mooring records are similar in character but with more gaps. The measurements show that the surface layer is ice-free at the beginning of the experiment. A steady rise in the ice thickness is observed from mid-October to mid-June along with a slower increase in the roughness. The maximum ice thickness exceeds 5 m after it solidifies and reaches up to 18 m during the thickening epoch. The percentage of ice keels is observed to be less than 5% of the overall ice layer during most of the year and did not significantly change during the IT epoch. More leads are observed during the epochs of the ice formation and ice melting, suggesting the highly variable surface layer structure from water to slush ice or vice versa. A linear regression coefficient (Rsq) is found to be 0.7 between the daily mean ice thickness and the roughness under the IT epoch, but as the ice cover is exposed to more solar radiation during the melting time (ice melting epoch), the roughness does not vary simply with the ice thickness (Rsq = 0.38). It is also important to note that the ice draft, h, is a positive definite quantity, and because the roughness is a large proportion of the mean (Fig. 6), the statistics are clearly not Gaussian. Figure 7 displays the satellite images of the surface cover from four separate periods of the year at the T1 mooring. The ice was fairly uniform on yearday 452 (IT) with a little more breakup on yearday 473 (IT). More leads are evident on yearday 534 (ice melting), and the ice is almost gone on yearday 594 (open water). The satellite images for the ice formation and ice solidification epochs could not be accessed.

In this section, the main results of the analysis are presented. The various statistics presented are generally computed over some small time window (usually 10 days) such that the variation of the statistic over the different ice epochs can be revealed. The first quantities presented are the arrival angle statistics, derived from an incoherent beamformer. Next, the statistics of the reflected intensity, including the mean intensity and SI are presented, followed by the PDFs of the reflected intensities. The last quantity treated is the time (τ) intensity covariance function, which gives information about the pulse spread.

An incoherent beamforming method using just intensity was applied to the observations to estimate the arrival angles (see the  Appendix). An incoherent approach was used because of the variability in the phone to phone response (Fig. 2). Figure 8 shows the time dependence of the arrival angle for the direct and surface-reflected paths, mooring T1, and the 11 kHz signal. The result here is quite typical of the other moorings. Variations in the arrival angle due to the mooring pull downs (Fig. 4) are evident and agree with the arrival angles estimated from the mooring motion. The time dependence of the RMS arrival angle for the surface path at the T1 and DVLA moorings is also shown in Fig. 8. The values change over the year from a maximum of 4.5° during the ice melting epoch to a minimum of 1° during the IT period. The RMS angle is interpreted as a metric of the wide-angled-ness of the surface scatter. The largest values are found during the open water, ice formation, and ice melting epochs, and the lowest values are found during the IT period.

The next observables to be treated are the mean and SI of the reflected intensity, $I r c ( n , m , t )$⁠, time series that show fluctuations on the time scales from hours to months. We therefore seek to interpret the slow modulation of I_rc as a changing mean while the fluctuations around that slowly changing mean will be interpreted in terms of the SI. The SI is a second moment of the intensity, defined as $( ⟨ I 2 ⟩ − ⟨ I ⟩ 2 ) / ⟨ I ⟩ 2$⁠, and quantifies, as the name implies, the degree of signal twinkling. The SI is strongly impacted by the signal interference, in this case, caused by reflections off of different regions of the rough ice surface. If the SI is much less than one, this implies a single dominant reflected path. For the SI greater than one, there are multiple paths, causing a complex interference pattern. As the number of paths goes to infinity, the central limit theorem dictates that the SI asymptotes to one (Colosi, 2016). The slow modulation component is obtained by low pass filtering, that is,

where the low pass cutoff frequency is 0.1 cycles per day (cpd). The fluctuations from which the SI are computed are given by

To get the mean reflected intensity, $I ¯ r c ( n , m , t )$ is normalized with the yearlong time mean from each phone. This gives a relative variation over the year and removes the phone to phone variation. The mean reflected intensity is then defined as the depth-frequency average of this normalized quantity,

A depth-frequency average is used because no significant variation was observed over the depth or frequency. The reflected intensity is, therefore, not a reflection coefficient; it just displays the relative variation of the reflected intensity over the year. Figure 9 shows the depth-frequency averaged reflected intensity as observed at the DVLA and T moorings, respectively: the values do not show much variation across all of the moorings. For the DVLA, the reflected intensity reaches a peak of roughly two under open water conditions. However, during ice formation, the values decrease sharply only to rise again during the ice solidification period. For the IT epoch, the reflected intensity slowly decreases from a value of roughly 1.5 to 0.5. For the early part of the ice melting phase, there continues to be some loss of the reflected intensity, reaching a minimum of 0.25, but then the reflected intensity rapidly increases as the ice disappears.

The SI also shows an interesting variation throughout the different epochs. Here, the SI is calculated along a sliding 10-day window (⁠ $δ t = 10$ days) using

where $t ̂$ is the center of the time window. Figure 10 shows the time series of the depth-frequency averaged SI as observed at the DVLA and T moorings and, again, the frequency averaging was performed because no significant frequency dependence was observed in our analysis. Figure 10 shows that no strong variation of the SI is observed across all of the moorings. The highest reflected intensity (Fig. 9) and lowest SI, with values well below one, were observed during the open water conditions. In the early phase of ice formation, the SI raises abruptly and then settles to a value near one during the IS and early IT phases. Over the IT, there is a slow increase in the SI up to a value of 1.5. The largest values of the SI are observed during the ice melting epoch, and then there is a return to the open water conditions with a low SI.

Unlike the reflected intensity, $S I ( n , m , t ̂ )$ shows a clear depth dependence. Figure 11 shows the frequency averaged SI as a function of the depth for three different ice epochs (OW, IT, and IM). There is a clear trend toward a lower scintillation at deeper depths. This effect is likely due to the depth dependence of the Fresnel zone for this geometry (Sec. IV), where deeper depths have a larger Fresnel zone, i.e., a larger region over which the diffraction averages over the rough surface. The predicted scaling is that the SI is inversely proportional to the Fresnel zone, hence, the depth scaling is found to be SI  $∝ z − 1 / 2$ (Sec. IV), which is close to what is observed.

In Sec. III B, during all of the ice conditions, it was found that the SI was between 1 and 2.5, implying that the propagation is in the partial and full saturation regimes. Here, the PDFs of the reflected intensity are examined to gain more insight into the propagation regime. Figure 12 shows the PDF of the $δ I r c ( n , m , t )$ for the OW, IT, and IM epochs on the DVLA and T1 moorings. The observed PDFs are compared to an exponential PDF, which would be expected in full saturation, and a modulated exponential (ME), which would be expected in partial saturation (Colosi, 2016). The ME distribution is defined as (Colosi et al., 2001)

where $x = I / ⟨ I ⟩$⁠, SI = 1 + 2b, $P 0 ( x ) = e − x / ⟨ I ⟩$⁠, and b, the modulation factor, should be less than about 0.5. There are two physical interpretations of the ME. One view holds that with a power-law spectrum, the small scales lead to saturated behavior (SI = 1), but the large scales in the rough interface modulate the interference making SI > 1. Another interpretation is that the modulation is a result of non-stationarity in the interface. In any case, the observed PDFs fit the ME PDF exceptionally well (b = 0.09 and 0.1 for the IT epochs on the DVLA and T1 moorings, respectively), even during the ice melting period (b = 0.46 and 0.47 on the DVLA and T1 moorings, respectively) where the SI is on the border of where the ME approximations should start to break down. For the open water cases, where SI < 1, the exponential is not a good fit, and it is estimated that the scattering regime is closer to the weak scattering case.

The pulse time spread is a statistical quantity derived from the time-lagged intensity covariance function, which can be considered physically as the inverse of the coherent bandwidth of the signal and does not depend on the signal time wander (Colosi, 2016). In the time domain, it can also be interpreted in terms of the time delayed paths relative to the specular path, which become ensonified by the effects of the rough surface. Figure 13 shows the time-lagged intensity covariance for the direct and reflected paths on yearday 321. Like the SI, an intensity covariance function is computed and then averaged over the depth, frequency, and some geophysical time window, δt, to produce a time spread estimate that is a function of the time of the year. The calculation is done as follows:

where the covariance is normalized by the zero lag value. The time-lagged intensity covariance for the reflected path is observed to be larger than that of the direct path during the entire year, as shown in Fig. 13. To get a time spread estimate, a Gaussian model is used (Colosi, 2016) of the form

where α is the time width of the unperturbed pulse and τ₀ is the time spread. The first term in this equation is simply the intensity covariance of the unperturbed signal, and the second term generates the time spread. The observations from Eq. (9) were used to estimate the unperturbed pulse time width, α, and the scattering induced time spread, τ₀ [Eq. (10)]. This is done in a two-step process. First, by forming the intensity covariance of the direct arrival [Eq. (9)], (the τ values are limited to values where there is a direct path energy), and α can be estimated by finding the e-folding scale. Knowing α, the observed intensity covariance from the reflected arrival can now be fit with Eq. (10) for the one unknown, τ₀.

Figure 14 shows a time series of the depth-frequency averaged pulse time spread, which is τ₀/α, as observed at the DVLA and T moorings. To demonstrate the depth-dependency in the pulse spread, the upper and lower DVLA HMs were analyzed independently. The spread appears to be greater at the lower depths. This impact is most likely caused by a larger Fresnel zone at the deeper depths, where a wider region is available further away from the specular point (Sec. IV). Seasonal differences in the pulse spread are seen in relation to the open water and different epochs of ice evolution. As compared to other epochs, the time spread is comparatively higher during the ice formation and then declines significantly to about 0.12 once the ice solidification phase is reached, where there is a smoother, more homogeneous surface. The spread increases slightly during the IT phase and reaches 0.2 as the ice melts. It returns to the open water values when the ice has almost disappeared.

Absent a precise ice scattering model, which will be the focus of future work, the goal of this section is to put the observations of the RMS arrival angle, reflected intensity, SI, and pulse spread into a conceptual physical context associated with the various phases of the annual ice development. This can be accomplished using simple models, existing in the literature. But before that is done, some acoustical concepts will be covered for use in the subsequent discussion.

The geometry of the scattering problem is shown in Fig. 1. The source is at the position $( 0 , 0 , − D )$ and the receiver is at $( R , 0 , − Z )$⁠. The specular point is located at $( x r , 0 , 0 )$⁠, and the rough surface is located at $( x , y , h ( x , y ) )$⁠. The specular point is given by $x r = R D / ( D + Z )$ and provides a total distance traveled by the surface interacting path as $L r = r 1 + r 2$⁠, where $r 1 = x r 2 + D 2$⁠, and $r 2 = ( R − x r ) 2 + Z 2$⁠.

The received field is due to the contributions from the entire rough surface. However, significant contributions arise only from the Fresnel zone, which is an ellipse with principle radii in the x and y directions, as discussed by Yang et al. (1992),

where $θ = sin − 1 ( D / R 2 + D 2 )$ is the approximate grazing angle of the specular ray. These radii result from a quadratic expansion of the path length around the specular point (Clay and Medwin, 1977; Yang et al., 1992; Yang and McDaniel, 1991). Here, it is seen that R_fx > R_fy and the Fresnel ellipse increases for the increasing Z and λ and decreasing θ. $R f x ≃ 1.4 R f y$ for the CANAPE geometry, where R_fy is between 3 and 10 m.

In a summary of the findings of the theoretical treatments of the rough surface scattering, Thorsos (1984) found that the effect of the roughness on the forward scatter when the incident pulse length is long and the beamwidths of the transmitter and receiver are wide was that there is essentially no alteration of the received intensity compared with the flat interface case. This was interpreted as a consequence of the energy conservation and the fact that long pulses and wide or near omnidirectional beams cannot resolve any broadening of the scattered pulse or change of the received angle. It may also be a consequence of using the quadratic expansion of the path length in the Fresnel approximation.

In the present case, short pulses were used to interrogate the rough surface, although the beamwidths of each transducer were still rather wide. With short pulses, one can observe a broadening of the received signal due to the acoustic field scattered from different parts of the Fresnel zone, which was noted in Sec. III D. The more diffuse the scattered field, the more spreading is observed. The further from the specular point, the larger the time delay and, thus, the more spread. The RMS angle from the incoherent beamformer should be linked to this increased time spread, which is explored later in the discussion. In the literature on wave propagation in random media, the pulse time spread depends on the cross-frequency coherence function (Colosi, 2016) and, in fact, the time spread parameter α can be related to the reciprocal of the signal coherent bandwidth. This interpretation may be related to the geometric interpretation presented above but requires more theoretical work to connect the two.

In Thorsos (1984), the high-frequency limit of the Kirchhoff approximation is used to model the ensemble average pulse for the forward scattering geometries. Here, the model is slightly modified to include a boundary condition that is appropriate for the water–ice interface by including the plane wave reflection coefficient, R, evaluated at the specular point. This simple multiplication by the reflection coefficient is the result of the high-frequency Kirchhoff approximation applied to the penetrable surfaces, as shown in Chap. 13 and Appendix L of Jackson and Richardson (2007). The ambiguities between the material properties and roughness exist for ice as well as the seafloor acoustic remote sensing applications. The modified version with adjustments in the notation is

where $τ = t − t s$ is the time delay from the specular time, α is the unperturbed pulse length (similar to the intensity covariance), and $erf [ z ]$ is the error function (Abramowitz and Stegun, 1972). The parameter A₀ is an area factor, which we take to be the Fresnel zone area, and I_t is the transmitted intensity at 1 m. The parameter τ₀ is an elongation time or time spread, which is due to the scattered field arriving from different places on the rough surface. Equation (13) does not take into account the time wander caused by changing the receiver locations or variations in the mean water–ice interface over the scale of a Fresnel zone. In the present analysis, the intensity covariance was used because it is insensitive to wander (Sec. III D). Be that as it may, Thorsos (1984) defines the time spread parameter τ₀ in terms of the geometry and RMS slope as

where $tan γ 0$ is the RMS slope of the rough interface, and γ₀ is the RMS slope expressed as an angle in radians. It should be noted, again, that the high-frequency limit of the Kirchhoff approximation was used to obtain Eq. (14).

The coherent reflection can also be taken into account in this model but is likely absent in these measurements. The coherent reflection coefficient is parameterized by the Rayleigh factor,

and is given by $exp ( − Φ 2 / 2 )$⁠. For the high frequencies used in this work, $Φ$ varies between 7 and 70 radians, corresponding to the RMS height of the ice interface of 0.1–1.0 m. Therefore, the exponential is extremely small, and the coherent reflection is unimportant in this experiment.

In the model detailed above, several of the parameters of the rough surface and geometry affect the observables studied in this work. The ambiguities between the material properties and roughness exist for ice as well as the seafloor acoustic remote sensing applications (Jackson and Richardson, 2007). First, the water-ice reflection coefficient, R, directly impacts the received intensity, as I is proportional to $| R | 2$⁠. This quantity is studied in Sec. IV B. Second, the rough interface does not change the peak intensity, only the pulse elongation through the RMS slope, γ₀. Therefore, in this model, the trends in the reflected intensity can be attributed only to the changing material properties of the ice, not the roughness. The effects of the time spreads can be attributed solely to the changing RMS roughness within the Fresnel zone and not the material properties. It is this separation of the effects of the material properties and ice roughness that makes the remote sensing of each parameter plausible.

However, there are limitations to the scattering model used here. The ice roughness often has an exponential roughness covariance function, which has a power-law wave number spectrum (Gavrilov and Mikhalevsky, 2006; Wadhams, 2012). This type of multiscale surface is often poorly modeled by the standard Kirchhoff approximation (Thorsos, 1990) with larger errors likely for the high-frequency limit. Therefore, at this time, it is not possible to perform remote sensing of the ice roughness using pulse spreads because a more appropriate model would need to be developed. However, as the reflected intensity is proportional to $| R | 2$ in the standard and high-frequency Kirchhoff models, remote sensing of the ice properties shows more promise, and some detailed modeling is performed in Sec. IV B, which shows a good match with the data.

Similar limitations exist with the existing models for the SI, which is a fourth moment of pressure. Assuming that the material properties of the ice are locally temporally stationary over the averaging window and spatially homogeneous over the Fresnel zone, the intensity fluctuations may be viewed as solely a function of the experiment geometry and interface roughness. Work by Yang and McDaniel (1991) and Yang et al. (1992) could be used to link the SI to the roughness covariance function and Fresnel radius. Key ideas from these two references may have an impact on any quantitative remote sensing technique using the SI.

These two models for the SI involve a double spatial integral over an exponential of the phase structure function, which is related to the roughness structure function. Here, a dimensionless number, Λ, takes into account the relationship between the Fresnel radius, roughness correlation scale, and wave number as in $Λ ∝ R f / ( k L f 2 )$⁠, and is also used in Colosi (2016) for the medium fluctuations. When Λ is small, a significant correlation exists between the contributions from the points within the Fresnel zone and can drive the SI to unity (saturation regime) or even greater than unity (Colosi, 2016; Yang et al., 1992). For small values of Λ, as the mean square phase, $Φ$⁠, increases, the SI increases from zero to slightly larger than unity and then falls back to unity as $Φ → ∞$⁠. A small value of Λ is likely for this experiment as the Fresnel radius varies between 3 and 14 m, and the ice roughness correlation length (excluding keels) ranges between 15 and 75 m (Gavrilov and Mikhalevsky, 2006). The integrals presented in Yang et al. (1992) are complicated, and specialized techniques were developed for the Gaussian roughness covariance. For the more realistic covariance functions, different integration methods must be developed, which is an opportunity for future research. Even given these limitations of the models of Yang et al. (1992) and Yang and McDaniel (1991), the overall trends in terms of the dimensionless parameters Λ and $Φ$ hold and are used for the power-law fluctuations in Colosi (2016).

For simplicity, the reflection coefficient is treated, here, using a flat three-layer system composed of a liquid ocean half-space, an elastic ice layer, and an air half-space (Brekhovskikh, 1960; Morozov and Colosi, 2017). This is called the elastic layer reflection coefficient (ELRC) model, and it is described enough in the literature that the material is not repeated here. The key ice parameters are the thickness, H, the density, ρ, the p-wave and s-wave speeds, c_p and c_s, respectively, and their attenuation factors, α_p and α_s (Hobæk and Sagen, 2016). From the measurements reported in the previous literature (Gavrilov and Mikhalevsky, 2006; Jensen et al., 2011; Kuperman and Schmidt, 1986; McCammon and McDaniel, 1985; Yang and Giellis, 1994), the historical p-wave speed and attenuation values range from 3000 to 3600 m/s and 0.07 to 0.76 dB/λ, respectively (Alexander et al., 2012; Hobæk and Sagen, 2016). For the shear wave speed and attenuation, the ranges are 1500–1800 m/s and 0.05–2.5 dB/λ, respectively. With the near disappearance of the multiyear ice, the conditions are expected to have changed considerably since these earlier measurements based on the ice-floe's history and the temperature and salinity of the ocean water at the time that the ice was formed (Hope et al., 2017). Recently, Duda et al. (2021) used c_p = 2500 m/s, c_s = 1200 m/s, $α p = 0.07$ dB/λ, $α s = 0.25$ dB/λ, and $ρ ice = 910$ kg/m³ to model the first year ice. The sound speed of the water is taken as c = 1430 m/s, and the density is ρ = 1020 kg/m³. For air, the values are c = 340 m/s and $ρ = 1.0$ kg/m³.

The CANAPE observations show that the statistics of the ice-scattered acoustic fields change considerably over the seasonal evolution of the ice and, thus, one set of ice values for the whole year is not a realistic model. Seasonally changing ice properties are, therefore, considered relative to the various epochs of the ice evolution. The first phase is ice formation and is when there can be a slushy water/ice mixture, pancake ice, and dendritic ice formations. During this phase, one would expect low shear speeds, p-wave speeds between the water and solid ice and high p-wave attenuation due to the porosity of the ice/water mixture and dendrites. The second phase, ice solidification, sees a gradual trend toward solid ice properties such as those described by Duda et al. (2021). The third and longest phase IT, conceivably, only has two parameters that change, namely, the ice thickness, H, and the roughness. The temperature and salinity structure of the ice varies over the year, also resulting in changes in the wave speeds (Laible and Rajan, 1996), but the ice porosity and thickness are reported to have a greater impact on the acoustic properties than the salinity and temperature variation within the ice (Alexander et al., 2013; McCammon and McDaniel, 1985; Yew and Weng, 1987). The last phase of ice melting is where H is diminishing rapidly, there are melt ponds and ice leads, and more slushy ice/water mixtures are expected. This period is expected to have a drop in the p-wave and s-wave speeds and an increase in the p-wave attenuation much like in the ice formation.

The ideas developed above are now used to qualitatively describe the possible processes driving the variability of the acoustic field statistics. Figure 15 shows the seasonal evolution of the reflected intensity, pulse time spread, and SI observed at mooring T1. The reader should also consider the evolution of the RMS receiver angle (Fig. 8), which is an indicator of the angular spread of the scattering. The four epochs of the ice formation, ice solidification, IT, and ice melting are marked in Fig. 15. Figure 16 shows the seasonal evolution of three quantities: (1) ice thickness, H, (2) the square of the ELRC (comparable to the reflected intensity in Fig. 15), and (3) the ice parameters, c_p, c_s, α_p, and α_s; the ice density is assumed to be constant at 910 kg/m³. The ice parameters used in the ELRC model are listed in Table I.

During the ice formation phase, the reflected intensity is seen to diminish rapidly, the spread slightly increases, and the SI is elevated relative to the pre-ice phase. The drop in the reflected intensity is attributed to the elevated p-wave loss due to the slushy ice/water mixture, and in the ELRC calculation, this behavior is mimicked qualitatively (Fig. 16). This hypothesis can be supported by previous studies, such as Williams et al. (1992), which found that the p-wave loss in ice is maximum within the skeletal layer. The p-wave critical angle is also likely to be a factor in the decreases in the reflected intensity. It is found to be equal to 51.55° for C_w = 1430 m/s and C_p = 2300 m/s. The spread and SI results could be attributed to the rough surface effects. Because the roughness is weak during this period, the increasing spread is potentially caused by the increased number of scattered paths away from the specular point. This hypothesis is supported by the slight increase in the RMS angle during this period (Fig. 8). Similarly the SI is elevated by a factor of 2 as the result of more interference from scattered paths.

At the start of the ice solidification phase, there is a notable inflection point for all three observables. The reflected intensity is seen to increase rapidly as the ice firms up. The hypothesis is that this increase is because of the reduction in the p-wave attenuation. Here, the loss due to the shear wave production is likely not significant as the ice is still thin and solidifying. Again, the ELRC captures this increase adequately (Fig. 16). The model yields a decent performance in this epoch when C_s = 1200 m/s; however, C_s is assumed to increase from 1200 to 1300 m/s and C_p is expected to increase from 2300 to 2500 m/s during this phase, which should improve the correlation between the model results and the observed reflected intensity. The formation of colder, less porous ice might be responsible for the rise in C_s and C_p during this phase. The spread and SI are both dropping considerably. During the solidification process the roughness is expected to increase, but a reduction is seen in the RMS arrival angle (Fig. 8). The reduction in the RMS angle is consistent with the reduction in the time spread as these variable angles could be associated with paths away from the specular path. A full-year time series of a 10-day average RMS reflected angle and spread with a sliding window of 1 day were evaluated using the linear regression analysis to assess the relationship between them. In line with our predictions, the spread was positively related to the RMS reflected angle with Rsq = 0.68, 0.66, and 0.57 for 12 kHz, 11.5 kHz, and 11 kHz, respectively. For the SI, there is a trend toward one, and because the PDFs are close to exponential (Fig. 12), the acoustic fields are apparently tending toward saturation. Saturation could occur as a result of the increased ice roughness, leading to more scattered paths, but those scattered paths would have to be close enough to the specular path so that they would not contribute to the time spread.

When the IT phase begins, there is, again, an inflection point. The reflected intensity is seen to slowly diminish, whereas the spread and SI slowly increase. The slow decrease in the reflected intensity is likely because of the slow thickening of the ice and, therefore, increasing the coupling into the ice shear waves (Fig. 15) and attenuation, which the wave experiences when it propagates through the ice layer and back. There is also likely some decrease in the reflected intensity due to the increasing ice roughness and ice keels. The reflected intensity starts at a value of about 1.5 and reaches a minimum shortly after the ice maximum of about 0.15. This pattern is mimicked in the computed reflection coefficient with no roughness contribution (Fig. 16). Interestingly, the time spread is rather insensitive to the thickening and the evolution of the ice over this phase, and the low values are, once more, strongly correlated with the low RMS angle values (Fig. 3; r = 0.75) and partially correlated with the percentage of ice keels during the IT epoch (Fig. 6; r = 0.44). The lossy out-of plane scattering and high angle scattering are possibly caused by ice keels during this epoch (Ballard, 2019; Simon et al., 2018). Apparently, the time delayed paths that give rise to the spread remain attenuated in energy over this time. The SI, on the other hand, shows a significant increase over the ice growth phase with a trend away from full saturation (SI = 1) to a more partial saturation, where correlated paths are constructively interfering to give high intensity glints. Again, because the spread is small, these paths, contributing to the scintillation, must be close to the specular point. The rise in the SI could also mean that there is more ice roughness non-stationarity because the PDFs match the ME fairly well. In this analysis, we cannot discriminate between the partial saturation behavior and non-stationarity.

The last epoch is the ice melting phase, which for the reflected intensity and spread, show a gradual rebound to the open water conditions, mostly tracking the diminishing ice thickness (Figs. 15 and 16). The reflected intensity continues to decrease as the ice melts in the early part of the ice melting phase. This decrease can be attributed to the decline in the p-wave and s-wave speeds and rise in the p-wave attenuation due to the melt ponds and slushy ice/water combinations, but it may also be related to the brine channel melting and ice becoming more porous as a result of warm water entering underneath the ice (Jackson, 1994). The behavior of the SI is much more interesting. Shortly after the ice maximum, the SI rises abruptly to a strong peak of 3.5–4.5, corresponding roughly to the minimum in the reflected intensity. From this maximum, the SI then declines steadily to the open water conditions. The cause of these high values for the SI is unknown, but it is hypothesized that keels, leads, and melt ponds may play an important role such that the signals are seeing a combination of ice and open water conditions (Fig. 6).

In this work, the sea ice scattering statistics of 11–12.5 kHz acoustic transmissions in the Beaufort Sea region of the Arctic Ocean are quantified over an annual cycle, and five significant surface scattering epochs are defined by shifts in the RMS arrival angle, reflected intensity, spread, and SI. The observations show apparent changes in the quantities with the varying time scales in each epoch. The reflected intensity is highest during the open water period and lowest during the early ice melting period. It varies significantly during the ice formation, ice solidification, and ice melting epochs but steadily decreases during the IT epoch. Unlike the reflected intensity, the SI is lowest during the open water period with values slightly below one and highest during the early ice melting period. Again, the SI varies dramatically during the ice formation, ice solidification, and ice melting epochs. It increases rapidly before settling to a value near one during the ice solidification and gradually evolved to a value of 1.5 in the later part of the IT phase. The SI and spread, unlike the reflected intensity, are observed to be both time and depth dependent.

Although the increased roughness results in a longer time delay under the Gaussian surface layer, the spread is found to be relatively lower under the IT than under the the ice formation and ice solidification. The hypothesis, here, is that two different mechanisms may have been responsible for the observed variations in the spread. First, increasing the square RMS phase and decreasing the isotropic correlation length reduce the contributions to the mean square pressure by driving the structure function's exponential to zero. Second, the increased roughness and ice keels lead to out-of plane scattering and high angle scattering, causing less energy to be scattered in the direction of the receiver. The lower values in the RMS angle during the IT phase support this hypothesis as only the area near the specular point may scatter the rays into the receiver. What to conclude from the pulse spread in the changing ice conditions remains an open question.

The ELRC model is used to estimate the ice parameters, including the p-wave and s-wave attenuations and speeds, and describe the possible processes controlling these variations qualitatively. There are two minimums in the ELRC: one in the ice formation and one in the ice melting. The elevated p-wave loss caused by the slushy ice/water mixture may be responsible for the decreases during the ice formation and ice melting. The results of the two processes in the squared ELRC are similar to the fluctuations observed in the data (r = 0.78). The reflected intensity is seen to increase rapidly as the ice solidifies. The present analysis suggests that the observed rise in the reflected intensity is caused by the reduction of the p-wave attenuation and the increase in the p-wave and s-wave speeds. The reflected intensity decreases steadily as the ice thickness increases. During the ice melting, the melt ponds and slushy ice/water mixtures are predicted to raise the p-wave attenuation and decrease the p-wave and s-wave speeds.

In summary, the sea ice scattering statistics in the Arctic Ocean were quantified in this study with the aim of using this knowledge to monitor the ice properties. Five distinct ice epochs were defined, and the ice parameters were predicted as a function of the time for each epoch. The conclusion is that the observed changes in the reflected intensity result from a combination of processes involving the ice composition, thickness, roughness, and ice keels at the ice/water interface. A future study will focus on developing a physics-based quantitative model, which links the scattering statistics to the specific ice and surface characteristics. The experimental arrangement needs to be improved to make the results less prone to error. Namely, better broadband sources using coded signals would reduce the intertransponder interference and allow the more frequent transmission to examine the temporal statistics. Furthermore, appropriate receivers designed to record mid-to high-frequency signals are necessary to quantify the signals' vertical statistics. These advancements in modeling/theory and observational capability will be prerequisite to fulling the promise for acoustically observing the changing ice conditions across the Arctic.

This research was supported by the Office of Naval Research (ONR) and M.K. was supported by an ONR Ocean Acoustics Graduate Student Fellowship under Award No. N00014-19-1-2203. The IPS ice draft measurements were funded by the ONR Arctic and Global Prediction Program (ONR 322AG) under Award No. N00014-15-1-2782 to the Woods Hole Oceanographic Institution. This material is based on work supported by the ONR under Award No. N00014-15-2068 to the Scripps Institution of Oceanography. We thank Richard Krishfield and Andrey Proshutinsky of the Woods Hole Oceanographic Institution for generously providing the essential IPS ice draft measurements. Any opinions, findings, and conclusions or recommendations expressed in this publication are those of the authors and do not necessarily reflect the views of the ONR.

The angular dependent time delay between the reference hydrophone and other hydrophones, $δ τ ( m , θ b )$⁠, is explained here (Fig. 17). The output of the beamformer for the reflected and direct arrivals is defined as

where N is the total number of hydrophones, $δ τ ( m , θ b )$ is the angular dependent time delay, $L d , r$ is the direct and reflected path distance, and θ_b is the bearing angle. The arrival angle, $θ d , r ( n , t )$⁠, relative to the mooring is then determined by the maximum of $I beam ( n , θ b , τ , t )$ for each n and t (Fig. 4). Last, $θ d , r ( n , t )$ is corrected for the mooring tilt.

It is assumed that the acoustic waves are propagating along the array with the sound speed at the reference hydrophone c₀. Using the reference hydrophone location, $r → m r ( t )$⁠, and other hydrophone location, $r → m ( t )$⁠, the distance and angle between the hydrophones are given, respectively, by

where $r m ( t ) = ( x m ( t ) , y m ( t ) )$ and $r m r ( t ) = ( x m r ( t ) , y m r ( t ) )$⁠. The direction of propagation is perpendicular to the wave front, and the distance between the same wave fronts on the different hydrophones, l_i, can be used to compute the time delay by $δ τ ( m , θ b )$= $l i ( m , θ b ) / c 0 ( m r , t )$⁠. Writing $l i ( m , θ b ) = R i ( m , t ) cos ( θ i − θ b )$⁠, it can further be written as

Figure 18 shows the beamformer output of the reflected path for the yearday 510.