The Lombard effect in singing humpback whales: Source levels increase as ambient ocean noise levels increase

Humpback whale (Megaptera novaeangliae) males sing an elaborate song throughout their range, but the song is most prevalent in their low-latitude wintering areas. The humpback whale song is composed of themes that can be broken into phrases which are made up of units (Payne and McVay, 1971). Each subsection within a song may be repeated multiple times, but the order remains the same (Payne and Payne, 1985; Payne and McVay, 1971). Most singing humpbacks are physically separated from other whales, and they tend to stop singing when they join another whale that is not singing (Tyack, 1981). Within a given year, every singing humpback in a population sings the same song, but this song changes throughout the breeding season and from year to year (Payne and Payne, 1985; Tyack, 1981; Winn et al., 1981).

Humpback whales produce song units with mean source levels (SLs) that vary depending on the unit type, the individual, and the occurrence (Au et al., 2006). Au et al. (2006) measured the root mean square (RMS) SLs for song units recorded from three humpback whales in the Auau Channel of the Hawaiian Islands and reported average RMS SLs ranging from 149 to 169 dB $re 1 μPa$ at 1 m, depending on the unit type and whale. RMS song SLs measured in Stellwagen Bank National Marine Sanctuary in Massachusetts Bay were an average of 170 dB $re 1 μPa$ at 1 m with a standard deviation of 3 dB (Cholewiak et al., 2018). Besides song, humpback whales of both sexes produce social sounds. These vocalizations may be used to coordinate whales joining together, signal competition, and provide information about the caller's identity (Dunlop et al., 2008). Unlike song, which is thought to be for long-distance communication, social sounds may be for communication between closer whales (Dunlop et al., 2008). Non-song social calls are reported to be produced at SLs ranging from 131 to 190 dB $re 1 μPa$ at 1 m, depending on the noise level (NL), call type, and context (Dunlop, 2016b,a; Dunlop et al., 2013; Dunlop et al., 2014; Fournet et al., 2018b; Fournet et al., 2018a; Thompson et al., 1986).

Humpback whales have been observed to change their social call SL in response to noise, an example of the Lombard effect (Dunlop, 2016b; Dunlop et al., 2014; Fournet et al., 2018a). Migrating humpback whales off Australia increased the RMS SL of their social calls 0.9–1.5 dB/1 dB increase in the background NL due to natural wind sources (Dunlop, 2016b; Dunlop et al., 2014). The humpbacks seemed to call at SLs about 60 dB above the 81–108 dB NL (Dunlop et al., 2014). Humpback whales in their summer foraging areas off Alaska also increased social call SLs in increased background noise (Fournet et al., 2018a). This background noise was due to both local vessels and natural sources, and SLs of humpback calls increased by 0.8 dB/1 dB increase in the background NL (Fournet et al., 2018a).

In addition to increasing their SLs, humpback whales respond to increased background NLs in other ways (Dunlop, 2016b; Dunlop et al., 2010; Risch et al., 2012). Migrating humpback whales were observed increasing their surface activity, including breaches and pectoral fin slaps, during periods of high natural background noise, but did not significantly alter their behavior during vessel noise (Dunlop, 2016b; Dunlop et al., 2010). Humpback whales singing in Stellwagen Bank National Marine Sanctuary decreased their time spent singing during times of anthropogenic frequency-modulated (FM) pulses between 400 and 1000 Hz, which were produced to detect groups of fish (Risch et al., 2012).

Although most investigations into the Lombard effect in cetaceans have reported nearly a 1 dB increase in SL for every dB increase in the background NL (e.g., Dunlop, 2016b; Dunlop et al., 2014; Fournet et al., 2018a; Holt et al., 2009; Parks et al., 2010), recent studies have not shown this same compensation ability, and most studies with other taxa have not reported full SL compensation for increases in background noise. Minke whales (Balaenoptera acutorostrata) increased the SL of their boing call an average of 0.24 dB/1 dB increase in the background NL (Helble et al., 2020). Bottlenose dolphins (Tursiops truncatus) increased the apparent output level of their whistles 0.1–0.3 dB/1 dB increase in background NLs (Kragh et al., 2019). Some frogs have shown no Lombard effect (e.g., Love and Bee, 2010), but male tùngara frogs (Physalaemus pustulosus) increased the amplitude of their in-air calls 0.1–0.3 dB/1 dB increase in natural background NLs in a laboratory (Halfwerk et al., 2016). Terrestrial animals, such as birds and primates, have exhibited a similar Lombard effect (reviewed in Kragh et al., 2019). Great tits (Parus major), for example, increased their song amplitude approximately 0.3 dB/1 dB increase in white NL (Zollinger et al., 2017). The first study on nonhuman primates showed that two macaques (Macaca fascicularis and Macaca nemestrina) increased the amplitude of their calls 0.2 dB/1 dB NL increase, which overlapped with their call bandwidth (Sinnott et al., 1975). When animals do not or cannot increase their SLs the same amount as increased NLs, their communication range will decrease. The impacts of a smaller communication range depend on the duration of the noise and the purpose of the call. These impacts are often difficult to assess when the purpose of the vocalizations is unclear.

Approximately 10 103 humpback whales are part of the central North Pacific stock and spend their winters around the Hawaiian Islands (Muto et al., 2019). These whales migrate north in the spring and spend their summers feeding primarily off northern British Columbia, southeast Alaska, and in the Gulf of Alaska (Barlow et al., 2011; Calambokidis et al., 2001). The humpback whale song has been recorded, and singing whales have been tracked in the area of the U.S. Navy's Pacific Missile Range Facility (PMRF) off Kauai, HI (Helble et al., 2015; Henderson et al., 2018). In this area, singing humpback whales swim with a mean speed of 3.5 km/h and show several different behavioral states, including directed travel, repeated stationary dives, milling, or a combination of the three (Henderson et al., 2018). Since it is not known how singing humpback whales in wintering areas respond to changes in background noise, these whales are the focus of this investigation.

