Snowmelt Onset and Caribou (Rangifer tarandus) Spring Migration

Abstract

Caribou (Rangifer tarandus) undergo exceptionally large, annual synchronized migrations of thousands of kilometers, triggered by their shared environmental stimuli. The proximate triggers of those migrations remain mysterious, though snow characteristics play an important role due to their influence on the mechanics of locomotion. We investigate whether the snow melt–refreeze status relates to caribou movement, using previously collected Global Positioning System (GPS) caribou collar data. We analyzed 117 individual female caribou with >30,000 observations between 2007 and 2016 from the Bathurst herd in Northern Canada. We used a hierarchical model to estimate the beginning, duration, and end of spring migration and compared these statistics against snow pack melt characteristics derived from 37 GHz vertically polarized (37V GHz) Calibrated Enhanced-Resolution Brightness Temperatures (CETB) at 3.125 km resolution. The timing of migration for Bathurst caribou generally tracked the snowmelt onset. The start of migration was closely linked to the main melt onset in the wintering areas, occurring on average 2.6 days later (range −1.9 to 8.4, se 0.28, n = 10). The weighted linear regression was also highly significant (p-value = 0.002, $R^2=0.717$). The relationship between migration arrival times and the main melt onset on the calving grounds ($R^2$ = 0.688, p-value = 0.003), however, had a considerably more variable lag (mean 13.3 d, se 0.67, range 3.1–20.4). No migrations ended before the main melt onset at the calving grounds. Thawing conditions may provide a trigger for migration or favorable conditions that increase animal mobility, and suggest that the snow properties are more important than snow presence. Further work is needed to understand how widespread this is and why there is such a relationship.

1. Introduction

The movement of an organism determines the fate of population dynamics, the structure of ecosystems, and the evolution and diversity of life on our planet. While scientists continue to disentangle the relationship between environmental stimuli and associated movements, the mechanisms and consequences of organismal movement remain ambiguous. The proximate and ultimate drivers of animal movement are among ecology’s greatest mysteries [1,2]. A comprehensive understanding of long-distance migration is essential in animal conservation and management.

Migratory tundra ecotypes of caribou and reindeer Rangifer tarandus sspp. (henceforth: caribou) exhibit the longest terrestrial migrations on Earth, typically stretching from the boreal forest to coastal Arctic tundra across the circumpolar north [3]. These arctic areas are undergoing some of the most rapid ecological changes due to global warming, including dramatic changes in snow and ice cover (e.g., [4,5]). In tandem, the majority of caribou populations have been falling, with serious impacts on both natural systems and human communities that rely on those populations [6,7]. Notably, the Bathurst caribou herd in Northern Canada has experienced among the most dramatic declines of any herd, a 98% drop from over 400,000 individuals to fewer than 8000 over the last 30 years [8]. Failure to recover and sustain these populations will have significant negative impacts on northern ecosystems, as this umbrella species is known for playing an essential role in contributing to nutrient cycling, the structure of plant communities, and the abundance of predators and scavengers in Arctic ecosystems [4,9], as well as being of absolutely critical material importance to northern Indigenous cultures [10,11]. Beyond the material importance of caribou for subsistence hunting, food security, and health, communities rely on caribou for cultural connections, identity, spiritual ceremonies, wellness, and livelihoods [12,13,14,15,16].

Although caribou populations are known to cycle dramatically over multi-decadal time frames, the current, more or less synchronous, decline in multiple populations is of high conservation concern, and understanding the role of changes in the biotic and abiotic environment is a major focus of research efforts [7,8]. These relationships are complex and difficult to identify for such a wide-ranging organism in such a seasonally and interannually dynamic environment, but they are central to our understanding of the effect of climate on the survival of these populations.

One mechanism by which climate change might influence caribou populations is by impacting mobility. Freshwater lakes and rivers, highly prominent features in the Arctic tundra [17], have been experiencing altered dynamics due to a rapidly changing climate. Bodies of water are also important potential barriers to the movements of terrestrial wildlife that use ice-covered lakes and rivers for migration and dispersal [18]. Arctic surface air temperatures have warmed at more than twice the global rate [19], resulting in unpredictable increases in precipitation as well as more frequent occurrences of extreme weather events (i.e., freezing rain) [6]. These impacts are amplified during winter, when precipitation (deep snow and/or freezing rain) causes unpredictable access to forage. Freezing rain events have been associated with high mortality in all age classes of caribou [6,20]. Furthermore, such conditions may negatively impact the herd’s ability to move efficiently through their migratory range. Landscape and snowpack conditions may hinder movement and cause caribou to overexert themselves, using more of their energy reserves that would otherwise be used for traveling to calving grounds. Given the requirement to make long-distance seasonal movements, the overall ease of movement across a large landscape is an important and understudied environmental variable [21] with potential consequences for demography and populations [22].

