A statistical acoustics approach for estimating population-scale bowhead whale migration speed and direction

The Bering-Chukchi-Beaufort (BCB) population of bowhead whales, Balaena mysticetus, is of great cultural importance to the Alaskan Natives of the North Slope Borough, the Northwest Arctic Borough, and the Siberian Yupik of Saint Lawrence Island. In the fall, the whales' westward migration in the Beaufort Sea supports a subsistence hunt by three villages (Utqiaġvik, Nuiqsut, and Kaktovik). The majority of the BCB population typically summers in the eastern Beaufort Sea (e.g., Moore and Reeves, 1993; see also Citta , 2021). Beginning in late August, whales travel westward along the North Slope of Alaska, heading for their nominal overwintering grounds in the Bering Sea. The fall migration in the Beaufort Sea is generally close and parallel to shore, mostly in water depths of 20–50 m. Certain aspects of the migration are believed to have changed over the past 20 years, with demonstrated shifts in migration timing (Szesciorka and Stafford, 2023) and in the whales' winter distribution (Citta , 2021; Citta , 2023; Insley , 2021).

Past studies of bowhead whale traveling speeds during migration have reported a wide range of average values, from 0.8 to 3 km/h (e.g., Richardson , 1987; Mate , 2000; Martin , 2023), to the more commonly reported 3–5 km/h (e.g., Rugh and Cubbage, 1980; Davis , 1983; Würsig , 1983; Wartzok , 1989), and beyond, up to 7 km/h (Zeh , 1993; Mate , 2000; Citta , 2015). This range of different estimates is related to differences in methodology, particularly the amount of time between consecutive samples used for the speed calculation, which in the above cited studies varied from minutes to week. Also, the whales' behavior, e.g., transiting versus feeding, will have a large effect on the mean traveling speed. Both aerial surveys and satellite telemetry provide much useful information about the timing and/or composition of the bowhead migration, but they are logistically expensive and obtain limited sample sizes. Aerial surveys also face safety concerns and are weather dependent. For this reason, this paper explores whether passive acoustic monitoring can provide additional information about migration speed and direction.

During their travels, bowheads produce a wide variety of sounds, which past work has roughly divided into “simple” frequency-modulated (FM) calls and “complex” calls (Clark and Johnson, 1984; Cummings and Holliday, 1987; Moore , 2006; Blackwell , 2007). During the autumn migration, bowhead whales produce these sounds at a rate between 0.5 and 5.4 (inter-quartile range) calls per whale per hour (Blackwell , 2021), which is too infrequent to permit tracking individual whales across a tracking site, as has been accomplished for more vocal species like humpback whales, Megaptera novaeangliae, or minke whales, Balaenoptera acutorostrata (Helble , 2015; Durbach , 2021). Thus, measuring the direction and the speed of individual migrating bowheads is infeasible using passive acoustic tracking.

This paper demonstrates how a statistical analysis of 7 years of passive acoustic localization data can provide insight into the direction and speed of a large-scale migration, even if individual whales call just a few times per hour. Section II describes the equipment, deployment geometry, data preprocessing, and localization strategies used to generate estimates of call rates as a function of location and time. It then discusses the statistical cross-covariance methods used to infer the seasonal migration direction and speed. Section III presents the estimated migration speeds, broken down by year and location along the migration route. Section IV examines evidence for a consistent migration speed across multiple years and sites, identifies evidence for a smaller-scale eastern movement at certain sites, and uses the data to flag the presence of an unusual calling whale. It then provides a preliminary examination into whether the effective migration speed changes over the course of a single season, when computed over weekly intervals.

