Effects of Pacific Summer Water layer variations and ice cover on Beaufort Sea underwater sound ducting^a)

A combination of environmental factors in the Beaufort Sea area of the Arctic Ocean (Fig. 1) creates an anticyclonic circulation feature called the Beaufort Gyre and also a unique ducted sound propagation situation. Beneath a ubiquitous surface sound duct, a subsurface duct often exists in these waters, having a sound-speed minimum (as a function of depth) at 75–250 m depth. This duct can trap sound so that it is not strongly attenuated by surface interactions (with waves or sea ice) that can absorb sound or direct it into the seabed (Freitag et al., 2015; Webster et al., 2015; Duda, 2017). When the duct is present in its common form, the sound-speed minimum lies in a layer of relatively fresh and cool Pacific Winter Water (PWW) originating in the Bering Sea area (Pickart et al., 2005; Itoh et al., 2012). The duct is below a similarly fresh but warmer layer of Pacific Summer Water (PSW) (Pickart, 2004; Spall et al., 2018). Warmer and saltier water of Atlantic origin is beneath the duct. Figures 2 and 3^,1 illustrate the geometry of the duct, showing some measured temperature, salinity, and sound-speed profiles from the area. With duct conditions in place (suitably warm PSW), trapped sound has little interaction with surface waves and/or surface sea ice. In this paper, sounds from Canada Basin Acoustic Propagation Experiment (CANAPE) deep-water transmitters received at a shelf-edge site (Fig. 1) are presented showing that the sound duct can function effectively in the area over long time periods and also that the ducting effect can weaken or disappear for a portion of the one-year record. At times, areas with good ducting can have large horizontal extent, and sound can propagate hundreds of kilometers with essentially only cylindrical spreading loss and ionic/molecular volume absorption, as previously examined (Duda, 2017). Using information garnered from a circulation model that is guided by local conditions, although limited in ability to converge to actual oceanic conditions and not intended to fully reproduce observations, this study suggests that the variable ducting behavior is caused by time-dependent lateral temperature variations at the depth of the PSW layer, which result from straining or mixing of the PSW. Stated succinctly, the PSW layer is not continuously warm enough, over time and space, to allow perpetual ducting.

An interruption of the ducted propagation in CANAPE in the late winter and spring of 2017 has been reported by Ballard et al. (2020), who analyzed sound from positions T1–T6 (Fig. 1) that they recorded near our receiver S1. They discussed the effect in terms of westward flowing waters of the Chukchi Slope Current (CSC) and attenuation from sea ice cover. Here, we provide evidence that ice changes are not needed to explain the interruption of ducted sound propagation and that eddy and jet flow features throughout the basin, not only in the CSC, may produce defects in the duct. These defects may appear intermittently at any time, but information presented here suggests that they would be more common in spring and early summer due to annual replenishment of warm water in the late summer and fall that strengthens the ducting effect in an average sense.

Motivated by the observations of ducting variability, this paper examines physical properties of the duct with established hydrodynamical (regional ocean) modeling tools and examines ducting physics with acoustical modeling tools. Aspects of two acoustic effects in the area are highlighted, with each demonstrated in Secs. III and IV using field observations. These are (1) duct intermittency caused by eddies, filaments, and jet flow and (2) variations of the interaction of sound that is typically ducted (normal modes that are ducted) with the ice cover, with weak interaction when the duct is strong. Temperature variations of water at the depth (density) of the PSW layer are identified to be an important cause of anomalous high attenuation with respect to an efficient ducting condition. These variations are caused by stirring of waters with greater and lesser heat content and straining of gradient features by the eddy field.

The regional ocean model used here provides four-dimensional (space and time) hydrographic fields of the Beaufort and Chukchi Seas to study the effects and to allow two- or three-dimensional sound propagation simulation at different times. The computational engine is the Massachusetts Institute of Technology General Circulation Model (MITgcm) (Marshall et al., 1997). The acoustic simulations are done with two-dimensional parabolic equation (2DPE) modeling that includes seawater volume absorption and a covering layer intended to simulate the effects of ice cover. The acoustic simulation method is similar to that used in Duda (2017) with one important difference: the range-independent water column constraint applied to the environment for the simulations in Duda (2017) is relaxed for this study. Initial results from this project were presented at a conference in 2019, and the conference paper provides a succinct preliminary explanation of our findings (Duda et al., 2019). Note that the ocean simulation results shown in that paper are taken from different model output than the results shown here. The results in this paper were produced using improved initial conditions that had higher temperature PSW.

Section II describes the data collection and modeling methods used in this work. Section III describes the Pacific Water duct in the Beaufort Gyre area, defines ducting parameters, and shows how they vary in a data set. Section IV shows propagation behavior in Beaufort Sea field data and further explains the basic properties of ducted propagation in the area. Section V shows acoustic simulation results for more than one year along a single 300-km-long track in the area, with simulated rough ice cover not allowed to vary. Both seasonal variations and day-to-day variations (to a lesser degree) are examined. Section VI is a summary.

The Beaufort Sea sound propagation findings shown here are primarily model-based, with the applicability of the environmental model justified by a brief comparison with recent observations. The tools used are a regional ocean model forced by surface fluxes that includes an evolving ice-cover module, and a 2DPE computational acoustic model. The environmental conditions used in many of the 2DPE simulations are taken from the ocean model.

The regional ocean and sea ice model is an implementation of the MITgcm (Marshall et al., 1997; Losch et al., 2010) for the Chukchi Sea and the Canada Basin area. The model simulates spatiotemporal evolution of the three-dimensional ocean and sea ice field in 2016. The 2DPE code for acoustic simulations includes an ice-simulating layer above the water. This 2DPE modeling method was used in a previous study of subsurface ducted sound propagation in the Beaufort Sea (Duda, 2017). Details about these simulation methods and the field data collection methods are given in Secs. II A–II C.

The ocean-sea ice coupled model of the Chukchi Sea and Canada Basin (including the Beaufort Sea) solves the hydrostatic primitive equations, including a nonlinear equation of state (Jackett and McDougall, 1995), and the sea ice momentum and thermodynamic equations using the MITgcm suite of codes. It uses rectangular Arakawa-C grids on the horizontal plane and a fixed z-level grid in the vertical direction in the discretized space. The sea ice component uses an elastic-viscous-plastic formulation of the ice rheology (Hunke and Dukowicz, 1997) and a three-layer scheme (two layers of ice and one layer of snow) of thermodynamics (Semtner, 1976; Winton, 2000) to simulate the sea ice concentration, thickness, and motion. The ice model result does not include features such as ridges and keels because they have horizontal length scales that are smaller than the grid spacing.

