To estimate SL we applied the sonar equation for an omnidirectional source and receiver with no processing gain (Urick, 1983; Lurton, 2010) such that:

here RL is the received sound pressure level of the sound measured at the sonobuoy, and TL is the transmission loss that occurs as the sound propagates through the environment from the whale to the sonobuoy. Our method for estimating SL was to make calibrated measurements of RL from calls originating at a known location, and to estimate TL for each detection of that call. In the sections below, we present a detailed description of the methods used to measure RL, triangulate call location, and estimate TL.

Call locations were obtained via triangulation of bearings from directional sonobuoys (detailed description below) which required detection of the same call on two or more different sonobuoys deployed concurrently. The mean SL for each triangulated call was obtained from all of the detections that comprised that triangulation. Thus, our sample size, n, for each call type was the number of calls, not the number of detections.

The detection of the same call at two or three hydrophones allows for the estimation of the variance due to various measurement errors separately from the variance due to whale depth and the natural variability in actual source level.

For a given call i detected at two hydrophones 1 and 2 we obtain two estimates of SL

where ∈ is an error term associated with the measurement error at each hydrophone (error in received level measurement plus error in transmission loss). To model the errors, we assume that ∈ are drawn from the same normal distribution for each measurement with mean 0 and standard deviation σ_∈. We then estimate the standard deviation of ∈ based on the n measurements as

Our estimate of σ^∈ can then be compared to the standard deviation of SL across all calls, σ. If measurement error is much smaller (i.e., σ^∈≪σ), then that would indicate that most of σ is due to natural variability. However, if σ^∈ is of similar magnitude to σ, then most of the variability in our estimate of SL would appear to arise from measurement error and little can be learned about natural variability.

Listening stations during ENRICH survey were conducted by deploying SSQ-955 HIDAR sonobuoys in DIFAR mode to monitor for, and measure bearings to, vocalising whales while the ship was underway (Gedamke and Robinson, 2010; Miller et al., 2015). For most of the survey, listening stations were conducted with a single sonobuoy for one hour every 30–60 nmi in water depths greater than 200 m following a protocol nearly identical to Gedamke and Robinson (2010). However, some portions of the voyage were dedicated to whale studies or fine-scale krill surveys. During 205 h of these periods, arrays of two or three sonobuoys were deployed and recorded concurrently in order to triangulate calls following the protocol described by Miller et al. (2015).

The sonobuoy recording hardware, software, and sound-processing steps were identical to other recent voyages in 2013, 2015, 2018, and 2020 (Miller et al., 2015, 2016; Calderan et al., 2020; Jackson et al., 2020), and are also described in detail in the ENRICH voyage project report and passive acoustics report (Double, 2019; Miller B. S. et al., 2019). In brief, radio signals were received on an omnidirectional or Yagi directional antenna connected to a VHF receiver. Audio from the VHF receiver was digitized and recorded via a sound-board connected to a desktop computer.

Compass calibration followed procedures described in Miller et al. (2015, 2016). Calibrated sound pressure levels, for each detection and for measurements of noise levels, were obtained following the “intensity calibration” method described by Rankin et al. (2019). This frequency-domain method involved multiplying the frequency calibration-curves for the omnidirectional hydrophone given in the DIFAR specification (e.g., see Maranda, 2001; Greene et al., 2004), the radio receiver output level (nominally flat in frequency and 1.0 ± 0.2 V rms @ 75 kHz for a WinRadio G39-WSBe), and sound board (flat in frequency and 8.4 volts peak-peak limits for the analog instrument input with gain set to 20 dB).

Aural monitoring and visual inspection of spectrograms was conducted by an analyst (SC, BM, EM, AŠ, or KS) during the voyage. Monitoring was conducted in real-time, and the intensity scale of the spectrogram was adjusted by the operator to suit the ambient noise conditions. The spectrogram used to visualize sounds from ABW and fin whales had a 250 Hz sample rate, 256 sample Fast Fourier Transform, and advanced 32 samples between time slices.

