Spring haul-out behavior of seals in the Bering and Chukchi Seas: implications for abundance estimation

Data collection

For this study we used haul-out behavior data and location estimates from bio-loggers deployed on bearded, ribbon, and spotted seals in the Bering, Chukchi, and western Beaufort seas by multiple organizations as part of collaborative investigations from 2005 through 2021. Seals were captured using nets and bio-loggers were attached during studies based in coastal communities or on research ships (Fig. 1). Ship-based capture events occurred during spring near the southern ice edge in the Bering Sea between 2005 and 2018. Land-based capture events occurred between 2005 and 2020 from May to October, generally between the Alaska coastal communities of Scammon Bay in the Bering Sea, Utqiaġvik in the Chukchi Sea, and Nuiqsut in the Beaufort Sea. Data from additional deployments along the Kamchatka peninsula in the western Bering Sea are also included. We refer readers to the primary literature for detailed capture and bio-logger attachment methods (see publications listed in Table S1). NOAA-led research was conducted under the authority of Marine Mammal Protection Act Research Permits Nos. 782-1676, 782-1765, 15126, and 19309 issued by the National Marine Fisheries Service, and Letters of Assurance of Compliance with Animal Welfare Act regulations, Nos. A/NW 2010-3 and A/NW 2016-1 from the Alaska Fisheries Science Center/Northwest Fisheries Science Center Institutional Animal Care and Use Committee (IACUC). ADF&G and NSB field work was covered by Research Permits Nos. 358-1585, 358-1787, 15324, and 20466 and by ADF&G IACUC permits Nos. 06-16, 09-21, 2014-03, 2015-25, 2016-23, 0027-2017-27, 0027-2018-29, 0027-2019-041.

We created a subset of haul-out behavior data from 250 bio-loggers deployed on 35 bearded, 110 ribbon, and 105 spotted seals to include only records from 1 March to 15 July between 2005 and 2022. Bio-loggers were of the ‘SPLASH’ or ‘SPOT’ family of tags developed by Wildlife Computers (Redmond, WA, USA). Deployments consisted of either a single ‘SPLASH’ device, a single ‘SPOT’ device, or both types. Devices were either adhered to the hair on the seal or attached through the rear flipper inter-digital webbing. The use of bio-loggers adhered to the back or head provides some benefits over flipper mounted devices (e.g., increased satellite transmittal rates, locations at sea) but these fall off during the following annual molt, which, depending on deployment date, limits the duration of haul-out data they provide especially during the focus months of our study. Additionally, bio-loggers attached to the head or dorsal region are often dry while the seal is floating at the surface, inducing a slight positive bias in the hourly percent-dry values reported by the bio-logger. For this study, in cases where both bio-logger types were deployed, we preferred hourly percent-dry observations from the flipper tag. All data were transmitted by the deployed instruments via the Argos satellite network and location data were either derived from Argos transmissions or transmitted FastLoc GPS data.

Sex and age class (non-dependent young-of-the-year, sexually immature subadults, and mature adults) were estimated at the time of deployment by various combinations of length, claw growth ridges (McLaren, 1958), and pelage characteristics. Seals determined to be less than one year were classified as young-of-the-year. For those bio-loggers deployed on young-of-the-year and transmitting into the next year (6 ribbon seals; 3 spotted seals), the age class was advanced to subadult on 1 March of the following year—the assumed anniversary of their birth. Subadults are those seals likely greater than one year of age but less than four years. Adults are individuals that are likely older than four years. Table 1 provides a summary of these deployments and data received from them.

Haul-out behavior data were recorded in a manner standard across Wildlife Computers bio-loggers and transmitted via the Argos satellite network as hourly percent-dry timelines. For each hour of a day, the wet/dry sensor was polled by the tag firmware every few seconds and the percent of the hour in the dry state was calculated (Fig. 2). On board the bio-logger, hourly percent-dry calculations were rounded to the nearest 10% inclusive of 0% and 100% along with additional values at 3% and 98%. This compression resulted in additional data transmission as each message consisted of two complete 24-hour records. Memory capacity allowed caching of percent-dry records for several weeks or months and each message was transmitted several times to ensure reception at the satellite. Bio-loggers were deployed and programmed in a manner to maximize data transmission during the spring pupping and molting period, though hourly percent-dry data were not always successfully transmitted. This is due to a variety of factors including satellite coverage, tag availability (e.g., tags mounted to the rear flipper often do not transmit while at sea), tag performance, duty cycling, and atmospheric interference. Fortunately, missing records do not substantially bias inference about haul-out probabilities (Conn et al., 2012).