The model configuration is based on the setup used for a study of the climatological flow patterns in the area (Spall et al., 2018) with some important modifications. The model has a rectangular domain of 1465 km × 2339 km with closed lateral boundaries covering a region larger than the Chukchi Sea and Canada Basin, with all derivatives computed in this domain. To insert coastlines and bathymetry into the calculation, the domain is “draped” over a spherical earth shape using a projection calculation (Fig. 4), as done by Spall et al. (2018). This domain is not aligned exactly with the Spall et al. domain. The projection stretches the domain at the edges, resulting in bathymetric feature mismatch; the errors are very small at the center of the domain where we study the conditions. Horizontal resolution of the model varies in space with the highest resolution of 1 km in both directions in a study region spanning 167°W to 147°W and 71°N to 76°N, covering the main area of interest: the shelf-edge region from Barrow Canyon to the Chukchi Plateau and waters to the north. The horizontal resolution then decreases gradually outward to 5 km on the boundaries.

There are 120 vertical layers, with the uppermost 40 layers being 4 m thick (z ∈ [−160, 0] m), the next 50 layers being 5 m thick (z ∈ [−410, −160] m), the next 16 layers being 10 m thick (z ∈ [−570, −410] m), and the layers below increasing thickness gradually. The model uses partial cells to fit the bathymetry, which is based on the International Bathymetric Chart of the Arctic Ocean (IBCAO) version 3.0 with a resolution of 30 arc seconds (Jakobsson et al., 2012). A quadratic bottom drag with a drag coefficient of 2 × 10⁻³ is used. The Smagorinsky horizontal viscosity with a nondimensional coefficient of 2.5 and the KPP vertical mixing parameterization (solving vertical viscosity and diffusivity) with background mixing coefficients of 10⁻⁵ m² s⁻¹ are used.

The model is initialized with three-dimensional temperature (T) and salinity (S) fields on 1 January 2016 from the Arctic subpolar gyre state estimate (ASTE) (Nguyen et al., 2017). The initial velocity is zero everywhere in the model. Following Spall et al. (2018), the lateral boundary conditions are no-slip and no normal flux. To preserve the subsurface water mass in the model, the T and S fields north of 77.4° N and below z = −500 m (far away from the region of interest) are restored toward the ASTE monthly mean field in 2016 with a time scale of 30 days.

To mitigate the impact of the closed eastern boundary and preserve the overall water volume in the model, a narrow channel of 60 m deep is used to connect the shelf at the eastern end to the Bering Strait, so an imposed influx at the Bering Strait is the same as the outward flux on the eastern end. The T, S, and flow imposed in the channel are designed to mimic the observed seasonal characteristics and influx of Pacific water (Weingartner et al., 2005; Woodgate et al., 2005) following Spall (2007) and Spall et al. (2018). The conditions in the channel are based on those used by Spall (2007), who applied analytical formulas of temperature and salinity spatiotemporal variation in the strait based on observations reported in Woodgate et al. (2005). To account for the recent warming of the Alaskan Coastal Current water in the Bering Strait as reported by Woodgate (2018), we have updated the temperature formula in Spall (2007) by increasing both the maximum temperature in the Alaskan Coastal Current water and the magnitude of its seasonal temperature variation from 3 °C to 4 °C. For further orientation on the latest conditions, Woodgate (2018) shows that the near-bottom Strait temperature ranges from −1.8 °C in winter months to +2.25 °C in September, and the salinity ranges from 32.95 g/kg in April to 32.05 g/kg in December. Volume transport ranges from 0.6 sverdrups (Sv) in January to 1.5 Sv in June. Both temperature and transport are up from the 1990s. The cited papers give additional suggestions for seasonally varying adjustments to these numbers to generate water column means. Effects of varying the Strait conditions on the Beaufort Sea conditions have not been examined with the exception of zero volume transport, a condition computed and reported by Spall et al. (2018).

The model is run for 3 years and forced on the surface by daily meteorological conditions in 2016 from the North American Regional Reanalysis (NARR) (Mesinger et al., 2006). The NARR 2016 condition time series is repeated in each model year. This annual cycle was the most recent available when computing time was available, and the computations were performed. The first 2 years are treated as a spin-up period over which the ocean field in the basin starts from rest and reaches a quasi-equilibrium state. The model output of temperature and salinity in the third year is examined in this study and used to provide the spatially and temporally varying environmental condition for the acoustic simulation (21 December in the second year to 15 December in the third year). Note that rather than trying to reproduce the real ocean, the hydrodynamic simulation is to produce representative spatiotemporal variations of the hydrographic condition over a year-long period to understand their impact on the acoustic propagation.

The acoustic simulations are made using a 2DPE code that includes an ice-simulating layer buoyed above the water. The code uses the Padé rational expansion method and algorithms used in the Range-dependent Acoustic Model (RAM) (Collins, 1993). RAM instructions were ported from fortran to matlab^® by Dr. Matthew Dzieciuch and then modified by us to include volume absorption and a variable-thickness layer at the surface intended to model the effects of an ice layer with a rough underside (Duda, 2017). The Padé expansion size of four is used. The fluid-type parameterized model seabed was designed to be weakly reflecting. Seabed sound speed was 1504 m/s at the seabed, increasing linearly to 1700 m/s at 100 m below the seabed and then increasing linearly to 3200 m/s at 300 m below the seabed. Seabed density was 1000 kg/m³. Attenuation was 3.2 dB/wavelength at the seabed, increasing linearly to 3.5 at 100 below the seabed and further increasing linearly to 8 dB/wavelength at 300 m below the seabed.