Curiously, previous studies have shown weak, if any, relationships between spring migration timing and snowmelt timing at the respective calving and wintering grounds [4,21,23]. These studies were limited by only considering the presence and disappearance of snow cover, not the immediate dynamics of snowmelt. Furthermore, there is strong indirect evidence that snow quality is of greater importance than snow presence [4,21,22,23]. In particular, deeper and wetter snow is thought to be significantly more energetically costly to move through. Here, we introduce a unique framework for identifying a potential mechanistic driver of caribou spring migration related to snow and ice phenology. We combine Calibrated Enhanced-Resolution Brightness Temperatures (CETB) [24] with Global Positioning System (GPS) tracking collar data of female Bathurst caribou to relate the timing of caribou migration with snowmelt events.

We estimate the onset, duration, and end of spring migration using the hierarchical range-shift analysis (HRSA) method [21]. Melt onset, brief early melt events, and melt duration are evaluated using the melt onset detection algorithm, which is a function of passive microwave brightness temperature and diurnal amplitude variations (DAVs) [25,26,27]. The passive microwave CETB data provide information on the properties of snow and ice at improved resolutions (3.125 km) compared to legacy data (25 km). In contrast to previous studies that examined only the presence and absence of snow and ice [4,21], we considered the physical state of the snow, as passive microwave radiation is highly sensitive to the frozen or liquid state of water in the snowpack [28]. We predict that snow phenology (i.e., the initiation and/or duration of melt and melt–refreeze cycles) is related to caribou migration timing such that the following hold:

• Caribou migrate in response to changes in snow phenology and begin their spring migration once melt onset has initiated.
• Earlier snowmelt leads to earlier caribou migration.
• Longer periods of melt–refreeze cycles delay caribou arrival at their calving grounds, causing longer migration durations.

2. Materials and Methods

We analyzed the spring migration timing of 117 female caribou with over 30,000 observations from 2007–2016 from the Bathurst herd in the Northwest Territories region of Canada (Figure 1) and compared these estimates against the snow property covariates.

We estimated the following population-level spring migration parameters:

• Timing of (individual and herd-level) female caribou departure from their wintering ranges.
• Timing of (individual and herd-level) female caribou arrival at their calving ranges.
• Duration of female caribou spring migration.
• Respective spatial coordinates and areas of wintering and calving ranges.

• These estimates were obtained using a hierarchical range-shift analysis model [21] applied to the female Bathurst herd individuals in this region.

We compared spring migration timing against the following snow property covariates derived from CETB data [24]:

• Timing of early melt and main melt events.
• Timing and duration of snow and ice melting and refreezing cycles.
• Timing of the end of melting and refreezing events.

• We estimated snow properties using a snowmelt detection algorithm [26,27,29] that was informed by brightness temperature data, the process of day/night thaw and refreeze, and refined by distinguishing between early, short-lived melt events and main melt events, defined by a specified number of consecutive days of melt.

We collected all estimates during the spring migration period of female Bathurst caribou, which ranges from 20 April to 1 June annually [30]. We extended this period to approximately two weeks prior and two weeks following migration (1 April–20 June) to account for any migratory activity that occurred outside of the typical spring migration period.

2.1. Caribou Location Data

The Bathurst caribou herd historically ranges over a vast area in Northern Canada (Figure 1). In spring, they migrate to traditional calving grounds in the tundra southwest of the Bathurst Inlet near Canada’s northern coast (66.8°N, 109.5°W) (Figure 1). Their wintering areas have historically included a large area along the tree line separating the tundra barren lands and the boreal forest in areas from Great Bear Lake in the northwest and as far as Northern Saskatchewan in the southeast. However, recent years have seen a significant contraction in the Bathurst range, as their population has fallen from over 400,000 in the 1980s to fewer than 8000 according to recent estimates [31]. Unlike most of the other seasonal movements of the Bathurst caribou, which are highly variable both spatially and temporally, the onset of spring migration is largely synchronous [21].

We used data from a total of 117 female Bathurst caribou with 31,405 associated observations (

Table 1. Summary of female Bathurst caribou herd location data across Northwest Canada 2007–2016 that were analyzed in this study. Data include the year, number of individuals included in the study and observations associated with those individuals. Collar types among animals analyzed include Telonics-GPS/Argos (55.4%), GPS (31.6%), Telonics-GPS/Iridium (10.2%), and Telonics Argos (2.8%).

Year	Number of Female Caribou Analyzed 1	Usable Observations 1
2007	20	557
2008	13	379
2009	15	1118
2010	20	1877
2011	18	3034
2012	21	4578
2013	14	3001
2014	18	4268
2015	32	8199
2016	25	4394
2007–2016	117 unique females	31,405

¹ The female caribou statistics reported here are post-filtering results, which include female caribou with observations that span at least 50% of the spring migration period between 1 April and 20 June of each year.

) in Northwest Canada from 2007 to 2016. The caribou were tracked with GPS collars by the Government of the Northwest Territories’ Department of Environment and Climate Change (GNWT ECC) as part of ongoing efforts to help recover the population [8].