The objective of this study was to calculate the SLs of humpback whale song units recorded on the PMRF off the island of Kauai and compare these SLs over a range of background NLs. Any change in SL as a function of the background NL was compared with previous publications on the observed Lombard effect in humpback whale non-song calls and other animal species. In this analysis, methods follow similar procedures to those used by Helble et al. (2020) for minke whale calls in the same area over the same time, so results are directly comparable between these studies. Knowledge about how humpback whale acoustic behavior changes with natural background noise fluctuations is necessary to put behavioral changes during anthropogenic disturbances into context.

The U.S. Navy's PMRF is located off the northwest coast of the island of Kauai in the Hawaiian Islands. Since 2011, an array of time-synchronized hydrophones from the underwater range has recorded at least two days per month, in addition to recording during U.S. Navy mid-frequency sonar training events. Although the number of hydrophones in the array and the sampling rate has changed over the years, from August 2012 to July 2017, the array configuration used in this study remained the same, containing 14 broadband hydrophones with a 96 kHz sampling rate. Starting in 2014, additional opportunistic recordings, spanning several weeks, were made at a 6 kHz sampling rate. The 14 offshore hydrophones used to localize humpback whales were at depths of 3150–4700 m and covered a rectangular-shaped grid approximately 20 km to the east/west and 60 km to the north/south (Fig. 1). These hydrophones were divided into four subarrays, containing a center hydrophone and four corner hydrophones. All data recorded at 96 kHz were down-sampled to 6 kHz before processing for sampling rate consistency. The system was designed to have a specific free field voltage sensitivity, and all hydrophones were validated to meet this specification within ±3 dB. Recordings from throughout the year were analyzed, although humpback vocalizations were only recorded in fall, winter, and spring. The recording effort (in hours) for each month was calculated, along with the number of acoustic localizations of humpback whale song units within 2–10 km of the center hydrophones of each array for each month (Fig. 2).

Several assumptions were made throughout this study. These assumptions are introduced here and their validity is discussed further in Sec. IV. The humpback whale vocalization directivity was assumed to be zero and both the animal source and receiving hydrophones were treated as omnidirectional. The main components of most humpback song units ranged from 150 to 1000 Hz, so the SL and NL measurements that were presented were limited to that band. The NLs were recorded on the bottom hydrophones and assumed to be a proxy for the noise experienced by the whale. The geometrical spreading transmission loss (TL) model used in this study was tested across different ranges and compared with other sound propagation models. Even though this model seems to perform well for the vocalizations used in this study, TL is more complex than what this model suggests and may be affected by other properties of the ocean like sea-surface roughness and internal waves.

Passive acoustic whale locations were estimated by detection and feature extraction, cross-correlation of those features to obtain time difference of arrivals (TDOAs) of the signal at each hydrophone, and comparison of these measured TDOAs with theoretical TDOAs across the search area. These steps are outlined in detail in other publications, using vocalizations from humpback whales (Helble et al., 2015), Bryde's whales (Balaenoptera edeni; Helble et al., 2016), and minke whales (Helble et al., 2020) and are, therefore, only summarized in this paper.

The generalized power-law (GPL) detector (Helble et al., 2012) was used to detect humpback whale song units (units). The GPL detector determined the start and end time of each unit and used a spectral “templating” procedure that subtracted the underlying noise in each frequency band from the detection, leaving only the spectral contents of the signal. These templates were later cross-correlated across hydrophones to obtain TDOAs.

To localize the singing humpback whales, the 14 hydrophones were divided into 4 subarrays [(A),(B),(C),(D)] of five hydrophones each (Fig. 1). If the unit was detected on the center hydrophone and at least three of the four corner hydrophones in a subarray, it was localized. The signal templates were cross-correlated to calculate the TDOA of the unit between the center hydrophone and the corner hydrophones. The subarray configuration was chosen so that a direct path solution always existed on the four hydrophone pairs across the monitored area, and so the maximum allowable TDOA between the center hydrophone and the corner hydrophones in the subarray was limited to the direct-path propagation time between them.

To calculate the location of the singing humpback whales, a minor modification was made to the methods outlined in Sec. II of Helble et al. (2015). Rather than using sequences of units, single humpback whale song unit templates were cross-correlated to estimate the TDOA of the unit between pairs of hydrophones. This modification was also used for minke whales (Helble et al., 2020) and Bryde's whales (Helble et al., 2016) and allowed for a precise location to be assigned to each unit produced. For single tonal humpback song units, timing delay errors were on the order of 40 ms, resulting in localization standard deviations of less than 60 m (Helble et al., 2015). As described in Sec. II D, localization accuracy was necessary for modeling the TL between the whale location and the recording hydrophone, which was used to calculate the SL.

Localized humpback whale song units were grouped into individual tracks using a semi-automatic tracker previously described by Klay et al. (2015). Localized units out to 20 km from the center hydrophone were considered for tracks in order to reduce the chance that a single whale's path was separated into multiple tracks. The units were recursively examined so that the elapsed time and distance between units aligned with reasonable assumptions about humpback whale swimming and singing behavior. A humpback whale track required a minimum of 12 localized units. Successive localizations were allowed to be separated by a maximum of approximately 1 km and 15 min. The number of tracks was a rough estimate of how many individual whale encounters were included in this study. These tracks allowed for efficient verification of localized detections as well as analysis of how presumed individuals responded to NL changes.

All tracks were validated to consist of humpback whale song units by an analyst. The analyst viewed a map of the track, the inter-unit intervals, and the corresponding spectrograms for a subset of units along each track to determine if that track was produced by a humpback whale. Tracks were consistent for singing humpback whales, and so detections here are assumed to be song units. However, it was not feasible to manually verify all individual units, and so it is possible that a small portion of these detections are social calls.