The ice-imitating surface layer is given a complex density, which can be adjusted to make the reflection properties of the lower side of the simulated ice layer equal to those of the lower side of an elastic layer with designated physical properties (Zhang and Tindle, 1995). The complex-density method reflection-matching method is not as accurate for this layering situation as for the half-space seabed model for which it was derived. This is because of unaccounted-for within-ice vertical interference patterns from the downward sound reflection at the ice upper surface, However, we expect no error or bias from this discrepancy in excess of the errors introduced by modeling heterogeneous ice as a single homogeneous layer. This homogeneous ice layer structure is an approximation because Arctic Ocean ice is highly variable and heterogeneous, for example, containing depth-dependent water-fraction and 3D brine channels. The ice thickness is not taken from the sea ice module of the ocean model because roughness of the ice underside (nonuniform depth at wavelengths less than a few hundred meters), which is not modeled in the module, often dominates sound attenuation by ice cover. The parameters used here for the ice model are ice compressional wave attenuation 0.5 dB/m, shear wave attenuation 2 dB/m, compressional wave speed 2500 m/s, shear wave speed 1200 m/s, and ice density of 0.910 g/ml. Using these parameters, the Zhang and Tindle (1995) expressions give a complex density of 0.356 + 0.994i g cm⁻³. The sound field within the ice layer is not properly computed with this method. For 200-Hz sound, the reflection coefficient is −4.4 dB at normal incidence, has a minimum of −12 dB at grazing angle 30°, and rises to 0 dB at grazing angle zero.

An ice surface with sloped boundaries scatters sound into the seabed and thus results in less sound in the water than a flat level surface with the same complex density. The depth of the rough underside of the ice as a function of distance x from the source, z_i(x), is calculated using either of two methods, a discreet keel method and a random short-wavelength method, or both together. These methods are described in the  Appendix. Simulation results for 200-Hz sound transmitted in a range-independent duct with various surface conditions are reported in Duda et al. (2019). Flat ice results in 1–3 dB less ducted sound with respect to open water conditions after 400 km, with the range produced by variation of PSW conditions. The addition of keels results in 6 dB loss with respect to open water for a weak duct and less loss for a strong duct.

Sound propagation data for a 10-month period from one 296-km length path in the study area are presented here. This is a subset of data collected in the CANAPE effort that took place October 2016 to November 2017. Sound was transmitted from six locations in deep water in the north. Figure 1 shows the locations of these transceiver stations T1–T6, which were moored by Scripps Institution of Oceanography. In addition to sound propagation data collected for transceiver-to-transceiver paths, sound was received by recording instruments labeled S1–S5 moored by Woods Hole Oceanographic Institution (WHOI) at the shelf edge northwest of Barrow, Alaska (southwest of the transceiver array), also for nearly 10 months (Fig. 1). The S1–S5 receivers were roughly 220–500 km from the transceiver stations. Chip-scale atomic clocks were utilized in the WHOI recorders to maintain precise time and synchronize acoustic data across the receivers.

Here, sound received level (RL) from T3 recorded at S1 is shown. Table I contains position, water depth, and path length information. RL has pressure units reported in dB re 1 μPa and is reported after match filtering to the source waveform, with normalization such that the filter-output peak level is close to the waveform root mean square (rms) level. The T3 sound source transmitted linear FM sweep (chirp) signals at 4-h intervals. The sweeps began at 225 Hz and ended at 325 Hz, with duration 135.0 s. Source level was 184 dB. The receivers were four-hydrophone Several Hydrophone Receiving Unit (SHRU) receivers sampling continuously at 3906.25 Hz. The source instruments were at 173 m depth except when high currents pulled them down. Receiver S1 had hydrophones at depths of 142.5–150 m, equally spaced.

The Beaufort Gyre circulation feature lies north of Alaska in the Canada Basin of the Arctic Ocean. The anticyclonic gyre motion arises from atmospheric forcing and impacts the pathways in the region of the relatively fresh water of Bering Strait origin (Spall et al., 2018). Observations show that a layer of warm and fresh PSW can lie above a very persistent layer of cold fresh PWW, which in turn overlies salty Atlantic Water. Above the PSW, a layer of variable cold temperature and low salinity with low sound speed is found, forming a thin surface duct. The PSW layer creates a PWW duct with a sound velocity minimum lying at 75–250 m depth. Below 350 m depth, and from the surface to the PSW, the sound speed increases with depth, which means that all sound that is not confined to the PWW duct is surface-interacting. A recent article presents a discussion of propagation in this duct (Duda, 2017). Figure 2 shows data from an ice-tethered profiler (ITP) device (unit ITP97; Toole et al., 2011) that illustrate PSW/PWW duct conditions over the course of 10 months. The instrument drifted mainly toward the southeast [Fig. 2(a)]. The recorded temperature and conductivity profiles show strong variability of the PSW [Fig. 2(b)], which affects the sound duct. Two example sound-speed profiles (profiles 82 and 259) are shown in Fig. 3. The figure includes a temperature-salinity (TS) diagram that shows how the salinity controls the stratification and how low temperature can occur at the potential density of the PSW (∼24.7 kg m⁻³) with only slight salinity changes so that the sound speed is reduced and the duct eliminated. The figure also depicts the calculation of a duct parameter that we utilize, P = DΔc, where D is height of the duct, defined as spanning from the depth of sound-speed maximum above the local minimum down to the lower depth of sound speed equal to that maximum, and Δc is the range of sound speed in the duct (upper-limit sound speed minus minimum sound speed). Figure 2 also shows a second duct parameter, P₂ = DΔc^1/2, which is motivated by P₂ being inversely proportional to the mode-n cutoff frequency of a mid-water duct computed using two turning-point conditions (no tunneling into surface duct) (Kucukosmanoglu et al., 2021). For the conditions shown here, the mode-n low frequency cutoff is f_c = K/P₂, where K = 3π (c_m/2)^3/2 (n − 1/2), which is 90 189 for c_m = 1440.5 and n = 1.

Acoustic normal modes for 200 and 950 Hz computed for ITP97 profile 82 (see Figs. 2 and 3) are shown in Fig. 5. About four modes are trapped in the duct for 200 Hz. Fifteen modes are trapped at 950 Hz. The plot shows the trapped modes to have group velocities that are all very similar and that are uniquely fast compared to the modes that occupy or intrude into the slow velocity surface duct.

The CANAPE T3 to S1 propagation observations show that sound can be ducted and travel great distances in the Beaufort Sea area with little loss from surface interaction and also that the ducted propagation path can disappear. The data confirming this are shown here, with some explanation of the data analysis methods.