Between 2007 and 2014, as part of their exploration activities in the Beaufort Sea, the Shell Exploration and Production Company implemented an acoustic monitoring program to study the effects of industrial activities, particularly seismic airgun surveys, on bowhead whales during the fall migration (August to October). Each year, during these months, arrays of directional autonomous seafloor acoustic recorders (DASARs) (Greene , 2004) were deployed at five sites in the Beaufort Sea between Kaktovik and Harrison Bay, a 280-km span (Fig. 1). Each site consisted of 3 to 13 DASARs deployed in a triangular grid with 7-km separation. This study analyzes data from sites 2 through 5 between 2008 and 2014, while ignoring site 1 due to the much smaller area that it covered. In some years, site 4 contained more than seven DASARs, but the analysis presented here restricts the DASAR participation to the same seven-DASAR configuration at each of the four sites, as shown in Fig. 1. During 2010, site 2 provided no useable data due to pack ice covering the site, so that particular site-year combination is absent from the analysis.

DASARs are autonomous recording packages equipped with an omnidirectional acoustic pressure sensor (sensitivity of −149 dB re V/1 μPa) and two horizontal directional sensors capable of measuring the north-south and east-west components of acoustic particle velocity (Greene , 2004). This arrangement permits the azimuth of received sounds, such as bowhead whale calls, to be measured from individual DASARs. Each time sequence is sampled at 1 kHz, with a maximum usable acoustic frequency of 450 Hz due to antialiasing filter roll-off. Coincident bearings to calls detected on different DASARs are combined via triangulation to yield two-dimensional call positions, from which the range of each call to every DASAR can be estimated (Greene , 2004). This ability to measure bearing from a single point allows a location to be estimated using only 2 DASARs (instead of 3–4 time-aligned nondirectional sensors).

The scale of the acoustic dataset motivated the development of methods for automatically detecting, classifying, and localizing bowhead sounds, exploiting the directional localization capabilities of the DASARs (Thode , 2012). A subset of these data from all years was also processed manually, to serve as a consistency check on the automated results. Both approaches have been extensively described and evaluated in other publications (Blackwell , 2015; Thode , 2016; Blackwell , 2017; Thode , 2017; Thode , 2020; Thode , 2021), but only the automated analysis is used in the current study, because it was the only approach that included processed data from all deployment days and thus permitted cumulative counts for an entire season.

Figure 2 illustrates the statistical analysis strategy at a particular site, with each site being analyzed independently. As stated in Sec. II A, each site has seven DASARs arranged in a triangular array. Figure 2 displays how four “western” DASARs are arranged along a north-south line roughly 3 km west of the vertical centerline of the array deployment, and three “eastern” DASARs are arranged along a north-south line roughly 3 km east of the array deployment vertical centerline.

The analysis begins by defining “eastern” and “western” rectangular regions, or “zones,” at the site, which can be envisioned as a “Start” line and “Finish” line for migrating whales. Each zone is W km wide along the Universal Transverse Mercator (UTM) easting coordinate and 2H km tall along the northing coordinate, with the zone positioned such that the vertical midpoint of the zone aligns with the horizontal line that passes through the geometric center of the array. The vertical boundaries of the zone thus lie ±H km above and below the horizontal centerline. Choosing a value of W is a balance between two requirements: It needs to be wide enough to ensure that an individual whale is likely to make at least one sound within the zone, but it needs to be narrow enough to ensure sufficient time resolution in the subsequent cross-covariance analysis. The baseline value of W is set to 1 km, based on the median cue rate estimates of 1.3 calls/whale/h by bowhead whales (Blackwell , 2021), but the results are also robust when the width is 2 km. The zone half-height value of H is set to 20 km for the baseline analysis; in practice, the exact value of H was found to matter little, even when H becomes infinite. The robustness of the analysis to these parameters eliminates the need to make specific assumptions about the typical calling rates of individual whales.

Both zones are separated by a distance D, as measured from the vertical centerlines of each zone (note that these zone centerlines are not the same as the array centerline). The baseline analysis uses a zone separation of 3.5 km (half the distance between DASAR packages), because previous studies found that over 80% of the calls generated within this range from a DASAR had a >99% probability of being localized by the system, assuming that calls are generated with a uniform spatial distribution surrounding a site (Thode , 2020). Larger zone separations of 7 and 15 km are also analyzed. See Table S1 in the supplementary material for the distance layout of 6 pairs of zones when W = 1 km and D = 3.5 km.