The Bathurst herd was equipped with monitoring devices, including high-frequency radiocollars and Argos satellite tags, which occasionally have large spatial errors of tens of kilometers. Argos collars typically have median intervals between transmissions of up to eight days, and more recent collars (GPS) have intervals as short as two hours [21]. Because we were interested in the migration to calving grounds, which male caribou generally avoid, we analyzed only movement data from adult females. Furthermore, individuals with fewer than 40 days of measurements in an 80-day period (1 April–20 June) were excluded from our analysis. Additionally, those days needed to span most of the period of observation (e.g., not days 1–40). Additional inconsistencies were addressed by implementing a migration analysis that was designed to deal with multi-day intervals between location fixes and missing data (more information on the migration analysis follows in Migration Timing Analysis, Section 2.2).

2.2. Migration Timing Analysis

Statistical methods to identify animal migrations have typically focused on analysis at the individual level. For instance, a movement ecologist might analyze individual caribou trajectories to understand how an individual responded to environmental cues while migrating between wintering and calving ranges. Nonetheless, many animal movement studies are currently focused on population-level inference, where multiple individual trajectories pooled at the population level would be analyzed for evidence of consistent behavioral response to environmental cues [32]. While analyzing migrations at the individual level might be useful, it tends to be a prohibitive task, given that a single movement dataset has hundreds of individuals with diverse sampling regimes across many decades [21]. Furthermore, despite individual variability in their wintering and calving ranges, caribou are social animals and exhibit largely collective migrations influenced by one another’s movements [33,34]. To obtain population-level inference, the typical approach is a hierarchical model with random effects for individuals sampled at the population level [32]. In this study, we used a population-level statistical method called the hierarchical range-shift analysis (HRSA) model, which was developed by Gurarie et al. [21] to characterize caribou spring migrations.

We also derived the displacement rates (speed) that individual female Bathurst caribou traveled at each time using the Calculate Motion Statistics tool in ArcGIS Pro version 2.9.4. We used unique collar IDs as track identifiers (i.e., unique entities) to calculate the geodesic distance (the angle between two location points multiplied by the circumference of the earth). Displacement rate is the displacement traveled from observation n to observation $n+1$ divided by the time interval (in hours) between subsequent locations.

2.3. Hierarchical Range-Shift Analysis (HRSA) Model

The movement of an individual animal in time $Z_i\left(t\right)$ is modeled as a possibly autocorrelated and fluctuating ranging process $r_i\left(t\right)$ around a mean process $m_i\left(t\right)$ [21]. The mean process describes the range shift and is specified by the central locations, times (i.e., beginning and duration of migration), and ranging area of movement events. The mean process $m_i\left(t\right)$ shifts linearly from the individual’s wintering grounds range $m_1$ to the individual’s calving grounds range $m_2$ and begins at an unknown time t. The time elapsed during the shift from $m_1$ to $m_2$ is specified as $\Delta t$, which is the duration of the movement event [21]. The HRSA returns estimates of the migration duration, timing (departure and arrival time), location, and variability estimates for a specified year at both the individual and population levels. Because the HRSA models the entire population’s seasonal migration, at a multi-month scale, it is generally highly robust to errors in individual locations and irregular sampling of location fixes. These estimates can be compared to other variables such as the timing of snowmelt, melt–refreeze status, and diurnal amplitude variation in brightness temperatures. The model is fitted to all the migration data in a given year using Bayesian Markov Chain Monte Carlo (MCMC) via the STAN package in R (rstan, [35]). Figure 2 shows an example for four females in 2011. Further details can be found in Gurarie et al. [21], and code to fit the HRSA is now available bundled in the package TuktuMigration at https://github.com/EliGurarie/TuktuMigration (Version 0.0.0.1; Accessed on 16 April 2022).

2.4. Brightness Temperature Data

The Calibrated Enhanced-Resolution Brightness Temperature (CETB) Earth System Data Record (ESDR) [24] uses the radiometer version of the Scatterometer Image Reconstruction (rSIR) method to produce enhanced-resolution gridded data at 3.125 km. The CETBs include twice-daily (separate day and night observations) calibrated brightness temperature ($T_B$) data for observations from all available passive microwave radiometers [36], including Special Sensor Microwave/Imager (SSM/I) instruments on a series of the Defense Meteorological Satellite Program (DMSP) satellites, which have near-global spatial coverage. The rSIR method incorporates all measurements from a given sensor that overlap a grid cell for the selected time period, weighted by the antenna Measurement Response Function [36].

