The phenology of the subnivium

The subnivium is a seasonal refuge that exists at the interface between the snowpack and the ground, and provides a haven for a diversity of species to survive extreme winter temperatures. Due to the fitness of many plants and animals being strongly influenced by winter conditions, much attention has been given to changes in the timing of snow cover extent and duration in seasonally snow-covered environments; however, these broad-scale characteristics do not capture the finer-scale dynamics of the subnivium. To study the factors associated with subnivium development, we quantified three critical phenophases of the subnivium: establishment, maintenance, and disintegration along a latitudinal and land cover gradient in the Great Lakes Region of North America. We hypothesized that subnivium phenophases would depend primarily on snow depth and air temperature, but that these would be mediated by latitude and land cover. We found that patterns in both establishment and disintegration were affected by latitude more than land cover, but that variability in the timing of early season snowfall events overrode the effects of both factors in subnivium establishment. In contrast, disintegration was predictably later in more northerly sites, regardless of interannual variation in weather patterns. We found that the subnivium was the result of a balance between ambient temperature, snow depth, and snow density, but that ambient temperatures constrained the system by contributing to the frequency of snowfall and inducing changes in snow depth and density. Areas in lake effect zones, characterized by high snow depths and persistent snow cover, may be the last refugia for subnivia-dependent species given the predicted shifting climate regimes of the 21st century.