Antarctic blue whale and fin whale sounds were detected and classified manually by the analysts, and the PAMGuard DIFAR module (Miller et al., 2016) was used to measure the received level of and bearing to detections (in dB re 1 μPa and degrees from true north, respectively). Fin whale calls were classified as either 20 Hz pulses or 40 Hz downsweeps as described in Miller et al. (2021a). Antarctic blue whales calls were classified into three types, stand-alone unit-A, Z-call, and D-calls, based on the description by Rankin et al. (2005) and using the naming conventions in Miller et al. (2021a). Calls designated stand-alone unit-A consisted of only the 25/26 Hz first tonal unit of ABW song with no other units detected, whereas Z-calls consisted of the full three-unit vocalization.

All detections used in this manuscript were those made manually in real-time during the voyage. However, precise start and stop times for all detections and upper and lower frequency bounds of 20, 40 Hz and D-calls were determined onshore post-voyage by visual inspection of the spectrogram by an analyst (BSM). Received levels, in dB re 1 μPa, were measured over the full duration of each detection in the frequency domain by applying Parseval’s theorem and DIFAR and recording-chain calibration factors following the method of Rankin et al. (2019). For stereotyped calls (Z-calls, unit-A, and 20 Hz), RLs were all measured in the same fixed frequency band for that call type. For D-calls and 40 Hz calls, RLs were measured over the band defined by the analyst for each call. For Z-calls, the detection start and stop times were restricted to the full duration of unit-A, and the frequency bounds were fixed to the band from 24 to 28 Hz. Frequency bounds of stand-alone unit-A were also fixed 24–28 Hz band. Though unit-A was narrowband, the use of a fixed 4-Hz band meant that our measurements of Antarctic blue whale unit-A should be directly comparable to those made in other studies (Širović et al., 2007; Bouffaut et al., 2021). RLs of fin whale 20 Hz pulses were measured over the band from 15 to 25 Hz, again for compatibility with estimates of SL of this call type from the Antarctic (Širović et al., 2007) and other oceans (Charif et al., 2002; Weirathmueller et al., 2013; Miksis-Olds et al., 2019). The RLs of D-calls and 40 Hz calls were measured over the full bandwidth of each call, as defined by the analyst’s frequency bounds. Figure 1 shows examples of each of these call types and the analyst detection boundaries.

Georeferenced locations of whale calls (latitude and longitude) were obtained by performing additional analysis on the bearings obtained during the voyage from the PAMGuard-DIFAR module (Miller et al., 2016). This post-voyage analysis was performed in Matlab 2019b (Mathworks Inc., Natick, MA, United States). A combined detection-matching and localisation algorithm was used to estimate the call location.

For each detection, all detections on other sonobuoys that could potentially correspond to the same call were identified to determine which, if any, were the best match. A detection was considered a potential match if it was of the same sound type, its arrival occurred within the maximum time-differences-of-arrival (TDOAmax) for that pair of sonobuoys, and that detection had not previously been determined to match a different detection. TDOAmax was calculated from the known deployment locations of sonobuoys using a constant speed of sound (1400 m/s based on hydrographic sampling during the voyage). For recordings from pairs of sonobuoys, potential matches included all of the detections that arrived within TDOAmax on the second sonobuoy. For triplets, potential matches included all of the possible combinations of pairs of matches that arrived within TDOAmax on all three sonobuoys. For each set of candidate matches, the triangulated location was calculated via Lenth’s maximum likelihood method (Lenth, 1981). To apply Lenth’s method, latitude and longitudes were converted to cartesian coordinates (universal transverse Mercator projection) and triangulated positions converted back into latitudes and longitudes using the Matlab mapping package M_Map (Pawlowicz, 2019).

The TDOAs for each candidate location were then modelled assuming a constant sound speed of 1400 m/s. The candidate location where the modelled TDOAs best matched those calculated from the measured arrival times of each detection was accepted as the best matching location. Locations were only considered potentially acceptable if the mismatch between modelled and measured TDOA was less than 30% of the total maximum TDOA [i.e., Σ(TDOAmodelled − TDOAmeasured) ≤ 0.3 ΣTDOAmax]. The threshold for exclusion of 30% of the total TDOAmax was chosen to reject grossly mismatched detections whilst accommodating the three main sources of errors in TDOA: (1) sonobuoys that had drifted away from their deployment locations, (2) variability in the analysts’ marking of the time bounds of each detection, and (3) lack of coherence between detections arising from differing propagation conditions. Since triplets of sonobuoys were always likely to have a worse total match than any one of the three pairs of sonobuoys they contained, potential matches from triplets of sonobuoys were all calculated first, and pairs were only investigated if no triplet candidates were above the exclusion threshold.