In this study, we used the DMSP-F15 SSM/I rSIR 37 GHz vertically polarized (37V GHz) channel observations from February 2000 to August 2021. The rSIR algorithm transforms swath data to enhanced-resolution grids by defining a group of multiple overlapping pixel measurements in the neighborhood of each grid cell [24,36]. The brightness temperature is estimated at each pixel from the overlapping observations, weighted by the Measurement Response Function (MRF) of the channel and rotation of the observation footprint. The MRF is estimated using the antenna gain pattern, antenna scan angle, and the integration period of the sensor channel [36]. The CETB data were combined as gridded spatio-temporal arrays (later referred to as data cubes) of brightness temperature, average time (of input values), and standard deviation of contributing measurements used to estimate the brightness temperature at each pixel. We created subsets of the CETB data cubes by selecting all pixels within specified latitudinal, longitudinal, and temporal ranges.

As reported by Hilburn [37], F15 SSM/I radiometer instrument retrievals were degraded after 14 August 2006 by activation of the radar calibration (RADCAL) beacon [37]. While strong RADCAL interference was detected on the SSM/I 22 GHz vertical polarization channel (22V), no effect was reported in the 37V GHz channel, which we used for this study.

Rather than using data from multiple sensors, we used the F15 SSM/I observations because the observations included the full period for which we had obtained caribou location data. The F15 orbital nodes have drifted several hours earlier in the day over the sensor lifecycle, from approximately 11 A.M./7 P.M. (morning/evening) local overpass times in 2000, to approximately 5 A.M./1 P.M. (morning/evening) in 2015 at 70°N locations [38]. The drift in observation time will have impacted brightness temperature magnitude for some dates, as does the shift in daylight throughout the season. There is a small shift in $DAV$ distribution for the second half of the period of observation that may have impacted the end of high $DAV$ but should not have impacted the melt onset date because at melt onset, $DAV$ is typically many times the threshold value.

2.5. Snowmelt Onset Detection Algorithm

We used twice daily CETBs to identify melt onset and intervals of melting and refreezing of the landscape. At a given frequency, the brightness temperature ($T_B$) is defined as a measure of received radiation intensity from a surface and behaves as a function of emissivity ($ϵ$) and the physical temperature of the surface ($T_s$) such that:

$$\begin{array}{ccc}\hfill T_B& \sim & ϵT_s.\hfill \end{array}$$

(1)

Snowmelt can be detected by passive microwave sensors due to a change in the signal between dry snow and wet snow, which emits energy close to that of blackbody energy levels [39,40]. The presence of water in a snowpack increases emissivity, resulting in higher brightness temperatures. Surface melt is typically captured by higher-frequency wavelengths (e.g., 37V GHz), as they are sensitive to shallow depths of the snowpack [39,40].

Diurnal amplitude variations ($DAV$) are calculated as the absolute value of the difference between consecutive morning and evening brightness temperature observations:

$$\begin{array}{ccc}\hfill DAV& =& |{T_B}_{morning}-{T_B}_{evening}|\hfill \end{array}$$

(2)

We calculated $DAV$ using CETB rSIR data from F15 SSM/I. If a $T_B$ observation is missing, due to coverage gaps, for example, at low latitudes, there will be gaps in $DAV$, but this was not an issue at the latitudes of the caribou observations. We used $T_B$ and diurnal amplitude variation ($DAV$) to estimate the onset of the main melt and duration of melt–refreeze cycles. To detect melt onset (the presence of liquid water), melt, and occurrences of melt–refreeze, we defined surface wetness of the snowpack as indicated by simultaneous increases in the $T_B$ and high $DAV$ values [41]. Snowmelt timing was identified as dates when both thresholds for $T_B$ and $DAV$ were met as follows:

$$\begin{array}{ccc}\hfill T_B& >& 247\mathrm{K}\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill DAV& >& 10\mathrm{K}\hfill \end{array}$$

(4)

These thresholds were developed using CETBs to detect melt in snow-covered high elevation SSM/I rSIR pixels in Colorado [27]. Both the melt and end of high $DAV$ (EHD) detection algorithms require discrete thresholds [25], a number of threshold exceedances, and defined observation windows [26,27,42] to determine the day of year and location(s) when melt occurred.

We identified two kinds of melt events: early melt onset and main melt onset (see Figure 3). Early melt is defined as an event when the snowmelt thresholds in a pixel (both Equations (3) and (4)) are met for one day (one out of two consecutive measurements), meaning that a brief period of snowmelt occurred in the morning or evening. Main melt is defined as an event when the snowmelt thresholds are met for five out of seven consecutive days (five out of fourteen consecutive measurements), meaning that a prolonged period of snowmelt occurred (melting during the day and refreezing at night or vice versa). The main melt onset date definition is robust enough to avoid the detection of earlier, brief melt periods that were followed by significant periods of accumulation and refreeze and allows for missing data in lower-latitude locations [27].

The melt detection algorithm determines the number of times the thresholds are concurrently met and identifies the annual melt onset date for each pixel when it has been sustained for the specified number of days.

