Improving Peary Caribou Presence Predictions in MaxEnt Using Spatialized Snow Simulations

The Arctic has warmed at twice the global average over recent decades, which has led to a reduction in the spatial extent and mass balance of snow. The increase in occurrence of winter extreme events such as rain-on-snow, blizzards, and heat waves has a significant impact on snow thickness and density. Dense snowpack conditions can decrease or completely prevent foraging by Peary caribou (Rangifer tarandus pearyi) by creating “locked pastures,” a situation where forage is present but not accessible under snow or ice. Prolonged and severe weather events have been linked to poor body condition, malnutrition, high adult mortality, calf losses, and major population die-offs in Peary caribou. Previous work has established the link between declines in Peary caribou numbers in the Canadian Arctic Archipelago and snow conditions, however these efforts have been limited by the quality and resolution of data describing snow physical properties in the Arctic. Here, we 1) investigate whether a snow model adapted for the Antarctic (SNOWPACK) can produce snow simulations relevant to Canadian High Arctic conditions, and 2) test snow model outputs to determine their utility in predicting Peary caribou occurrence with MaxEnt modelling software. We model Peary caribou occurrence across three seasons: July – October (summer forage and rut), November – March (fall movement and winter forage), and April – June (spring movement and calving). Results of snow simulations using the Antarctic SNOWPACK model demonstrated that both top and bottom density values were greatly improved when compared to simulations using the original version developed for alpine conditions. Results were also more consistent with field measurements using the Antarctic model, though it underestimated the top layer density compared to on-site measurements. Modelled outputs including snow depth and CT350 (cumulative thickness of snow layers surpassing the critical density value of 350 kg·m-3; a density threshold relevant to caribou) proved to be important predictors of Peary caribou space use in each of the top seasonal models along with vegetation and elevation. All seasonal models were robust in that they were able to predict reasonably well the occurrence of Peary caribou outside the period used to develop the models. This work highlights the need for continued monitoring of snow conditions with climate change to understand potential impacts to the species’ distribution.


