Understanding the spatiotemporal distribution of snow refugia in the rain-snow transition zone of north-central Idaho

Knowledge of snow cover distribution and disappearance dates over a wide range of scales is imperative for understanding hydrological dynamics and for habitat management of wildlife species that rely on snow cover. Identification of snow refugia, or places with relatively late snow disappearance dates (SDDs) compared to surrounding areas, is especially important as climate change alters snow cover timing and duration. The purpose of this study was to increase understanding of snow refugia in complex terrain spanning the rain-snow transition zone at fine spatial and temporal scales. To accomplish this objective, we used remote cameras to provide relatively high temporal and spatial resolution measurements on snowpack conditions. We built linear models to relate SDDs at the monitoring sites to topoclimatic and canopy cover metrics. One model to quantify SDDs included elevation, aspect, and an interaction between canopy cover and cold-air pooling potential. High-elevation, north-facing sites in cold-air pools (CAPs) had the latest SDDs, but isolated lower-elevation points also exhibited relatively late potential SDDs. Importantly, canopy cover had a much stronger effect on SDDs in CAPs than in non-CAPs, indicating that best practices in forest management for snow refugia could vary across microtopography. A second model that included in situ hydroclimate observations (December–February (DJF) temperature and March 1 snow depth) indicated that March 1 snow depth had little impact on SDD at the coldest winter temperatures, and that DJF temperatures had a stronger effect on SDD at lower snow depths, implying that the relative importance of snowfall and temperature could vary across hydroclimatic contexts in their impact on snow refugia. This new understanding of factors influencing snow refugia can guide forest management actions to increase snow retention and inform management of snow-dependent wildlife species in complex terrain.