We also identified the end of high $DAV$, which is the end of the cyclical phase changes between the liquid and solid states of water within the snowpack and is followed by a transition from an ice- or snow-covered landscape to a consistently wet-snow or snow-free landscape. End of high $DAV$, which occurs after melt onset events, was detected using the similar workflow as described for melt onset detection, with the end of high $DAV$ conditions considered satisfied when $T_B>247\mathrm{K}$ and $DAV<10\mathrm{K}$ for seven out of ten consecutive days (seven out of twenty consecutive measurements). Previous studies defined the end of high $DAV$ as the last date when $DAV$ thresholds were met for more than three out of five consecutive days [43]; we found that a requirement of seven out of ten consecutive days of exceeded thresholds was more effective to capture the end of high $DAV$ in this region.

While the melt onset date was assigned to a pixel on the first date of a series of exceedance conditions, the end of high $DAV$ date was assigned to the last date in the window when the end of high $DAV$ conditions were satisfied. Therefore, those conditions would have been sustained for seven out of ten consecutive days, and the last occurrence event would be the end date assigned to a given pixel. This algorithm is based on the observation that melt is followed by a period of high $DAV$, where the snowpack is melting during the day and refreezing in the evening (i.e., a period of melt–refreeze cycles). At the end of these cycles (end of high $DAV$), the snowpack may be actively melting during both the morning and evening until snow disappearance [29], which usually follows quite quickly. An example of the estimation of melt onset and end of high $DAV$ can be seen in Figure 3. In Figure 4, we compare the timing and magnitude of $T_B$ (A), $DAV$ (B), and the latitudinal component of the migration (C) to show the overall synchroneity of the big shifts. Relationships between the passive microwave and caribou datasets were analyzed by comparing the output from each dataset.

2.6. Moderate Resolution Imaging Spectroradiometer (MODIS)

The Moderate Resolution Imaging Spectroradiometer (MODIS) Terra MOD09A1 Version 6.1 products provide estimates of surface spectral reflectance at 500 m resolution [44]. MODIS bands 1 through 7, two quality layers, and four observation bands were used to select a value for each pixel in an eight-day composite of products. This is performed to correct for atmospheric conditions, such as gases and aerosols [44]. These data were used in the visualization of snow and ice extent during the caribou migration and melt onset period.

We used MODIS true color imagery at 500 m resolution to capture landscape conditions at the beginning of each month within the Bathurst caribou migratory period (1 April–1 June). We selected 8-day composites for three defined time intervals (30 March–6 April, 1–8 May, and 25 May–1 June) each year and overlaid the overlapping caribou location data during those time intervals. This analysis was performed in ArcGIS Pro version 2.9.4. In Figure 5, we plot the location of each caribou for the corresponding time period as a function of snow and ice extent annually.

2.7. Comparing Caribou Spring Migration to Snowmelt Timing

Using 37V GHz CETB data, we identified the following:

• Dates of early melt events, which were defined as a brief period of snowmelt. Snowmelt thresholds ($T_B>247\mathrm{K}$ and $DAV>10\mathrm{K}$) were exceeded for just one day.
• Main melt events, defined as prolonged, robust start of snowmelt. Snowmelt thresholds ($T_B>247\mathrm{K}$ and $DAV>10\mathrm{K}$) were exceeded for five out of seven consecutive days. The beginning of the event is the melt onset day.
• The end of high $DAV$ (EHD), defined as the end of a period of frequent fluctuations in snowmelt and refreeze. End of high $DAV$ thresholds ($T_B>247\mathrm{K}$ and $DAV<10\mathrm{K}$) conditions were met for seven out of ten consecutive days.

Using a corresponding decade of multiple daily observations of female caribou position, we identified the onset of movement, interpreted as migration start, and arrival at calving grounds, interpreted as migration end, for female caribou in the Bathurst herd between 2007 and 2016. We then compared these snow and migration observations to each other as well as to corresponding observations of the snow cover extent and disappearance from daily MODIS data.

To capture the landscape conditions experienced by the herd each year, we overlaid female caribou departure and arrival locations onto melt onset and end of high $DAV$ maps in their respective locations. Example female caribou departure locations relative to melt onset are shown for 2011 in Figure 6, where the departure location is used to extract the melt onset date for each year. For each year, we then extracted time series pixels of brightness temperature and $DAV$ at each location and plotted each event onto time series graphs. We used these statistics to determine the timing of female caribou movement events with respect to melt events.

We quantified the relationship between the start of migration against the early and main melt events in the wintering range of the caribou (locations of departure), and the end of migration against the early and main melt events in the calving grounds (locations of arrivals) by performing linear models weighted by the reciprocal of the herd-level standard deviation. We report regression coefficients, R² and p-values. We also report summary statistics of the lag in days between the corresponding migration timings and snowmelt phenology. All statistical analyses were performed in R v4.3.0 [45].