Although only data from one path are shown here, sounds from the deep sources recorded at all five WHOI receivers showed ducted propagation occurring for a portion of the experiment. In addition, recordings of the deep source signals by other systems in the WHOI receiver area have been reported (Collins et al., 2019, Ballard et al., 2020). Figure 6 shows RL over the year for the T3 to S1 path, with a matched filter applied to the frequency sweep signals. The pattern of arriving energy in the vertical plot direction shows diminishing energy from November through April and increasing energy from May to August. A 50 dB or greater signal level variation over the year was seen, with the signal-to-noise ratio dropping to near unity or below in April, compromising signal detectability.

The horizontal plot direction shows RL versus time for each arriving pulse. The changes in this pattern that occur over time provide useful information about the environmental effects on propagation. The pulse arrival temporal structure can be understood by examining simulated pulses for this environment. Figure 7 shows simulated pulses made for conditions typical of the region. These typical range-dependent sound-speed conditions are taken from the ocean model described in Sec. II. The top panel shows a pulse resulting from 240-km propagation in water of 3000 m or greater depth, with no ice cover, with the sound source at 173 m depth in the PWW duct. The center panel shows the form that this pulse evolves into after propagating an additional 80 km into shallow water of 300 m depth, imitating the T3 to S1 propagation conditions. The late surface-interacting pulse section and the subsurface ducted pulse section maintain their form to a large degree, with weak or nonexistent seabed interaction. Note that the slowest mode, surface confined mode 1, is not energized [see Fig. 5(b)]. The main ducted mode arrives at ∼221.9 s. After this time, higher modes first lag as a function of increasing mode index, and then dispersion reverses at mode 7 or 8, and they move earlier in time with respect to mode index. Figure 8 shows the sound-speed conditions used for the simulations of Fig. 7. The upper 20 m contains cold water with salinity less than 29 g/kg lying immediately above a commonly found near-surface temperature maximum layer, neither of which affects sound in the way the PSW-related fast (warm) features near z = 100 m do.

With the exception of sound trapped in the PSW/PWW duct arriving at 166.4 and 221.9 s, the depicted sound energy is influenced by surface interaction. The early-arriving sound corresponds to deep-cycling ray paths or, equivalently, to high-order normal modes with turning depths (Munk et al., 1995) that are deep, below the duct. The late-arriving sound corresponds to shallow rays (or modes with shallow turning depths). The ducted arrivals that are clearly seen in the simulated pulses shown in Fig. 7 are difficult to identify in the field data presentation of Fig. 6, but they are visible in Fig. 9, which shows line plots of some of the Fig. 6 data, showing that the measured and simulated pulses are very similar in form.

Based on knowledge of the modeled pulse behavior for ducted propagation (Fig. 7), the experiment waveforms (Figs. 6 and 9) can be parsed into arriving energy for ducted sound and surface-interacting sound, respectively. This was done by finding the peak of the arriving energy and assigning it to the ducted sound. This will be more thoroughly justified by analyzing sound simulations in Sec. IV B. First, the noise level at the time of each pulse arrival was found by computing the mean RL in a 0.8-s window centered at 200.5 s. Next, the power of the ducted sound is extracted. To do this, the peak RL in the record is identified, and the mean RL is calculated in a 20-ms window centered at the peak. This is assigned as the ducted peak signal level, if sufficiently above the noise level and verified to not be an artifact at an unreasonable travel time with respect to other identified ducted peaks. The third sound level to be calculated is the power of (late) surface-interacting sound arriving after the ducted sound. This is given by the mean RL in an 80-ms window centered 220 ms after the time of the peak that was identified as the ducted arrival. Figure 10 shows time series of the noise level, the ducted peak levels, and the RL for late surface-interacting sound for transmissions with adequate signal-to-noise ratio. Travel-time data are also shown for completeness.

Two notable characteristics of the time series in Fig. 10 are that the ducted sound level and the extracted surface-interacting sound level had different behavior as the year progressed and that the ducted sound always had higher intensity than the surface-interacting sound. Both time series have short-time scale variations at periods of hours to days, which result from waves, small-scale medium variations, and advected medium variations that can each create interference pattern variations (from phase variations) and amplitude variations (Colosi, 2016). Also, changes to the surface conditions as the ice advects over the path will cause changes to the acoustic field. Looking at the remaining long-time scale variations, the ducted sound had a very steady level from 1 November to 1 February. The level then steadily dropped down to meet the noise level in early April. From mid-April to mid-July, the level rose again and then became steady at the signal strength found in the early months. The surface-interacting sound behaved differently. It did not stay steady over the first few months like the ducted sound, but instead diminished over time as one would expect as the ice cover forms and grows. After a brief period where this sound did not exceed the noise level, the level of this sound steadily rose, but more slowly than the ducted sound. The RL values of the ducted sound and the non-ducted surface-interacting sound each deviate by more than 40 dB in an apparently annual cycle. The rapid variations of the transmitted pulses are lower than typically seen in temperate zone experiments. Figure 10(c) shows rms fluctuations in 6-day-long windows for both types of sound, computed as the standard deviation of RL time series measured in dB. The RL time series entries are themselves average values computed over short travel-time windows (20 and 80 ms) described earlier. Interestingly, these are below 4 dB for lengthy periods, often below 3 dB. [5.6 dB is the full saturation benchmark, unity scintillation index (Colosi, 2016).] Of equal interest, the fluctuations rise April through June, with highest values reported for the ducted sound (actually, sound at the time of the ducted peak, which includes surface-interacting sound transiting across the duct at the receiver location). Note that averages over shorter travel-time windows are likely to have higher rms values.

The asynchronicity of the year-long pattern of the two types of energy suggests that different factors were controlling the two sound field components. The uniform level of the ducted sound from November to February suggests that the growth of patchy, thin, or weakly ridged first-year ice had no effect on this sound (the path was ice-free at mooring deployment). However, increasing sound loss due to developing ice is a sensible explanation for the diminishing non-ducted sound energy. With ducted sound level being insensitive to ice cover, the attenuation of the ducted sound after February 1 must have another cause. It is shown in Sec. IV B that intermittent alteration of the duct, specifically the cooling of the PSW layer along the track, allows sound to leave the duct by time-variable amounts, thus causing fluctuation of the ducted sound level. This is further supported by many simulation results using many different ice geometries, not all shown here. Importantly, these show that ice has no effect on the fully ducted sound energy, even a population of large keels 20–30 m in height. Sound in a range-dependent duct that has zones of weak ducting, i.e., zones with a relatively cool upper cap of PSW, may be more properly referred to as partially ducted.