Tags that fell off due to molt, attachment failure, or seal mortality and remained on ice or land may have continued to send data to satellites; i.e., a detached bio-logger that is dry (either on ice or land) will record and transmit data suggesting the seal is hauled out. Therefore, end times of each deployment were identified by examining bio-logger locations, percent-dry records, and dive behavior (if available) to determine when bio-loggers ceased providing data consistent with seal behavior. For example, a data record that ended with several consecutive days (∼10+ days) of 100% dry observations and with locations indicating the tag was on land were truncated to the final stretch of 100% dry observations. The vast majority of deployments ended with the device detaching in the water and the deployment end date was obvious. There is no perfect algorithm for identifying deployment end dates and each deployment in question was considered separately. While not perfect, we are confident our reliance on expert opinion and examination of multiple data streams provided the best option. Data outside of the deployment start and end times were discarded prior to analysis.

Of key interest in this study was the relationship between haul-out behavior and weather covariates that vary with time and seal location. The use of modern bio-loggers that record and transmit behavioral data while simultaneously providing location estimates was key to this objective. We explored the use of a continuous-time correlated random walk (Johnson et al., 2008) movement model to predict locations at specific times. However, the sparse nature of data from some bio-loggers, especially those mounted to the rear flipper, resulted in poor modeling performance or convergence issues. For this study, we calculated a weighted average daily location where the inverse of the estimated Argos or FastLoc GPS location error was used for the weight. Each Argos location estimate was assigned an error radius based on either the categorical location quality (3 = 250 m, 2 = 500 m, 1 = 1,500 m, 0 = 2,500 m (Lopez et al., 2013); we chose 2,500 m for location classes A and B) or, when available, the estimated error radius from the Argos Kalman filter algorithm. Location estimates from FastLoc GPS were all assigned an error radius of 50 m. On days when haul-out observations were present but location data were missing we used the seal’s last calculated weighted average daily location; days when the location intersected with land were removed from the seal’s record. We recognize that bearded and spotted seals haul out on land. However, assessing the relationship between haul-out behavior and weather covariates and seals’ availability for aerial surveys on land was outside the scope of this study. Additionally, any daily locations on land were likely more reflective of coordinate averaging and measurement error, rather than actual use of coastal features.

Explanatory variables

In addition to sex and age class, we analyzed variables that might help explain variation in haul-out probabilities. These included day-of-year (for seasonal effects) and local solar hour (for diurnal effects). We calculated local solar hour using the {solaR} package (Perpiñán, 2012) within the R statistical environment (R Core Team, 2021) based on the weighted daily average locations. We also linked the weighted average daily locations to weather values from the North American Regional Reanalysis (NARR) model produced by the National Centers for Environmental Prediction (Mesinger et al., 2006). The NARR model assimilates observational data to produce a long-term picture of weather over North America and portions of the surrounding seas. Weather variables are made available across the region 8 times daily. For this study, NARR weather values were limited to the extent of our study area over the Bering, Chukchi, and Beaufort seas at 3-hr intervals based on the original grid resolution of 32 km (1,024 km²). The following weather variables are known to affect haul-out behavior in other Arctic pinnipeds (Reder et al., 2003; Udevitz et al., 2009; Perry, Stenson & Buren, 2017) and were interpolated and assigned to daily seal locations using a bilinear method: (1) air temperature at 2 m above the Earth’s surface, (2) wind consisting of northerly and easterly vector components converted to wind speed using the Euclidean norm, (3) barometric pressure at sea level, and (4) precipitation(Table 2).