3. Results

A decade of observations of snowmelt timing and caribou spring migration timing have shown strong and consistent patterns. We focused on the 2007–2016 decade, for which there were sufficient migration data and snowmelt timing data. Migration typically began northeast of Yellowknife, Canada (64°N, 114°W), and ended southwest of the Bathurst Inlet in Canada (66°N, 109°W). Migration started between day of year 107 (17 April) and 126 (6 May) and ended between day of year 141 (31 May) and 159 (8 June). Migration durations ranged from approximately 17 to 45 days (see

Table 2. Results of melt detection and the female Bathurst caribou herd movement analysis for 2007–2016. Included here are the number of early melt flags, melt onset date (early and main), end of high $DAV$ (EHD), start of migration, migration duration, and end of migration.

Year	Number of Early Melt Flags	Early MOD 1 (DOY ± Uncertainty)	Main MOD 1 (DOY ± Uncertainty)	EHD 1,3 (DOY ± Uncertainty)	Departure Date 2 (DOY ± s.e.)	Migration Duration 2 (Days ± s.e.)	Arrival Date 2 (DOY ± s.e.)
2007	5	98 ± 0.5	114 ± 0.5	140 ± 0.5	122 ± 1.0	29 ± 1.2	152 ± 2.2
2008	5	100 ± 0.5	119 ± 0.5	127 ± 0.5	120 ± 1.4	36 ± 2.2	156 ± 3.6
2009	4	101 ± 0.5	118 ± 0.5	152 ± 0.5	120 ± 2.5	45 ± 3.5	165 ± 6.0
2010	3	96 ± 0.5	103 ± 0.5	124 ± 0.5	108 ± 1.2	42 ± 2.3	151 ± 3.5
2011	1	117 ± 0.5	123 ± 0.5	128 ± 0.5	126 ± 1.0	17 ± 1.1	144 ± 2.1
2012	4	111 ± 0.5	117 ± 0.5	132 ± 0.5	121 ± 1.0	23 ± 4.8	145 ± 5.8
2013	4	100 ± 0.5	128 ± 0.5	142 ± 0.5	128 ± 1.8	26 ± 4.1	154 ± 5.9
2014	4	109 ± 0.5	116 ± 0.5	141 ± 0.5	118 ± 1.5	26 ± 2.1	145 ± 3.6
2015	6	100 ± 0.5	118 ± 0.5	134 ± 0.5	116 ± 0.8	29 ± 1.1	145 ± 1.9
2016	0	119 ± 0.5	119 ± 0.5	138 ± 0.5	118 ± 0.9	17 ± 1.0	136 ± 1.9

¹ MOD = Melt onset date; EHD = End of High $DAV$; DOY = Day of Year. Uncertainty in days. ² Caribou data in this table represent female animals with sufficient observations during April–June. ³ There is a possibility that a change in observation time due to F15 drift caused an increase in EHD duration, although we did not notice a step change or trend in these values. EHD vs. migration timing was not significant.

).

Time series plots of caribou motion statistics and snow characteristics were created to test whether there was a temporal relationship between weather variables, in particular, snow characteristics such as melt and refreeze, and female caribou herd-level migration timing. Migration (shown here as latitudinal movement, northward migration) to the Bathurst calving grounds showed strikingly similar timing to that of brightness temperature as depicted in Figure 4. The dramatic increase in the female Bathurst herd’s latitudinal movement is a signature of their spring migration to the Bathurst Inlet, where they return each year just before calving season. There was a consistent pattern in their migration across all years, which can be seen in the sudden peaks in latitudinal (and longitudinal, though it is not shown) movement with respect to time. Increases in movement (i.e., spring migratory behavior) appeared to be coincident with increases in both brightness temperature and large $DAV$ increases. The initiation of spring migration aligned relatively well with the period of peak values in $DAV$ across all years. There were some years in particular (i.e., 2010 and 2016) when female caribou were widely dispersed with respect to the latitudinal ranges of their wintering grounds, and were analyzed as two groups, but the initiation of movement appeared to be mostly synchronous for individuals across all years.

Female caribou displacement rates (speeds) were relatively low prior to spring migration, rapidly increased following the start of migration, and peaked as they approach their calving grounds. The motion statistics were used in plotting the time series of brightness temperature and movement data for 2007–2016 (Figure 7, Figure A1, Figure A2, Figure A3, Figure A4, Figure A5, Figure A6, Figure A7, Figure A8, Figure A9 and Figure A10). We selected two years to illustrate in more temporal detail the relationships and different types of melt and migration patterns. The year 2007 (Figure 7A–E) illustrates a year in which there were multiple early melt events and a prolonged period of high $DAV$, which is typically indicative of continuous melting and refreezing occurrences. The year 2011 (Figure 7F–J) illustrates a year in which there was one brief early melt event followed by a main melt onset event and a relatively brief period of melting and refreezing preceding continuous morning and evening snowmelt and transition to open continuously melting (or snow free).