Introduction
Snow conditions vary widely across a range of spatial and temporal scales in complex terrain due to interactions between topography, vegetation, and hydrological processes that control snow deposition and ablation. Fine-scale variations in snow cover and properties can cause differences in local hydrology, which has cascading effects on vegetation and wildlife. When topographic, vegetation, and hydrometeorological factors create optimal conditions for the retention of snow, this area can be thought of as a 'snow refugium. ' We define 'snow refugia' as areas with generally later snow cover relative to surrounding areas due to mechanisms including radiation shading, cold air retention, low canopy snow interception, and wind redistribution of snow. Snow refugia may be critical for maintaining biological processes and biodiversity in complex terrain and enhancing the resilience of biota to climate changes and variability (Dobrowski 2011, Curtis et al 2014, McLaughlin et al 2017, Niittynen et al 2018, Pastore et al 2022. For example, evergreen forests may be expected to persist in snow refugia because the cold temperatures and more frequent frost favor evergreen species over deciduous species (Pastore et al 2022). Similarly, snow refugia may support aspen stands in mixed sagebrush-aspen environments (Soderquist et al 2018, Marshall et al 2019. Snow refugia will also maintain the snow cover, depth, and extent that snow-dependent wildlife species such as bumblebees (Bombus sp.), wolverines (Gulo gulo), and snowshoe hares (Lepus americanus) need for thermal cover and phenological matching (Kudo and Ida 2013, Curtis et al 2014, Mills et al 2018. Therefore, there is a need to understand and predict where snow refugia occur to promote the conservation or formation of these areas. Terrain and hydrometeorological processes affect snow conditions in ways that vary across scales (Clark et al 2011). Orographic effects and rain shadows control precipitation variations in complex terrain (Leung et al 2003). Increasing elevation is typically associated with lower air temperature, although this varies with latitude, pressure, season, and hydrometeorological conditions (Stone and Carlson 1979). Cold-air pooling can create temperature inversions in low-elevation areas and topographic concavities (Lundquist et al 2008, Daly et al 2010, Curtis et al 2014. Additionally, canopy cover and air temperature affect the efficiency of canopy snow interception (Roth and Nolin 2019), canopy unloading, and turbulent (sensible and latent) energy transport (Strasser et al 2008). While snow may unload from the canopy, some sublimates directly from the canopy or melts and passes though the snowpack, causing forests to have lower snow accumulation than nearby unforested areas (Dickerson-Lange et al 2017, Jeníček et al 2020. Tree canopies shade direct shortwave radiation (0.28-3.5 µm) from reaching the snow surface while simultaneously emitting longwave radiation (3.5-100 µm), causing snow around tree trunks to melt more quickly (Dickerson-Lange et al 2017). Debris on the snow surface decreases snow albedo and subsequently increases net snow cover radiation (Hardy et al 2004, Dickerson-Lange et al 2017. Slopes and concavities also reduce net shortwave radiation at certain times of the day due to terrain shading (Curtis et al 2014). Finally, snow redistribution by wind and gravity changes the spatial distribution and depth of snow (Bernhardt et al 2012, Marshall et al 2019. These processes can work with or against each other to produce snow conditions that are highly variable even at very fine scales (Jost et al 2007).
Current snow data products are limited by their spatial or temporal resolution, cost, and effort, making them insufficient for identifying snow refugia. Remote cameras are useful for monitoring snow cover in complex terrain (Sirén et al 2018). Cameras collect data at a finer spatial resolution (i.e. the viewshed of the camera) than is captured by airborne or satellite remote sensing methods but can be deployed relatively inexpensively with a denser spatial distribution than in situ observation methods such as snow surveys or snow telemetry sites. Cameras can be programmed to take timelapse images, allowing for regular monitoring of snow within the camera viewshed. A temperature measurement is also recorded with each image, which is useful for identification of local freezethaw cycles and cold-air pools (CAPs). Cameras can be programmed to take both motion-triggered and timelapse images, allowing wildlife researchers to observe wildlife detections and snow conditions with one instrument. Snow depth can be determined when a snow stake is installed (Dickerson-Lange et al 2017, Sirén et al 2018 or even in the absence of one (Strickfaden et al submitted), and are extremely reliable for determining snow presence. Thus, remote cameras are a cost-effective, reliable tool for estimating snow conditions and identifying snow refugia, particularly in complex forested terrain.
We sought to improve our understanding of how biophysical conditions control the locations and dynamics of snow refugia at fine scales (∼30 m 2 ) in complex forested terrain spanning the rainsnow transition zone. Our specific objectives were to quantify the relative timing and sensitivity of snow disappearance dates (SDDs) across the study area and to understand how combinations of topography, vegetation, and hydrometeorological conditions interact to produce snow refugia. We deployed a stratified network of remote cameras in a complex forested terrain to monitor snow conditions at fine scales. Empirical models based on these data were developed to (a) determine the general locations of snow refugia based on static (terrain and vegetation) variables, and (b) assess the sensitivity of actual SDDs to differences in dynamic (hydrometeorological) variables. In these analyses, areas with later-than-average SDDs are considered to function as snow refugia. The resultant models provide insight for forest management to conserve and create lateseason snow cover, map potential locations of snow refugia, and identify priority areas of conservation and management of snow-dependent wildlife species.

Study area
The study area for this project is Moscow Mountain in Latah County, ID (figure 1). The study area has a wide variety of topography, canopy cover, and hydroclimatic conditions that result in large variations in snowpack distribution and properties. Moscow Mountain spans approximately 800-1500 m a.s.l. Forest cover is primarily coniferous and is comprised of ponderosa pine (Pinus ponderosa) and Douglas-fir (Pseudotsuga menziesii) forests at low elevations and western red cedar (Thuja plicata) and fir (Abies sp.) forests at high elevations (Falkowski et al 2009). Winter precipitation is generally rain-dominated at the lower elevations and snowdominated at the higher elevations, which is representative of mid-elevation landscapes in the interior Pacific Northwest spanning the rain-snow transition zone. The Moscow Mountain Snow Telemetry station (SNOTEL 989, 1430 m a.s.l.) receives an average of approximately 1060 mm of precipitation annually. The mean peak snow-water equivalent (SWE) for Moscow Mountain from water years 2001-2021 was 529 mm, and the peak SWE in 2021 was 594 mm. The mean December-February (DJF) temperature on Moscow Mountain from water years 2001-2021 was −0.55 • C, and the mean DJF temperature in 2021 was −2.0 • C (USDA 2022).

Snow data collection
We deployed 134 camera stations from October 2020 to May 2021 to collect snow presence and depth data. Study sites were stratified by elevation, aspect, and canopy cover and quasi-randomly selected to permit reasonable access via existing trail and forest road networks. Elevation was classified from <925 m, 925-1050 m, 1050-1175 m, 1175-1300 m, and >1300 m and aspect was classified as N, S, E, or W. Canopy cover was classified as 'sparse' for <35%, 'moderate' for 35%-75%, and 'dense' for >75% canopy cover based on handheld densiometer measurements. Reconyx® Hyperfire or Hyperfire II cameras were installed at each site and programmed to record one image each hour of the day, including nighttime hours. Cameras were retrieved following complete snow cover ablation from April to May 2021.
On the deployment date, we recorded the latitude, longitude, and elevation of the camera. We mounted and locked cameras onto trees at a height of 2-3 m to prevent snow from blocking the cameras. During camera deployment, we took reference images for superimposing a 'virtual' snow stake (VSS) onto images to measure snow depth (Strickfaden et al submitted). We used VSSs instead of physical snow stakes to reduce material and personnel costs. We placed a reference snow stake with 2 and 10 cm gradations at 5, 10, and 15 m within the viewshed of the camera and allowed the camera to take motiontriggered images. We took an additional set of reference images on camera retrieval to account for potential changes in the camera's viewshed during its deployment. VSSs were then superimposed onto images using the 'edger_multi' function in the edger package (Strickfaden et al submitted) and used to estimate snow depths. We maintained cameras every few weeks to change batteries, ensure proper function, and remove obstructing vegetation.

Air temperature
We placed an external LogTag® TRIX 8 temperature recorder (hereafter 'LogTag') at each camera station. We programmed the LogTags to record temperatures at 45 min intervals. The LogTags were housed inside plastic protective coverings and radiation shields. Each radiation shield was comprised of 15 cm of PVC pipe covered in aluminum foil tape to reflect shortwave and longwave radiation. We drilled holes into the radiation shields to allow for increased air flow (Terando et al 2017). We suspended the LogTag in its housing at the same height as the camera in a sheltered area to minimize the potential impacts of direct shortwave radiation.

Hemispherical photography
To quantify total vegetation cover and estimate incoming shortwave radiation, we took hemispherical photographs at each camera site the following summer. We used a Canon ™ EOS 70D SLR camera with a Sigma 8 mm circular fisheye lens. We took hemispherical photographs on days with little to no wind or precipitation and under diffuse light (e.g. at very low sun angles or during overcast conditions) to maximize the contrast between sky and vegetation.
At each camera site, we placed the DSLR camera with attached lens on a tripod 5 m into the camera viewshed and oriented it vertically. We took photographs at several exposures and manually selected the photograph with the best exposure for each camera site. We analyzed the best image from each location using Hemisfer software (Thimonier et al 2010). We manually selected threshold values for distinguishing sky pixels from vegetation pixels for each image. We calculated potential and effective shortwave radiation assuming 50% cloud cover given that winter conditions are a mix of sunny and overcast. Hemisfer also outputs a count of sky and vegetation pixels in the image which we used to calculate a percent cover value. Because the 180 • image captures some understory vegetation in addition to the overstory vegetation, this percent cover value from the hemispherical photographs (hereafter 'total cover') is different from the canopy cover value obtained using the densiometer (hereafter 'canopy cover').

Cold air pooling model
To identify potential CAPs, we developed a simple procedure following methods used by Ashcroft and Gollan (2012). We used a 1 m resolution DEM of Moscow Mountain (Bright et al 2019) to calculate flow accumulation and relative elevation of each pixel. First, we calculated flow accumulation using the Hydrology tools in ArcGIS (ESRI 2011) as a proxy for dense air flow accumulation and took the natural log of this value. Second, we calculated relative elevation by determining the elevation of the lowestelevation pixel within a 25, 50, and 100 m radius from each pixel. We subtracted the elevation of a particular pixel from the elevation of the lowest pixel in the focal radius to determine that pixel's relative elevation. We rescaled both the log flow accumulation and relative elevation rasters from 0 to 1. Finally, we subtracted the log flow accumulation raster from each relative elevation raster calculated at different focal radii. Negative differences indicate pixels had both high log flow accumulation and low relative elevation for a given focal radius, which together contribute to potential accumulation of cold air. Therefore, if the difference between the two values was negative, we assigned the pixel a value of 1 (CAP pixel), and if it was positive, we assigned it a value of 0 (non-CAP pixel). We merged rasters for the three focal radii to determine if a pixel was identified as a potential CAP in any of the three neighborhoods.
We assessed the ability of this approach to identify true CAP locations by analyzing the minimum and maximum temperatures from the LogTags on certain dates. We selected dates by visually examining the incoming shortwave radiation flux, net longwave flux, and wind speed curves on the SNOTEL website. We found six dates with smooth shortwave flux curves, negative longwave flux, and low wind speeds, which indicate conditions that would create cold-air pooling (Pastore et al 2022). We built suites of linear models (LMs) to estimate minimum and maximum temperatures on the six selected dates. After accounting for date, elevation, aspect, and canopy cover, LogTags in CAP pixels had lower average minimum temperatures by 0.46 • C and lower average maximum temperatures by 0.94 • C than LogTags in non-CAP pixels. Camera stations were not always placed in the bottoms of CAPs which have the greatest temperature difference (Clements et al 2003), which may explain why we found smaller temperature differences than other CAP studies. Despite this, we found sufficient evidence that this layer adequately identified CAP pixels because temperatures of LogTags in CAPs were suppressed. Thus, this layer was used to indicate whether a camera site was in a CAP.

Image processing
We manually recorded snow presence and snow depth in each image. We measured snow depth using any VSSs present in images. If snow depth >150 cm (the height of the VSS), it was recorded as '150+' . The SDD for a camera was defined as the first day there was no snow cover in the immediate viewshed of the camera. Though some cameras had larger viewsheds than 15 m, we only assessed snow presence within the distance of the farthest VSS, which was ⩽15 m from the camera. Snow events after seasonal snow cover at the camera site disappeared were not counted towards the SDD following Dickerson-Lange et al (2017). We processed images using Timelapse2 software (Greenberg 2020).

Data analysis
We used data from the 2020-2021 season to build two suites of LMs with SDD in Julian days as the response. Our objective was to build two models: a static model which could be used to determine the general locations of snow refugia based on terrain and vegetation variables collected by satellites or estimated with models, and a dynamic model which could be used to estimate the actual SDDs by incorporating interannually fluctuating variables collected at the camera stations (i.e. snow depth and temperature). For the static models, we included elevation (classification or value in meters); slope (degrees); aspect (cardinal direction or cosine of radians); canopy cover (classification or percent); CAP potential (indicator); and mean effective shortwave radiation (W m −2 ). For the dynamic models, we also included total cover from the hemispherical photographs, mean DJF temperature from the LogTag ( • C); and March 1st snow depth (cm), which approximately coincides with peak SWE at this location. March 1st snow depth may have been biased low if snow depth exceeded 150 cm at a site (17% of sites). Combinations of the predictors and interactions between them were tested, and the best model for each model set was selected using Akaike's information criterion for small sample sizes (AICc; Hurvich and Tsai 1989) and the MuMIn package (Barton 2015) in Program R (R Core Team 2021).

Results
We collected 562 000 timelapse images at 134 camera sites in 2020-2021. SDDs ranged from 18 March to 26 May. Although snow depth in the images could only be measured up to 150 cm, depths of up to 203 cm were recorded during snow density and hardness sampling (Strickfaden et al submitted).
The static model with the most support contained elevation, cardinal aspect, and an interaction between percent canopy cover and CAP potential (figure 2, table 1). The coefficient of determination (R 2 ) of this model was 0.71. SDD was predicted to increase by one day for each 20 m increase in elevation (95% CI: [17.1, 22.6]). Sites on west-facing aspects had the earliest SDDs, with south-facing aspects 2.5 days later (95% CI: [−1.3, 6.3]), east-facing aspects 6.9 days later (95% CI: [3.0, 10.9]), and north-facing aspects 11.4 days later (95% CI: [7.5, 15.3]) than west-facing aspects. For sites not in CAP pixels, SDD was predicted to increase by one day for each 10% increase in canopy cover (95% CI: [6.6, 18.8]), while for sites in CAP pixels, only a 4% increase in canopy cover was needed (95% CI: [4.0, 72.9]; table 2).
The dynamic model to predict SDDs with the lowest AIC included cardinal aspect, percent total cover, and an interaction between mean DJF temperature and March 1st snow depth ( figure 3, table 3). The R 2 value of this model was 0.79. Sites on westfacing aspects had the earliest SDDs, with southfacing aspects 1.5 days later (95% CI: [−1.7, 4.6]),      (table 4). SDD was predicted to increase by one day for each 0.07 • C decrease in mean DJF temperature (95% CI: [0.05, 0.10]), each 5 cm increase in March 1st snow depth (95% CI: [3.7, 8.5]), and each 10% increase in total cover (95% CI: [5.5, 15.6]). The interaction indicated that DJF temperatures had a stronger relationship with SDDs at locations with shallow snowpack than in areas with deeper snowpack. In addition, SDDs were relatively insensitive to snow depths in colder areas (figure 3). A considerable advantage of the static model is that the data required by this model are readily available as spatially distributed datasets, which allowed for the extension of this model beyond the camera sites to map SDDs across our entire study area (figure 5). Mapping of SDDs identified numerous isolated areas with later SDDs than their surroundings that likely function as snow refugia.

Discussion
We built two models to predict the locations of snow refugia at fine spatial scales (∼30 m 2 ). Each model elucidated a different effect topoclimates had on snow retention in our study area and year. The static model, which did not rely on snow or temperature data, suggested that SDDs in CAPs were more sensitive to canopy cover than SDDs in non-CAPs. The dynamic model, which relied on snow and temperature data but is more transferable to other years, suggested that SDDs were less sensitive to snow depth at locations with colder DJF temperatures relative to warm DJF temperatures. At very cold DJF temperatures, SDD was expected to be almost the same regardless of March 1st snow depth, similar to other observational (Hubbart et al 2015) and modeling studies (Sun et al 2022). Evaluation of the specific mechanisms that produced these two key findings is beyond the scope of this investigation but is likely due to differing SWE and importance of net radiation and turbulent Figure 5. Maps of (a) relative snow disappearance date (SDD) predicted with the static model, (b) canopy cover, (c) cold-air pooling potential, and (d) aspect across Moscow Mountain in Latah County, ID. In panel (a), the color of the camera sites (triangles) corresponds with the observed snow disappearance date at that site using the same color scale as the predicted SDDs. Pixels are 30 m 2 resolution except for panel (c) in which pixels are 1 m 2 resolution. Each 0.1 increase in relative SDD represents an approximate seven-day increase in actual SDD. Panel (e) shows the virtual snow stakes used for snow depth measurements and shows how SDD was determined from camera data. All snow was melted below the yellow dashed line (i.e. within 10 m of the camera). Note that the snow further than 10 m was not counted towards the SDD. energy transfer in these environments. These mechanisms could be elucidated with a combination of detailed observational and physically based modeling investigations.
Both models tended to underestimate SDDs at the camera sites with the latest actual SDDs by approximately 4.5 days, meaning these sites are serving as later-season snow refugia than are being predicted by the models. Thus, these are conservative estimates of snow refugia locations in our study area and year. We examined the locations of the most extreme outliers to assess potential causes of modeling errors but found no clear and consistent explanations.
Specific SDD timing will vary between years based on specific and intra-annual variations in the hydrometeorological conditions that affect snow deposition and ablation processes. Normalizing the predictions of the static model from 0 to 1 (i.e. interpreting the predictions as relative earliest to latest SDDs) may make the predictions more robust to annual temperature and precipitation differences. This is reasonable because the dominant physical processes underlying general differences in SDD with elevation, aspect, and cold-air pooling are unlikely to change dramatically due to climate change (i.e. high elevations, north aspects, and CAPs are still expected to have the latest SDDs regardless of climate change) but are likely to vary as the relative importance of different snow energy balance components are altered by climate variations. Additionally, we did not have enough data from prior years to determine the sensitivity of either model to interannual variability in temperature and precipitation. Further study of how hydrometeorological variations affect the relative timing of SDDs is warranted to further advance understanding of snow refugia dynamics.
Associations between canopy cover and SDD may change as certain temperature thresholds are exceeded because of climate change (Lundquist et al 2013). In colder climates, seasonal melt is dominated by shortwave radiation that is strongly reduced by forest canopies, whereas in warmer climates where midwinter melt is more common and shortwave radiation comprises a smaller portion of the net snow cover energy, canopy cover enhances net longwave radiation and midwinter melt rates (Lundquist et al 2013). Additionally, warmer air temperatures can cause snow to adhere to the forest canopy more efficiently. Assuming snow does not immediately slough off the canopy, this reduces accumulation under the canopy and allows for more sublimation from the canopy (Lundquist et al 2013, Roth andNolin 2019). Thus, as mean temperatures increase, denser vegetation cover may be detrimental to snow retention. Both best-fit models predicted that denser vegetation cover promoted snow retention in our study site. However, the mean DJF temperature at the Moscow Mountain SNOTEL station was −2.0 • C in 2020-2021 (USDA 2022), which is approaching the threshold of −1 • C at which canopy cover enhances snowmelt by increasing longwave radiation (Lundquist et al 2013). As climate change continues to increase mean temperatures in the rain-snow transition zone, SDDs may need to be re-evaluated to determine if patterns of SDD with canopy cover have shifted or if persistent locally decoupled CAPs (Daly et al 2010) reduce the sensitivity of SDD in these areas. Our results highlight the importance of fine-scale topographic variation in mediating temperature increases by creating CAPs where local temperatures may be colder than surrounding average temperatures.
Importantly, we did not quantify the effects of small canopy gaps (1-3 tree heights in diameter) or forest edges on snow retention. Small canopy gaps have a unique canopy structure that creates a distinct radiative environment compared to open or contiguous forest (Lawler and Link 2011). Edge environments also have unique snow processes and dynamics, and the effects of edges can extend up to 44 m into the forest (Webb et al 2020). Further investigation is needed to understand the role of small canopy gaps and edges in promoting snow refugia. We also did not explicitly account for wind drifting of snow in these models. However, prevailing winds may explain why the west-facing aspects exhibited the earliest SDDs, because winds are typically from the west during major storm events with observed preferential deposition on leeward aspects.

Conclusions
Snow data are often only available at relatively broad spatial resolutions (e.g. ⩾500 m 2 ) or at very fine but spatially sparse scales in canopy gaps, leading to potentially spurious estimates of local snow conditions in complex forested terrain and a poor understanding of the availability of snow refugia for snowdependent biota. We were able to quantify SDDs at fine spatial and temporal scales with certainty because we have literal snapshots of snow conditions in diverse and complex terrain from the remote cameras. The camera data, along with site data which are commonly recorded during camera deployment, allowed us to map SDDs across our study area with relatively little field effort and determine likely locations of snow refugia which were not captured by cameras.
These models elucidated important effects of topoclimates on snow retention in complex terrain, namely that SDDs in CAPs are sensitive to canopy cover while SDDs in very cold locations are less sensitive to snow depth. These insights are useful for management of forests and wildlife. For example, CAPs on high-and mid-elevation north-facing slopes can be managed for the optimal amount of forest cover dependent on the dominant radiative processes (i.e. denser cover in colder environments and sparser cover in warmer environments). Management for conditions promoting snow retention will ensure these vital refugia are available as climate change advances SDDs. Findings from these models are most applicable to complex forested terrain spanning the rain-snow transition zone of the interior Pacific Northwest. These findings may not be appropriate in areas with very simple terrain, considerable wind or gravitational snow redistribution, or ephemeral snow cover. Future research is needed to evaluate the transferability of these models across environments; develop and advance statistical and processbased models to adequately predict snow refugia; and evaluate the inter-annual persistence and climate sensitivity of the refugia identified here.

Data availability statement
The data that support the findings of this study can be found at https://doi.org/10.21429/bma6-xn17 as well as on the ScienceBase website at https://www.sciencebase.gov/catalog/item/63f7ab26 d34e4f7eda4565e7.