Introduction
Snow covers roughly 24 million square kilometers of the Northern Hemisphere annually, but its ephemeral nature causes wide fluctuations in extent and persistence (Lemke et al 2007). Snow cover extent has declined in the Northern Hemisphere over the past 90 years, with the most significant reductions beginning in the 1980s (Vaughn et al 2013). Since this time, snow onset has advanced slightly, but the most significant changes have been in snowmelt date, which has shifted earlier by approximately two days per decade (Chen et al 2016). These changes exhibit wide spatial heterogeneity, however, as the timing of snowmelt has remained relatively stable in North America, but advanced drastically in Eurasia (Peng et al 2013). Although much attention has been given to these broad-scale changes in snow cover (Chen et al 2016, Wang et al 2017, they do not shed light on the conditions of below-the-snow environments. The subnivium is the seasonal refuge that exists at the interface between the snowpack and the ground (Pauli et al 2013), which provides a haven for a diversity of plants (Bjork and Molau 2007), mammals (Korslund and Steen 2006), amphibians (O'Connor and Rittenhouse 2016), birds (Whitaker and Stauffer 2003), and arthropods (Hagvar 2010) to survive extreme winter temperatures. The formation of the subnivium results from continual sublimation and condensation occurring within a snowpack and the upward migration of water vapor from areas of high vapor density (closest to the ground) to low vapor density (near the snow surface) (Marchand 2013, Pinzer andSchneebeli 2009). This movement of vapor reduces the size of the ice crystals in the bottommost snow layer, creating a network of loosely-connected crystals whose low density traps heat released from the soil (Marchand 2013). When snow depths are sufficiently high, the low thermal conductivity of the snowpack insulates the subnivium, creating a warmer and more stable microclimate compared to external air temperatures (Pruitt 2005).
While climate limits the distribution of the subnivium at regional scales, other factors such as vegetation, microtopography, and wind can influence the formation of the subnivium at local scales (Heffernan et al 2014, Petty et al 2015. Snow accumulation and density, both of which directly control the formation and persistence of the subnivium, are affected by ambient temperature, wind, snowfall, and radiation fluxes (Harestad and Bunnell 1981). Deep, low density snow is most effective in maintaining the thermal stability of the subnivium (Ge and Gong 2010). With warmer ambient temperatures, ablation increases, reducing depth and increasing snow density throughout the entire snowpack (Zhang et al 2016). Under colder conditions, the temperature gradient between the bottom layer of snow and the air temperature increases, which increases snow density at the surface of the snowpack (Kim andJamieson 2014, Pfeffer andMrugala 2002). Despite increased density at the surface, air temperatures at or below 0 • C prevent the occurrence of surface snow melt and support the retention of snow depth (Marchand 2013). Wind influences snow depth and density by reducing snow accumulation through erosion and increased sublimation, but can also increase accumulation and compaction through redistribution (e.g. snowdrifts) (Gascoin et al 2013, Pomeroy andGoodison 1997). Finally, the energy balance between net long-wave and incoming short-wave radiation influences snow characteristics as absorption of radiation raises snow temperatures and increases melt rates (Essery et al 2008, Varhola et al 2010. Many local-scale factors influencing snow cover characteristics are modulated by land cover. Canopy cover intercepts up to 40%-60% of snowfall (Hedstrom andPomeroy 1998, Pomeroy et al 1998). Intercepted snowfall is exposed to wind, unprotected from solar radiation, and subsequently vulnerable to sublimation back to the atmosphere (Martin et al 2013). The effects of forest cover on accumulation are particularly important in coniferous forests, where snow can be retained in the canopy for days to months (Pomeroy et al 2002, Pomeroy andSchmidt 1993). Although interception may reduce snow depth, forest cover also promotes snow retention once a snowpack has been developed by providing a buffer against wind and blocking the snow from incoming shortwave radiation (Essery et al 2008, Lundquist et al 2013, Varhola et al 2010. Since the amount of longwave radiation emitted by forests generally does not exceed incoming shortwave radiation (Essery et al 2008), the thermal protections provided by trees can promote longerlasting snow conditions.
We investigated the phenology of the subnivium across three different cover types (coniferous forests, deciduous forests, and open sites) and three latitudinal ranges in the North American Great Lakes Region. While previous work on the subnivium has focused on the role of threshold depth (Pruitt 2005) and density (Marchand 2013), the life cycle of this ephemeral habitat has not been fully explored. Our goal was to quantify, for the first time, three critical phenophases of the subnivium: establishment, maintenance, and disintegration, and to evaluate how latitude and land cover mediate these phenophases. We predicted that subnivium conditions would depend primarily on snow depth and air temperature, with delayed establishment in coniferous sites due to increased canopy interception of snowfall, and in more southerly latitudes due to less snowfall. Once established, however, we expected that the duration of the subnivium would be longest in northerly, forested sites due to later snowmelt dates. Conversely, the duration of the subnivium would be reduced in southerly open sites characterized by shallower snow depths and greater exposure to shortwave radiation. Similarly, we hypothesized that the thermal stability of the subnivium would be lowest in southerly open areas, while northerly areas that have less ambient temperature variation during the winter season would have the most stable subnivia. Understanding the phenology of this seasonal refuge-the timing of its establishment, maintenance, and disintegration-is important for characterizing the current conditions to which species are exposed, and for predicting future exposure.

Study area
The climate of the Great Lakes Region is typified by cold winters and warm-to-hot summers (Peel et al 2007), but varies according to latitude and proximity to the Great Lakes (Andresen et al 2014). The greatest range of temperatures occurs from December to February and the coldest overall temperatures occur in northern interior areas, away from the Great Lakes (Andresen et al 2014). Areas in the lake effect zone (i.e. downwind of the lakes) are typically wetter with more moderate climates and large amounts of snowfall (Burnett et al 2003, Changnon and Jones 1972, Scott and Huff 1996. To capture the range of variation in subnivium phenology across the northern Great Lakes Region, we selected three latitudinal ranges (lower latitude: 42 • -44 • , middle latitude: 44 • -46 • , and upper latitude: 46 • +), and within each examined three cover types (deciduous, coniferous, and open) for a total of nine sites (figure 1, table S1 available at stacks.iop.org/ERL/13/064037/mmedia).

Data collection
During the winters of 2015-2016 and 2016-2017, we deployed three temperature strings at each site, with each string containing four individual temperature probes for a total of 12 probes per site. Each probe was separated by 0.3 meters, and we positioned strings with roughly 5 meters spacing. We staked the strings to ensure that the probes were flush with the ground. Probes recorded ground temperatures at five-minute intervals from November through April. As a result of the microtopography at several sites, 25 individual probes in winter 2015/16 and four probes in winter 2016/17 became encased in ice or submerged in standing water during all or part of the winter season and provided faulty temperature data; data from these probes were excluded from analyses.
We established three weather stations at each site to which we affixed a heated rain gauge for measuring the amount of liquid precipitation (Davis Instruments Corp Rain Collector and Rain Collector Heater), an anemometer for measuring wind speeds (Davis Instruments Corp Anemometer, mounted at mean height of 1.8 meters), and a temperature probe housed in a radiation shield (Davis Instruments Corp External Temperature Sensor, mounted at mean height of 1.6 meters) for measuring ambient temperature. Each rain gauge possessed a wetness sensor (Davis Instruments Corp Leaf Wetness Sensor) to identify the beginning and end of individual precipitation events with greater precision. Similar to the subnivium temperature probes, all weather station instruments recorded data from November through April, but at one-minute intervals. Between winter 2015/16 and 2016/17, we also installed a snow depth sensor (HRXL-Max Sonar WRS Series Ultrasonic Snow Depth Sensor, typical accuracy of 1%) at each site to obtain finer scale snow depth data, which was measured at one-minute intervals from November through April in winter 2016/17. For winter 2016/17, we used our field-measured snow depth values; in winter 2015/16, we used the midpoint of the daily snow depth range from the National Operational Hydrologic Remote Sensing Center's (NOHRSC) Interactive Snow Information map (www.nohrsc.noaa.gov/ interactive/html/map.html). NOHRSC uses a physically-based snow model to produce spatially explicit estimates of daily snow depth, which are reported as ranges of accumulation (e.g. trace to 5.1 cm, 5.1-10 cm, 10-20 cm, 20-30 cm, etc.). We confirmed that depths obtained from NOHRSC were Figure 2. Phenophases of the subnivium applied to a temperature probe at one of the forested study sites during the winter of 2016/17. Subnivium establishment occurs on the day when daily temperature stabilizes (i.e. the slope of the line segment thereafter ≈ 0). The subnivium is then maintained for the duration of time starting from the day of the first occurrence of a zero slope to the day of the last occurrence of a zero slope. The subnivium disintegration date is the day following the last occurrence of a stable slope. appropriate to use for winter 2015/16 by investigating the correlations between our field-measured depth data in winter 2016/17 and the midpoints of the daily snow depth range modeled by NOHRSC for the same time period (table S2). From these data we determined the date of the first snow event that resulted in at least 10 cm of snow accumulation for each site in each season, as well as a snow-off date. We obtained daily snow density values for each site for both winters from NOHRSC data to have an estimate that covered the entire depth of the existing snowpack.

Establishment, maintenance, and disintegration of the subnivium
We aggregated all ground temperature probe data to daily mean ground temperatures for each probe and fit piecewise linear models to the ground temperatures of each probe (Zeileis et al 2003) to determine the number of breakpoints for the entire winter season (Ryan and Porth 2007, Zeileis et al 2003, Zeileis et al 2002. Breakpoints represent the day where the slope of the modeled temperature changed significantly (Ryan and Porth 2007) and indicate shifts in the temperature profile of the subnivium (Zeileis et al 2002). For any number of breakpoints p, there are p + 1 segments in which the regression coefficient is constant; therefore, the segmented model takes the form where i is the day, j is the segment index, and k is the probe identifier (Zeileis et al 2002). Thus, y ik is the subnivium temperature at day i for probe k, x is the date, jk is the regression coefficient and u jk is the intercept for the segment. We trialed multiple piecewise regression models for each probe by varying the number of breakpoints and selected the model which minimized the total residual sum of squares for each segment ( and had the lowest Bayesian Information Criterion (BIC) value (Zeileis et al 2003). We used BIC because it is appropriate for selection when describing all possible models (i.e. breakpoints) (Zeileis et al 2003) and the Akaike information criterion has been shown to overestimate the number of breakpoints (Bai and Perron 2003). Once the number of breakpoints was identified, we obtained the slopes for each segment and used these to define three phenophases of the subnivium: establishment, maintenance, and disintegration. We considered the subnivium established on the day when daily temperature stabilized (i.e. the slope demonstrated stability for at least three days: 0.15 ≤ slope ≤ 0.15), given that the first snowfall with accumulation of at least 10 cm had occurred on or prior to that day (figure 2). Although Pruitt (2005) identified 20 cm as the minimum depth required for subnivium establishment, we found evidence for stabilized temperatures at depths less than this amount when analyzing ground temperatures at hourly time scales. Therefore, we used a threshold depth of 10 cm to include these shorter-term stabilizations and their potential effects on subnivium establishment. Additionally, piecewise regression requires a minimum of three observations, or in this case days, to fit an appropriate line segment.
We defined subnivium maintenance as the duration of time from the day of the first occurrence of a stable slope to the day of the last occurrence of a stable slope (figure 2). Since stable temperatures can occur when snow is not present, we only considered slope segments that fell between the day of the first major snow event and the snow-off date to ensure that we were capturing subnivium conditions. To assess the stability of the subnivium, we calculated the range between maximum and minimum ground temperatures during the maintenance period. Finally, we defined subnivium disintegration as the day that directly followed the last stable slope segment (figure 2). We then created pseudo-Julian dates with November 1st corresponding to 1, and represented each phenophase as a continuous date variable.
The responses of establishment and disintegration date, maintenance period, and subnivium temperature range, were analyzed using generalized linear mixed effects models (GLMM) with Gaussian distributions and temperature string as a random effect (Bolker et al 2009). All models were fit using the lme4 package in R (Bates et al 2015, R Core Team 2015. For maintenance period and subnivium temperature range, we ran separate mixed models for each year of data collection to account for interannual variability.

Predictors of subnivium phenology
Predictors for all responses (establishment and disintegration date, maintenance period, and subnivium temperature range) included latitude, land cover, and year. Covariates of minimum and maximum snow depth (cm), minimum and maximum air temperature ( • C), mean snow density (g cm −3 ), and mean wind speed (m s −1 ) were also included for maintenance periods and subnivium temperature range. Prior to inclusion in the models, we standardized all continuous predictors and checked for multicollinearity. We sequentially removed covariates from inclusion based on their variance inflation factors (VIF) until all covariates had VIFs less than 4 (Zuur et al 2010). For each response, we developed a set of candidate models that included additive combinations of the remaining explanatory variables. We included interactions between continuous covariates that we had hypothesized to be relevant, but did not test all combinations to avoid overparameterization. For this reason, and because we were missing one latitude and land cover combination in Year 1 (upper latitude, deciduous), we did not include an interaction of latitude and land cover in models for maintenance period or subnivium temperature range. We assessed all models with the Akaike information criteria adjusted for small sample size (AIC c ). We considered models that were within 2 ΔAIC c of the top model as competitive; however, models that contained uninformative parameters were excluded from the model selection (Arnold 2010). We averaged predictions for responses that had more than one competitive model based on AIC c weights (w ), but avoided model averaging parameter estimates given recent concerns (Cade 2015).

Timing of subnivium establishment and disintegration
Subnivia established at 79 out of 83 temperature probes in winter 2015/16 (Year 1) and 91 out of 104 probes in winter 2016/17 (Year 2) in all cover types and latitudes. To investigate the effects of latitude, cover type, and interannual variability on subnivium establishment and disintegration dates (represented as days after November 1st), we tested a set of 13 models (tables S3 and S4). For sites and probes where subnivia developed, the best model for establishment date included an interaction between latitude and winter season (w = 0.61) ( Interannual variability produced differences in the timing of establishment in lower and middle latitudes, but the early establishment of subnivia in upper latitudes was consistent across years ( figure 3(a)).
The best model for subnivium disintegration date included an interaction between latitude and season and between latitude and cover type (w = 0.99) (  figure 3(c)). Disintegration dates in middle latitudes were consistent across years and cover types, while in lower latitudes  disintegration varied by about a month between years. The disintegration of subnivia in the northerly sites displayed variability in both years and cover types, with a three-week difference in timing between years and a two-month difference between cover types. Upper-latitude deciduous sites had the latest disintegration dates (157.25 ± 3.72), while subnivia disintegrated two weeks earlier in upper-latitude open sites (140.52 ± 2.51) and two months earlier in upperlatitude conifer sites (101.10 ± 2.87, figure 3(d)).

Duration of subnivium maintenance
Out of 22 candidate models (table S5), the model that best explained maintenance period in Year 1 included latitude, an interaction between minimum air temperature and maximum snow depth, an interaction between maximum snow depth and mean snow density, minimum snow depth, and mean wind speed (w = 0.40) (table 1, figure S1). The top model for Year 2 included cover type, an interaction between minimum air temperature and maximum snow depth, an interaction between maximum snow depth and mean snow density, and mean wind speed (w = 0.85) (  figure 4(a)). Higher minimum temperatures ( • C) produced shorter maintenance periods in both years (Year 1: −2.77 ± 0.82, Year 2: −16.04 ± 3.12, figure 4(b)). Higher snow density (g cm −3 ) also resulted in shorter maintenance periods, although the magnitude of this As maximum snow depth increases, the maintenance period also increases. As minimum air temperature increases, maintenance period decreases. The magnitude of this effect is much stronger in year 2. effect was much stronger in the second year (Year 1: −2.53 ± 2.69, Year 2: −18.86 ± 3.51).
The importance of minimum air temperature, maximum snow depth, and mean snow density became clearer when investigating their interactions. Warmer air temperatures produced shorter maintenance periods, but this effect was lessened by higher snow depths, and was most extreme at colder air temperatures ( figure 5(a)). At minimum air temperatures of −30 • C, a wide range of snow depths were able to maintain subnivia, but the length of this maintenance period fluctuated by depth from around 75 days at low depths of 5-15 cm to around 150 days at high depths of 80-90 cm. The differences in maintenance periods produced by depth were most extreme at low air temperatures, but as temperature increased the variation between maintenance periods produced by differing depths decreased.
Depth and temperature alone did not sufficiently describe the conditions necessary for subnivium maintenance because their effect was mediated by snow density. Although the interaction of temperature and depth showed that longer maintenance periods were possible with higher snow depths and colder air temperatures, high snow densities counteracted this effect ( figure 5(b)). At low snow depths, all density values had similar predicted maintenance periods of around 30 days, but as snow depth increased, maintenance periods rapidly increased for snowpacks with low densities (0.10-0.20 g cm −3 ) ( figure 5(b)). At intermediate densities (0.22-0.26 g cm −3 ), the increase in maintenance periods under increasing snow depths reduced in magnitude, and was almost flat at a density of 0.26 g cm −3 . Finally, high snow densities (0.28-0.30 g cm −3 ) overrode the positive effect of snow depth on subnivium duration ( figure 5(b)).

Subnivium temperature range
Out of 15 possible models (table S6), two candidate models of subnivium temperature range emerged for each year. The top model for Year 1 included cover type, latitude, minimum air temperature and minimum snow depth (w = 0.48), while the second-best model included latitude, an interaction between minimum air temperature and minimum depth, mean wind speed, and mean density (w = 0.20) (table 1, figure  S2). The top model for Year 2 included latitude, an interaction between minimum snow depth and mean snow density, minimum air temperature, and mean wind speed (w = 0.64), while the second-best model included latitude, an interaction between minimum air temperature and minimum snow depth, mean wind speed, and mean density (w = 0.26) (table 1, figure  S2).
Subnivium temperature range represented a measure of stability, with smaller ranges corresponding to more stable subnivia. Consequently, since subnivia in middle latitudes had the smallest temperature ranges in both years, they could be considered the most stable (Year 1: 3.96 ± 0.78, Year 2: 2.98 ± 0.82). Temperatures in lower latitude subnivia were the least stable (Year 1: 7.01 ± 0.64, Year 2: 5.74 ± 0.56). In Year 1, upper latitudes displayed ranges similar to lower latitudes (6.24 ± 0.83), while in Year 2 upper latitude ranges were similar to middle latitudes (3.35 ± 0.64). Thus, while subnivia in upper latitudes were stable and comparable to those in middle latitudes in Year 2, overall both upper and lower latitudes were more susceptible to interannual variability in factors such as snow depth, snow density, and air temperature.
The predictors of subnivium temperature stability were similar between years. Increases in minimum air temperature, mean wind speed, and minimum Colder air temperatures increased temperature variability in the subnivium, but since snow depth insulated the subnivium it minimized the effect of ambient temperature (figure 6(a)). Density further mediated the magnitude of this effect; higher densities (0.26-0.30 g cm −3 ) produced increases in temperature variability even with increasing snow depth ( figure 6(b)).

Discussion
In describing three important phenophases of the subnivium-establishment, maintenance, and disintegration-we found that while latitude, land cover, and interannual variability influenced subnivium phenology, optimal subnivium conditions depended on the relationship between air temperature, snow depth, and snow density. Although the relative impact of these three factors varied across the different subnivium phenophases, air temperature was a universal constraint through its effect on snowpack depth and density, and varied according to phenophase. Above-freezing temperatures promoted rapid subnivium establishment, while below-freezing temperatures reduced snowmelt and supported longer subnivium maintenance periods. Concurrently, cold air temperatures increased temperature variability in the subnivium during the maintenance phase. Thus, the creation of conditions for a stable and longer subnivium period required that air temperatures fell within a certain range.

Establishment and disintegration date
In subnivium establishment, interannual variability was most important in lower and middle latitudes; whereas, in upper latitudes establishment date was similar between years. While the timing of establishment progressed gradually in winter 2015/16 (Year 1), a region-wide snow storm occurred in winter 2016/17 (Year 2), causing spatial synchronization of establishment dates across all latitudes. Additionally, contrary to our expectations, there was no difference in establishment between cover types within or between latitudes. This suggests that latitude, as a proxy for the amount of snowfall an area typically receives, plays a larger role than cover type in subnivium establishment, but that variability in early season snow storms can override the onset of subnivium conditions. Disintegration was not as prone to variability in weather patterns and was more synchronized between years and cover types within the same latitudinal range. As predicted, subnivium disintegration exhibited a gradual pattern of later disintegration at more northerly sites. Middle latitudes exhibited the highest consistency in disintegration dates, which unlike lower and upper latitudes, fell within a similar range when comparing across years and cover types. This suggests that conditions in middle latitudes during the maintenance period leading up to disintegration were such that neither interannual variability nor differences among cover types were sufficient to override the subnivium conditions created by the balance of air temperature, snow depth, and snow density.

Duration of subnivium maintenance
We predicted that subnivium conditions would depend primarily on snow depth and air temperature. While this was partially supported, we found Figure 6. Predictions of subnivium temperature range highlighting the interaction between (a) minimum air temperature and minimum snow depth and (b) minimum snow depth and mean density. As minimum air temperature increases, subnivium temperature range decreases, which means that the subnivium becomes more stable. The most stable subnivia (i.e. temperature range = 0 • C) occur when the minimum air temperature is at or greater than −20 • C. As minimum snow depth increases, the variability in subnivium temperatures decreases at most densities. Higher densities promote an increase in temperature variability at higher snow depths; however, this increase is relatively small in magnitude. that air temperature provided the overarching constraint on subnivium existence, and each subnivium phenophase resulted from a balance between air temperature, snow depth, and snow density. Consequently, describing the subnivium in terms of a threshold of one of these factors alone is insufficient due to their interdependency.
Colder air temperatures were required to promote subnivium maintenance. Warmer air temperatures increase ablation, which reduces depth, increases density, and reduces the overall insulative capacity of snow cover (Ge and Gong 2010). Warmer air temperatures led to shorter subnivium maintenance periods, but this effect was reduced at higher snow depths (figure 5(a)). Since higher snow depths mitigated the effect of warmer air temperatures, they promoted longer maintenance periods, but this outcome was offset by high snow density ( figure 5(b)). The potential for snow density to counteract the effect of snow depth became stronger in the later part of the winter season since snow density increases over time due to changes in the snowpack induced by air temperature, like melting and refreezing (Judson and Doesken 2000). Our results show that at high densities (e.g. >0.25 g cm −3 ) maintenance periods decreased even with increasing snow depth, reducing subnivium maintenance from approximately 30 days at shallow depths (5-25 cm) to 10-15 days at midrange depths (25-45 cm, figure 5(b)). In fact, the range of snow densities at our surveyed sites (0.10-0.30 g cm −3 ) represents a low to moderate range compared to those found in other studies at similar latitudes (Dawson et al 2017, Derksen et al 2014.

Subnivium stability
We measured thermal stability during the maintenance period by quantifying overall temperature variability. We found that warmer air temperatures reduced subnivium temperature variability by reducing the gradient between the bottom layer of snow and the air; however, without sufficient snow depth, the subnivium lacked insulation and was vulnerable to swings in air temperature outside of the optimal range for minimizing thermal variability (figure 6(a)).
By exploring differences in subnivium temperature ranges between latitudes, we found that middle latitude sites (44 • -46 • ) had optimal ranges of air temperature, snow depth, and density relative to upper and lower latitudes. Subnivia in these sites had the narrowest temperature range and this result was consistent in both years, demonstrating that conditions in these areas buffered the subnivia from the effects of interannual variability and differences in land cover. It is, therefore, possible that sites with more stable subnivium temperatures represented the trailing range boundary of optimal conditions for subnivium maintenance; however, it is likely this range boundary will shift northward given the projected northward shift from snow-to rain-dominated regions (Ning and Bradley 2015).

Caveats
Although our study represents a novel exploration of the subnivium and its phenology, we recognize several potential limitations. First, given the number of parameters in our models of maintenance period and subnivium stability, we were unable to model explicit interactions between latitude and land cover due to issues of model overparameterization. For example, collapsing all sites of a particular cover type throughout our entire latitudinal gradient (42 • −46 • +) into a single unit avoided problems with model convergence. Second, we were unable to take direct measures of solar radiation and were not able to explicitly account for this variable, however, many of the abiotic conditions we captured (e.g. ambient temperature, snow density) change in tandem with solar radiation and likely account for these effects.

Biological implications
Understanding the phenology of this critical, yet understudied, seasonal refuge is crucial for predicting how global change will affect the subnivium and overwintering species that depend on the subnivium. Although overwintering species have different physiological tolerances and coping strategies, longer maintenance periods with consistent and predictable temperature stability will provide thermal protection and fitness benefits to a wide diversity of subnivium-dependent species, while also sustaining other ecological processes. For example, the subnivium ecosystem influences patterns and frequencies of soil freeze-thaw cycles (Groffman et al 2001), which in turn regulate soil biogeochemistry (Kreyling et al 2012). Disturbances to the subnivium in either extent, duration, or thermal stability can therefore disrupt these prevailing regimes and result in phenological mismatches (Wheeler et al 2015), leading to root damage in plants (Sanders-DeMott et al 2017) and microbial mortality (Schimel et al 2004). Subnivium disturbances also produce phenological mismatches in wildlife, with decreased survival rates in freezetolerant wood frogs (O'Connor and Rittenhouse 2016), arthropod communities (Bokhorst et al 2012), and small rodents (Korslund and Steen 2006).
Future climate change scenarios predict warmer winter temperatures, which are often accompanied by an increase in precipitation falling as rain rather than snow and an overall increase in air temperature variability (e.g. cold air outbreaks) (Vaughn et al 2013). This means that most of the Great Lakes Region will experience decreased snow depth, increased snow density, and temperature variability, which will reduce the extent, duration, and thermal stability of the subnivium. In this future scenario, areas that receive high amounts of lake-effect snow may be important refugia for subnivium-dependent species. Future research needs to quantify how climate change will affect subnivium phenology to understand the impacts of these changes on the myriad of species that depend on this seasonal habitat.