INTRODUCTION
The Arctic has warmed at twice the global average over recent decades because of a number of processes and feedbacks (Serreze and Barry, 2011;Pithan and Mauritsen, 2014;Davy et al., 2018). A direct consequence is the increased occurrence of winter extreme events, such as rain-on-snow (ROS) (Liston and Hiemstra, 2011; and blizzards (Day et al., 2018;Dolant et al., 2018), which have significant impacts on how the cryosphere responds to climate change (Fountain et al., 2012). More specifically, the spatial extent of snow (Derksen et al., 2016), sea ice (Stroeve et al., 2014;Serreze and Stoeve, 2015), and ice caps (Gardner et al., 2013;Papasodoro et al., 2015) has decreased. Negative snow anomalies could impact many Arctic ecosystem processes as seasonal snow cover plays an important role in the surface energy balance (Lund et al., 2017) through its high albedo and low thermal conductivity (Riche and Schneebeli, 2013;Domine et al., 2018). It is also a key hydrological variable, acting as an important freshwater reservoir (Barnett et al., 2005), and has major implications for permafrost thermal regimes (Romanovsky et al., 2010).
Snow cover and physical properties also influence the survival and reproductive success of wildlife in Arctic ecosystems. Unfavourable conditions created by the densification of the snowpack though the formation of ice crusts or wind slabs may affect the survival of ungulates (e.g., caribou) by blocking their access to food (Putkonen and Roe, 2003;Hansen et al., 2019). Both Vikhamar-Schuler et al. (2013) and Ouellet et al. (2017) found that a threshold of 350 kg·m -3 represents a critical snow density value associated with population declines of both Svalbard reindeer (Rangifer tarandus platyrhynchus) and Peary caribou (Rangifer tarandus pearyi), respectively. Dolant et al. (2018) further demonstrated that sustained heavy winds could contribute to snow densities exceeding this threshold, which resulted in a die-off event for barren-ground caribou (Rangifer tarandus groenlandicus) on Prince Charles Island, Nunavut, in 2015 -16. With dense fur, short muzzles, and wide hooves, Peary caribou are particularly well adapted to harsh Arctic conditions (COSEWIC, 2004). Unlike adjacent caribou designatable units, which assemble in herds, Peary caribou live in small groups of five to ten individuals and move within and between islands in search of suitable forage (Miller et al., 1977). They select the most nutritious parts of a great variety of plants according to their seasonal availability, such as flower heads in the summer for their high energy content (COSEWIC, 2004) or legumes in the winter for their high digestibility and protein content (Larter et al., 2002). Overall, their habitat consists of Arctic tundra and a highly variable topography from flat surfaces in the southwest area of the Canadian Arctic Archipelago (e.g., Victoria Island) to mountainous terrain in the northeast (e.g., Ellesmere Island). Their distribution occupies a territory larger than 1.9 million km 2 across northern Nunavut and the Northwest Territories ( Fig. 1) (COSEWIC, 2015).
Peary caribou were listed as endangered under the federal Species at Risk Act (SARA) in Canada in 2011. Winter extreme events are one of the main threats associated with their recent decline; understanding how changing snow conditions influence their distribution and survival will be important to the species' recovery (COSEWIC, 2015). While the literature clearly demonstrates that snow conditions influence caribou survival (Hansen et al., 2019), the lack of high-quality data characterizing snow conditions across the Arctic Archipelago has been identified as a critical knowledge gap towards better understanding space use (Johnson et al., 2016;Boelman et al., 2019). Both Johnson et al. (2016) and Jenkins et al. (2020) developed predictive species distribution models based on environmental predictor variables including snow. Jenkins et al. (2020) used monthly snow depth means at 24 km resolution (from the Canadian Meteorological Centre daily snow depth analysis database) in late winter models across the eastern Canadian range of Peary caribou but found other predictor variables were more strongly linked to Peary caribou presence. Johnson et al. (2016) used snow data from climate models at a 25 km resolution (Canadian Centre for Climate Modelling and Analysis) and found both snow depth and surface snowmelt contributed to their top seasonal models developed for the entire Canadian Arctic Archipelago. There remain, however, large uncertainties associated with the ability of both these climate models to adequately represent snow conditions in the Arctic and to identify the key snow properties known to influence FIG. 1. Peary caribou distribution defined using a standard convex polygon methodology encompassing aerial counts, telemetry, and Aboriginal traditional knowledge (1970 -2015) (from Johnson et al., 2016, permission granted). The red circle shows the Bathurst Island complex study site. caribou foraging patterns and fitness (e.g., snow density, stratification, and ice layers).
Despite recent progress in the identification of ROS events from space-based observation in the Arctic (Grenfell and Putkonen, 2008;Dolant et al., 2016), little is known about the spatial and temporal patterns of ROS or the occurrence of other winter extreme events such as blizzards and the cumulative impact of these extreme events on surface energy balance and snow conditions. The SNOWPACK model (Lehning et al., 2002a, b), originally developed for use in alpine areas, provides one tool to address this gap by integrating snow and atmospheric data to characterize snow conditions. Ouellet et al. (2017) created a platform to spatialize SNOWPACK simulations using meteorological forcing data to enable the production of snow layers at spatial and temporal scales relevant to wildlife research. The platform includes the Antarctic version of the SNOWPACK model (Groot Zwaaftink et al., 2013) adapted to recognize differences in snow deposition and accumulation processes in polar regions with potential for application to the Canadian Arctic environment.
The objectives of this study are to 1) test our ability to characterize snow conditions in a study area of the Canadian High Arctic using a snow model developed for the Antarctic, and 2) evaluate the utility of the modelled snow data in combination with other environmental predictor variables (including vegetation, elevation, and human disturbance as in Johnson et al., 2016) at predicting Peary caribou presence using MaxEnt modelling. Our approach is novel in that it includes one of the first attempts to include snow physical properties generated from a snow model to directly evaluate key climate threats to foraging accessibility (e.g., densification from ROS, ice crusts, wind slabs) to predict Peary caribou space use.