The spectral density was calculated and integrated over the frequency bandwidth of interest to estimate the sound pressure spectral level of the received humpback whale song units and the background noise. The spectral level, or mean square received level (RL), measured as $μPa2$⁠, is

where f_s is the sampling frequency, nFFT is the number of samples used in each fast Fourier transform (FFT) window, and $Sp(fi)$ is the spectral density and is summed over n frequency bins. The spectral density has units of $μPa2/Hz$ and is calculated by incoherently averaging nT time segments of the squared-magnitude of the fast Fourier transformed signal (⁠$|Xj(fi)|2$⁠) as in

The factor of 2 starting the right-side of the equation accounts for energy at negative frequencies. In the denominator, the ratio of the sampling frequency and the FFT length (⁠$fs/nFFT$⁠) normalizes by the bin width. The sum of $wi2$ is the sum of the square of the window function that is multiplied by each of the j time series segments before Fourier transforming. An nFFT of 1280, an overlap of 75%, a sampling frequency f_s of 6 kHz, and a Hamming window were used in this analysis.

To estimate the NL at the time of a song unit, a time series sample was selected from just before and after that unit. A 1 s buffer was used between the signal and the noise samples so that any residual signal not included in the detection was not included in the noise sample. $Sp(fi)$ [Eq. (1)] was summed from f₁ = 150 Hz to f_n = 1000 Hz. This band matches the frequencies used for the unit templates and covers the dominant frequencies of humpback song units. Approximately 2.5 s of noise was used in the noise sample, which corresponds to nT = 60 time segments. The noise samples taken before and after the unit were similar, indicating that there was no signal present in the noise measurements.

The noise in this 150–1000 Hz band was primarily due to ambient noise sources in contrast to local point sources. In deep water, the primary contributors to noise in this band are wind, waves, and shipping traffic (Wenz, 1962). No Navy exercises were taking place during the recordings used for this study and because the study area is a restricted area, any vessel noise would have been from distant shipping traffic and not local shipping activity. Although no times with detections were used for noise calculations, overall ambient noise may also include distant humpback whale singing activity.

To estimate the RL of humpback whale song units, $Sp(fi)$ was again summed from f₁ = 150 Hz to f_n = 1000 Hz over the duration of the signal as determined by the GPL detector. A song unit detection contains both the signal from the unit and the signal from the background noise, so the unit templates as described in Sec. II B of Helble et al. (2015) were used in place of $Xj(fi)$⁠. The unit templates isolate the spectral contributions of the song units from the background noise.

The NL and song unit RL measurements were converted into decibels, using $10 log10(RL)≡RLdB,$ where RL is in units of $μPa2$ and RL_dB is in units of dB $re 1 μPa$⁠. This method calculates the RMS RL (RL_dB) and NL (NL_dB), which is the method used for the remainder of this paper.

Monte Carlo simulations were used to estimate the accuracy of RL measurements for humpback whale song units in all likely NLs. Five song units were chosen for these simulations to cover the diverse spectral and temporal characteristics of song units (Fig. 3). The signals were then reduced in amplitude and added to 160 min of randomly selected ocean noise recorded at the PMRF over all likely signal-to-noise ratios (SNRs), defined as

where RL_dB is the RL and NL_dB is the NL, both in units of dB $re 1 μPa$⁠. The GPL detector was used to detect these song units inserted into the noise, and the RLs were estimated. The RL measurements were compared against the known RLs for each of the SNRs tested.

To estimate the SL of humpback whale song units, the RL of the units must be added to the TL between the source and receiver as described by

where TL_dB is in units of dB (RL at range r relative to RL at 1 m from the source), RL_dB is in units of dB $re 1 μPa$⁠, and SL_dB is in units of dB $re 1 μPa$ at 1 m. This equation assumes an omnidirectional source and receiver.

TL was estimated using two methods. The first method was to use the range-dependent acoustic model (RAM) to estimate the TL between the location where the song unit was produced and the hydrophone location where the unit was recorded (Collins, 1995). Peregrine is a C-language interface to the split-step Pade parabolic equation acoustic propagation code Seahawk (Heaney and Campbell, 2016; Heaney et al., 2017), which is based on the RAM and was used for this TL estimation in the same way as it was used by Helble et al. (2020). The TL over the humpback singing bandwidth was calculated by incoherently averaging the TL between the source and receiver in the 150–1000 Hz band in 5 Hz increments. TL values were interpolated in 60 deg radial increments from each hydrophone. Since the TL varied less than 1 dB as a function of the azimuth over the range used for this study, this azimuthal interpolation was justified. TL was calculated for whale depths between 5 and 100 m, covering the depths at which most of the units were expected to be produced (Henderson et al., 2018). Bathymetry data were retrieved from the National Oceanic and Atmospheric Administration (NOAA) National Geophysical Data Center U.S. Coastal Relief Model with 3 arcsecond resolution (National Geophysical Data Center, 2005). Past seasonal sound speed profiles (SSPs) were calculated from the 2018 World Ocean Atlas (Locarnini et al., 2018; Zweng et al., 2018). The sediment was treated as an acoustically thick halfspace (implemented as 20 wavelengths at the given frequency, containing an exponential absorptive sponge along the bottom of the sediment layer). Various grain sizes on the Krumbein phi (⁠$ϕ$⁠) scale (Krumbein and Sloss, 1951; Wentworth, 1922) were tested in the TL model. The TL was calculated over all likely bottom compositions and SSPs, and combinations of SSPs and bottom types that resulted in the highest and lowest TL values were noted.

The second method for estimating TL was to use the geometrical spreading and attenuation loss equations described by Urick (1983). For slant ranges from the source to the hydrophone greater than the seafloor depth at the source location, the SL was estimated from the RL and the location of the source by adding the losses from both spherical and cylindrical spreading,