Additional filters were applied to exclude poor triangulations and facilitate accurate estimation of source levels. Triangulations where whales or sonobuoys were in waters shallower than 1000 m were excluded due to our inability to model these TL accurately (see next section). Furthermore, we excluded triangulations that were greater than 50 km from the sonobuoy as these were likely to be inaccurate since they were more than three times further than our typical sonobuoy separation of 15–20 km. On the other end of the distance scale, we also excluded triangulations that were within 2 km of a sonobuoy, since these were more likely produce larger errors in TL and SL resulting from sonobuoys drifting away from their deployment locations.

Lastly, we considered the maximum angle of intersection, θ, of the bearings of each triangulation. Assuming that errors in bearings are symmetric around the true mean, then distant triangulations from bearings that are closer to parallel are more likely to be biased further away than those that are at right angles. Thus, SL that include these near-parallel bearings could be biased high. To check whether this occurred in our dataset, we investigate SL for an additional subset of calls where the angle of intersection of bearings, θ, was within 45° of a right angle (i.e., sin θ < 0.71). This subset of triangulations should be comparatively free of any bias that may arise from the limited angular precision of sonobuoys.

To reduce uncertainty in our modelling of TL, we considered only the most stable sound propagation regime: deep-water propagation between a shallow source (whale) and shallow receiver (sonobuoy). Specifically, we included only data when the sonobuoy and call location were in waters deeper than 1000 m. At shallower seafloor depths, TL are expected to be more sensitive to the bathymetry and geoacoustic model of the seabed, so would be expected to generate less accurate TL and SL.

Transmission loss modelling was performed using the parabolic equation method via the software RAMGeo (Collins, 2002) as distributed in AcTUP (Duncan, 2005), a wrapper and interface to the Acoustics Toolbox (Porter, 2020). To estimate TL at different potential whale depths, we applied the principle of reciprocity (Jensen et al., 2011). Here our working definition of reciprocity is that the TL between source A and receiver B will be identical even if the locations of A and B are switched. For all propagation models, the source was modelled at the known depth of the sonobuoy, 140 m. The grid of receivers (whales) was modelled with 5 m depth resolution, and 10 m of range resolution. Antarctic blue whales are known to produce sounds at depths typically between 25 and 30 m when they are in the sub-Antarctic (Bouffaut et al., 2021), and other subspecies of blue whales are also known to produce sounds at similar depths, though singular (stand-alone) version of calls tends to be produced slightly shallower (Lewis et al., 2018). There is more variability in the estimates of the depths at which fin whales produce sound. These range from 10–15 m (Stimpert et al., 2015) to 50 m (Watkins et al., 1987). Here we considered a range of whale depths from 10 to 50 m to assess the sensitivity of our method to depth. We also consider in more detail TL models with whale depths of 20 m for stand-alone blue whale unit-A and D-calls; 25 m for blue whale unit-A in Z-calls, and 15 m for both fin whale call types as these depths represent our best guess of whale calling depth for each call type from these prior studies. For these models we present detailed histograms and residuals for the spatially mixed geoacoustic model (described below). TL for ABW unit-A, Z-calls, and fin whale 20 Hz pulses were modelled at a frequency of 25 Hz, while D-calls and 40 Hz pulses were modelled at a frequency of 50 Hz.

Transmission loss was modelled independently for each detection (i.e., measured RL). For each detection, seafloor depth was extracted from the etopo1 database (Amante and Eakins, 2009) along the bearing connecting the sonobuoy and call location. The distribution of depth ranges along each of these bathymetric profiles is summarized in Figure 2. The sound speed profile for that detection was then calculated from the ENRICH voyage CTD data that were nearest in time, usually within 24 h (Figure 2). Sound speeds were calculated from pressure (depth), salinity (conductivity), and temperature using the formula of Chen and Millero (1977). Since the CTD was deployed only to 1000 m, sound speeds extracted from the world ocean atlas (Boyer et al., 2018) were used to fill the remaining depths between 1000 m and the seafloor.