Over non-overlapping time intervals in time increments (t_inc) of 15 min, the number of calls detected within each zone was counted, creating “eastern” and “western” time sequences of counts. All calls are accepted regardless of call bandwidth, provided that the precision of the localization estimate is less than the width of the zone. Figures 3(a) and 3(b) show one realization of the eastern and western time sequences at site 2 in 2008, in 15-min increments, measured from the specific zones labeled by the asterisks in Fig. 2(c) (i.e., “Start 6” and “Finish 6”). Because of the episodic and sporadic nature of the migration, both time sequences show considerable fluctuations with time, as different groups of animals pass across the site. The similarities between the two time sequences are apparent, when viewed across the time scale of an entire season, although the absolute call counts vary between the sites. What is not visible across this time scale is that the pattern of the eastern time sequence [Fig. 3(a)] is shifted slightly earlier in time relative to the western time sequence [Fig. 3(b)], as would be expected for calling whales that are migrating east to west.

Figure 3(c) displays the cross-covariance function (thick solid green line) between the two sequences in panels (a) and (b). One sees that the time delay t_peak associated with the highest correlation is +0.75 h, which means that sequences of calls detected within the eastern zone statistically occur 0.75 h earlier than the western sequences, implying a western migration. If t_peak is interpreted as the most likely time required for a calling group of whales to traverse 3.5 km from the centerlines of the eastern Start to the western Finish zones, then one obtains a migration speed of 3.5 km/0.75 h = 4.67 km/h to the west.

To incorporate more localization data at a site, this analysis can be repeated by selecting the five additional zone pairs shown in Figs. 2(b) and 2(c) and computing additional cross correlation functions, which are displayed as five additional thin solid lines in Fig. 3(c). Since all these functions were computed using the same separation D between each zone pair, their ensemble average can be computed to enhance the migration signature, which is displayed in Fig. 3(c) as a thick dashed line. This final ensemble-averaged cross-covariance function (EACCF) will be the only quantity presented in future figures, unless stated otherwise. The use of all six zones allows all location data within 4.75 km of the vertical array centerline to be incorporated into the analysis.

The analysis can be repeated with larger zone separations D, which increases the number of zone pairs available for averaging and increases the resolution of the migration speed signal; however, at very large zone separations (D > 15 km) the correlation peak vanishes as call counts become uncorrelated, at least for a zone strip width of 1 km.

We have made several assumptions about the acoustic behavior and detectability of bowhead whale calls. We assume that the acoustic behavior of the animals is spatially and temporally consistent over the time scale t_peak and spatial scale D. We also assume that the acoustic propagation and ambient noise conditions are strongly correlated between zones over these scales. Note this assumption is weaker than assuming that every call generated within a zone must be successfully detected and localized; instead, this approach only requires that the ratio of the detection probabilities between the two zones remains constant over the time scale t_peak. The reason why this method is robust to detectability concerns is that the cross-covariance analysis depends not on the absolute numbers of calls detected within each zone, but on the changes in call numbers between time intervals. As long as the propagation conditions, ambient noise levels, and animal behavioral states in both zones remain unchanged over the t_peak time scale, then the relative temporal structure of both zones will be preserved, even if the absolute propagation and ambient noise conditions differ between the zones. As a concrete hypothetical example, if shipping noise is present in the western zone but not in the eastern zone, so that some calls in the western zone are masked, the cross correlation signal will still appear, because as long as the masking effects of the noise remain constant over time, the structure of the cross correlation function (e.g., time delay associated with peak correlation) will remain the same, although the magnitude of the cross correlation will decrease because a fixed proportion of western zone calls remain undetected. However, if the shipping noise is so intense that all calls become undetectable, the technique will fail.

As the spatial separation between the zones (D) reaches very large values, so will the required t_peak interval, and the more vulnerable the analysis becomes to a decorrelation in these factors between zones. Furthermore, at larger zone separations slight variations in swimming speed between individuals crossing the start zone yield larger variances in positions that eventually exceed the strip width in the end zone.