Study Site
The study focused on the Bathurst Island complex (BIC), located in the Queen Elizabeth Islands in Nunavut, Canada (Fig. 2). This region was selected because it is one of the few places in the Arctic where Peary caribou have been fitted with GPS collars, and where animal locations are available on a weekly basis (or more) throughout the year. The study area approximates 19,000 km 2 that includes Bathurst Island and five large islands (Cameron, Vanier, Alexander, Massey, and Helena) and incorporates Qausuittuq National Park and Polar Bear Pass National Wildlife Area (Miller, 1995). There are no permanent settlements on the BIC, but Inuit from the community of Resolute Bay hunt on the island complex.
The main plant functional types present are cushion forbs, dry graminoids, bryophytes, prostrate dwarf shrubs, and cryptogams (Walker, 2000;Gould et al., 2003). The nearest meteorological station, which is in Resolute on Cornwallis Island (74.72 N,94.97 W), reports average temperatures of −32.0°C in January and +4.5°C in July. Annual precipitation rarely cumulates above 160 mm, making this one of the driest regions on Earth (Environment and Climate Change Canada, 2019).

Telemetry Data
The telemetry data (GPS collars) used for this study were acquired from the Government of Nunavut (Jenkins and Lecomte, 2012). Argos satellite collars were put on seven female Peary caribou in different areas of the BIC. Location data from the collars spanned from April 2003 to May 2006. GPS collars were programmed to collect one location every two days during the months of April through June and one location every five days for the rest of the year.
We used two of the three years of available telemetry data for our analysis. We developed our MaxEnt models using data from July 2003 to June 2004 (Fig. 2, Table 1). While we recognize that using only one year of data may limit the application of our models to Peary caribou management, we were focused on enhancing our ability to model snow in Arctic environments and providing a cursory evaluation of the potential utility of the model to wildlife management. We restricted our model development to one year of data to increase temporal matching between snow conditions and caribou presence data. The inability to capture spatial and temporal changes in resources by averaging conditions over multiple years can result in poor model fit and weak inferences about factors governing animal distributions (see Boyce et al., 2002). Nevertheless, we used data from 2005 to validate the robustness or transferability of our models to other times. We examined the relationship of area-adjusted frequencies of Peary caribou telemetry data across 10 equal interval bins of predicted presence probabilities (use versus expected) following the external validation procedure described by Boyce et al. (2002). A good model should have a strong positive correlation with more animal locations in higher ranked bins.