SL_dB is the SL (dB $re 1 μPa$ at 1 m), RL_dB is the RL (dB $re 1 μPa$⁠), r_T is the transition range in meters at which geometrical spreading transitions from spherical to cylindrical, α is the attenuation loss coefficient in dB/km, and r is the slant range from the whale to the hydrophone in meters (Urick, 1983). For slant ranges from the source to the hydrophone less than r_T, the SL was calculated using spherical spreading only,

The transition range, r_T, was estimated to be one water depth since the whales vocalize near the surface and the hydrophones are raised just above the seafloor (Henderson et al., 2018). To confirm, a variety of values were tested for r_T, and the binned average SLs were plotted as a function of range. Additionally, TL as a function of range with various r_T values was compared against the Peregrine model. The attenuation loss coefficient α was primarily influenced by frequency dependent absorption by the water for the relatively short ranges and deep water used in this study. Using the method described by Ainslie and McColm (1998) and a median song unit frequency of 600 Hz, the attenuation coefficient α was calculated to be $α=0.03 dB/km$⁠. If instead values at the lowest (150 Hz) and highest (1000 Hz) frequencies of the detected song units were used to calculate α, the estimated TL shifted by less than 0.3 dB. This small difference across the full frequency range justified the choice of selecting the median value of 600 Hz. Attenuation is minimal compared to geometrical spreading, so “geometrical spreading and attenuation loss equation” is shortened to “geometrical spreading equation” for the remainder of this paper.

The TL estimates from both methods (the Peregrine model and the geometrical spreading equation) were added to the measured RLs of humpback whale song units to calculate SLs. The average SL should be independent of range and, therefore, constant over all detectable distances. However, at farther ranges, lower SL units may be masked from detection, limiting the detections to units with greater SLs. This masking would cause average SLs to trend higher at farther distances (discussed in more detail in Sec. II E).

To determine the best TL model for this study, the humpback whale song unit RLs from ranges of 0–20 km were measured. These RL values were averaged over 10 m horizontal range increments and plotted as a function of range. Using both the Peregrine model and the geometrical spreading equation, TL was calculated for each of the units. As with RL, TL was averaged in 10 m horizontal range increments and plotted as a function of range. Peregrine was unable to estimate the TL from song units at horizontal ranges between 0 and 2.5 km due to the high-angle propagation and the inherent limitations of the parabolic equation from a source near the surface and a receiver near the bottom of the deep ocean. The TLs and resulting SLs estimated from the two methods were compared over the ranges available.

Passive acoustic detection and localization of marine mammals is affected by the acoustic environment. Both the TL as sound travels through the environment and the background NL are factors in the received SNR and affect the probabilities of detection and localization [see Eqs. (3) and (4)]. Masking occurs when a signal is unable to be detected because of the addition of other sound from the environment. Masking is a primary concern when measuring changes in vocal behavior to ensure that changes are real and not artifacts of the inability to detect the signals of interest.

The probabilities of detection and localization of humpback whale song units were modeled by simulating animal source locations randomly distributed within a 20 km radius from the center hydrophone of each subarray and vocalizing over the range of estimated SLs (see Sec. II F) in the observed background NLs. The estimated probability of detection at each hydrophone, $P̂D$⁠, was calculated by

where $ρ(r,θ)$ is the probability density function (PDF) of whale singing locations in the horizontal plane, and $g(r,θ)$ is the detection function (Buckland et al., 2001). A homogeneous random distribution of animals over the area of interest, $π(w22−w12)$⁠, was assumed and, therefore, $∫∫ρ(r,θ)=1/(π(w22−w12))$⁠. Since the humpback whales were assumed to sing near the surface, the detection function, $g(r,θ)$⁠, was assumed to be a function of range and azimuth only. The detection range of w₁ to w₂ is measured from the recording hydrophone to the source. Often, w₁ is set to zero with the animal directly above the hydrophone but is included here as a variable for reasons explained subsequently. The detection function also depends on the azimuth due to differences in bathymetry across the search area. The probability of localization depends on the probability of detection at the center hydrophone in the subarray and at least three of the corner hydrophones.

Monte Carlo simulations were used to calculate the probability of detection (⁠$P̂D$⁠) at each hydrophone. The humpback whale song unit detection function was dependent on the TL, SL, and NL such that $g(r,θ;TL,SL,NL)$⁠. To characterize the GPL detector performance, five humpback whale song units (Fig. 3) were randomly distributed throughout the search area in all background noise conditions, and their probabilities of detections were calculated. First, the noise was removed from these five high SNR units following the methods described in Helble et al. (2012). Next, the amplitudes of the song units were adjusted so that the SLs ranged from 140 to 190 dB $re 1 μPa$ at 1 m in 0.5 dB increments. To simulate RL, these signals were reduced in amplitude based on the estimated TL from each simulated whale position. The simulated received signal was randomly added into noise taken from 160 min of the PMRF noise samples and processed with the GPL detector. This process was repeated for all five units across the search area, and $P̂D$ was calculated for song units in all likely combinations of SL and NL at each hydrophone. In these simulations, only the amplitude of the signals was reduced by the TL and no distortions (such as multipath) that may have affected the detectability of units were simulated. However, since w₂ was limited to ranges of primarily direct-path propagation, signals recorded on the PMRF were minimally distorted by the environment.

To calculate the probability of localization (⁠$P̂L$⁠) or the probability that a song unit was detected on the center hydrophone of an array and at least three corner hydrophones, $g(r,θ)$ at the center hydrophone was multiplied by the highest three of four $g(r,θ)$ probabilities from the corner hydrophones, where r and θ from each adjacent hydrophone were adjusted to reference the same source position as defined by the center hydrophone. The resulting probability of localization function, $gL(r,θ)$⁠, at the center hydrophone was inherently less than the detection function, $g(r,θ)$⁠, due to the requirement of the song unit being detected on the center hydrophone and at least three of the supporting hydrophones. As with the probability of detection calculation [Eq. (7)], $gL(r,θ)$ was multiplied by the PDF and then integrated over the range and azimuth to get $P̂L$⁠. The probability of localization for each of the four subarrays was similar since the depths and bathymetries did not vary much across the range.