For all seal species, we considered the following variables when modeling the hourly haul-out behavior: day-of-year, solar hour, temperature, wind speed, barometric pressure, precipitation, and wind chill (represented by a wind:temperature interaction (Udevitz et al., 2009). Ribbon and spotted seal models included age-sex class and interactions between day-of-year and age-sex class, but we omitted these from bearded seal models due to poor representation of age-sex classes (Table 1). Bearded seal models included a latitudinal effect (and an interaction with day-of-year), because bearded seals occupy a substantial range during the spring and we were interested in possible differences in the timing of haul-out behavior along a latitudinal gradient. We omitted the latitudinal effect from ribbon and spotted seal models because, during the spring, these species are most prevalent near the southern ice edge in the Bering Sea (Conn et al., 2014).

Haul-out modeling

Haul-out records for seals are often characterized by sequential hours spent hauled out on ice alternating with long periods in the water (Fig. 2). Commonly used statistical models for binary data (e.g., logistic regression) assume independence among responses, an assumption that is clearly violated if hourly responses are modeled. Any analysis that ignores temporal autocorrelation in responses will thus have overstated precision (Betts et al., 2006).

To properly account for temporal dependence and to take advantage of computational efficiency, we used generalized linear mixed pseudo-models (GLMPMs; Ver Hoef, London & Boveng (2010)) to model variation in haul-out behavior as a function of (1) covariate predictors, (2) temporally autocorrelated random effects, and (3) individual random effects representing heterogeneity in individual behavior. We used the glmmLDTS package (Ver Hoef, London & Boveng, 2010) to implement GLMPMs. We explored two different model formulations for our data and we fit separate models to bearded, ribbon, and spotted seal data sets as we expected differing behavior by species. Separate models for each species were also needed because a single, very large data set proved computationally intractable. In our first model formulation and for each species, we fitted a year-independent model that predicted average haul-out behavior as a function of demographic, weather, seasonal, and diurnal effects. Second, for ribbon and spotted seals (which had considerably more data than bearded seals), we fitted models that included all the effects from the first model, but also permitted annual variation in haul-out timing. This second set of models was used to examine whether haul-out patterns varied by year and to determine the annual timing of apparent peaks in haul-out behavior. For both models, we assumed an hourly Bernoulli response (a binary 0-1 response dependent upon whether the tag was greater than 50% dry for a given hour) where the linear predictor was modeled on the logit scale. This is consistent with previous approaches (London et al., 2012; Ver Hoef et al., 2014) and only 7.0% of our observations fell between 10% and 90% hourly percent-dry.

We followed Ver Hoef et al. (2014) in using linear, quadratic, and cubic effects of day-of-year to represent temporal changes in behavior during the season. However, unlike previous models for harbor seals (London et al., 2012) and ice-associated seals (Ver Hoef et al., 2014), which treated hour-of-day as a 24-level categorical variable to capture diurnal cycles, we adopted a continuous formulation based on Fourier series that provides a flexible model while preserving the inherent circularity needed for time-of-day effects (i.e., hour 0 should be equal to hour 24). It also represents hour-of-day with 6 parameters, which is a considerable reduction when compared to a 24-parameter variable. According to this approach, we used the following specification for hour-of-day effects: ${\displaystyle H_t={\alpha }_{1}cos\left(\frac{\pi t}{4}\right)+{\alpha }_{2}sin\left(\frac{\pi t}{4}\right)+{\alpha }_{3}cos\left(\frac{\pi t}{6}\right)+{\alpha }_{4}sin\left(\frac{\pi t}{6}\right)+{\alpha }_{5}cos\left(\frac{\pi t}{12}\right)+{\alpha }_{6}sin\left(\frac{\pi t}{12}\right)}$ where H_t gives the effect for solar hour t and α_i are estimated parameters (regression coefficients).

For the second set of models examining inter-annual variation in sea-ice use, we fitted models with year by day-of-year interactions. However, in this case we only included year:day and year:day², omitting the main effects of year as well as year:day³ interactions because models with the latter effects were numerically unstable. However, the modeled interactions were sufficient to allow shifts in haul-out distribution, as one can show mathematically that a simple horizontal shift in timing of haul-out distributions does not affect the main effects or cubic terms in a polynomial regression model. Bearded seals were not included in this examination of inter-annual variation because of limited data across many years in the study.

A typical model fitting exercise would also include a model selection process. However, AIC (and similar criteria) is not suitable when using pseudo-likelihoods, because pseudo-data generated in the model fitting process (Ver Hoef, London & Boveng, 2010) differ between models (Ten Eyck & Cavanaugh, 2018). After fitting GLMPM models, we instead used ‘type III’ F-tests to calculate p-values (Ver Hoef, London & Boveng, 2010) to evaluate model performance and important terms. We also produced predictions of haul-out behavior as a function of three influential predictors (e.g., solar hour, day-of-year, age-sex). Weather covariates for these predictions were based on daily or hourly smoothed weather covariate values across the study region. Such predictions were then used to develop haul-out probability surfaces, explore conditional effects of weather covariates, and determine annual peaks in haul-out activity. The timing of peak haul-out behavior was further used to regress against the annual maximum sea-ice extent in the study region. Predictions before 15 March and after 30 June were not included in visualizations or other evaluations to avoid spurious model predictions at the edge of the data range.

Visualizing the marginal or conditional effect of an individual weather covariate (where all other weather covariates are being held at mean values) on haul-out probability was difficult in this analysis because of the collinearity between covariates as well as the spatial and temporal variation across such a large region. The relationship of each weather covariate with haul-out probability, averaged over the other weather conditions, was more variable than model coefficients would imply. That said, important insights can be gained from plots of marginal effects. To create these plots, we predicted haul-out probability across the full range of each weather covariate while fixing hour of the day at local solar noon and day-of-year at 15 May. For the other weather covariate values, we chose not to use a fixed mean value because we expect weather to vary within day over the season (e.g., the temperature at solar noon will gradually increase from March through June). To account for this, we fit a simple generalized additive model for each weather covariate with smooth terms for day-of-year and solar hour. We used predicted values from the generalized additive model in lieu of holding other weather covariates at a fixed mean value which would not capture seasonal change. The visualizations also include vertical lines representing 95% confidence intervals around the predicted haul-out probability to better communicate the variation in model uncertainty.

We also assessed whether the annual variation in maximum spring sea-ice extent in the Bering Sea influenced the seasonal peak of seal haul-out probability. In particular, we used sea-ice concentration data from the Nimbus-7 SMMR and DMSP SSM/I-SSMIS Passive Microwave Dataset, Version 1 (Cavalieri et al., 1996) to calculate maximum sea-ice extent. All sea-ice concentration grid cells (25 km²) in the study area with greater than 15% concentration were counted daily to get the total sea ice extent for each day between 15 February and 15 July across all years. Maximum spring sea-ice extent was simply the largest daily count of grid cells with greater than 15% concentration for each year. A separate regression model, built on the results of the haul-out model, was used to evaluate the relationship between the annual computed peak haul-out day (as the response) with the maximum sea-ice extent (as the predictor).

Bearded seals

Age and sex class were not included in the model for bearded seals due to our lower sample size for adult and young-of-year age classes. As such, results are shown for all ages (Fig. 5; see also Fig. S1). Additionally, after approximately 9 months, 7 devices deployed on the rear flipper of bearded seals reported implausible hourly percent-dry data (100% dry for several weeks but indicative of movement and increasing transmission rates (see Boveng & Cameron, 2013)). All data after the first instance of unrealistic values were censored from this analysis. In addition to a peak around local solar noon, the bearded seal model predicted additional haul-out activity around local midnight. In concert with the lower magnitude of haul-out probability, bearded seal haul-out behavior was also more protracted throughout the spring season compared to ribbon and spotted seals (see below). Overall, bearded seals were less likely to haul out and had a bi-modal distribution of haul-out probability across the day.

When exploring the influence of weather, bearded seal haul-out probability was strongly affected by wind (F_1,42728 = 130.468; p = <0.001) and temperature (F_1,42728 = 19.5; p = <0.001) with much higher haul-out probability during periods of higher temperatures and low wind speeds (Fig. 6). Not surprisingly, wind chill (F_1,42728 = 14.54; p = <0.001) was also important. Barometric pressure (F_1,42728 = 7.779; p = 0.005) was also significant factor although less apparent (Fig. 6). Any effect of precipitation was not a significant influence on haul-out probability (F_1,42728 = 0.519; p = 0.471).

Ribbon seals

Ribbon seals exhibited a pattern of gradually increasing haul-out probability in April that peaks in late May for subadults and in early June for adults (Figs. 7; see also Fig. S2). There is an apparent seasonal progression with subadults hauling out earlier in the season followed by adult males and, then, adult females. Haul-out behavior was clearly centered around local solar noon and expanded to other hours later in the spring as seals entered their molting period. Subadults showed an earlier start and more intense haul-out activity in April and May. The young-of-the-year records began after weaning and the model predictions seemed to indicate development of in-water activities (e.g., swimming, foraging) in May and, like adults, haul-out behavior was centered around solar noon. Adult females had a more protracted haul-out season compared to males, and more time was spent hauled out in June compared to adult males and subadults.

The haul-out probability for ribbon seals was mostly influenced by temperature (F_1,99540 = 6.87; p = 0.009) and wind (F_1,99540 = 49.314; p = <0.001) with barometric pressure (F_1,99540 = 3.446; p = 0.063) having a milder impact. Ribbon seals were more likely to haul out when temperatures were relatively warm and less likely to haul out at higher winds and lower pressure values (Fig. 8). Neither wind chill (F_1,99540 = 1.83; p = 0.176) nor precipitation (F_1,99540 = 0; p = 0.989) were a significant influence on haul-out probability. Compared with bearded seals, the effect of weather covariates on the predicted haul-out probability for ribbon seals was less striking. Because our ribbon seal model included age and sex class, we can visualize the different influences of weather covariates on those classes and see that subadults differ from adult males and females (Fig. 8).

Spotted seals

Compared to ribbon seals, spotted seals showed a longer spring haul-out season that was less intensely centered on solar noon (Figs. 9; see also Fig. S3). Adults of both sexes spent considerable time in April, May and June hauled out. Adult spotted seal males had a more protracted haul-out season compared to females, and more time out of the water in June (Fig. 9). This likely reflects the triad behavior in spotted seals when suitor males haul out with nursing females. As with ribbon seals, the young-of-the-year records began after weaning and the model predictions reflected development of in-water activities (e.g., swimming, foraging) in May and haul-out behavior centered around solar noon.

Spotted seal haul-out behavior was less affected by the weather covariates compared to ribbon and bearded seals but their influence on the model was still significant in some cases. Temperature (F_1,115189 = 5.384; p = 0.020), wind (F_1,115189 = 45.718; p = <0.001), and barometric pressure (F_1,115189 = 9.445; p = 0.002) were all significant. Spotted seals were more likely to be on the ice when temperatures were relatively warm and less likely to haul out at higher wind speeds. Wind chill (F_1,115189 = 0.72; p = 0.396) and precipitation (F_1,115189 = 0.773; p = 0.379) were not as influential as the other covariates. Differences in the magnitude of response between the age-sex classes were present and consistent across each of the weather covariates (Fig. 10). There was a consistent ranking of adult males being the most likely to haul out, followed by adult females, and, then, subadults. This differs from ribbon seals which showed more overlap between adult males and adult females and that subadults were most likely to haul out across the presented range of weather covariates.

Annual variation in haul-out timing

The second set of models, which included annual variation in haul-out patterns, uncovered significant contributions for linear and quadratic interactions between day and year for only spotted seals (day:year, F_15,115144 = 4.445; p = <0.001; day²:year, F_15,115144 = 5.854; p = <0.001). Ribbon seals showed no significant contribution for interactions between day and year (day:year, F_10,99510 = 0.516; p = 0.880; day²:year, F_10,99510 = 0.549; p = 0.856). Predicted distributions of haul-out activity were largely unimodal, but varied some among and within years with respect to the timing and magnitude of haul-out peaks (Fig. 11). It is important to note that predicted variation in annual haul-out patterns likely reflected both process error and sampling variability. While we did remove any years where only one deployment in a species + age:sex group was present, there were still some years where the pattern shown was informed by a small number of individuals that may not represent population-level patterns.

The annual timing of peak haul-out probability for ribbon seals and adult male spotted seals appeared to have no relationship with the amount of yearly maximum sea-ice extent in the study area. For ribbon seals and adult male spotted seals, p-values were substantially larger than 0.05 (ribbon seal adult females; R² = 0.004, p = 0.896; ribbon seal adult males: R² = 0.059, p = 0.693; ribbon seal subadults: R² = 0.007, p = 0.828; spotted seals adult males: R² = 0.004, p = 0.889). Adult female and subadult spotted seals both showed a negative trend (peak haul-out occurs later in years with less sea ice) but neither with a significant relationship for adult female spotted seals(spotted seal adult female: R² = 0.456, p = 0.141; spotted seal subadults: R² = 0.369, p = 0.062).