The validity of the approach can be checked by comparing the resulting t_peak,D time delays for different zone separations D, under circumstances where D is small enough that the eastern and western time series remain correlated. If these estimates are related to migration speed, then the quantity D/t_peak,D,which represents the mean population migration speed, should be similar regardless of the particular value of D selected. The values of D, W, and t_inc are interrelated and must be chosen to maximize sample size yet preserve sufficient time resolution in the time sequence to measure a relatively precise time lag. Fortunately, the analysis is robust to the specific values chosen, and so no specific assumptions about animal calling rates are required.

In principle, this entire analysis can be repeated with north-south zones to determine the north-south component of the migration speed, yielding not only a migration speed but also a compass direction. Here, the analysis is restricted to only the east-west coordinate as this migration is known to consist of mostly westerly travel.

In summary, the baseline statistical analysis uses a zone separation D = 3.5 km, a zone width W of 1 km, and a t_inc of 15 min, with additional analyses conducted for D = 7 and 15 km. Supplemental analyses using zone widths between 0.5 and 2 km and time increments up to 30 min differed little from these baseline results.

The next few plots are all computed for a zone separation of 3.5 km; equivalent plots derived from separations D = 7 and 15 km are included in the supplementary material.

Figure 4 displays the EACCF computed at each site [i.e., dashed line in Fig. 3(c)], broken down by year, for a zone separation of D = 3.5 km. Each plotted EACCF has been normalized such that the local maximum is 1. The arithmetic mean (i.e., equally weighted average) of the normalized EACCFs across all years is also displayed as a dashed line. For reasons that become apparent shortly, the average excluding 2009 and 2011 is also displayed. These plots help reveal whether a single site exhibits unusual characteristics relative to other sites consistently across seasons.

Figure 5 reorganizes the normalized EACCFs by year, with the dashed line now representing an arithmetic average across sites. This rearrangement permits one to investigate whether a particular year has unusual features relative to other years. For example, the years 2009 and 2011 exhibit features different than the other years, which is why they have been removed from the averages in Fig. 4.

Figures 6 and 7 consolidate all the previous results by taking the site-averaged or year-averaged curves (black dashed lines from each subplot in Figs. 4 and 5) and replotting the x axis in terms of effective migration speed using the transformation speed = D/t_delay. Analysis results are shown for zone separations of 3.5, 7, and 15 km.

If the peak time delays between eastern and western zones arise from animal movements, then the resulting peaks should be independent of the value of D analyzed, when displayed as a function of speed. Negative speed values indicate travel to the east; positive values travel to the west. The relative heights of two peaks on the same plot indicate the relative fraction of the population contributing to that migration speed.

Figures 4 and 5 show consistent local maxima at +0.75 h for the 3.5-km zone separations, and the supplementary material show maxima of +1.5 to +1.75 h for the 7-km separation (Figs. S1 and S2), and +3 h for the 15-km separation (Figs. S3 and S4). These peaks emerge regardless of whether averaged across years for a given site (Sec. III A) or across sites for a given year (Sec. III B). The positive values for time delay indicate a westward migration, as expected from historical visual surveys, satellite telemetry, and traditional knowledge. When replotted in terms of effective migration speed, all three zone separations line up their peaks at values of 4–5 km/h westward (Figs. 6 and 7). Our estimates are based on observations over a relatively small distance (3.5 km), and a short time interval (a few hours). In addition, our DASAR arrays were in areas where satellite telemetry studies, which largely overlapped in time with our data collection (Olnes , 2020) have shown bowhead whales mainly to be transiting. In view of the variability in travel speed estimates from past studies (see Sec. I), it therefore seems reasonable to compare our values to a study with similar criteria. Rugh and Cubbage (1980), for example, observed migrating whales from a promontory at Cape Lisburne, where whales were followed over a distance of <10 km and <2 h of time. They report an average travel speed for “undisturbed migrating bowhead whales” of 4.7 ± 0.6 km/h, a value very similar to our own results.