To keep the probabilities of detection and localization across the study area close to one over the observed NLs, the maximum allowable radius (w₂) from the center hydrophone was set to 10 km. To reduce errors in estimated SL related to the uncertainty of position or directionality of the singing whale, the minimum radius (w₁) was set to 2 km. Any errors in depth, location, or unknown directionality of the song units would have a greater impact on the TL at these closer ranges. All SLs included in this analysis were from within this 2–10 km range.

SLs were calculated by adding the measured RL of each humpback whale song unit to the expected TL from the animal's position [Eq. (4)]. The NL in the 150–1000 Hz band before and after each unit was also measured and saved.

The relationship between the humpback whale SL and the ocean NL was modeled with a generalized additive model (GAM) and the “mgcv” package in R (Wood, 2017). An identity link function of the SL response variable was used and the error terms were assumed to have a Gaussian distribution. The NL predictor variable was modeled with a cubic regression spline smoothing term with five knots (k = 5) to capture the nonlinearities in the relationship between the predictor and response variable but not over-fit the data. To ensure that the number of knots were not over-specified, the effective degrees of freedom were used as a guide (Wood, 2017).

If song units were masked, the SL results could be overestimated. If a greater proportion of units were missed in higher NLs than lower NLs, then the estimated average SL would be biased high in increased NLs, artificially inflating any Lombard effect. The impact of masking was minimized in this study by limiting the detection range to 10 km from the center hydrophone based on the results from the probabilities of detection and localization calculations. Because it was impossible to ensure all units were detected, the sensitivity of the relationship between SL and NL was investigated by increasing the weighting of the left tails of the SL distributions. The number of song units in the masked region was simulated as the function $f(x)=axb$⁠, where f(x) is the number of units simulated at each SL interval (x), b is a constant controlling the rate of decay in the tail, and a was chosen so that the tail distribution generated the detected number of units just above the masked region and reached a value of zero at SLs of 145 dB RMS $re 1 μPa$ at 1 m, the lowest assumed SL based on these data. Values of $b=(1,2,3,4)$ were all tested where b = 1 results in the most extreme left-tail distribution (straight line). GAMs were fitted to these altered distributions in the same way as before.

Humpback whale song unit SLs were also analyzed in 5 dB NL bins. The average SL and variance in each bin were calculated and compared using the nonparametric one-sided Wilcoxon's rank sum test and one-sided Ansari-Bradley test. Histograms of SLs were plotted, and the shape and character of the distributions as a function of noise were visually examined. The histograms were fit to the data using nonparametric kernel smoothing distributions evaluated at 100 evenly spaced points over the range of SL values for each NL bin.

To investigate the response of individual whales to changing background noise, the relationship between the estimated SL and NL were examined on a per-track basis. Tracks that contained at least 200 song units and spanned a range of ocean NLs of 10 dB or more in the 150–1000 Hz band were considered for this analysis. These more strict conditions eliminated short tracks with low sample sizes and tracks that were in similar noise conditions for the full duration. The estimated SLs for the units in each track were plotted against the corresponding measured NLs. The slope of the linear fit of each track indicates the average individual response to the noise. The distribution of these slopes was analyzed. The GAM fits used for the aggregated data described previously were not applied on the individual tracks due to lower sample sizes and the high degree of variability of the song unit SLs.

Opportunistic recordings totaling 604 days from 2012 to 2017 were used to detect and track 83 974 humpback whale song units. These units were all produced at ranges of 2–10 km from the center hydrophone of each subarray and formed 524 tracks through the PMRF. Most tracked units occurred between December and April each year even though there was recording effort throughout the year (Fig. 2). Song unit RLs (Sec. III A) were added to their estimated TLs (Sec. III B) to calculate the SLs of the units (Sec. III D). TL was verified as accurate using a search area of 0–20 km and 289 467 song unit RLs.

To estimate RL measurement accuracy, example humpback whale song units were inserted into background noise to compare known RLs with RLs estimated using the GPL detector. Five example humpback song units were tested for SNRs from −15 to 30 dB (Fig. 4). The mean RL measurement error was −0.1 dB $re 1 μPa$ and the maximum error was 3.2 dB. A lower limit of −15 dB SNR was chosen to ensure accurate RL estimates while reducing the number of missed units and, therefore, also reducing the effects of masking. A negative SNR is possible because NLs were integrated across the full 150–1000 Hz band, but each unit did not include components from this full band [Eqs. (1) and (3)]. Probabilities of detection and localization values accounted for this SNR limit in their calculations.

RLs and associated NLs were measured for all detected song units. For the 83 974 humpback whale song units produced 2–10 km from the measuring hydrophone, the 25th, 50th, and 75th percentiles of the RL measurements were 93, 97, and 101 dB $re 1 μPa$⁠, respectively. The NLs associated with these units and averaged across hydrophones had 25th, 50th, and 75th percentiles of 90, 93, and 95 dB $re 1 μPa$⁠, respectively. Again, RL and NL values were similar because the spectral densities were integrated over a large band and, therefore, the frequencies of peaks in intensity were not necessarily the same.