SNOWPACK Snow Model
A spatialized platform of the SNOWPACK (Lehning et al., 2002a, b) model developed by Ouellet et al. (2017) provides an opportunity to overcome the lack of in-situ snow data in the Arctic. We extracted snow depth and the cumulative thickness of snow layers surpassing the critical density value of 350 kg·m -3 (CT350; Ouellet et al., 2017). CT350 was calculated by adding the thickness of the layers of snow in which density is greater than the threshold of 350 kg·m -3 at three-hour increments over the entire season. The number of layers depends on precipitation events and the predicted snow depth, with the transition between solid and liquid precipitation set at +1.2ºC. Cumulating layer thickness of density values higher than 350 kg•m -3 provides an indication of the thickness and persistence of potentially unfavourable caribou foraging conditions. A hard crust that  is several centimetres thick will not increase the average density of the whole snowpack dramatically but can be problematic for caribou (Tyler, 2010;Vikhamar-Schuler et al., 2013;Langlois et al., 2017;Ouellet et al., 2017). As such, cumulating the layers of density provides additional information about foraging conditions not captured by snow density. We extracted mean values of daily air (2 m) and surface temperatures (˚C), relative humidity (%), wind speed (m•s -1 ), incoming/reflected shortwave and incoming longwave radiation (W•m -2 ), and cumulative precipitation over the three-hour period (kg•m -2 or mm) from NARR (North American Regional Reanalysis product from the National Centers for Environmental Prediction Environmental Modeling Center) for input into the SNOWPACK simulations. While NARR data are available at a 32 km resolution, snow simulations were rescaled to 1 × 1 km using land cover, topography, and soil albedo information to better match the spatial resolution of the Peary caribou location data. Beaudoin-Galaise (2016) suggested that local information on the soil configuration in SNOWPACK (i.e., topography, albedo, land cover) introduces variability in snow simulations despite the pixels sharing the same meteorological forcing information. Their results suggested that slopes greater than 3º will affect the CT350 values, depending on aspect with respect to main wind direction. We visually inspected the SNOWPACK outputs to evaluate whether the use of 1 × 1 km land cover, topography, and albedo information introduced finer-scale variation in snow conditions compared to a 32 × 32 km resolution of the meteorological data.
Several studies have highlighted problems with using SNOWPACK in open tundra environments. For example, Langlois et al. (2012) discussed problems in microstructure, such as a large overestimation of snow grain size in strong kinetic metamorphism conditions. Domine et al. (2019) further investigated how energy transfer and thermal conductivity created problems in the snow profile, reporting a tendency to overestimate density values near the bottom of the snow cover. We evaluated whether the use of the Groot Zwaaftink et al. (2013) Antarctic version of SNOWPACK (hereafter referred to as the Antarctic version) could overcome some of the effects of potential problems in microstructure (Gouttevin et al., 2018) on the snow density profile (Fig. 3). This version accounts for the harsh Antarctic wind conditions, which contribute to hardening of the top layer of the snow to densities as high as 600 kg·m -3 . This version also only cumulates snow once a given wind speed threshold (4 m·s -1 ) is surpassed to account for very important wind redistribution processes. In the model, each precipitation event is considered as a new layer with initial volume fractions of air, ice, and water. A densification scheme takes place as a function of time along with wind speed (see below) influencing snow accumulation in the Antarctic version of the model. We conducted two simulations using 1) the original SNOWPACK version, and 2) using the Antarctic version. Both were forced with NARR data and rescaled at 1 km using local topography, albedo, and soil type information.
To evaluate the performance of SNOWPACK over open tundra areas and justify the use of the Groot Zwaaftink et al.
(2013) approach, we extracted wind speed values from our own meteorological station near Cambridge Bay (69.06 N, −104.76 W) (https://grimp.ca/data/cambridge-bay-1/) for the 2014 -15 winter season. No snow data are available for this type of evaluation on BIC, however the snow type ("tundra snow") is the same as in Cambridge Bay following the classification from Sturm et al. (1995) based on snow geophysical properties such as thickness and density and the processes governing these properties such as wind speed, precipitation, and temperature. Results suggest that the 4 m·s -1 threshold is reached for 63.3% of the period, which justifies the use of the Antarctic version in the region.
We compared snow simulations to in-situ data collected at our Cambridge Bay site where our group has collected snow information since 2015 (Vargel et al., 2020) and where snow conditions are expected to be similar to BIC. The data were collected along six transects established in the Greiner watershed spanning over 600 km 2 encompassing various ecotypes, as described in Ponomarenko et al. (2019). Each transect consisted of snow depth measurements taken every 100 m with a full snow pit dug every 1 km. Around each snow pit, a series of snow depth measurements were recorded in circle increments in a 30 -40 m radius using a magnaprobe (Sturm and Holmgren, 2018), allowing a collection of over 1000 depth measurements at each site. Further details on snow measurements can be found in Langlois et al. (2012Langlois et al. ( , 2020. We compared the density simulated by SNOWPACK to snow measurements for the two layers analyzed: 1) the top layer mainly consisting of drifted snow, and 2) the lower layer consisting of depth hoar, that is, large snow grains forming a low-density snow layer (Fig. 4). Results show that the original version of SNOWPACK (circles) largely underestimates the top layer density values while the density of bottom layers is strongly overestimated, both outside of the range from field measurements (boxplots). When using the Antarctic version (triangles), both top and bottom density values greatly improve and become consistent with the observed snow density standard deviation in these layers.

Environmental Predictor Variables
Snow depth and the cumulative thickness of snow layers surpassing the critical density value of 350 kg·m -3 (CT350) were extracted from the SNOWPACK model as key predictor variables of forage condition (as described above). Other environmental predictor variables were selected based on available literature and Inuit Qaujimajatuqangit ( Table 2).
The Brazel (2006) dataset is a characterization of the BIC vegetation. This product is a vector digitization of seven maps of the vegetation communities in the Queen Elizabeth Islands produced by the Geological Survey of Canada. Because this dataset was digitized, it has an estimated error of up to 250 m. The maps were created from field surveys conducted in 1986 and include 43 vegetation classes that provide detailed information to target important areas for Peary caribou. The classes describe the habitat by the first and second most common group of species in the area (Luzula, lichen, purple saxifrage, etc.) or by type of cover (lake, wetland, unvegetated area, etc.). We did some exploratory analyses investigating whether merging vegetation classes influenced modelling outputs. Merging classes appeared to have no effect so we used the original 43 vegetation classes.
The potential impacts of increasing anthropogenic development in the Arctic region on Peary caribou are not well understood (COSEWIC, 2015). Current disturbance data consisting of GPS points of permanent human structures located on the BIC were included in the model. Since the BIC is uninhabited, the structures include about 20 oil exploration wells on Alexander Island abandoned in the 1990s. While there are very few of these locations, we integrated them because, according to local Indigenous knowledge, Peary caribou are sensitive to anthropogenic features (Miller et al., 1977;Taylor, 2005; Canadian Wildlife Service, 2013).

Configuration of the MaxEnt Model
MaxEnt is a distribution modelling software package that is one of the most widely used for species' distribution modelling in the fields of biogeography and ecology (Ahmed et al., 2015). MaxEnt is very popular because it offers users many different methods and settings to create species' distributions using maximum entropy (Phillips and Dudík, 2008;Phillips et al., 2017). Phillips et al. (2006) describe the model as a general approach for the characterization of probability distributions and for deriving predictions from incomplete information (presence-only data). Environmental variables are used to predict species' distributions, in our case, using GPS locations.
In this study, we created MaxEnt models for each Peary caribou biological season defined to reflect seasonal changes in forage availability according to different indices available in the literature (e.g., changes in body conditions, changes in snow cover, and temperature) and refined based on Inuit Qaujimajatuqangit (see Johnson et al., 2016 for more details). The seasons were defined as 1) summer (from July to October) when individuals increase fat reserves for winter and enter the rut (i.e., mating), 2) winter (from November to March) when winter extreme events impact forage availability, which affects both adult survival and pregnancy rates, and 3) spring (from April to June) when the caribou move to calving grounds and increase energy reserves for calving and lactation. We modelled the seasons separately to capture variation in behavior in response to seasonal conditions (Johnson et al., 2016). All of the data (Peary caribou presence and environmental data) were transformed into raster format with pixels of 1 × 1 km. Environmental conditions in each pixel were averaged over the entire biological season.  in models for all biological seasons since snow patches can be found throughout the year at a latitude of more than 75°. Predictor variables were tested for collinearity, and no relationship was found to be significant between any of the variables (Table 3). We used forward selection of predictor variables for inclusion in MaxEnt models. We ran single-predictor variable models first. Only variables that produced models with an area under the curve (AUC) value above 0.6 were used to develop multivariate models. This allowed us to minimize the number of models tested while ensuring the retention of variables that explained more than that expected from random chance (AUC = 0.5, Elith et al., 2006). Each model was tested with 10,000 background points for optimal convergence using 1000 iterations and 10-fold cross-validation (Kohavi, 1995). Response curves were set to linear, quadratic, and hinge features. Hinge features allow for the fitting of more complex, nonlinear functions similar to general additive models (Elith et al., 2011). We left the regularization parameter, a penalization term that prevents fitting overly complex models, set at the default of 1. Instead, we depended on Akaike information criterion (AIC) to assess overall model performance and the trade-off between model fit and complexity (Anderson et al., 2000). The top candidate models are those with the lowest AIC scores. We report the percent permutation importance as a measure of variable importance for our top models (Smith and Santos, 2020). The percentage reflects the loss in a model's ability to predict occurrences when the variable of interest is permuted randomly. Values are standardized so they sum to 1, with high percentages indicating high variable importance.
Relative probabilities of presence are expressed on a logistic scale, with zero representing the lowest rank and one the highest. We chose the logistic output over the raw exponential output for ease of interpretation and comparison with other models (Elith et al., 2011). We set τ to 0.5 to represent a typical site where the animal is present on the logistical scale; in other words, to indicate that there is a 50% chance that Peary caribou will occupy a suitable site (Elith et al., 2011). The latter seems reasonable given that 1) Peary caribou are capable of moving to any area in our study site (no dispersal barriers), 2) biases in our "presence" sampling are unlikely due to GPS collar locations (see above), and 3) the Peary caribou species is a habitat generalist, using a variety of habitats within its range to satisfy its life history requirements (see Jenkins et al., 2020).