Figure 5 also shows that 2009 and 2011 have unusual EACCF patterns. In 2011, no consistent migration signal is flagged, and the cross-covariance structure changes completely for different zone separations [as seen by the different migration speeds displayed in Figs. 7(b) and 7(d)]. The EACCFs these years are highly variable [Figs. 5(b) and 5(d)], and the ensemble averages have reduced values compared to other years [Figs. 7(b) and 7(d)]. The most straightforward explanation of these breakdowns is that 2009 and 2011 had the lowest number of calls detected over the entire 7-year program (Fig. 8). For example, in 2009 and 2011 sites 2 to 5 recorded a total of roughly 8000 calls, vs the ∼19 000 and 24 000 calls detected during 2010 and 2012. The relative paucity of calls may have degraded the cross-covariance analysis for these years [Figs. 7(b) and 7(d)].

A surprise that emerges when examining Figs. 6 and 7 is that a local correlation peak often appears for eastward travel at speeds between 4 and 5 km/h in 2008 and 2009 [Figs. 7(a) and 7(b)], and across sites 3 and 5 [Figs. 6(b) and 6(d)], as indicated by the local maxima visible for negative time delays and swimming speeds. The relative amplitudes (and thus relative population) of the eastward peaks are typically smaller than those of the western peaks, but they are persistent. The eastern migration signal is even stronger when data are averaged across years at sites: Both sites 3 and 5 show eastern travel delay times with correlation peaks nearly equal to those of the western migration [Figs. 6(b) and 6(d)].

Figures 4(b) and 5(b) suggest that the bowhead acoustic data are displaying an unusual pattern at site 3 in 2009. Figure 9, which displays the 2009 site 3 data in a format like that of Fig. 3, reveals that over 900 calls were detected over just 4 h between 18:00 and 22:00 on September 14, which is 41% of the 2200 calls detected at the site across the entire 2009 season.

The sequence of calls in the eastern and western time series [Figs. 9(d) and 9(e)] clearly shows signs of eastward travel, as the calls peak in the western zone [Fig. 9(e)] first. This single event dominates the EACCF for the entire season [Fig. 9(c)].

This situation prompted a review of the spectrograms and associated localizations at site 3 in mid-September. This review revealed that over 29 h, between 17:00 on September 13 and 22:00 on September 14 (when the ambient noise level rose to mask the sounds), a single animal was discovered to be generating around 7 short highly harmonic calls per minute (Fig. 10), or about 420/h. Directional data from DASAR D at site 3 (DASAR 3D) showed that the animal was initially west of site 3 and traveling west on the evening of February 13, but then reversed course sometime during the morning of the 14th and traveled east back toward site 3. It passed nearly south of DASAR 3D at 18:05, which implies it must have passed through zone Finish 6 before 18:00 (see Fig. 2 for relative location of zone finish 6 in relation to DASAR 3D). It reversed westward, presumably passing into zone Finish 6 again, and then resumed its eastward travel, passing south of DASAR 3D a second time (during this bout of calling) at around 19:05 on September 14. The whale's switchback to the west means that it crossed the western zone twice but the eastern zone once, explaining why two peaks are visible in Fig. 9(e) but not Fig. 9(d).

Such call sequences have been documented elsewhere (Stafford , 2012), but are not a common feature of the bowhead whale repertoire during the autumn migration (e.g., Blackwell , 2007). The statistical analysis was thus able to flag an unusual circumstance across years and sites, which, due to low call counts in 2009, had an overbearing effect on that year's results.

The overall results suggest that during the fall migration, a minority of animals temporarily head east along the migration path. Instead of a continuous westward flow of animals, the actual migration is like a river with eddies—with locations or times where some animals mull or reverse course, possibly for opportunistic feeding. This loitering behavior may be becoming more widespread; there is evidence that bowhead whales may be present in the Alaskan Beaufort Sea throughout the summer (Clarke , 2018; Olnes , 2020; Citta , 2021).

All the results presented previously were obtained by cross-correlating time sequences constructed across an entire season [e.g., Figs. 3(a) and 3(b)], which provided enough data samples to yield a clear average migration speed estimate across the season. A reasonable follow-on question is whether sufficient resolution exists in the time sequence estimates to detect changes in the mean migration speed within a season. Such changes could be due, for example, to the successive migration of different age classes of whales (Koski and Miller, 2009; George , 2021) throughout the season.