Three different models of the geoacoustic properties of the seafloor were used, following the seabed sediment properties as described by Lurton (2010). The first model was meant to approximate highly porous sediments with large grain size (e.g., ooze, clay, mud) and here the seabed was assumed homogeneous with a compressional sound speed of 1500 m/s and density of 1200 kg/m3. The second model was meant to approximate less porous sediments with small grain size such as sand, and the seafloor was assumed homogenous with a compressional sound speed of 1749 m/s and density of 1941 kg/m3. The acoustic parameters for these models were based on the only geoacoustic data available for the study area: densities and sound speeds measured from sediment cores taken in the survey area during a geoscience voyage in 2000 (Brancolini and Harris, 2000). A third model was also considered that was a spatial mix of the first two models with the sediment type determined from the Census of Seafloor Sediments (Dutkiewicz et al., 2015). Here, the small-grain model was used when the Census of Seafloor Sediments indicated the sediments were “sand” at the whale location, and the large grain model was used everywhere else (Figure 3). Use of the parabolic equation and geoacoustic model meant that we did not need to calculate eigenrays, or make assumptions about whether a detection was a direct arrival or from a multipath surface or bottom reflection. As a reference solely for visual comparison purposes, we also plotted alongside each of our geoacoustic TL models, transmission loss resulting from simple spherical spreading (TL = 20 log₁₀r) and cylindrical spreading beyond a “transition range” of our mean depth of 2250 m: TL = 10 log₁₀ 2250 + 10 log₁₀r (Urick, 1983; Hall, 2005; Lurton, 2010).

We illustrate our method for estimating SL with real examples of each of the key steps for a pair of detections for a single call Figure 4. Different coloured panels in this figure illustrate features of the detection matching and localisation algorithm, environmental data, and propagation modelling for a single D-call. Here the TDOA measured between detection #1 on sonobuoy A and detection #2 on sonobuoy B is 7.70 s, and the TDOA estimated from the localisation from sonobuoy bearings is 7.00 s. Our tolerance (0.3 * TDOAmax) is 2.18, so these detections are deemed to be a match (reception of the same call). The RL and TL for this example illustrate the complexity of shallow-water high-latitude sound propagation. RLs were 108.9 dB, and 109.3 dB at detection #1 and detection #2, respectively. This is an interesting quirk because detection #1, which was only 3 km from the whale, had RL that was 0.4 dB lower than that from detection #2 which was nearly 13 km away. However, the larger TL of 86.0 vs. 84.4, respectively, at 20 m depth potentially provides some explanation. The addition of RL and TL yielded SL of 195.0 and 193.7 for detection #1 and detection #2, respectively. Thus, our estimate of SL for the call is 194.3 and our estimate of measurement error, σ_∈, is 0.92 for this call. In this example, detection #2 appears to consist of obvious multipath arrivals. Looking closely at the spectrogram, one can see a faint copy of the call preceding the detection. This faint preceding call is most likely the direct arrival, but could also be a weaker multi-path. None of these seeming irregularities pose any problem for our methods, however, they are worth noting because they could cause problems for algorithms that require common assumptions about arrival order or TL (e.g., assumption that all detections are direct arrivals, or assumption that there is a monotonic decreasing relationship between TL and distance).

During 205 h of concurrent acoustic recordings from multiple sonobuoys (150 h of 2-audio-channel and 55 h of 3-audio-channel recordings), 1915 calls met all of our criteria for SL estimation, including 532 unit-A calls, 350 complete Z-calls, 850 D-calls, 62 20 Hz pulses, and 121 40 Hz downsweeps measured over 19, 14, 18, 7, and 12 recording-days, respectively (Table 1). The locations of these calls were predominantly south of 65°S, over the continental slope, and adjacent deeper waters (Figure 5). RLs of triangulated calls ranged from 93 to 125 dB re 1 μPa. Unlike the results of Širović et al. (2007), our RLs did not show a clear trend with log of distance as would be expected for a source at fixed depth with constant SL and with a TL that follows a simple geometric spreading law of the form X log r (Figure 5). Nor did our RL appear to closely follow a simple analytical model of cylindrical spreading in a surface duct with a transition range of r₀, such that TL = 10 log10 r₀ + 10 log10 r (Urick, 1983; Hall, 2005; Lurton, 2010).