Characterization of Snowpack Conditions
Using the Antarctic version over the Cambridge Bay meteorological study site greatly improved snow thickness and density simulations. Figure 3 shows SNOWPACK simulations at Cambridge Bay where snow data were collected during winter 2016. The figure (bottom graphs) clearly illustrates the problematic density inversion described in the literature for Arctic snow cover (Gouttevin et al., 2018;Domine et al., 2019). Figure 4 shows the results of the comparison of simulated snow values to snow measurements. The original version of SNOWPACK (circles) largely underestimated the top layer density values while the bottom layer density is strongly overestimated, with both falling outside the range of field measurements (boxplots). Both top and bottom density values were greatly improved and more consistent with the observed snow density in the layers with the Antarctic version (triangles), although the Antarctic model still underestimated the top layer density compared to on-site measurements.
Given the observed improvement in density simulations using the Groot Zwaaftink et al. (2013) Antarctic version of SNOWPACK, the methodology was applied over BIC from July 2003 to June 2004. Figure 5 shows the results of the rescaling of the spatialized snow simulations using the Antarctic model to 1 × 1 km using land cover, topography, and soil albedo information. Soil albedo has an impact on snow onset timing, however the small albedo variability over our study sites (between 13% and 20%) does not suggest that it would have a large spatial impact on snow cover duration, so the gained variability within the 32 km pixels is most likely linked to topography.
The modelled snow conditions over BIC from July 2003 to June 2004 are described in Table 4. As expected, maximal snow depth increased over this period, as did maximal CT350. The maximum value of CT350 occurred in April to June (40051.0 cm). The mean CT350 values were also very high in the November to March and April to June biological seasons (respectively 4134.6 and 4351.8 cm). Similar to Ouellet et al. (2017), we found many areas exceeding the CT350 threshold; these areas represent