An example of this possibility is shown in Figs. 11 and 12, which display D = 3.5-km EACCF estimates for sites 2 and 4, computed over 1-week intervals with 50% overlap. Thus, each row in each subplot advances 3.5 days on the calendar and has been normalized so that the local maximum across the row scales to 1.0, to enhance the visibility of any long-term trends, regardless of the exact number of calls detected during that week.

If the time delay between the eastern and western time sequences remained constant throughout the season, then the images in Figs. 11 and 12 should show bright vertical lines that mark the unchanging migration delay throughout the season. This static migration delay is visible at site 2 for 2008, 2012, and 2014 [Figs. 11(a), 11(d), and 11(f)] and the first half of 2013 [Fig. 11(e)]. Figure 12 also displays times of fixed migration delay in panels (a), (e), and (g). However, some images suggest that large interruptions in the normal migration pattern do occur. For example, during the last half of the 2013 season [Fig. 11(e)], the most probable time lag jumped from 0.5–0.75 h to 1.5–2.5 h over the course of 2 weeks from September 1 to 14, before reverting to more typical values after September 14. A similar shift is also observable at site 4 over a slightly earlier time window [Fig. 12(f)]. Oddly enough, this interruption is not clearly visible for site 3 (see Fig. S5 in the supplemental material).

If this shift is not an artifact, multiple changes in the migration could have produced it: The actual swimming speed of the population may have changed, the direction of the migration may have shifted from primarily east-west to north-south (thus slowing the effective east-west migration speed), or the animals may have switched into a different behavioral mode, such as foraging, and thus have swum more irregular, even backtracking, paths within the site.

A statistical approach for estimating the mean fall migration speed of the BCB bowhead whale population has been demonstrated to yield consistent results across sites and most years. The observed typical speed of 4–5 km/h to the west is consistent with at least some previous studies, in which travel speeds of transiting whales were measured over relatively short distance and time intervals. The analysis also suggests eastward travel occasionally occurs during the migration at similar speeds, a conclusion verified for one particularly unusual situation at site 3 in 2009. Finally, some evidence exists that this analysis might detect deviations from the normal migration flow over time scales of a week, a topic that merits future investigation.

While the results here have been restricted to analyzing travel in the east-west direction, this analysis could be extended to incorporate simultaneous measurements of north-south effective speed, yielding the compass direction of the migration in addition to its speed.

The ability to infer the mean travel speed and two-dimensional direction of a migrating population in weekly increments has several interesting potential applications. The first is that it would provide an ability to conduct much more rigorous statistical correlations between whale movement and remotely sensed oceanographic and biological features (e.g., wind, ice coverage, phytoplankton), techniques that can currently only be crudely deduced from presence/absence acoustic detections on a single hydrophone. The brief discussion about the anomalous shift in migration speed over 2 weeks at site 4 in August 2013 provides a glimpse of this potential. The second application is the observation of long-term shifts in migration behavior across multiple years, a topic of great interest in a region undergoing rapid climate change. A final application, also illustrated in this paper, is the ability of the technique to flag unusual situations in large datasets, like the atypical calling behavior of a single animal at site 3 in 2009 or the site 4 migration slowdown in 2013.

While this approach was developed for bowhead whales, the fundamental principles may be applicable to other migrating species where individuals are occasionally calling while moving.

See the supplementary material for additional tables on the locations of zones in Fig. 2 and alternate versions of Figs. 4 and 5 that use zone separations D of 7 and 15 km. Versions of Figs. 11 and 12 for sites 3 and 5 are also included.

We dedicate this paper to J. Craighead “Craig” George, a friend and colleague who served as a guiding light and advisor through more than two decades of bowhead research in Alaska. We also thank Nikki Diogou for quickly searching NOAA satellite databases for possible environmental perturbations associated with Figs. 11 and 12. This work was supported by the North Pacific Research Board, Project 2117.

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.