Two methods, the Peregrine model and the geometrical spreading equation, were compared for estimating the TL. The two TLs were calculated for 289 467 RLs at distances of 0–20 km from the measuring hydrophone. The TL was averaged in 10 m range bins and plotted as a function of range [see the lower part of Fig. 5 for Peregrine (purple) and the geometrical spreading equation (black)]. The Peregrine model required environmental inputs to estimate the TL. A sediment grain size of $ϕ=5$ was assumed, and SSPs were calculated from temperature and salinity data that most closely matched the date of the song unit. Changing the assumed grain size and SSP to other likely values had negligible effects on the TL within the 10 km range that was used when calculating the SL. At farther ranges, the TL became more sensitive to these environmental values, and at 20 km, the TL varied by up to 7 dB using $ϕ$ = 4–8 and all likely SSPs. The sediment grain size value $ϕ=5$ was chosen based on past TL experiments that used sonobuoys to measure the TL in the region. The details of these experiments are not given here because sediment grain size did not noticeably affect the TL in the 2–10 km range of interest for this study. In addition, the Peregrine model used an assumed humpback whale singing depth of 30 m. Changing the assumed animal singing depth to the most likely values of 5–100 m (Henderson et al., 2018) also had negligible effects on the TL. The variability of the Peregrine TL for any given range was due to differences in bathymetry along the song unit transmission path that were dependent on the azimuth. The only adjustable parameter in the geometrical spreading model was the transition depth, r_T. This transition depth was assumed to be the water depth in this study. Choosing other values of r_T resulted in less agreement between the Peregrine model and the geometrical spreading model and caused SL to vary as a function of range (more about this in the next paragraph). Navy surface assets in the area with known SLs and similar frequency ranges were also used to verify r_T, and these estimated SLs closely matched the known SLs. In the 2–10 km range of interest, the greatest differences between the geometrical spreading model and Peregrine occurred at the closest ranges (3.4 dB at a range of 2 km). Complex surface-bottom interactions occur at close ranges that are not accounted for in the geometrical spreading model, and these interactions may explain some of the model discrepancies. Although close agreement between these two models does not make them good approximations of TL, their agreement is reassuring since it is less likely for both models to be just as erroneous.

The TL estimates from the two models were added to the song unit RLs that had been averaged in 10 m range bins (Sec. III A) to estimate the SL as a function of range (upper part of Fig. 5). The average SL values fluctuated by less than 5.0 dB from 2 to 10 km with a slightly negative slope (0.11 dB/km using a linear fit) when the geometrical spreading equation was used to calculate the TL. The average SL values varied by 5.9 dB or less with a slightly positive slope (0.12 dB/km using a linear fit) when the Peregrine model was used to calculate the TL. To verify if masking could be biasing the SL estimates upward at farther ranges, this process was repeated using only song units that occurred in noise backgrounds of 90 dB $re 1 μPa$ or less (minimal to no masking expected), which resulted in no appreciable change in Fig. 5. The SL was only estimated for units produced between 2 and 10 km from the center hydrophone (indicated by the red dotted vertical line in Fig. 5). The localization and depth uncertainty at close ranges resulted in proportionally more uncertainty in SL than those at farther ranges. The range was limited to 10 km to minimize the effects of masking (discussed further in Sec. III C) and because of greater TL uncertainty at farther ranges.

Because both models gave similar results, either model would be suitable for estimating the TL in the study area. For this study, the geometrical spreading equation was used. The SL did not change appreciably as a function of range when this model was used to calculate the TL, and the computation time was significantly faster for this model when compared with the Peregrine model.

To assess the impact of masking on the SL results, the probabilities of detection and localization were calculated for all observed humpback whale SLs across all background NLs. The estimated probability of localization, $P̂L$⁠, was determined for all likely combinations of SL and NL over a range of w₁ = 2 km and w₂ = 10 km from the center hydrophone of each subarray. Assuming a random spatial distribution of animals on the range, the average of $P̂L$'s over the four subarrays gave the average probability of localization as a function of the SL and NL. The probabilities of detection and localization were close to 100% over the 2–10 km search area from each center hydrophone for all observed NLs. Example probabilities of detection and localization versus range are shown for subarray (D) (Fig. 6). The assumed song unit SLs for this example were distributed with a mean of 168 dB $re 1 μPa$ at 1 m and a variance of 30 dB, which was the observed SL distribution in the lowest NL bin of 80–85 dB. The NLs used for these simulations were 85 and 100 dB $re 1 μPa$⁠. These maps, therefore, represent the worst-case scenario for masking because they simulate the case where the SL distribution does not change as the NL increases (no Lombard effect). The probabilities of localization in these examples are $P̂L$=99.9% for NL = 85 dB $re 1 μPa$ and $P̂L$ = 80.4% for NL = 100 dB $re 1 μPa$ over the 2–10 km search area. The average probability of localization as a function of the SL and NL over all the subarrays was calculated and is shaded in Fig. 7. Areas of black background indicate $P̂L$ = 0%, and areas of white background indicate $P̂L$ = 100%.

Humpback whale song units were estimated to have a median RMS SL of 173 dB $re 1 μPa$ at 1 m measured over 150–1000 Hz and averaged over all NLs (Table I). The 25th and 75th percentiles of the SLs were 169 and 176 dB $re 1 μPa$ at 1 m, respectively. Humpback whale song unit SLs increased as background NLs increased (Fig. 7). The average RMS SL (black, linear fit in Fig. 7) increased 0.53 dB/1 dB increase in the NL (95% confidence interval, 0.52–0.54 dB/dB). The main distribution of the SLs was well above the masking zone.