Mean depth along the 4155 bathymetric profiles used in TL modelling was 2250 m, and sound-speeds ranged from 1441 to 1508 m/s across all CTD casts (Figure 2). The depth of the sound-speed minimum ranged from 40 to 160 m across all CTD casts, with a median of 110 m (Figure 2). According to the census of marine sediments (Dutkiewicz et al., 2015) approximately 25% of the triangulated detections were located over sand (Figure 3).

The different geoacoustic models yielded substantially different estimates of TL, with the clay model leading to highest TL, the sand model leading to lowest TL, and the mixed geoacoustic model yielding intermediate results, but closer to the clay model (Figure 6).

The geoacoustic model and modelled depth of whale had a strong influence on SL. The higher TL of the “clay” geoacoustic model yielded higher SL, and vice-versa for the lower TL of the “sand” geoacoustic model. Shallower whale-depths yielded larger SL across all geoacoustic models (Figure 7). However, the differences in SL among the three geoacoustic models were largest at shallow depths and decreased as whale depth increased.

The assumptions about whale calling depth had the largest effect on estimated TL and consequently SL. For all call types SL estimates were 15–20 dB higher if the whale was assumed to be calling at a depth of 10 m compared to 50 m (Figure 7). For our best guess of whale depth and geoacoustic model, SL across all call types ranged from 188 to 192 dB, with standard deviations of calls ranging between 5.5 and 7.8 dB and measurement error, σ_∈, ranging from 5.7 to 6.7 dB (Figure 8). SL for the subset of triangulations with intersection angles near right-angles were 2–3 dB lower than the full datasets for respective call types, except for unit-A within a Z-call, which was 1 dB higher (Figure 8).

For models where whale depths for each call type corresponded closely to those reported in previous studies (20 m depth for stand-alone unit-A and D-calls, 25 m for unit-A in Z-calls, 15 m for fin whale 20 and 40 Hz calls), mean SL measured in our study were all within 4 dB of each other. Our modelled SLs at our best estimate of whale depth for each call type were broadly similar to previous estimates in the Antarctic (Figure 8 and Table 1). Our estimates of SL of unit-A within a Z-call was 2 dB higher than previous estimates which both had a mean of 189 dB (Širović et al., 2007; Bouffaut et al., 2021), however, they are all within standard deviations of these other measurements. Our estimates of stand-alone unit-A of 188 dB is the first estimate of SL for this call. Compared to SL of unit-A produced within a Z-call it is 3 dB lower than our estimate from this study, and 1 dB below previous studies (Širović et al., 2007; Bouffaut et al., 2021). Our estimate of a mean SL of 192 dB for fin 20 Hz pulses was 3 dB higher than previous estimates of 189 dB made in the Antarctic and elsewhere (Širović et al., 2007; Weirathmueller et al., 2013), but again within standard deviations. Our estimate of SL of 190 for D-calls was considerably higher than calls that were measured in the North Atlantic (Berchok et al., 2006; Akamatsu et al., 2014). Thode et al. (2000) presented source spectrograms of 3 D-calls with spectral densities between 160 and 170 dB re 1 μPa Hz^–1 from 40 to 80 Hz. Converting these to comparable units by adding 10log₁₀(B), where bandwidth, B = 40 Hz, we obtain an SL of 176–186 dB. Our estimate of 190 is 4 dB above the high end of that range. Our SL of 189 dB for fin whale 40 Hz downsweeps are the first published estimate of SL for this call type in the Southern Hemisphere. These appear to be similar to the higher of the two measurements reported in the Northern Hemisphere (Wiggins and Hildebrand, 2020).