Individual female caribou in the Bathurst herd showed relatively low activity with respect to displacement rate and distance prior to the start of their spring migration. The displacement rate increased rapidly following main melt onset, and the initiation of this increase in activity overlapped with peak values in $DAV$. The female caribou departed shortly after main melt onset events in both cases, but the duration, displacement rate, and total daily distances of their migrations differed significantly. In 2007, the migration duration was approximately 33 days, while the migration duration in 2011 was approximately 16 days. Female caribou appeared to move relatively longer distances at a slower rate during spring migration in 2007 when compared to 2011. Following the calving season, there were multiple instances in which the brightness temperature and $DAV$ signals fluctuated dramatically and appeared to be in synchrony with the distance and displacement rate traveled by the herd during these intervals. This was particularly the case in 2011, which appeared to have multiple dramatic fluctuations in brightness temperature and $DAV$ starting at day 180 (29 June) and ending at day 225 (13 August). During this time, there was a sudden increase in female caribou distance and speed. This was also the case in 2007 on day 210 (29 July), when there was a sudden brief increase in DAV, as well as caribou distance and displacement rate. The snow coverage seemed to be much larger and prolonged during the migration period in 2011 than it was in 2007 as suggested by the frequent early melt events in 2007 compared to 2011.

The timing of migration for Bathurst caribou generally tracked snowmelt onset (Figure 8). In particular, the start of migration was closely linked to the main melt onset in the wintering areas, occurring on average 2.6 days later (range −1.9 to 8.4, se 0.28, n = 10). The weighted linear regression was also highly significant (p-value = 0.002, $R^2=0.717$,

Table 3. Weighted regression statistics ($R^2$, regression coefficient (standard error), and p-value) of Bathurst caribou herd-level migration timing (departure/start and arrival/end) against snowmelt phenology (onset and main melt). The snowmelt was estimated in the wintering ranges for departure, and on the calving ground for the arrivals. Bold-facing indicates highly statistically significant relationships.

	Departure			Arrival
Snowmelt	${\mathit{R}}^{\mathbf{2}}$	$\mathit{\beta }$ (s.e.)	p-Value	${\mathit{R}}^{\mathbf{2}}$	$\mathit{\beta }$ (s.e.)	p-Value
Onset	0.165	0.25 (0.20)	0.244	0.364	0.40 (0.19)	0.064
Main	0.717	0.76 (0.17)	0.002	0.688	0.59 (0.14)	0.003

). The relationship between the migration arrival times and the main melt onset on the calving grounds ($R^2$ = 0.688, p-value = 0.003), had, however, considerably more variable lag (mean 13.3 d, se 0.67, range 3.1–20.4). No migrations ended before the main melt onset at the calving grounds.

Relationships between migration timing and early onset of melting were not statistically significant (Table 3); however, it bears noting that no migration (at the herd-level) began before the early melt onset.

The relationships between melt onset events and the female Bathurst herd spring migration timing are shown in Figure 9A–C, which show melt onset versus the start of spring migration, the end of spring migration, and the duration of spring migration at the herd level. Figure 9A shows the relationship between early melt events (dashed line) and main melt events (solid line) and the female Bathurst herd migration onset in the wintering region.

Across all years, the female Bathurst herd seemed to migrate along the edge of snow and ice as opposed to following a wave of vegetation growth, which can be seen in Figure 5. The entire migration occurs well before any flush of green, in strong contrast to the “green wave” surfing which characterizes the spring migrations of temperate ungulates [46].

There was a dramatic difference in the timing of melt for 2010 when compared to other years. The female Bathurst herd movement appeared to be contingent on the presence of snow and ice, as they migrate further north to their calving grounds.

We also analyzed the spring migration timing of individual female caribou in the Bathurst herd against snowmelt timing between 2007 and 2016. Individual animals did not appear to have as strong of a relationship to early or main snowmelt timing (Appendix B Figure A11).

4. Discussion

The mass spring migration of barren-ground caribou occurs over large distances during a time period that largely overlaps with the dramatic transformation of the sub-Arctic and Arctic environment from being entirely covered with snow to being nearly snow-free (i.e., snow presence). We explored that linkage between the spring migration and the receding snow by analyzing the movements of over a hundred female caribou over a ten year span against the remotely sensed observations of snow phenology (i.e., snow quality or snow melt–refreeze status). We found strong evidence that snow pack melt timing is an important predictor—perhaps even a trigger—of caribou migration. These results are notable, as previous efforts to identify the snow-related triggers of migration largely showed little or no effect of snow cover, depth, or snow water equivalent [4,21,23].