Another way of looking at the SL data is to group them in different NL bins. SLs were grouped into 5 dB NL bins (Fig. 8, Table I). The song units produced in each NL bin had median SLs that were significantly greater than the median SLs of units in the NL bin centered 5 dB lower (one-sided Wilcoxon's rank sum test, $p≪0.001$ for all four comparisons). Additionally, the variance of the SLs in each NL bin significantly decreased as the NL increased for all but the greatest NL bin, which also had the smallest sample size (one-sided Ansari-Bradley test, $p≪0.001$ for the first three comparisons). Song units produced during the highest NLs were most likely to be affected by masking. The dashed portions of the tails of the histograms in Fig. 8 indicate the SLs that may be artificially suppressed due to masking (i.e., $P̂L<100%$⁠).

Although masking only occurs well into the tail of the SL distributions (Figs. 7 and 8), the potential effects of masking on the measured SLs of the humpback whale song units were investigated. A range of heavier left-tailed SL distributions were simulated to determine how more units in the masking zone might affect the results. The GAM fits to these differently weighted distributions are plotted in red in Fig. 7. The solid red line indicates the GAM fit if no units were missed due to masking, while the dashed red line indicates b = 2. The most likely values are between these two lines since b = 1 is a linear model and produces an improbable elbow in the SL distribution.

The slopes of the SLs as a function of the NL are shown in the lower portion of Fig. 7. As stated previously, the linear fit suggests that humpback whales increase their RMS SLs by 0.53 dB/1 dB increase in background noise in the 150–1000 Hz band. The GAM had the greatest slope at background NLs around 91 dB $re 1 μPa$ with a SL increase of 0.68 dB/1 dB NL. The uncertainty of the GAM fit increases above 95 dB RMS $re 1 μPa$ at which point masking plays a proportionally larger role and the number of song units decreases. But even with the increased effect of masking, the GAM slopes decreased in these highest NLs. The GAM was able to explain approximately 16% of the variability in humpback whale song unit SLs (deviance explained). In addition, the spread of the predicted SLs was smaller than the spread of the observed values. The residuals were approximately normally distributed except for in low NLs when the model had a tendency to overpredict the SLs. No evidence of heteroskedasticity or unmodeled relationships between residuals and SLs was present.

On a per-track basis, the individual humpback whale SL generally increased with increasing background NL. In total, 88 tracks met the criteria of consisting of over 200 song units and spanning a noise range of 10 dB or more. The shortest track had 216 units, and the longest track had 3000 units. The total number of units used for this analysis was 75 157. An example track is mapped in Fig. 9(a) with color indicating the recorded NL along the track. This track had a duration of 9.5 h and contained 864 song units. In this track, the SL increased by approximately 0.19 dB/1 dB increase in background noise [Fig. 9(b)]. The majority of tracks that met these criteria (75 out of the 88 tracks) had increased song unit SLs during periods of increased background noise [Figs. 9(c) and 9(d)]. The median of ΔSL/ΔNL for these tracks was 0.34 dB/1 dB (25th percentile = 0.17, 75th percentile = 0.56), which is less than the aggregate SL response.

Humpback whales changed the intensity of their song units as the background NL changed. When background NLs increased, the intensity of the song units increased. These changes occurred in real time on hourly and daily time scales. The Lombard response was stronger than what was observed in minke whales in this same study area (Helble et al., 2020), which may be because of the larger size of humpback whales and/or the diversity of unit types that they produce. Humpback whales have previously only been studied exhibiting the Lombard effect while producing their social sounds and have been reported to increase their SLs by 0.8–1.5 dB/1 dB increase in the NL (Dunlop, 2016b; Dunlop et al., 2014; Fournet et al., 2018a). The observed average response of a 0.53 dB increase in the SL per 1 dB increase in the NL at the PMRF is substantially less than what has been reported previously. Humpback whales were most responsive to (i.e., greatest slope at) background NLs around 91 dB $re 1 μPa$ with a SL increase of 0.68 dB/1 dB NL. It also appears that humpback whales had a decreased sensitivity to the highest NLs (i.e., decreasing slope). The median song unit SL at the PMRF was 173 dB $re 1 μPa$ at 1 m. The average SLs of social calls reported in the previous studies about the Lombard effect in humpback whales were 131–165 dB (Dunlop, 2016b; Dunlop et al., 2014; Fournet et al., 2018a). The maximum intensity of vocalizations produced by an animal is likely limited by physiology, therefore, since the social calls are less intense than the song units, perhaps the whales have a greater dynamic range to increase their social call SLs than to increase their song unit SLs.

This measured Lombard effect for the humpback whale song is greater than the response for many other species (e.g., Kragh et al., 2019; Sinnott et al., 1975; Zollinger et al., 2017) but less than the reported response for most marine mammals (e.g., Dunlop, 2016b; Dunlop et al., 2014; Fournet et al., 2018a; Holt et al., 2009; Parks et al., 2010). The location of the PMRF array in deep water and covering approximately 1200 km² allowed for a sample size that included over 500 encounters with singing humpback whales and almost 84 000 song units. This song unit sample size is 2 orders of magnitude greater than sample sizes used by previous researchers (e.g., Dunlop, 2016b; Dunlop et al., 2014; Fournet et al., 2018a; Holt et al., 2009; Parks et al., 2010). With small sample sizes, it is possible to have skewed results. The effects of small sample sizes in both low numbers of vocalizations and a small sample of whales were evident when analyzing individual responses to noise. For example, individuals in this study had Lombard responses of up to a 0.99 dB increase in SL per 1 dB increase in NL, and some showed no Lombard response or a negative response (Fig. 9). In addition, shallower water environments that are used in some studies will result in higher uncertainties in the TL estimations due to reflection, refraction, and absorption by the seafloor that are more difficult to model. The deep water in the PMRF study area allowed for selecting song units that had little interaction with the seafloor. Finally, masking can skew results. If lower SL vocalizations were unable to be detected in high noise conditions and masking was not accounted for, then the Lombard response could appear greater than it really was. The search area at the PMRF was chosen so that the probability of localization would be high for all expected SLs and NLs. Further, the effects of masking on the Lombard response were estimated using GAM modeling in case some low SL units were still masked.

SLs of animal vocalizations can be highly variable and dependent on a multitude of factors, including behavior, vocalization type, and noise. For example, behavioral state can impact the SLs of vocalizations. In this study, the humpback whales were singing in contrast to producing social calls. Social calls are thought to be for communication between closely spaced animals (e.g., group coordination, male competition) and, hence, have a lower SL, whereas song may be for longer-distance communication (Dunlop et al., 2008). Different vocalization types within both the song and social call repertoires are likely more difficult to increase in intensity (Au et al., 2006). Also, these results suggest that individuals may respond less when NLs change within a singing bout than when the NL changes between their singing bouts. In addition, location can affect the SLs of vocalizations. These humpback whales were in deep water off of Hawaii and may sing at different intensities than if they were in shallow water (e.g., Au et al., 2006), along their migration route (e.g., Dunlop et al., 2013), or in feeding areas (e.g., Fournet et al., 2018b). The presence of predators may also impact the SLs of vocalizations. Besides these factors, this study and many other studies have found relationships between SLs and NLs. It was not surprising to find that that NLs were only able to explain a small portion of the variability in the SLs because so many other factors also likely influence the intensity of singing. But, because these NLs and SLs are so intricately related and easy to objectively measure (compared to behavior), the SLs of animal vocalizations should always be reported along with the associated NLs.

Since humpback whales in this study did not fully compensate for increasing noise, their communication space was reduced during periods of high noise. To demonstrate this reduction in communication space, an assumption was made that humpback whales need −5 dB of SNR in the 150–1000 Hz band to effectively communicate (assuming a greater SNR is needed for information transmission than for simple detection). The TL was modeled by using the geometrical spreading equation at close ranges [Eq. (5)] and the Peregrine model for ranges beyond 20 km when the propagation of sound is highly dependent on bathymetry, depth, and sediment type and thickness. When using the geometrical spreading equation, the assumed transition depth was 4000 m (average water depth), and when using the Peregrine model, the assumed environmental conditions were the same as those stated in Sec. II D. With these assumptions in mind, if a humpback whale produced a song unit at 168 dB $re 1 μPa$ at 1 m in 82.5 dB NL (averaging the 80–85 dB NL bin limits), the maximum allowable TL for information transmission would be 90.5 dB [using Eqs. (3) and (4)]. The communication range in these low NLs would depend on the direction of the intended communication, but, for example, for deep ocean propagation to the northward direction of the PMRF, communication range predicted by the Peregrine model would be approximately 90 km (using the environmental assumptions from Sec. II D). In contrast, if a humpback whale produced a song unit at 177 dB $re 1 μPa$ at 1 m in 102.5 dB NL, the maximum allowable TL for information transmission would be 79.5 dB, and the communication range would be approximately 20 km (using either the Peregrine model or the geometrical spreading equation). Therefore, an increase in the NL of 20 dB is likely to decrease the communication range (distance) in the deep ocean by 4.5 times and communication space (distance-squared) by 20 times. Since the humpback song is a male vocalization, attraction of mates might be reduced in increased background noise conditions. However, the range of SLs recorded suggests that humpback whales should be able to maintain their communication space at least in the lower NLs of the range recorded at the PMRF but they do not. Perhaps the maximum communication space is not necessary for the purpose of their song.

Several assumptions were made in this study. The NL at the bottom was considered to be a proxy of the NL at the surface. It is likely that the NL recorded by the PMRF bottom-mounted hydrophones is less than the NL experienced by the whale. The noise sources within the frequency band of humpback whale song units are wind and waves, shipping, and biological sources (Wenz, 1964). Since all of these sources originate near the surface, the NL near the surface is likely greater than the NL near the bottom and, therefore, the communication range may be less than predicted. To calculate the SLs of song units, it was assumed that both the source and receiver were omnidirectional. The hydrophones on the PMRF were designed and tested to be omnidirectional, so the omnidirectionality of the receivers is supported. It is possible that humpback whales are not omnidirectional sources. Au et al. (2006) plotted the SLs estimated using a vertical hydrophone array and found that higher frequency harmonics may have vertical directionality. However, for the lower frequencies studied here, directionality is assumed to be minimal. If the song units are directional, the measured RLs and estimated SLs from bottom-mounted hydrophones would be lower than if the recordings were made on-axis. The SL and NL values reported in this study were limited to the frequency band containing the main components of the humpback whale song. This band was selected because experiments with other taxa have shown that the noise band that influences the SL the most is the band that covers the main frequencies of the animal vocalization (e.g., Hage et al., 2013; Halfwerk et al., 2016; Manabe et al., 1998). However, these controlled exposure experiments have been done with terrestrial taxa, and more work is needed to understand how marine mammals respond to noise in bands outside of their primary vocalization band.

Singing humpback whales responded to increasing noise by increasing the SLs of their song units, but they did not fully compensate for increasing background noise. The humpback whales studied off of Kauai sang at a greater intensity when the background NL increased. Increases in ambient noise could reduce humpback whale male communication space in the important breeding area around the Hawaiian Islands. It is unknown how this decrease in communication space could affect the population. These observed effects will help contextualize effects of anthropogenic noise sources, which in this area includes U.S. Navy training activities.

The Lombard responses of whales should continue to be studied. More work is required to determine the factors that are influencing the range of Lombard responses observed in other locations. Species differences, the behavior that the whale is engaged in, the location of the whale, and the noise source may influence a whale's response to noise. In addition, the purpose and function of marine mammal vocalizations need to be better understood to predict the effects of these behavioral changes. This work is necessary to fully assess the impacts of noise on marine mammals.

This research was supported by the Office of Naval Research (Code 322, Marine Mammals and Biology), Commander, U.S. Pacific Fleet (Code N465JR), and the Naval Facilities Engineering Command Living Marine Resources Program. This paper was presented along with results from Helble et al. (2020) at the Ocean Sciences Meeting in San Diego, CA, in February 2020. The authors would like to thank Glenn Ierley for his helpful feedback, which improved this manuscript.