Caribou Seasonal Habitat Use
The results associated with the model adopted for each caribou biological season are given in Table 5. All top seasonal models among the 34 tested included variables describing critical snow conditions (snow depth or CT350) for Peary caribou. It is worth noting that the disturbance variable was never found to improve model performance.
The top model for the July to October season had three variables: snow depth (41.3%), elevation (34.9%), and CT350 (23.8%) ( Table 6). Caribou occurred in low to mid-elevation habitats, as well as low snow depths and density (Fig. 6). The mean value of pixels in this model, which represents the probability of Peary caribou being present, is 0.66 ± 0.24. This value is greater than that expected from random chance.
The November to March model had four variables: elevation (44.2%), snow depth (29.6%), vegetation (18.1%), and CT350 (8.1%) ( Table 6). The vegetation classes most used by caribou were Luzula, lichens, and herb barrens. The dwarf shrub vegetation classes were least used. Peary caribou occurrence over large areas of the BIC was mostly low (red) in this season (Fig. 6). The areas with the highest probability of presence (in blue) were near the coasts and on the most northeastern point of the BIC (Cameron Island) where the elevation is low and snow depth is average (about 50 cm). The mean value of pixels in this model was 0.134 ± 0.218, which is substantially lower than the July to October season mean value (0.66 ± 0.24).
Peary caribou did not appear to respond to either disturbed areas or CT350 during the April to June season. In fact, these univariate models explained little more about Peary caribou space use than what would be expected by random chance (i.e., AUC = 0.5). Instead, Peary caribou habitat use was best explained by vegetation (63.9%), elevation (20.6%), and snow depth (15.5%) (Tables 5 and 6). Caribou used vegetation classes including Luzula, purple saxifrage, lakes, and wetlands. Their presence was low in unvegetated areas and herb barrens. Peary caribou occurred  in areas of low elevation and high snow depths. The mean value of all pixels for this model is 0.408 ± 0.223 (Fig. 6). All three top seasonal models were reasonably well validated. The April to June and July to October models showed moderate levels of transferability to other time periods with respect to predicting Peary caribou occurrences (r s = 0.74). While the robustness of the November to March model was slightly lower (r s = 0.69), the model was able to discern areas with a low probability of Peary caribou presence from areas with a higher probability of Peary caribou presence.