The ducted propagation in the range-independent situation is straightforward. Using ray and modal sound propagation descriptions, the vertical oscillation of rays in the duct explains the lack of surface interaction, as does the alternate view of modes trapped in the duct (Fig. 5). Rays or modes not fully trapped in the duct will interact with a rough sea surface or rough ice cover, and their energy could scatter to other ray angles or other modes, including bottom-interacting modes or steeper rays. In this scenario, energy could also be lost by excitation of ice vibrations. No such loss occurs with a well-developed duct, known from our many simulations.

However, a range-dependent duct with a gap in the high-velocity PSW upper cap can allow otherwise trapped rays to escape, a process that can also be understood in terms of coupled-mode propagation for abrupt PSW structure or adiabatic mode propagation for smoother PSW structure. Energy transfer from fully ducted to surface-interacting modes would reduce the ducted energy level in situations where the surface interaction induces sound loss. A flat air-water surface would not produce loss, but other conditions are expected to send energy into the seabed or the ice, so that an interrupted duct would introduce loss of otherwise ducted signal.

The insensitivity to surface conditions of representative range-independent ducting can be demonstrated with 2DPE simulations. Simulations were run for 200-Hz center-frequency pulse propagation with 40-Hz bandwidth for a source at 173 m depth. Results are shown in Fig. 11. Three surface conditions are modeled: flat water surface (no ice), a uniform-thickness complex-density ice-like layer (flat ice), and variable-thickness ice-like layer (ice keel model 3).

Figure 11 shows that the ice severely attenuates the amount of surface-interacting sound that is present 300 km from the source. The coupled-mode propagation (scattering) caused by the range-dependent keel ice geometry fills in the quietest sections of the flat ice simulation (near 204.8 and 207.5 s). On the other hand, the peak RL of sound at 208.2 s confined in the continuous PSW/PWW duct is unchanged. This indicates that another process must have been responsible for the temporary absence of loud ducted sound (Fig. 10). The candidate attenuation-inducing process introduced earlier in this section is the escape of sound from the duct when PSW is absent or cooler than average. Under those conditions, sound can transition from ducted to surface-interacting and would then be subject to loss from two processes: reflection with substantial energy loss from the ice and conversion to high-angle sound by reflection from a sloped portion of ice that would subsequently penetrate into the seabed. Conversion by slope reflection to low-angle sound would cause rapid ice loss, process one.

The duct formed by PSW lying above PWW is not uniform over time and space, as shown by ITP data (Fig. 2 shows data from one ITP) and also suggested by our modeling (Fig. 8). The PSW layer is notably variable over time and space in the model, with consequences of sound-speed profile (duct) variability and alteration of sound ducting effects.

A time series of 200-Hz narrowband acoustic simulations along a single 300-km path was made using the ocean model output. The path is shown in Fig. 1. The path is in the region of the CANAPE deep-water to shallow-water propagation paths and has direction and length similar to some of those paths, but the acoustic simulations are not intended to reproduce the sound measured in the experiment. This is because the highly nonlinear ocean model does not incorporate data during the modeled time period (year 2016) and does not contain the same eddy field as the ocean. Therefore, the acoustic modeling shown here for conditions computed by the ocean model is intended to show expected variability of sound and typical propagation conditions and to provide explanations for the major effects. Thus, it is implicitly assumed that processes at work are steady over time so that statistics of duct interruption, as examined here, are steady from year to year. Under this assumption, observation statistics of 2017 (for example, mean TL in July) can be explained by studying the conditions in the model that has 2016 forcing. Oceanic and acoustic conditions may remain relatively unchanged in the future, but there are hints of climatic alterations in the Arctic that would change things. Note that an important difference between the model and the reality is the weakening of the Beaufort High during winter (January through March) of 2017, during the field experiment (Moore et al., 2018). That atmospheric high forces the Beaufort Gyre, a major element of the transport of Pacific water past the T3 to S1 sound path and the model sound path analyzed here.

Four snapshots of modeled sound amplitude |ψ | are shown in Fig. 12, where acoustic pressure p(x,z) = ψ(x,z) x^−1/2 exp(ik₀x −π/4) and k₀ is the background acoustic wavenumber (Jensen et al., 2011). Also, shown are the sound-speed conditions along the path and the sound level at one depth. The keel model 3 ice layer is used in all simulations, with different random small horizontal scale ice-depth variations added to the keel geometry at each time step. The sound in the duct at 300-km distance from the source varies by over 20 dB in the four realizations. The primary source of sound variability at the depth of the duct is variation of the sound speed at the top of the duct (the PSW); the sound-speed minimum in the duct (the PWW) and the water beneath the duct are stable in comparison to variations of the sound-speed maximum value above the duct. Figure 12 shows that sound escapes the duct and that the ducted sound level drops where the duct structure is anomalous, e.g., x ∼ 180 km on 18 May [Figs. 12(d)–12(f)], x ∼ 100 km on 16 July [Figs. 12(g)–12(i)], x ∼ 20, 100, and 250 km on 9 August [Figs. 12(j)–12(l)].

Both loss mechanisms mentioned previously contribute to sound energy loss in the simulations: loss via reflection from the ice with a low reflection coefficient and redirection of sound into steep angles, and thus into the seabed, by reflection from a sloping ice surface (no ice and large waves can also do this). When the warmest class of PSW is present, the ducting is strong, and when the PSW temperature is lower and the duct therefore weak or absent, sound can be strongly attenuated by these processes.