Our results point strongly towards snow properties rather than presence as being the most important for the caribou. Specifically, the bulk of the caribou migration occurs while snow is present on the landscape, including—often—on the calving grounds when the caribou arrive in late May. With their large hooves against a relatively low mass leading to exceptionally low foot loading among other unique physiological adaptations [47], caribou are unsurprisingly supremely adapted for moving in the snow. However, the freeze–thaw dynamics which correspond to the main melt onset correspond to snow that can be both deep and heavy, and particularly energetically costly to move through. By beginning the migration just as that transition occurs, the caribou likely stay ahead of the most difficult snow to move through. An additional explanation is that the snowmelt timing may be a signal of imminent ice break-up. The taiga shield barren lands that cover the bulk of the Bathurst caribou range are dominated by surface waters—lakes and rivers cover roughly 35% of the land area. While caribou are excellent swimmers, it is much more efficient to move on ice. Furthermore, they are generally hesitant to cross water where the melt is incomplete (e.g., pancake ice conditions). The specifics of surface water conditions on caribou movements remain to be studied in more detail.

There was less evidence that the arrival at calving grounds was influenced by the timing of snowmelt on calving grounds and the end of high fluctuations in brightness temperature (DAV) on both the wintering and calving grounds. As demonstrated in earlier studies [21,22], the arrival time at calving grounds appears to be driven largely by the physiological need to arrive in time to calve and is less influenced by weather conditions en route.

We did not find significant relationships between snowmelt characteristics and female caribou spring migration timing when we analyzed the female caribou at the individual level (Figure A11). This may be due to the idea that, as social animals and collective decision makers, caribou do not respond to environmental stimuli alone and may be influenced by the behavior of nearby individuals. The herd-level relationships were tested using the mean location of all individuals, which was based on the density and aggregation of individual female caribou in their wintering and calving regions.

There are many potential energy trade-offs that may occur during migration with respect to landscape conditions, in particular, the pace of migration is influenced by both access to forage and ease of mobility. Landscapes covered by harder snowpacks may restrict access to forage [48], but are mechanically easier to travel through [49], which is consistent with our observation that the duration of migration was shorter when melt onset occurred later on wintering grounds from which they departed. Upon arrival at the calving grounds, however, the need to move at an efficient pace is lessened, as they have reached their destination and may take advantage of regions where forage is more accessible. This is one possible explanation for why caribou continued to move, and the duration of their migration increased when melt occurred later on calving grounds.

MODIS imagery aided in our visualization of movement events with respect to landscape conditions over time. We observed that caribou not only departed shortly after each main melt event had been triggered, but also appeared to move along the edge of the snow and ice-covered landscape as they migrated to their calving grounds. We also observed a slight latitudinal and longitudinal shift in wintering grounds. It appears the herd’s wintering range shifted in the direction (northeast) of the melt onset gradient. In more recent years, the herd has shifted to locations where melt onset typically occurs later, in which landscape conditions would remain advantageous, as they migrate across conditions that are mechanically favorable (i.e., a harder snowpack during spring migration).

These observations provide some evidence of the relationship between snowmelt timing and caribou movement and highlights the need for more comprehensive analysis of these variables in testing these relationships. The adaptive behavior (i.e., migration timing, duration, and locations) of caribou across all years of this study raises the question of whether caribou make the necessary adaptations when faced with environmental variability and unpredictable snow and ice dynamics. These adaptive mechanisms may explain our observed population-level responses to the ongoing consequences of rapid climate change.

5. Conclusions

The start of the female Bathurst herd spring migration as defined by a rapid increase in distance and displacement rate is highly correlated with robust melt events (defined for this region using a main melt event requirement of five out of seven consecutive days).

Migration start is highly correlated with main melt onset date, with an $R^2$ = 0.717 for main melt events, and is not correlated with early melt onset date, with an $R^2$ = 0.165 for early melt events. Overall, earlier snowmelt was correlated with earlier onset of migration, although individual caribou did not always follow this. Migration also took longer when melt was earlier, with possible implications for animal energy levels and health at the calving grounds.

However, the relationship did not prove to be significant when we analyzed individual female caribou responses to snowmelt timing (Appendix B Figure A11). Although the close temporal association does not imply causation, the transition from frozen to wet, and melting and refreezing provide either triggers for migration or favorable conditions that increase mobility. Our approach is novel as an initial attempt to use an enhanced-resolution passive microwave product to derive environmental variables that yield the necessary snow and ice property brightness temperatures to analyze such a relationship. This relationship should be further explored with additional controls for data gaps (i.e., number of caribou observations and inconsistency of readings across years) and compounding evidence of local conditions upon caribou departure, arrival at calving grounds, and during their migratory journey. A causal relationship between snowmelt timing and caribou migration would allow for the anticipation of caribou migratory behavior and potential shifts in herd ranges. Understanding links between these phenomena greatly enhances our ability to sustainably manage and support populations of caribou when considering key regions in the spatial extent of future management plans.