DISCUSSION
This study highlights several new and ongoing improvements to snow simulations in the Arctic and demonstrates how these simulations have value in further understanding the effects of snow properties on wildlife in species distribution models such as MaxEnt. We made various modifications to the Antarctic SNOWPACK model related to snow distribution and wind compaction  to improve modelled outputs. Refinement in spatial resolution of meteorological forcing in the future may result in better Peary caribou occurrence models for the Arctic because it will allow for a better representation of the variation in conditions within the 32 km NARR pixels. We did not improve the snow simulations by refining the spatial resolution in this study. Instead, we increased the sensitivity of the simulations to local parameters affecting grazing conditions. This process included the application of the Antarctic version of SNOWPACK and the incorporation of information on soil albedo, topography, and vegetation at a scale of 1 × 1 km. This last modification did introduce some finer-scale variability in snow conditions across our study area (Beaudoin-Galaise, 2016). The MaxEnt results clearly demonstrate the value of the snow characteristics in predicting Peary caribou presence. At least one of the modelled snow variables was selected in each season. This result is an improvement from previous attempts to model the distribution of Peary caribou (e.g., Johnson et al., 2016;Jenkins et al., 2020), where data associated with snow had high uncertainty and were limited in terms of its value in describing suitability of foraging conditions. The dynamics of snow properties that can be modelled using SNOWPACK have been identified as a key limitation to wildlife research in the Arctic, particularly under current threats from climate change (Boelman et al., 2019). For this reason, we will continue to develop these tools at spatially relevant scales based on this proof-of-concept approach. Future work will focus on using global environmental multiscale -local area model (GEM-LAM) forecasts (2.5 km resolution) to force the SNOWPACK model and will investigate their ability to improve Peary caribou distribution modelling. However, the temporal application of GEM-LAM forecasts will be limited to 2017 onwards based on available data, which highlights some of the potential trade-offs with future model improvements. We also plan to develop a snow depth distribution model, using topography indices proposed by Winstral et al. (2014). We will explore different processes affecting snow depth by describing explanatory variables using the random forest algorithm and using generalized linear and additive models, which will allow for a more defined spatial assessment of preferable locations for Peary caribou.
Both snow variables (snow depth and CT350) were important in explaining Peary caribou presence in the July to October summer season when weather conditions are mild but snow is still present on portions of the landscape (Johnson et al., 2016). Peary caribou may avoid snow to maximize vegetation accessibility in this season to increase fat reserves as quickly as possible; fat reserves are linked to improved pregnancy rates (Thomas, 1982). Vegetation was not found to be an important predictor in this model perhaps because greening at higher elevations during the summer opens more areas for grazing so that individuals select habitats based on availability of vegetation rather than specificity.
In the November to March season, caribou use of areas with average snow depth (about 50 cm) could be linked hypothetically to the trapping effect from vegetation that increases snow depth but reduces snow density (Busseau et al., 2017;Ponomarenko et al., 2019). Caribou are generally able to dig through low-density snowpacks of less than 0.5 m (Thomas and Edmonds, 1983). Caribou were not present in high snow density areas (e.g., > 350 kg·m -3 ) because it impedes foraging (Vikhamar-Schuler et al., 2013;Ouellet et al., 2017). The effects of vegetation in the model confirm the findings of Parker and Ross (1976) and Thomas et al. (1976) who reported winter feeding craters that frequently contained mosses, rushes, and lichens. Jenkins et al. (2020) also found dense and barren lichen and moss were important in their winter model. From April to June, Peary caribou used habitats at low elevations with high snow depths. This finding may be an artifact of Peary caribou moving to coastal areas where snow melts and plants emerge first (Miller et al., 1977). Feeding on high-energy plants like rushes (Luzula) and purple saxifrage is critical in this season (Miller et al., 1977;Gunn and Fournier, 2000), as females need to increase their energy reserves for calving and lactation after a very physically challenging winter period (Miller and Gunn, 2003). In such habitats, vegetation and micro-topography are likely the main determinants of snow depth and rate of snowmelt. Vegetated areas will have higher snow depths because of the trapping effect of the vegetation on snow FIG. 6. Graphical representation of the habitat model for each Peary caribou biological season. The red to yellow colours indicate habitats that will be avoided by the caribou during the season shown (presence probability 0.00 to 0.50), and the green to blue colours indicate habitats that will be selected by caribou (presence probability 0.50 to 1.00). (Sturm et al., 2001;Gouttevin et al., 2018). Despite higher snow depths, snow is unlikely an impediment in these areas at this time of year. Higher ground temperature under bare rocks in spring, in combination with a low snow depth due to wind exposure and a low albedo when exposed, will create a patchwork of areas with accelerated snowmelt.
Elevation was also an important predictor variable throughout the year with individuals tending to use lowto mid-elevation areas and avoid higher elevation areas, presumably because of the lack of shelter. The need for shelter can be critical during storms as demonstrated in a study by Dolant et al. (2018) who linked caribou mortality to sustained wind speeds that densified snow beyond the 350 kg·m -3 threshold. Interestingly, human disturbance was not an important predictor variable in any season. The negative impacts of human disturbance on caribou have been documented in a range of studies (Nelleman and Cameron, 1998;Festa-Bianchet et al., 2011;Johnson and Russell, 2014;Plante et al., 2018). The data limitations in our study only allow us to conclude that the individuals of the BIC Peary caribou subpopulation do not seem to be deterred by abandoned or inactive oil exploration wells.
The three models appear to be consistent with smallerscale foraging studies on Peary caribou. Our ability to show that the patterns persist at larger scales suggests that food or, more specifically, forage accessibility is an important limitation for Peary caribou unlike other types of caribou (Rettie and Messier, 2000). Peary caribou are habitat and forage generalists, tracking vegetation phenology across the Arctic tundra according to seasonal availability (Miller et al., 1977). It is clear from our models and others developed by Johnson et al. (2016) and Jenkins et al. (2020) that while Peary caribou are broadly distributed across Arctic landscapes, they do occupy a specialized niche and respond to both biotic and abiotic factors influencing their environment. Improving geospatial tools that can be used spatially and temporally to describe physical snow properties such as density will be important to further understand this niche and, more broadly, the role that snow plays in the ecology of other wildlife species of the Arctic (e.g., muskoxen).

CONCLUSIONS
We conducted the first spatialized snow simulations of 1 km using local topography, albedo, and soil type information over the Bathurst Island complex. The snow simulations allowed us to build on the work of Johnson et al. (2016) to produce improved models predicting the distribution of Peary caribou. Our work adds a new dimension to the findings of Vikhamar-Schuler et al. (2013), Ouellet et al. (2017), Langlois et al. (2017), and Dolant et al. (2018) by highlighting that caribou behaviour and space use are largely dependent on favourable snow conditions throughout the year. Work to replace the NARR meteorological data with GEM-LAM in the snow simulations is currently underway, as are efforts to develop the detection of rain-on-snow events and ice layers from passive microwaves (Dolant et al., 2016(Dolant et al., , 2018 for use in Arctic species distribution and habitat modelling. Work to include a vapour flux between the soil and snow to improve snow microstructure developed under temperature gradient metamorphism would pave the way to microwave coupling and snow state variable assimilation. We will assess the potential value of these future refinements to further our understanding of the possible impacts of climate change on Peary caribou.