A sequence of maps showing the evolution of the PSW layer for one year in the model is shown in Fig. 13. The colors show T_u, which is the temperature at the sound-speed maximum (zero-crossing of the sound-speed vertical derivative) at the top of the duct, determined using a sound-speed gradient zero-crossing analysis algorithm. The heat content of this duct-capping feature falls to a minimum in the June to July time period. Warm water enters starting in August by moving northward from the shelf along the north shore of Alaska. This flow continues over the months that follow, with the flow diminishing in January and February and being absent or nearly so from March through July. The heat carried by this incoming PSW is continuously mixed into the surrounding waters primarily by isopycnal stirring and mixing effects, with additional influences of diapycnal mixing. The overall effect of the PSW inflow, the stirring effects of the swirling eddy field, and the isopycnal diffusion processes is the continual evolution of the patchy duct-top PSW layer that is clearly evident in Fig. 2 (data) and Figs. 8 and 12 (model output) as well as in Fig. 13, with an annual cycle of mean duct-top sound speed due almost exclusively to T_u variation. This imparts an annual cycle to the ducted propagation. Note that westward movement into the T3-S1 path of relatively cold water at the density level which marks the duct top is evident on 8 April and 26 April, consistent with the Ballard et al. (2020) hypothesis that the CSC weakens the duct, as this plotted colder duct-top water would do.

Cold filamentary anomalies as well as cold eddies would each interrupt sound ducting, enabling surface attenuation. Here, we analyze the length scales and the eddy behavior to help determine the relative importance of the two feature types and quantify certain aspects of the features. The relationship between feature size of anomalies and coupled-mode propagation in a surface duct has recently been examined (Colosi and Zinicola-Lapin, 2021), and related behavior is expected in this duct. The character of the eddies and strained filaments that lead to PSW layer coolness can be analyzed in terms of covariance computed along lines (propagation tracks). The intent is to establish correlation length scales of the duct-modifying temperature fluctuations of the PSW layer. D'Addezio et al. (2020) use standard measures to evaluate submesoscale feature vertical correlation properties, which are applied here in the horizontal. These are the covariance matrix R, standard deviation vector σ, and the correlation C, where

and

The matrix X is the model variable of interest arranged into an N_x × N_t matrix, where N_x is the number of samples in the along-track dimension and number of samples N_t in the time dimension. The term $X¯$ is the row (time) average. The standard deviation vector is the diagonal matrix of standard deviations in the horizontal. The statistics R and C were computed for a few duct measures, including sound-speed range within the duct Δc, duct parameter P, and temperature at the sound-speed maximum at the top of the duct (T_u). These measures are highly correlated with each other and give similar analysis results due to the relative uniformity of the PWW and Atlantic layers with respect to the PSW layer. Here, we present statistics of Δc for one year along the 300-km model propagation track (Figs. 1 and 12). Figure 14 shows the analyzed Δc data, R, and C. An estimate of the duct interruption correlation length is 36 km. The correlation length is larger in the south, 80 km or greater, is minimum at x ∼160 km where the very large Δc values fall away over a short distance in the early months, and relatively short in the north. The duct Δc data in Fig. 14(a) show the relative lack of a warm water layer, resulting in poor ducting, over the time of period day 150 to day 250 (19 May to 27 August).

The anomalies of duct Δc seen in Fig. 14(a) are transient in nature and are not inconsistent with a field of eddies. Eddies are known to populate the area (Spall et al., 2008). However, the maps of the related and correlated quantity T_u show jet flow and filament features as well as eddies (T_u and Δc have very similar R and C). To evaluate qualitatively the relative importance of these two feature types, the OceanEddies eddyscan software has been used to identify eddies with closed T_u contours (Faghmous et al., 2015; Frenger, 2020). The reason for identifying the feature type is that the lifetimes and behavior of eddies, and the eventual acoustic predictability, may differ from the turbulent flow seen in this area. The software is written to identify and track eddies in ocean dynamic height fields that satisfy a specific criterion, such as having a single extremum within the closed contour. The algorithm uses a geometric eddy definition as opposed to a physical (flow-based) one, for instance, use of the Okubo-Weiss parameter (Weiss, 1991), and is well suited to tabulate eddy-like anomalies in any scalar field, such as the warm layer that caps the duct. The geometric-based eddy scanning algorithms are not an infallible method for identifying and quantifying eddies, and performance information can be found in Faghmous et al. (2015) and papers cited therein.

Figure 15 shows identified eddies at two times, separated into cold and warm eddy groups. Many notable features are not identified as eddies. Only a few identified eddies are larger than 113 km² (12 km diameter if circular). On 30 December, early in the analyzed portion of the simulation, warm PSW above 0.5 °C is plentiful in the region. On 10 July, there is virtually no PSW remaining in the simulations that is above 0.5 °C; by this time, the heat in the PSW has been carried out of the study region or diffused away through along-isopycnal mixing with cooler but fresher water of the PSW density class or through vertical mixing. Determining the heat budget of PSW in the entire modeled region, and in subregions, is beyond the scope of this paper. The stirring process that precedes mixing causes spice variations (variation of T and S on isopycnal surfaces) that are tied to sound-speed variations on isopycnal surfaces (Colosi et al., 2012). There are many filamentary temperature anomaly features that are not eddies. The largest eddies are at 73°–74° N, which is possibly the penetration latitude for jet features moving northward off the shelf at ∼150° W, for both warm jets in the fall and cold in the spring. At this point, the flow interacts with the Beaufort Gyre and evidently wobbles or bifurcates into eastward and westward branches, from a transport standpoint (Spall et al., 2018). Mainly, the eddies are quite small and do not dominate the temperature (i.e., sound-speed) patterns. Note that the constraints of the ad hoc geometric method do not select important anomaly features that are eddy-like but not round (e.g., at 75.5° N, 152° W, 30 December).

