Caves and other subterranean features provide unique environments for many species. The importance of cave microclimate is particularly relevant at temperate latitudes where bats make seasonal use of caves for hibernation. White-nose syndrome (WNS), a fungal disease that has devastated populations of hibernating bats across eastern and central North America, has brought renewed interest in bat hibernation and hibernaculum conditions. A recent review synthesized current understanding of cave climatology, exploring the qualitative relationship between cave and surface climate with implications for hibernaculum suitability. However, a more quantitative understanding of the conditions in which bats hibernate and how they may promote or mediate WNS impacts is required. We compiled subterranean temperatures from caves and mines across the western United States and Canada to: a) quantify the hypothesized relationship between mean annual surface temperature (MAST) and subterranean temperature and how it is influenced by measurable site attributes, and b) use readily available gridded data to predict and continuously map the range of temperatures that may be available in caves and mines. Our analysis supports qualitative predictions that subterranean winter temperatures are correlated with MAST, that temperatures are warmer and less variable farther from the surface, and that even deep within caves temperatures tend to be lower than MAST. Effects of other site attributes (e.g., topography, vegetation, precipitation) on subterranean temperatures were not detected. We then assessed the plausibility of model-predicted temperatures using knowledge of winter bat distributions and preferred hibernaculum temperatures. Our model unavoidably simplifies complex subterranean environments, and is not intended to explain all variability in subterranean temperatures. Rather, our results offer researchers and managers improved broad-scale estimates of the geographic distribution of potential hibernaculum conditions compared to reliance on MAST alone. We expect this information to better support range-scale estimation of winter bat distributions and projection of likely WNS impacts across the West. We suggest that our model predictions should serve as hypotheses to be further tested and refined as additional data become available. 

Subterranean microclimate data

We compiled temperature data collected from caves and mines (139 unique sites, including 75 caves and 64 mines) across western North America (9 U.S. states, 2 Canadian provinces) between 2006 and 2019 from several data providers and one published source (see publication Acknowledgements and Cited Literature). We only included datasets for which site location was provided with positional error < 10 km. Data were collected using temperature loggers (iButton, Onset HOBO loggers), which recorded temperature at subdaily intervals that varied by site. Many sites contained multiple loggers, placed at a range of distances from the site entrance (479 total loggers, mean of 3.4 per site). Some datasets were accompanied by metadata describing key attributes that may influence temperature at the logger placement site (e.g., distance from site entrance, description of site size or shape, site geology), or other logger placement details (e.g., wall, crevice, or ceiling). Although many loggers recorded humidity as well as temperature, logger saturation, a phenomenon in which the logger fails to accurately record humidity after reaching readings of 100% humidity, precluded the use of data from the majority of loggers. Therefore, we focused only on temperature here. 

We first screened all logger temperature data for deployment errors, particularly evident recording of surface temperature prior to and following placement of the logger at the site. We restricted our analysis to data recorded during the core winter months of December - February to restrict analysis to a consistent time frame that can be considered winter at all sites across the broad latitudinal range of our study. For loggers deployed over multiple winters, we summarized each winter period separately. Winter periods during which fewer than 14 days of readings were recorded were excluded from analysis (resulting in 131 sites, 427 loggers, and 699 logger-winters). We further restricted our dataset to loggers for which distance from the site entrance was provided, or for which cave maps were provided that could be used to generate a logger placement distance (54 sites, 156 loggers, 202 logger-winters). Recorded distances ranged from 0 to 324 m from the site entrance (median = 61 m). 

Because the compiled and filtered raw logger dataset contained >1.2 million observations, which were highly variable in their distribution across loggers as well as their frequency and regularity of recording, we conducted all analyses using a stratified random sample from the raw logger temperature dataset. We selected 250 records from each unique logger so as to retain all sites and loggers; for loggers with fewer than 250 records (19 loggers; 5%), all records were selected. Model results were insensitive to this choice of sample size. All logger data processing and analysis was conducted in R (version 3.4.1).

Surface climate and landscape data

We derived MAST and predictors representing key landscape attributes using Google Earth Engine, a cloud-based computing platform supporting large-scale analysis on an extensive catalog of remotely sensed, climatological, and other geospatial datasets. We estimated MAST at each cave site using the DayMet Version 3 dataset, which provides gridded daily surface temperature at 1-km resolution (1980-2018). Daily mean temperature was first calculated as the mean of daily minimum and maximum temperature. We then estimated MAST as the 20-year (1998 - 2018) mean of daily mean temperatures. 

We also derived site-level landscape predictors that are believed to impact winter cave temperatures (Perry 2012, Environmental Reviews 21:28). We extracted elevation of each site from the Shuttle Radar Topography Mission (SRTM) digital elevation model at 30-m resolution. Based on elevation, we derived a multiscale topographic position index (TPI), in which canyon and valley bottoms have low position and peaks and ridges have high position. Using a moving-window approach, TPI was calculated as the elevation of a focal raster cell minus the mean elevation within a given neighborhood surrounding the focal cell. We calculated TPI for three neighborhood sizes (500 m, 5 km, and 25 km squares), then averaged these to produce a multiscale index. We extracted Continuous Heat-Insolation Load Index (CHILI), a surrogate for effects of solar insolation and topographic shading (e.g., due to slope aspect) on evapotranspiration, from the Global ALOS CHILI product at 90-m resolution. Percent tree cover was extracted from the Terra MODIS Vegetation Continuous Fields product, which estimates percent tree cover within 250-m resolution pixels. We estimated the mean annual number of snow days at 500-m resolution from the MODIS Snow Cover dataset (V6) by averaging the number of days per year with at least 10% snow cover over the most recent 5-year period available (2014-2018), as well as mean annual snowpack (snow water equivalent (SWE) on April 1) for the same time period at 1-km resolution from DayMet (V3). We sampled groundwater table depth estimated from compiled point observations gap-filled with a mechanistic groundwater flow model (Fan et al. 2013, Science 339:940) as a proxy for the likelihood of water flow within caves or mines, assuming that flow was more likely at sites with shallower (i.e., closer to the surface) groundwater.

See uploaded ReadMe file.