Though we have focused on the calls of Antarctic blue and fin whales on their Antarctic feeding grounds, our method for estimating SL from sonobuoy detections should be generally applicable to other call types and environments – provided that there is sufficient environmental information to model TL with fidelity. Due to only having coarse information of seabed geoacoustics from a survey two decades prior, we opted not to include shallow water localisation in our results. Our inclusion of only deep-water detections could potentially result in censoring bias if whales change SL in response to water depth. However, Miksis-Olds et al. (2019) provide an illustrative example of the difficulties of SL estimation in shallow water without sufficient knowledge of seabed geoacoustics, so SL of blue and fin calls in shallow water remains a knowledge gap for future studies to address.

The only component of our method that is not general-purpose is our localisation and detection matching algorithm, which makes use of direction of arrival via the sonobuoys acoustic vector sensors. However, use of bearings in our localisation algorithm only served to simplify the logistics of deploying our arrays by allowing for localisations from as few as two sensors and by simplifying time-synchronisation since signals from all sonobuoys were recorded simultaneously on the same analog-to-digital converter using the same clock. In comparison TDOA methods typically require 3, if not more, omnidirectional sensors for effective localisation (Møhl et al., 2001; Spiesberger, 2001; Thode et al., 2006; Baumgartner et al., 2008; Miller and Dawson, 2009), and also often require additional hardware or software-processing for accurate time synchronisation (Møhl et al., 2001; Thode et al., 2006). Assuming that there is sufficient hardware to obtain a localisation (via bearings, TDOA, or some combination of the two), then our remaining methods, e.g., for measuring RL and estimating TL, can all be applied more broadly to any relatively narrow-band low-frequency sounds in almost any marine environment without modification. For calls that have larger bandwidth than those of blue and fin whales, broadband propagation modelling may need to be considered (see Helble et al., 2013 for an example of how such methods might be developed and applied). Furthermore, the general-purpose methods we presented are not mutually exclusive with other localisation algorithms. Thus, future studies could investigate the use of multipath, either independently (Bouffaut et al., 2021) or alongside bearings and TDOAs and propagation modelling (Nosal and Frazer, 2007) in order to better estimate range, depth, bottom composition, and measurement error.

Regardless of geoacoustic model, and regardless of whale depth, the variability in our SL was larger than previous estimates of ABW unit-A, similar to variability reported for blue whale D-calls, in the middle of the pack for fin whale 20 Hz pulses, and higher than the two estimates of SL for 40 Hz pulses (Table 1). This suggests that the precision of our SL estimates is either lower than some of the other methods, or there is more natural variation in SL (or whale calling depth) than has been captured in previous studies of Southern Hemisphere whales. The Gaussian distribution of SL for the signals measured in this study suggests that there could be a natural range of SL produced by blue and fin whales and further that this range may naturally vary with animal depth. This makes sense from a physiological perspective since every 10 m of depth results in an increase of 1 atm (100 kPa) of pressure, potentially constraining the SL of an air-based sound source as has been proposed for blue whales (Aroyan et al., 2000; Dziak et al., 2017).

The variability in our SL estimates comprises natural variability in the source level as produced by the animals, as well as errors in our measurements and models. Quantifying the amount of natural variability both within individuals and among different whales would be biologically interesting. However, within our results, we are unable to separate natural variability from errors in our measurements and models because variability from measurement error, σ_∈ was similar in magnitude to the total variability across all calls σ. The main sources of measurement and modelling error in our study include: errors in measuring received level, errors in localisation, and mismatch of TL models with reality.

While there are clearly limitations with the precision of SL estimates made using sonobuoys, a key advantage is that sonobuoys deployed as part of a dedicated cetacean research program provide more opportunities to estimate SL than tags on individual animals (small sample size), dipping hydrophone recordings (which rely on assumptions about the underwater behaviour of animals seen at the surface) or even single fixed hydrophones (where precise locations of calling whales are generally unavailable). Our SL for Antarctic blue whale sounds were estimated from a much larger number of calls than previous studies, and a similar number of fin whale calls as previous Southern Hemisphere studies, even though fin whales were not the primary research focus during the voyage. In addition to a higher number of calls, our estimates came from a larger number of recording sessions/different days, a typical proxy for number of individual animals.