Another measure of duct variability is the probability distribution function (pdf) of Δc, the duct “sound-speed range,” the difference between the sound speed at the top of the duct and the sound-speed minimum in the duct (see Fig. 3). A basic question is how steady is the ducting, which a pdf can help answer but cannot completely answer because the spatial interval between low Δc events also guides duct leakage. The pdf is computed for the Δc data along the line denoted in Fig. 15 (model grid line 550), for one year, with these data shown in Fig. 14(a) and in logarithm form in Fig. 16(a). Figure 16(b) shows the distribution of log₁₀Δc plotted in cumulative normal probability plot format. The pdf follows a lognormal distribution with two exceptions: there is a lack of extremely large values over the top 2% of samples, and the lowest 10% of values are smaller than would agree with lognormal (i.e., there is an excess of very small values, less than 1 m/s, possibly indicating a non-normal process occurring, such as mixing of two water masses (one showing no ducting), a bimodal process. The pdf of sound speed at the top of the duct, as opposed to Δc, is a closely related parameter that also fits closely to lognormal because the duct sound-speed minimum is quite stable. From a theoretical standpoint, a lognormal pdf is consistent with multiplicative processes acting together, from the central limit theorem. Candidate processes acting together in this situation are jets and eddies in multiple size classes, for example, small, medium, and large eddies, working in unison to stain the PSW layer temperature anomaly structure. The lognormal pdf has not been rigorously associated with statistics of oceanic scalar fields such as the PSW conditions examined here [Figs. 13, 15(a), and 15(b)], to our knowledge, but it has been associated with statistics of longitudinal velocity difference, a strain component (Jung and Swinney, 2005). An examination of the relationship between velocity statistics and resultant deformed scalar field statistics is well beyond the scope of this paper, although near-lognormal scalar turbulence fluctuations have been observed before (Burlaga and Lazarus, 2000).

Section IV C explained how the year-long cycle of sound ducting in the data is consistent with an interruption of the high sound speed (high heat content) in the density layer at the top of the duct, while Sec. IV D examined in more detail the model representation of this layer. To quantify the ultimate effect of this layer variability on sound using the model, this section presents a one-year-long time series simulation of 200-Hz propagation along the 300-km path (Fig. 15) where the ducted propagation has been examined (Fig. 12). The intent is to quantify how the temperature and sound-speed variations in the model (Figs. 12–14) would affect sound propagation, with ice properties not allowed to vary. The simulation is not intended to replicate the data because the MITgcm modeled fields are not constrained to be close to the observations using assimilated data; there are insufficient data available to allow such a practice, and the model is forced with 2016 surface fluxes, not the fluxes for the acoustic experiment time period. Instead, under the assumption that statistics of sound-speed conditions are steady from year to year, the effects on sound ducting of temporal changes in the local heat content and/or of the modeled stirring and mixing processes are examined.

The results are shown in Fig. 17. The ice cover is the keel ice model (keel model 3), kept constant over time, although very small ice-depth fluctuations at short wavelength are added randomly to impart expected natural variations. These are generated from an uncorrelated Gaussian random variable with a standard deviation of 0.1 m, which is bandpass filtered to retain only wavelengths between 100 and 4 m (see  Appendix). The top panel shows the time and space structure of the sound speed at the sound-speed maximum defining the top of the duct. Below that, the time-space structure of range-averaged (filtered) ducted narrowband 200-Hz sound with cylindrical spreading removed, at 180-m depth, is shown. The averaging is done with two-way application of a 16-km top-hat filter. The averaging is done to eliminate mode interference effects, and the result is representative of a broadband sound field (Harrison and Harrison, 1995). The averaged result exhibits a decline in the summer at ranges of 240–300 km from the source. The time series of the range-averaged sound level at a few distances from the source are also provided, clearly showing the behavior over one year. Specific weak-duct features are clearly responsible for large reduction in distant sound levels, on the order of 30 dB. The overall pattern has many similarities with the pattern seen in the data (Fig. 10).

Examining the time series, the simulated sound level at the depth of the duct is high and relatively constant along this track from December until April and then begins to drop off. In this winter period, variations of many dB are seen at longer distances, with less variation at short distance. Beginning in May, the sound at 240 and 300 km distance drops off and then stays consistently low, similar to what the ducted sound data show (Fig. 10, center panel, topmost line). The decline is delayed in the model by a few months compared to the data. This is not unexpected, because the ocean model used here does not assimilate data from this region. Also, a timing discrepancy is consistent with typical interannual regional fluctuations, for example, the 2017 weakening of the Beaufort High mentioned earlier (Moore et al., 2018). More plentiful and profound spatial gaps in the PSW during the low-sound period, visible in Fig. 17(a), are a major factor in the sound level decline. By August, the sound at 180 km is also reduced from its high winter level. During October, the ducting effectiveness returns to the high levels seen in the winter, although in December, a duct interruption reduces sound at distances greater than 200 km from the source.

The characteristics of order-200 Hz sound propagation in the duct created by the presence of PSW above PWW in the Beaufort Sea area of the Canada Basin have been presented using a combination of ITP data analysis, analysis of long-distance propagated sound, regional ocean hydrodynamic modeling, and underwater sound simulation. The basics of how sound behaves in the duct, which generally lies between 70 and 240 m depth, computed with range-dependent and range-independent PE modeling, were presented first. The ducting can be understood based on already published work and agrees with normal-mode analysis, as expected. A fundamental result is that sound in a continuous duct suffers little attenuation from a rough surface or from surface sea ice, even for rough ice and ice with large keels. This is caused by the character of the capping PSW layer, including its distance below the surface. Another fundamental result is that sound recorded in the Canada Basin area, although typically showing strong ducted arrivals, was not always effectively ducted and was attenuated for a period lasting many months.

To illustrate how duct variations cause changes in the energy of ducted sound, a joint study of the modeled PSW layer properties and sound propagation over a one-year period was made. The model representation of the layer shows a seasonal late-summer/fall input of anomalously warm PSW water at the duct-cap density moving north from Alaska and a late-winter/spring input of cool water in these same density and depth layers. The circulation of these waters northward and mostly westward in the Beaufort Gyre circulation pattern, combined with filamentary and eddy-filled turbulent flow patterns, gives a seasonal nature to the overall strength of the duct, measured by either duct-cap temperature T_u, duct sound-speed range Δc, or the simple duct parameter P. The details of the evolving flow pattern also impart variations to these quantities at horizontal scales of 10–100 km, giving the duct a spatially and temporally intermittent form.

The time series of acoustic simulations along a 300-km track illustrates the effect of duct gaps, and of the seasonal cycle of heat input, on long-range sound propagation in the region. The sound level time series in the simulation has a strong resemblance to a measured time series collected in the CANAPE program. The details of these two time series naturally differ because the model ocean fields differ in detail from the real ocean. However, they are similar enough to believe that the reduction of sound in the model caused by gaps of duct-capping warm PSW are a credible explanation of the variations seen in the data. A question remains about the similarity of the model ocean fields to the real ocean, necessary to believe that the processes observed in the simulation truly explain the sound measurements. To evaluate this, Fig. 18 shows the probability distributions of P and Δc for two areas of the model and ITP97 data set. The duct is more uniform in the ITP data, while the model conditions in the acoustic simulation area show a broad distribution, with a larger quantity of larger values (strong duct) and also a larger quantity of small values (weak duct). The model condition distributions for the ITP area have small peaks that are at lower P and Δc values than the dominant peaks seen for the data set, but these also show broad tails of very low values not seen in the data. It can be argued that the data set is small and is likely to depart from the more plentiful and possibly more representative model output set. On the other hand, both sets show important duct anomalies (gaps, or holes).

Work remains to more fully explore the simulations and the data set and to examine in detail (1) the escape of sound from the duct at a weak PSW anomaly, possibly using the theoretical framework of mode coupling, such as transport theory (Colosi, 2008) or a deterministic coupling examination (Preisig and Duda, 1997); (2) the surface processes that alter the sound, also possibly in terms of mode coupling or ray tracing; and (3) re-entry of sound back into the duct after passing a weak PSW anomaly, a process made non-symmetric by probable dissipation due to boundary interaction of the leaked sound.

Multiple processes contribute to the building up of the duct in this region and to its decay. The decay can be overall or concentrated at small anomalies (duct holes) that compromise acoustic ducting. Processes include time-varying transport from the Bering Sea, summer heating of northern Alaska shelf waters, northward transport of very warm PSW, advection of warm PSW within the pattern of circulation, eddy formation and strength, geostrophic shearing and filament formation isopycnal diffusion of PSW, and diapycnal diffusion PSW heat. Details for each of these may be altered by changing Arctic climate conditions, including changes in weather and ice. Future systematic changes to sound propagation in this region are difficult to predict with our current knowledge.

This work was funded by the United States Office of Naval Research (ONR) Ocean Acoustics Program, Grant Nos. N000141712624 and N000141512196. The regional ocean model was run at a United States Department of Defense High Performance Computing Facility, with computing access arranged through ONR. We acknowledge Peter Worcester and Matthew Dzieciuch of Scripps-University of California San Diego (UCSD) for generously providing essential technical information about their acoustic systems. We thank James Dunn, Peter Koski, and the captain and crew of R/V Sikuliaq and cruise Chief Scientist Mohsen Badiey for their efforts on behalf of this project. The Beaufort Gyre upward-look sonar data used to design the ice keel model were collected and made available by the Beaufort Gyre Exploration Program based at WHOI in collaboration with researchers from Fisheries and Oceans Canada at the Institute of Ocean Sciences. The parabolic model code was adapted from a matlab version of the acoustic code RAM, with the matlab version developed by M. Dzieciuch of UCSD.

The underside profile of the ice as a function of range along a transect, z_i(x) is constructed with two methods: (1) a spectral method to build stochastic fluctuations with length scales of tens of meters, between an outer scale and an inner scale, called z_si(x), and (2) a method to build a profile z_ki(x) that is composed of keels that extend beneath a flat ice underside in the space between keels. The profile z_i(x) is a linear combination of the two profile types.

In method (1), the profile is built from a normally distributed random variable with prescribed variance that is bandpass filtered to yield a smoothed z_si(x). In method (2), what is built is a collection of symmetric triangular keel objects, k_j, with j = 1:K being the keel index and K the keel quantity and with the maximum keel depth of each object d(k_j) and each keel position x(k_j) being drawn from random processes.

Implementing method (1) is straightforward. A normally distributed random variable is built, with one sample per step of the 2DPE. This is then bandpass filtered to retain the white spectrum between 4 and 100 m wavelength, with the result multiplied by the ice standard deviation (0.1 m for this work) to give the underside anomaly spatial series. A mean depth, chosen to be 1.2 m for this work, is added to the series to complete the process.

In method (2), a set of parameters controls the behavior of d(k_j) and x(k_j) and the produced keel ice underside depth z_ki(x). The aspect ratio of the keels is fixed, with keels being twice the half-width in total extent, with the half-width to depth ratio being 4.85, from the findings of Kuuliala et al. (2017). The model parameters are chosen to yield a keel field with an exponential keel depth pdf (Wadhams, 1980; Kuuliala et al., 2017) and an exponential interkeel distance pdf (Wadhams, 1980), with each pdf roughly consistent with upward-looking sea ice sensing sonar data from a mooring of the Beaufort Gyre Exploration Project (BGEP) (Woods Hole Oceanographic Institution, 2020); BGEP Mooring One data of May 2016 are used for this purpose. The time-sampled keels are converted to a distance measure by assuming an ice speed of 0.25 m/s, giving an estimated 4.3 keels/km. The method first makes an exponentially distributed collection of keel spacings, added together to yield a transect length. Then depths (thus widths) are assigned to each keel, also from an exponential distribution.

The BGEP-data keel depth pdf is not exactly exponential, and a model keel depth pdf made by fitting parameters of an exponential model to the data contains more large keels than the data. This model made in this way is called keel model 1 and is a fixed set of keels along a 300-km trackline. First, a sequence of interkeel distance is generated, δx_n, n = 1:N for N keels, chosen randomly from an exponential distribution P(δx) = μ⁻¹exp(−δx/μ), μ = 98.7, with minimum δx of 0.3 m. Next, depths d_n for each keel are chosen from a distribution estimated from the BGEP data, P(d) = μ⁻¹exp(−d/μ), μ = 3.4, and the surface is constructed, with keels not allowed to overlap, for consistency with the BGEP processing (one is deleted).

Because keel model 1 produces excess large keels, compared to BGEP data, despite the use of BGEP ice keel data, acoustic simulations reported here use keel model 3, an ice underside geometry that results when the keel model 1 keel depths are multiplied by 0.3 before the triangular keels are built. The simulations also have small stochastic ice underside variations of tens of cm added to keel model 3 geometry (method 1 is used), so that the ice underside z_i(x) = z_ki(x) + z_si(x). The major purpose of including the variable-thickness ice layer is to evaluate the sound ducting efficiency in the presence or absence of sound attenuation from rough ice bottom, less so to quantify the attenuation, which is ice parameter dependent and adjustable in the model. Figure 19 shows the keel model 3 ice underside depth. This model has keel sizes that are consistent with the Scripps and BGEP ice underside observations reported in Ballard et al. (2020). Note that the fixed keel depth-to-width ratio used here can be improved on to produce a better transect of an ice-depth profile with keels at all angles, with some keels appearing wider in an improved model. However, the model used here is sufficient for this study of ducted sound, where ice interaction and sound attenuation is a secondary effect that comes into play when the duct is interrupted.