The role of spring ecosystems as climate refugia in a semi-arid environment

Wet and cool microenvironments often serve as climate refugia in semi-arid regions. However, springs—locations where groundwater reaches the Earth’s surface - remain underexplored as climate refugia.This study investigated the potential of spring ecosystems as climate refugia in a semi-arid mountainous region of central Idaho, U.S.A. Using high-resolution PlanetScope imagery (2017–2024), we derived seasonal phenophases from a Normalized Difference Vegetation Index (NDVI) time series to assess ecological stability at 40 springs and surrounding non-spring areas. We fit a linear mixed effects model with phenophase as the dependent variable, spring and water year as random effects, climatic water balance (CWB), snow disappearance date (SDD), heat load index (HLI), topographic wetness index (TWI), and their interactions with site type (spring or non-spring) as predictors. We found that springs exhibited significantly lower interannual variability in end of growing season (EOS) timing (24 days less than non-springs). Higher annual CWB, reflecting greater precipitation relative to potential evapotranspiration, corresponded with later EOS timing for both springs and non-springs, but springs were less sensitive to annual CWB as shown by lower effect sizes. Springs phenology showed weak associations with TWI and HLI, underscoring their independence from topographically driven refugia. Our findings highlight springs as climate refugia due to their buffering of water limitations that stabilize late season phenology. Under climate change, water deficits will become more severe, making climate refugia like springs increasingly important. Future research should examine spring recharge processes and incorporate additional snowpack variables to monitor stability across a range of climate conditions. Title The role of spring ecosystems as climate refugia in a semi-arid environment List of Authors Grace Peven 1, Jan U.H. Eitel 1, Timothy E. Link 1, Eli W. Estey 1, Mary Engels 1 Institutional affiliations 1 College of Natural Resources, University of Idaho, Moscow, Idaho, USA Contact Information Corresponding author: Grace Peven, [email protected], [email protected]

Wet and cool microenvironments often serve as climate refugia in semi-arid regions. However, springs—locations where groundwater reaches the Earth’s surface - remain underexplored as climate refugia.This study investigated the potential of spring ecosystems as climate refugia in a semi-arid mountainous region of central Idaho, U.S.A. Using high-resolution PlanetScope imagery (2017–2024), we derived seasonal phenophases from a Normalized Difference Vegetation Index (NDVI) time series to assess ecological stability at 40 springs and surrounding non-spring areas. We fit a linear mixed effects model with phenophase as the dependent variable, spring and water year as random effects, climatic water balance (CWB), snow disappearance date (SDD), heat load index (HLI), topographic wetness index (TWI), and their interactions with site type (spring or non-spring) as predictors. We found that springs exhibited significantly lower interannual variability in end of growing season (EOS) timing (24 days less than non-springs). Higher annual CWB, reflecting greater precipitation relative to potential evapotranspiration, corresponded with later EOS timing for both springs and non-springs, but springs were less sensitive to annual CWB as shown by lower effect sizes. Springs phenology showed weak associations with TWI and HLI, underscoring their independence from topographically driven refugia. Our findings highlight springs as climate refugia due to their buffering of water limitations that stabilize late season phenology. Under climate change, water deficits will become more severe, making climate refugia like springs increasingly important. Future research should examine spring recharge processes and incorporate additional snowpack variables to monitor stability across a range of climate conditions.

Spring ecosystems, groundwater, groundwater dependent ecosystems, phenology, climate change, climate refugia, remote sensing

Climate refugia are areas that remain relatively stable amid climate variability and change, making them increasingly prioritized worldwide for safeguarding biodiversity and ecosystem services (Morelli et al., 2020; Keppel et al., 2024). The distribution of climate refugia (CR) is typically influenced by topographic and edaphic conditions that vary depending on the environmental constraints of a particular location (Dobrowski, 2011; Lawler et al., 2015). As such, what constitutes CR for one species or ecosystem may differ greatly for another (Stralberg et al., 2018; Michalak et al., 2020). For example, topographic factors like slope, aspect, and exposure (e.g., ridgelines vs. valley bottoms) form microclimate conditions which mediate physiological processes like transpiration and photosynthesis in plants (Ackerly et al., 2020). In water-limited environments, these CR often align with areas of greater water availability and reduced energy loading, such as shaded aspects, depressions, or cold air drainages that reduce evaporative demand (McLaughlin et al., 2017). Identification of CR primarily through topography risks overlooking CR shaped by hydrogeologic features, which can exert a disproportionate influence on surrounding ecosystems (Perla and Stevens, 2008; Millar and Westfall, 2019; Fey et al., 2019). Hydrogeologic features such as fissure springs, perched shallow aquifers, and areas of contrasting permeability create distinctive mesic microenvironments that provide critical water resources in otherwise arid landscapes (McLaughlin et al., 2017; Freed et al., 2019; Cartwright et al., 2020). These features may be even more critical CR than topographically driven mesic environments because they connect the surface to deep groundwater, creating potentially stable locations of water availability (Weissinger et al., 2016). Spring ecosystems are an example of a distinctive hydrogeologic feature and may serve as CR in water-limited environments (McLaughlin et al., 2017; Cartwright and Johnson, 2018; Cartwright et al., 2020). Springs are formed where the aquifer meets the Earth’s surface, discharging cold groundwater that produces a cooler, wetter microclimate compared to the surrounding landscape, particularly in semi-arid and arid regions (Stevens et al., 2021). Although spring hydrology varies, meandering subsurface groundwater flow paths create a temporal decoupling from recharge events (e.g., snowmelt) to the surface expression of groundwater at spring heads (Whiting and Godsey, 2016; Weissinger et al., 2016; McLaughlin et al., 2017). Such decoupling between regional hydroclimatic conditions and the microclimate at springs buffer the effects of interannual climate variability. Despite springs high potential as CR, their sensitivity to climate change and annual climate variability has been minimally studied (Weissinger et al., 2016), especially from an ecological perspective (Cartwright and Johnson, 2018), and primarily discussed in theoretical terms (McLaughlin et al., 2017; Cartwright et al., 2020). To better understand the potential of springs as CR, robust and accessible monitoring approaches are needed (Tang et al., 2016). Plant phenology, which tracks the timing of cyclical biological events such as leaf-out, flowering, and senescence, is a widely utilized approach for studying ecological responses to climate variability (Ford et al., 2016; Piao et al., 2019; Brooks et al., 2020; Dronova and Taddeo, 2022), making it well suited for assessing the role of springs as CR (Cartwright and Johnson, 2018). As such, elevated soil moisture at springs relative to the surrounding landscape is associated with drought resilience (Cartwright and Johnson, 2018; Fuchs et al., 2019) and faster post-fire vegetation regeneration at springs in water limited landscapes (Tsinnajinnie et al., 2021; Peven et al., 2024). While previous studies have used interannual phenological stability as an indicator of refugial potential at springs (Cartwright and Johnson, 2018), they have not explored the seasonal phenological phases, henceforth phenophases, most closely tied to water limitations, such as the end of the growing season (i.e., senescence) (Wu et al., 2022; Zhou et al., 2023). Furthermore, the interaction between topo-climatic factors, such as aspect, slope, and cold air drainages, and phenological stability at springs remains untested, to our knowledge. To address these knowledge gaps, we examined the interannual phenological variability of springs relative to the surrounding landscape and assessed how springs phenology correlates to interannual climate variability and topographic factors. Our overarching question was, do springs provide climate refugia in a semi-arid ecosystem of central Idaho ? Our study aimed to expand the understanding of springs as CR and the potential importance of incorporating hydrogeologic features in the process of CR identification.

Study area We conducted our study in the lower Big Creek watershed in central Idaho (Figure 1A). The Big Creek watershed is part of the headwaters of the Columbia River basin and sits within a federally designated wilderness area. Our study location is representative of many semi-arid mountainous environments in the western United States which are experiencing novel and intensified drought (i.e., snow, meteorological, and ecological) (Marshall et al., 2019; Moss et al., 2024). Additionally, wilderness protection provides an opportunity to study springs within a landscape effectively unaltered by modern human activities. The lower Big Creek watershed receives approximately 380 mm of annual precipitation with most of the precipitation occurring as snow in the winter and rain in the spring months, with minimal summer precipitation (gridMET; Abatzoglou, 2013). Big Creek, at University of Idaho’s Taylor Wilderness Research Station (TWRS) (Figure 1A), is characterized by an average January minimum temperature of approximately −10 °C, while average July maximum temperature is 28 °C. However, the large elevational gradient (>1200 meters) across the study area causes considerable climatic variation. The bedrock geology consists mainly of Diorite, Quartzite, and Granodiorite (Stewart et al., 2013). Dominant soil types include sandy and silty loam, but in steeply sloped areas (>30°) minimal or thin soil cover exists. Our study area is in the Idaho Batholith Hot Dry Canyons (16j; Level III/IV) EPA ecoregion (USEPA, 2013). At lower elevations there are patchy forests of ponderosa pine ( Pinus ponderosa ) and Douglas fir ( Pseudotsuga menziesii ) with sagebrush ( Artemisia ), bunchgrasses, and wildflowers. Higher elevations consist of denser forested patches of intermixed lodgepole pine ( Pinus contorta ) and Douglas fir with a shrubby understory (e.g., Ceanothus ). Spring sites generally have higher biomass than the surrounding landscape and harbor a mix of riparian and upland plant species such as Rocky Mountain maple ( Acer glabrum ), aspen ( Populus tremuloides ), red alder ( Alnus rubra ), willow ( Salix ), red osier dogwood ( Cornus sericea ), various forbs, and Douglas fir (Figures 1B & 1C). Spring selection and delineation We selected 40 hillslope and helocrene springs (Stevens et al., 2021) with full satellite visibility (i.e. no overhanging or obscuring landscape features) that were field verified as perennial spring sources (Peven et al., 2024). The distribution of springs is likely driven by a combination of geological faults, exposed bedrock, and steep topographic breaks. Elevations of spring sites range from 1183 - 2383 meters. Initially, we digitally delineated each spring’s footprint (i.e., area of elevated soil moisture influence) based on field observations and visual determination of spring-dependent vegetation from the World Imagery basemap in ArcGIS Pro (Esri et al., 2021) (0.3 meter pixel resolution). A 3-meter buffer was added around each spring footprint to ensure the inclusion of all surrounding spring-dependent vegetation. Each spring footprint polygon was duplicated to a randomly chosen area directly adjacent to the spring site with the same aspect, slope, and elevation to serve as a paired control (hereafter, non-spring) site that represented non-spring dependent vegetation (Figure 2A). To exclude non-vegetated pixels (e.g., rocks or bare soil) from the delineated spring and non-spring footprints, we determined a vegetation threshold for each site type based on seasonal maximum Normalized Difference Vegetation Index (NDVI) values. Using NDVI derived from PlanetScope (Planet Labs PBC, 2024) imagery (3x3 meter pixel resolution), we extracted the day of year (DOY) with the maximum NDVI between April and October 2023 from the preliminary spring and non-spring footprints (Figure 2A). Generally, the maximum NDVI occurred around mid-July across the entire study area, so we selected the highest quality image (zero cloud cover and full study area coverage) around this timeframe (July 22 nd, 2023). For all spring and non-spring sites we ran a sensitivity analysis to determine the NDVI threshold that produced the best footprint delineation. We considered the best delineation to be one that did not over- or under-exclude pixels within our initial footprint and matched our field observations of spring footprint area and size. Our sensitivity analysis resulted in a ≥0.6 NDVI threshold for spring sites and a ≥0.4 NDVI threshold for non-spring sites. Finally, the footprints were revised to only include pixels that met the described threshold criteria (Figure 2B). Satellite imagery processing PlanetScope multi-spectral surface reflectance imagery harmonized to Sentinel-2 collected between 2017 – 2024 was downloaded from Planet Explorer (Planet Labs PBC, 2024). All years where PlanetScope imagery was available were included. Each image was atmospherically corrected and orthorectified prior to downloading and contains four bands (blue, green, red, near infrared (NIR)) with a 3-meter spatial resolution and daily temporal resolution. PlanetScope imagery was selected for this study due to its high spatial and temporal resolution, thus its ability to detect small spring footprints and to capture fine temporal changes in phenology (Zhao et al., 2022; Eitel et al., 2023). All available images from PlanetScope satellite sensors (Dove Classic, Dove-R, and SuperDove) were obtained and subsequently filtered for high quality using the Planet-provided Usable Data Mask (UDM) GeoTIFFs to remove pixels containing snow, haze, and/or clouds. To quantify phenophases based on PlanetScope data, we extracted the mean NDVI (Eq. 1) time series across each spring and non-spring footprint from DOY 100 – 300 (~April– October) representing the growing season across all years using the terra package in R (Hijmans, 2024). NDVI time series have been widely used for remotely monitoring phenology (Zeng et al., 2022) including in some of the seminal work by Cartwright and Johnson (2018) that specifically focused on using NDVI to assess the drought resilience of springs.\(\mathbf{\text{NDVI}}=\ \frac{\rho\text{NIR}\ –\ \rho\text{Red}}{\rho\text{NIR}+\rho\text{Red}}\) (Eq. 1) While PlanetScope nominally has daily temporal resolution, cloud cover often limits the availability of daily imagery. To address gaps in available imagery, we applied a generalized additive model (GAM) to the NDVI dataset for each spring, site type (spring or non-spring), and year, using the mgcv package (Wood, 2011) in R (R Core Team, 2024). We iteratively fit GAM models with varying smoothing dimensions (5 through 10) and selected the model with the lowest Akaike information criterion (AIC) value (Akaike, 1974). We then used our final models to predict NDVI values for DOYs 100 - 300 for each group resulting in our final greenness time series (Figure 3). Refugial indicators and analysis Similar to Cartwright and Johnson (2018), we interpreted the following as indicators of climate refugia: 1) less interannual phenological variability (lower interquartile range (IQR)) at springs relative to non-spring sites, and 2) less sensitivity (lower effect size) to interannual climate conditions at springs relative to non-spring sites. Using our final NDVI (i.e., greenness) time series, we extracted a suite of phenophases commonly utilized in land surface phenology (Table 1; Figure 3) (Zeng et al., 2020). Each of the phenophases listed in Table 1 were calculated per spring, site type, and year (2017 – 2024). We used a repeated measures analysis of variance (ANOVA) to test whether there were significant differences in interannual phenophases and variability across all years between springs and their paired non-spring sites (Table 1) using the lme4 package (Bates et al., 2015). We checked for model assumptions using the Diagnostics for Hierarchical Regression Models (DHARMa) package (Hartig, 2022). When the ANOVA test resulted in significant differences (α <0.05) we ran a post-hoc pairwise test using the emmeans package with a Tukey adjustment to estimate the average differences between phenophases at springs and non-springs (Lenth, 2024). We also conducted paired t-tests to compare the mean interquartile range (IQR) between springs and non-springs for each phenophase. Interannual climate effect on phenology To assess how phenology is affected by interannual climate variability and topographic factors at spring and non-spring sites, we included the following variables in our analysis: climatic water balance (CWB), snow disappearance date (SDD), topographic wetness index (TWI), and the heat load index (HLI). The CWB represents a functionally integrative measure of water and energy availability and was calculated by subtracting the potential evapotranspiration (PET) from precipitation for each spring. PET was estimated using the Thornthwaite method (Thornthwaite, 1984) with daily temperature data from a weather station at TWRS. Daily air temperature was the only climate variable available from the TWRS weather station, making the Thornthwaite method the most practical option for calculating PET. When there were gaps in daily temperature, we used the linear relationship between maximum and minimum daily temperature from the closest weather station in Yellow Pine, Idaho and TWRS to build a complete daily temperature dataset across the study timeframe (Max temp R 2 = 0.81, RMSE = 5.2; Min temp: R 2 = 0.86, RMSE = 3.1). We then used the mean environmental lapse rate (-6.5 deg C/1000 meters) to estimate daily temperatures at each spring site. For precipitation, we used the 800-meter Parameters-elevation Regressions on Independent Slopes Model (PRISM) (PRISM Climate Group, 2014) data to calculate total annual, monthly, and seasonal (winter, spring, summer, and fall) precipitation at each site. We also hypothesized that SDD would influence phenology (Slatyer et al., 2022). We estimated the SDD at each site using the Normalized Difference Snow Index (NDSI) derived from the Moderate Resolution Imagery Spectroradiometer (MODIS) Terra satellite sensor (500 m pixel resolution /daily resolution) downloaded from NASA’s Application for Extracting and Exploring Analysis Ready Samples’ point request service (AppEEARS; AppEEARS Team, 2024). MODIS was selected for SDD since it contains the shortwave infrared (SWIR) band necessary for NDSI calculation. Our study area generally has high cloud cover and patchy image availability in the spring months, so we fit a GAM spline to all NDSI values at each spring site. From the fitted spline, we identified the annual SDD as the day of year when NDSI reached 0.4, which is a common SDD threshold (Nolin, 2010). To account for topographically-mediated microclimate influence on phenology, we calculated the TWI (Beven and Kirkby 1979) and the HLI (McCune and Keon, 2002) using a 10-meter digital elevation model (DEM) (USGS, 2024) and the Topography toolbox in ArcGIS Pro (Dilts, 2023) (Table 2). To statistically examine climate and topographic effects on each given phenophase, we fit the following linear mixed effects model: Phenophase ~ Climate variables*site type + Topographic variables*site type + spring name + water year with climate and topographic variables as fixed effects, including interactions between these variables and site type (spring or non-spring), and spring name and water year as random effects. All predictor variables were tested for collinearity using Spearman’s rank tests and variables with a Spearman’s ρ > |0.4| were not included in the same model. SDD and CWB were collinear (Spearman’s ρ = 0.46), so we fit two different models for these variables. We compared the standardized relative effect sizes (β) (i.e., z-scores) of each predictor between spring and non-spring sites. Model fit was assessed using the conditional and marginal R 2 values, root mean squared error (RMSE), and AIC. We again used the DHARMa package to test model assumptions.

Interannual climate variability While our analyses were confined to the 8 years where Planet satellite data was available, we observed notable variability in interannual climate conditions within the study timeframe (Table 3; Figure 4). For example, 177% of average (30-year) winter precipitation occurred in 2017 (149 mm) and 43% of average spring precipitation occurred in 2021 (59 mm) (Table 3), the latter aligning with the spring drought observed across the western U.S. (Affram et al., 2023). The study area is generally in a climatic water deficit, meaning that more PET occurs than precipitation each year, but this varies across the study area with elevation and terrain sheltering that produces distinct meso- and microclimates. Phenological differences between springs and non-springs We recorded a total of 636 phenological observations for each of our phenophases (40 springs and 40 non-springs over 8 years excluding 4 outlier observations). Our results indicate significant (α < 0.05) differences in all phenophases between springs and non-springs (Figure 5; Table 4). Notably, EOS timing was on average 23 days later for springs compared to non-springs (Table 4). Start of season (SOS) timing was later at springs, but only by 6 days on average. Across the study area, springs had an average growing season length (GSL) that was 17 days longer than non-springs. Unlike the magnitude, the variability as quantified by the interquartile range (IQR) did not significantly differ between springs and non-springs for all phenophases (Table 5). The IQR values for EOS timing and GSL differed significantly (α <0.05) between site types, but GSL variability is likely driven by EOS timing since SOS timing is similar between sites. The variability in EOS timing is 24 days less at springs relative to non-springs. On an annual basis, EOS variability was consistently smaller at springs except for in years 2017 and 2023 (Figure 6). In 2022, springs had a later EOS compared to other years, while non-springs experienced the largest variability compared to all other years. Climatic drivers of EOS phenology Interannual climate effect Since EOS timing and variability diverged most significantly between springs and non-spring sites we focused our climate sensitivity analysis on EOS timing. Additionally, EOS timing (i.e., senescence) is generally more limited by water availability in our study region and thus serves as a better indicator of the role of springs as climate refugia. Higher annual CWB values, indicating higher precipitation totals relative to PET, corresponded with significantly later EOS timing for both springs (standardized β = 0.49, p = <0.001) and non-springs (standardized β = 0.83, p = <0.001) (Table 6; Figure 7a). SDD had opposite directional effects on EOS timing between springs and non-springs (Figure 7b). SDD had a significant negative effect on EOS timing at springs (standardized β = - 0.11, p = 0.003), suggesting that a later SDD corresponds to earlier EOS at springs, while at non-springs SDD corresponds to later EOS timing (standardized β = 0.17, p = <0.001). Both the CWB and SDD marginal standardized effect sizes were higher for non-springs, suggesting a stronger coupling to interannual climate conditions relative to springs. We iterated through several model versions with annual, seasonal, and monthly CWB totals, but annual CWB totals consistently produced the best fit based on ANOVA model comparisons and ΔAIC (Akaike, 1974). Our model that included annual CWB, HLI, and TWI as predictors explained 58.5% (marginal R 2 ) of EOS timing variability and the full model with the random effects (water year and spring name) explained 78.7% (conditional R 2 ) of EOS timing variability (Appendix S1; Table S1). Our SDD model also explained a high proportion of variability in EOS timing with the full random effects structure (conditional R 2 = 77.5%), but SDD alone did not explain much variability (marginal R 2 = 17%). Both models had relatively high RMSE values (13.3 and 13.5 days), suggesting that they have limited accuracy for predictive purposes. Topography effect We observed a median TWI value of 3.6 with values ranging from 0.79 to 6.73 across sites. The median HLI value was 0.6 with values ranging from 0.28 to 1.04. Wetter sites, as indicated by increasing TWI, had a slightly negative yet non-significant relationship to EOS timing at springs (standardized β = -0.07, p = 0.097) and a positive and significant relationship to EOS timing at non-springs (standardized β = 0.11, p = 0.001) (Table 6; Figure 8a). Hotter and drier topographic sites (e.g., southerly aspects), as indicated by increasing HLI, had a small non-significant and positive association with EOS timing at springs (standardized β = 0.04, p = 0.306) and a significant and larger negative effect on EOS timing at non-springs (standardized β = -0.14, p = <0.001) (Table 6; Figure 8b),

While previous studies have explored the refugial potential of springs using NDVI variability (Cartwright and Johnson, 2018), they haven’t used high spatial (3 meters) and temporal (~daily) resolution imagery and phenophases more closely tied to regional climate constraints. Here, we used EOS timing derived from high-resolution PlanetScope satellite imagery to demonstrate that springs provided CR based on lower EOS timing variability across a range of climatic conditions (including the 2021 spring drought; Figure 4) and sensitivity to interannual climate relative to the surrounding landscape. Additionally, the lack of a significant association between spring phenology and topographic variables (TWI and HLI) and opposing directional relationships relative to non-springs, emphasizes the independence of these hydrogeologic features from topographic-based CR identification. Although hydrogeologic features such as springs are increasingly recognized for their role in providing CR, their incorporation into CR identification remains largely underexplored (Miller and Westfall, 2019; Cartwright et al., 2020; Morelli et al., 2020; Ishiyama et al., 2023). Springs have relatively stable interannual EOS timing Monitoring the interannual variability in NDVI is a widely used approach to measure stability across ecosystems (White and Lewis, 2011; Liu et al., 2016; White et al., 2022). While few studies have focused on spring ecosystems, Cartwright and Johnson (2018) utilized the mean and standard deviation of July NDVI values to quantify spring resilience to drought across several decades in a semi-arid region of Oregon. However, focusing on phenophases that are more mechanistically linked to the climatic constraints of semi-arid environments offers more nuanced insights into whether springs function as CR (Zhang et al., 2020). The lack of difference in SOS IQR between springs and non-springs suggests that groundwater from springs exerts during this phenophase and is likely driven by other constraints such as photoperiod and temperature regulation (Meng et al., 2020). By contrast, EOS timing had the largest and most significant difference in variability (IQR) between springs and non-springs suggesting groundwater from springs strongly influences this phenophase. Thus, when using phenology to monitor or identify CR, conclusions may differ depending on whether early-season, mid-season, or late-season dynamics are examined, given their distinct climatic drivers. Springs buffer interannual climate variations Climatic water balance Our results suggest that larger CWB values (i.e., more annual precipitation than annual PET) delay EOS timing in springs and non-springs. However, the relative effect of CWB is larger for non-springs than springs, implying that springs are less sensitive to interannual fluctuations in water availability driven by the combined effects of precipitation and evapotranspiration, making them an important CR in semi-arid environments with little summer precipitation. Further, EOS timing at springs during the severe drought year of 2021 showed little deviation from average values (Figure 6), indicating a buffering effect against drought conditions that are likely to become more frequent in the future (Moss et al., 2024). Consistent soil moisture from groundwater at springs likely dampens the effect of CWB on EOS timing and allows for longer growing seasons and later EOS timing compared to non-springs. This finding aligns with previous studies in semi-arid/arid settings that found water-related conditions (e.g., soil moisture, precipitation) as a main driver of EOS timing across vegetation types (Tercek et al., 2023; Wu et al., 2022; Zhou et al., 2023; Kloos et al., 2024). Importantly, the drivers of EOS timing are highly contextual to the limiting environmental conditions of a site (Dobrowski, 2011; Lawler et al., 2015). There is growing agreement that humid environments have a negative correlation between precipitation and EOS timing, while the opposite has been observed in arid environments (Liu et al., 2016; Zhou et al., 2023). For example, Kloos et al. (2024) found that warmer and drier conditions advanced EOS timing in a German deciduous forest. Conversely, and similar to our results, in a northern hemisphere-wide phenology study, Liu et al. (2016) demonstrated that in arid/semi-arid environments, total precipitation had a positive and larger climatic effect on EOS timing compared to temperature. Snowmelt timing (SDD) Snow disappearance date (SDD) had a negative effect on EOS timing at springs, suggesting that later SDD advances EOS timing. Contrary to this finding, at non-springs, the opposite phenomenon was observed. In both cases, however, it is important to note that SDD effect sizes were relatively small compared to those of CWB. There is ample evidence that earlier snowmelt timing is associated with advanced springtime phenology in plants (Winkler et al., 2018; Slatyer et al., 2022; Kelsey et al., 2020), however less evidence exists for an association between snowmelt and EOS timing (Potter, 2020), particularly in groundwater dependent ecosystems. Later snowmelt timing contributes to the persistence of shallow summer soil moisture (Blankinship et al., 2014; Harpold and Molotch, 2015) and thus could explain the positive effect on EOS timing at non-spring sites. Whereas at springs, earlier snowmelt timing may correspond to delayed EOS timing due to deeper snowmelt infiltration into the soil profile, recharging soil moisture at springs through interflow, and sustaining higher water availability at springs into the summer (Carroll et al., 2019). However, it is challenging to disentangle these climatic drivers since snowmelt, temperature, and precipitation dynamics interact with topography to shape summer soil moisture conditions. Topography does not predict phenological stability at springs Both TWI and HLI had insignificant effects on EOS timing at springs, confirming our expectation that spring phenology is largely independent of topographically-mediated microclimate. TWI had a small negative effect on EOS timing, suggesting that springs located within wetter and cooler landscape positions have earlier EOS timing. Springs typically have ample soil moisture, thus are more limited by temperature, so if they are located in ravines, depressions, and/or cold air pools, temperature limitations could be exacerbated by additional topographic sheltering (Pastore et al., 2022). Conversely, delayed EOS timing occurred at non-spring sites with increased TWI and decreased HLI, supporting the idea that in drier landscapes topographic sheltering from high evaporative demand provides microrefugia (Dobrowski, 2011). Our results diverge from Cartwright and Johnson’s (2018) study which found that springs at higher elevations with reduced topographic heat load (steep slopes, north-facing aspects) were associated with greater drought resilience compared to lower elevation and higher heat loaded springs. However, their study used the average and standard deviation of July NDVI across a 30-year period which may explain why our results for EOS timing differed.

Since our goal was to quantify generalized ecological stability at springs, we did not incorporate species-specific data or climate change velocity modeling as is common in climate change refugia studies (Barrows et al., 2020; Morelli et al., 2016; Keppel et al., 2024). We recognize that springs and the surrounding landscape possess different plant functional types, thus different physiological and evolutionary adaptations to environmental conditions which may affect EOS timing (Xie et al., 2018; Thoma et al., 2019). Availability of empirical and modeled climate data was a limiting factor in our climate effects analysis due to the ungauged and remote location of our study in a wilderness area. To more fully understand the hydrometeorological-spring dynamics in snow-dominated environments, more snowpack metrics are needed. Potential snowpack metrics for future research could include the fraction of liquid to solid precipitation, annual snow water equivalent (SWE) totals, and rate of snowmelt. We explored using the snow data assimilation system (SNODAS) model (NOHRSC, 2004), but when comparing to several regional snowpack telemetry network (SNOTEL) sites we found poor agreement in annual trends thus had low confidence in using these data (data not shown). Future work should explore additional high-resolution snow datasets (e.g., SnowModel, Liston and Elder, 2006) to capture snow-phenology relationships. We hypothesize that the high EOS timing variability at springs in years 2017 and 2023 (Figure 6) may be due to snowpack dynamics but we were unable to assess the relationships with the available data. Additionally, snowmelt lag effects from previous water years may impact phenology at springs (Weissinger et al., 2016). We explored the correlation between 1–4-year lags in SDD and EOS timing but found no strong evidence for lag effects (data not shown). However, future work should continue exploring lag effects with SWE-based metrics and consider a dynamic snow disappearance NDSI threshold to better capture the variability in snowmelt timing in mountainous terrain (Ma and Zhang, 2022). Future research considerations Global climate change models forecast increasingly hot and dry summers, increased interannual snowpack variability, and shifts to more rain than snow in the western United States (Klos et al., 2014; Marshall et al., 2019; Overpeck and Udall, 2020; Siirila-Woodburn et al., 2021). Increases in the rain to snow ratio may reduce groundwater recharge (Condon et al., 2020), thus the water availability (i.e., refugial stability) of springs. Since our models had low predictive power, it is unclear whether springs will continue to provide climate refugia in the future with continued baseline shifts in climate. Weissinger et al. (2016) found a clear relationship between diminished spring discharge and high PET in an arid location of the southwestern United States, but the strength of the discharge and PET relationship varied between springs. Springs that experience a reduction in water availability could become more suitable for upland species and transition into refugia for new plant communities (Morelli et al., 2020; Halofsky and Hibbs, 2009). This possible transition to new plant composition at springs may lead to positive climate feedback loops, reinforcing reduced soil moisture at springs (Penuelas et al. 2009). Springs fed by deeper and longer subsurface flow paths may have delayed responses to changing groundwater recharge dynamics and thus continue to provide CR. More investigation into the specific characteristics (e.g., flowpath modeling, groundwater aging, recharge area physiography) of individual springs in connection to hydrological and phenological stability is needed (Freed et al., 2019).

We demonstrated that springs provide CR in a semi-arid mountainous environment using indicators of temporal phenological stability and interannual climate sensitivity. Both end of growing season (EOS) variability and coupling to interannual climate conditions (i.e., climatic water balance and snowmelt timing) were significantly less at springs compared to the surrounding landscape. Importantly, topographic factors were not significant predictors of EOS timing variability, highlighting the necessity of including hydrogeologic features independent of topography in CR identification. Our findings suggest that spring ecosystems provide CR by offering a buffer to the increasing water limitations that impact the stability of late growing season phenology. Reduced phenological variability at springs could limit phenological mismatch and impacts to biodiversity in arid or semi-arid environments (Renner and Zohner, 2018). Future research should explore spring recharge characteristics and additional climate variables (e.g., snowpack dynamics) in monitoring the refugial stability of springs. Rising annual air temperatures will exacerbate climatic water deficits, accelerate the depletion of summer soil moisture, and consequently advance EOS timing in arid and semi-arid environments making CR like springs increasingly important. Author Contributions Grace Peven : Conceptualization, Methodology, Software, Formal Analysis, Data Curation, Writing – Original Draft, Visualization, Funding acquisition. Jan UH Eitel : Conceptualization, Methodology, Software, Writing –Review & Editing, Supervision. Timothy Link : Conceptualization, Writing – Review & Editing. Eli Estey : Software, Writing – Review & Editing. Mary Engels : Supervision, Writing – Review & Editing, Funding acquisition.

The authors would like to thank Trey Dejong for statistical help and Andrew Armstrong for helping retrieve air temperature data from TWRS. Funding Information Funding for this study was provided in part by the DeVlieg Foundation and the Curt and Adele Berklund Foundation. Data Availability Statement The data and code that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.14783217 Tables Table 1. Phenophases derived from an NDVI time series at each spring and non-spring site across the study timeframe (2017-2024). | Max NDVI timing (i.e., timing of peak greenness) | Maximum NDVI DOY | The timing of peak greenness indicates maximum annual plant growth. | | Start of season (SOS) and end of season (EOS) timing | DOY pre- (SOS) and post- (EOS) seasonal NDVI maximum where NDVI reaches 50% of seasonal amplitude | Start and end of the growing season is typically triggered by climate conditions (e.g., temperature, photoperiod, water availability) and indicates sensitivity to interannual conditions. | | Duration of growing season | Number of days between SOS and EOS | The length of growing season is important for plant growth and food/shelter availability for surrounding ecosystem and indicates response to climate conditions. | Table 2. Climatic and topographic variables used in the linear mixed effects regression model to assess their effects on phenophases at springs and non-spring sites. | Climatic Water Balance (CWB) | TWRS weather station (air temperature) and PRISM (precipitation) | Precipitation - PET Lower CWB values (mm) indicate drier and warmer conditions; higher CWB values indicate wetter and cooler conditions. | | Snow disappearance date (SDD) | 500-meter MODIS-Terra NDSI | SDD determined from MODIS NDSI threshold of 0.4. | | Topographic Wetness Index (TWI) | 10-meter digital elevation model (DEM) | Based on slope and flow accumulation, low TWI values (3) values indicate wetter sites. TWI can also serve as a proxy for cold air pools. | | Heat Load Index (HLI) | 10-meter digital elevation model (DEM) | Incident solar radiation based on aspect, slope, and latitude. Values closer to 0 indicate wetter and cooler topo-climate conditions (e.g., northwest-facing slopes), while values closer to 1 indicate warmer and drier topo-climate conditions (e.g., southeast-facing slopes). | Table 3 . Annual climate conditions from 2017 – 2024. Seasonal cumulative precipitation was categorized into seasons: fall (October – December), winter (January - March), spring (April - June), summer (July - September). All data are from gridMET centered at Taylor Wilderness Research Station (TWRS) (Abatzaglou, 2014). | 30-year average (1991-2020) | 381 | 97 | 84 | 136 | 65 | -134 | | 2017 | 490 | 127 | 149 | 147 | 69 | -92 | | 2018 | 419 | 126 | 116 | 155 | 23 | -129 | | 2019 | 360 | 77 | 80 | 124 | 80 | -168 | | 2020 | 396 | 76 | 123 | 164 | 33 | -118 | | 2021 | 270 | 93 | 69 | 59 | 50 | -307 | | 2022 | 373 | 98 | 68.3 | 174 | 34 | -180 | | 2023 | 397 | 82 | 70 | 159 | 86 | -155 | | 2024 | 323 | 68 | 79 | 118 | 58 | -234 | Table 4. Results from repeated measures ANOVA and post-hoc pairwise tests. The estimate is the difference between the spring and non-spring site with the +/- sign indicating the directional difference relative to the non-spring site. DOY = day of year, df = degrees of freedom, SE = standard error. | Max NDVI DOY | 595 | +9.46 | 0.86 | 10.9 | <0.0001 | | SOS DOY | +5.96 | 0.52 | 11.5 | <0.0001 | | | EOS DOY | +23.0 | 1.24 | 18.6 | <0.0001 | | | GSL | +17 | 1.4 | 12.2 | <0.0001 | Table 5. Results from a paired t-test between springs and non-springs IQR values for each phenophase across all years. The estimate is the difference between the spring and non-spring site with the +/- sign indicating the directional difference relative to the non-spring site. DOY = day of year, IQR = interquartile range, df = degrees of freedom. | Max NDVI DOY IQR | 7 | +1.78 | -4.5, 7.9 | 0.68 | 0.51 | | SOS DOY IQR | -1.19 | -4.3, 1.9 | -0.89 | 0.4 | | | EOS DOY IQR | -24.06 | -39.3, -8.8 | -3.74 | 0.007 | | | GSL IQR | -17.69 | -30, -5.4 | -3.39 | 0.011 | Table 6. Results of the linear mixed effect models. Significant results are bolded. β = standardized coefficient estimate, CI = confidence interval | Predictor | β | p-value | 95 % CI | β | p-value | 95 % CI | | Annual CWB (mm) | 0.49 | <0.001 | 0.37– 0.61 | 0.83 | <0.001 | 0.76 – 0.91 | | SDD* | -0.11 | 0.003 | -0.19 – -0.04 | 0.17 | <0.001 | 0.1 – 0.24 | | TWI | -0.07 | 0.097 | -0.15 – 0.01 | 0.11 | 0.001 | 0 – 0.21 | | HLI | 0.04 | 0.306 | -0.04 – 0.12 | -0.14 | <0.001 | -0.23 – -0.06 | | SDD effect was estimated in a separate model than the other listed predictors but included here for comparison | |||||| | * the non-spring estimates are the sum of the main (spring) and interaction effects | Figure Legends Figure 1. (A) Locations of 40 springs included in our study located in the lower Big Creek watershed in central Idaho, U.S.A. (B) Landscape view of an example spring in our study area and (C) a zoomed in view of the spring with a distinctive strip of shrubby spring-dependent vegetation. Figure 2 . (A) Initial spring and non-spring footprint delineations with July 2023 Normalized Difference Vegetation Index (NDVI) values and (B) final spring and non-spring footprints based on respective NDVI thresholds. Figure 3 . Example of a Normalized Difference Vegetation Index (NDVI) time series (green line) with phenophases marked by dashed black lines. See Table 1 for detailed explanation of each phenophase. Figure 4 . Seasonal precipitation difference in each water year (October 1 to September 30) relative to the 30-year average (1991-2020). Each color corresponds to seasons as defined in Table 3. Data from gridMET centered at Taylor Wilderness Research Station (Abatzoglou, 2013). Figure 5. Phenophase differences between all spring and non-spring sites summarized across all years (2017-2024). Figure 6 . Interannual EOS timing at springs and non-springs from 2017 to 2024. Figure 7. Interannual climate marginal fixed effects (solid lines) and their 95% confidence intervals (shaded area around each line) on spring and non-spring end of growing season day of year for (a) total annual climatic water balance (CWB) in millimeters (mm) and (b) SDD (day of water year). The coefficient (coeff) estimate from the linear mixed effects model indicates the standardized relative effect size of each predictor. Figure 8 . Topographic marginal fixed effects (solid lines) and their 95% confidence intervals (shaded area around each line) on spring and non-spring end of growing season day of year for the (a) topographic wetness index and (b) heat load index. The coefficient (coeff) estimate (β) from linear mixed effects model indicates the standardized relative effect size of each predictor.

1. Abatzoglou, J. (2013). Development of gridded surface meteorological data for ecological applications and modeling. International Journal of Climatology, 31, 121–131. https://doi.org/10.1002/joc.3413 Ackerly, D. D., Kling, M. M., Clark, M. L., Papper, P., Oldfather, M. F., Flint, A. L., & Flint, L. E. (2020). Topoclimates, refugia, and biotic responses to climate change. Frontiers in Ecology and the Environment, 18 (5), 288–297. https://doi.org/10.1002/fee.2204 Affram, G., Zhang, W., Hipps, L., & Ratterman, C. (2023). Characterizing the development and drivers of 2021 Western US drought. Environmental Research Letters, 18 (4), 044040. https://doi.org/10.1088/1748-9326/acc95d Akaike, H. (1974). A new look at the statistical model identification. IEEE Transactions on Automatic Control, 19 (6), 716–723. IEEE Transactions on Automatic Control. https://doi.org/10.1109/TAC.1974.1100705 AppEEARS Team. (2024). Application for Extracting and Exploring Analysis Ready Samples (AppEEARS). Ver. 3.69. NASA EOSDIS Land Processes Distributed Active Archive Center (LP DAAC), USGS/Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, USA. Accessed November 10, 2024. https://appeears.earthdatacloud.nasa.gov Barrows, C. W., Ramirez, A. R., Sweet, L. C., Morelli, T. L., Millar, C. I., Frakes, N., Rodgers, J., & Mahalovich, M. F. (2020). Validating climate-change refugia: Empirical bottom-up approaches to support management actions. Frontiers in Ecology and the Environment, 18 (5), 298–306. https://doi.org/10.1002/fee.2205 Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. Beven, K. J., & Kirkby, M. J. (1979). A physically based, variable contributing area model of basin hydrology / Un modèle à base physique de zone d’appel variable de l’hydrologie du bassin versant. Hydrological Sciences Bulletin, 24 (1), 43–69. https://doi.org/10.1080/02626667909491834 Blankinship, J. C., Meadows, M. W., Lucas, R. G., & Hart, S. C. (2014). Snowmelt timing alters shallow but not deep soil moisture in the Sierra Nevada. Water Resources Research, 50 (2), 1448–1456. https://doi.org/10.1002/2013WR014541 Brooks, B.-G. J., Lee, D. C., Pomara, L. Y., & Hargrove, W. W. (2020). Monitoring broadscale vegetational diversity and change across North American landscapes using land surface phenology. Forests, 11, 606. https://doi.org/10.3390/f11060606 Carroll, R. W. H., Deems, J. S., Niswonger, R., Schumer, R., & Williams, K. H. (2019). The Importance of Interflow to Groundwater Recharge in a Snowmelt-Dominated Headwater Basin. Geophysical Research Letters, 46 (11), 5899–5908. https://doi.org/10.1029/2019GL082447 Cartwright, J., & Johnson, H. M. (2018). Springs as hydrologic refugia in a changing climate? A remote-sensing approach. Ecosphere, 9 (3), e02155. https://doi.org/10.1002/ecs2.2155 Cartwright, J. M., Dwire, K. A., Freed, Z., Hammer, S. J., McLaughlin, B., Misztal, L. W., Schenk, E. R., Spence, J. R., Springer, A. E., & Stevens, L. E. (2020). Oases of the future? Springs as potential hydrologic refugia in drying climates. Frontiers in Ecology and the Environment, 18 (5), 245–253. https://doi.org/10.1002/fee.2191 Condon, L. E., Atchley, A. L., & Maxwell, R. M. (2020). Evapotranspiration depletes groundwater under warming over the contiguous United States. Nature Communications, 11 (1), Article 1. https://doi.org/10.1038/s41467-020-14688-0 Dilts (2023) Topography Toolbox for ArcGIS Pro. University of Nevada Reno. Available at: https://www.arcgis.com/home/item.html?id=247fbe56c7ff48229c9b1fe132d1b5e9 Dobrowski, S. Z. (2011). A climatic basis for microrefugia: The influence of terrain on climate. Global Change Biology, 17 (2), 1022–1035. https://doi.org/10.1111/j.1365-2486.2010.02263.x Dronova, I., & Taddeo, S. (2022). Remote sensing of phenology: Towards the comprehensive indicators of plant community dynamics from species to regional scales. Journal of Ecology, 110 (7), 1460–1484. https://doi.org/10.1111/1365-2745.13897 Eitel, J. U. H., Basler, D., Braun, S., Buchmann, N., D’Odorico, P., Etzold, S., Gessler, A., Griffin, K. L., Krejza, J., Luo, Y., Maguire, A. J., Rao, M. P., Vitasse, Y., Walthert, L., & Zweifel, R. (2023). Towards monitoring stem growth phenology from space with high resolution satellite data. Agricultural and Forest Meteorology, 339, 109549. https://doi.org/10.1016/j.agrformet.2023.109549 Esri, Maxar, Earthstar Geographics, and the GIS User Community. ”World Imagery” [basemap]. Scale Not Given. ”World Imagery”. June, 17, 2021. https://services.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer. Fey, S. B., Vasseur, D. A., Alujević, K., Kroeker, K. J., Logan, M. L., O’Connor, M. I., Rudolf, V. H. W., DeLong, J. P., Peacor, S., Selden, R. L., Sih, A., & Clusella-Trullas, S. (2019). Opportunities for behavioral rescue under rapid environmental change. Global Change Biology, 25 (9), 3110–3120. https://doi.org/10.1111/gcb.14712 Freed, Z., Aldous, A., & Gannett, M. W. (2019). Landscape controls on the distribution and ecohydrology of central Oregon springs. Ecohydrology, 12 (2). https://doi.org/10.1002/eco.2065 Ford, K. R., Harrington, C. A., Bansal, S., Gould, P. J., & St. Clair, J. B. (2016). Will changes in phenology track climate change? A study of growth initiation timing in coast Douglas-fir. Global Change Biology, 22 (11), 3712–3723. https://doi.org/10.1111/gcb.13328 Fuchs, L., Stevens, L. E., & Fulé, P. Z. (2019). Dendrochronological assessment of springs effects on ponderosa pine growth, Arizona, USA. Forest Ecology and Management, 435, 89–96. https://doi.org/10.1016/j.foreco.2018.12.049 Halofsky, J. E., & Hibbs, D. E. (2009). Controls on early post-fire woody plant colonization in riparian areas. Forest Ecology and Management, 258 (7), 1350–1358. https://doi.org/10.1016/j.foreco.2009.06.038 Harpold, A. A., & Molotch, N. P. (2015). Sensitivity of soil water availability to changing snowmelt timing in the western U.S. Geophysical Research Letters, 42 (19), 8011–8020. https://doi.org/10.1002/2015GL065855 Hartig F (2024). _DHARMa: Residual Diagnostics for Hierarchical (Multi-Level / Mixed) Regression Models_. R package version 0.4.7, . Hijmans R (2024). _terra: Spatial Data Analysis_. R package version 1.7-71, https://CRAN.R-project.org/package=terra Ishiyama, N., Sueyoshi, M., García Molinos, J., Iwasaki, K., Negishi, J. N., Koizumi, I., Nagayama, S., Nagasaka, A., Nagasaka, Y., & Nakamura, F. (2023). Underlying geology and climate interactively shape climate change refugia in mountain streams. Ecological Monographs, 93 (2), e1566. https://doi.org/10.1002/ecm.1566 Kelsey, K. C., Pedersen, S. H., Leffler, A. J., Sexton, J. O., Feng, M., & Welker, J. M. (2021). Winter snow and spring temperature have differential effects on vegetation phenology and productivity across Arctic plant communities. Global Change Biology, 27 (8), 1572–1586. https://doi.org/10.1111/gcb.15505 Keppel, G., Stralberg, D., Morelli, T. L., & Bátori, Z. (2024). Managing climate-change refugia to prevent extinctions. Trends in Ecology & Evolution, 39 (9), 800–808. https://doi.org/10.1016/j.tree.2024.05.002 Kloos, S., Klosterhalfen, A., Knohl, A., & Menzel, A. (2024). Decoding autumn phenology: Unraveling the link between observation methods and detected environmental cues. Global Change Biology, 30 (3), e17231. https://doi.org/10.1111/gcb.17231 Klos, P. Z., Link, T. E., & Abatzoglou, J. T. (2014). Extent of the rain-snow transition zone in the western U.S. under historic and projected climate: Climatic rain-snow transition zone. Geophysical Research Letters, 41 (13), 4560–4568. https://doi.org/10.1002/2014GL060500 Lawler, J. J., Ackerly, D. D., Albano, C. M., Anderson, M. G., Dobrowski, S. Z., Gill, J. L., Heller, N. E., Pressey, R. L., Sanderson, E. W., & Weiss, S. B. (2015). The theory behind, and the challenges of, conserving nature’s stage in a time of rapid change. Conservation Biology, 29 (3), 618–629. https://doi.org/10.1111/cobi.12505 Lenth R (2024). _emmeans: Estimated Marginal Means, aka Least-Squares Means_. R package version 1.10.3, . Liston, G. E., & Elder, K. (2022). A distributed snow-evolution modeling system (SnowModel). Journal of Hydrometeorology. 7(6): 1259-1276. https://research.fs.usda.gov/treesearch/26319 Liu, Q., Fu, Y. H., Zhu, Z., Liu, Y., Liu, Z., Huang, M., Janssens, I. A., & Piao, S. (2016). Delayed autumn phenology in the Northern Hemisphere is related to change in both climate and spring phenology. Global Change Biology, 22 (11), 3702–3711. https://doi.org/10.1111/gcb.13311 Ma, Y., & Zhang, Y. (2022). IMPROVED ON SNOW COVER EXTRACTION IN MOUNTAINOUS AREAS BASED ON MULTI-FACTOR NDSI DYNAMIC THRESHOLD. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XLIII-B3-2022, 771–778. XXIV ISPRS Congress “Imaging today, foreseeing tomorrow”, Commission III - 2022 edition, 6–11 June 2022, Nice, France. https://doi.org/10.5194/isprs-archives-XLIII-B3-2022-771-2022 Marshall, A. M., Abatzoglou, J. T., Link, T. E., & Tennant, C. J. (2019). Projected Changes in Interannual Variability of Peak Snowpack Amount and Timing in the Western United States. Geophysical Research Letters, 46 (15), 8882–8892. https://doi.org/10.1029/2019GL083770 McCune, B., & Keon, D. (2002). Equations for Potential Annual Direct Incident Radiation and Heat Load. Journal of Vegetation Science, 13 (4), 603–606. McLaughlin, B. C., Ackerly, D. D., Klos, P. Z., Natali, J., Dawson, T. E., & Thompson, S. E. (2017). Hydrologic refugia, plants, and climate change. Global Change Biology, 23 (8), 2941–2961. https://doi.org/10.1111/gcb.13629 Meng, L., Zhou, Y., Gu, L., Richardson, A. D., Peñuelas, J., Fu, Y., Wang, Y., Asrar, G. R., De Boeck, H. J., Mao, J., Zhang, Y., & Wang, Z. (2021). Photoperiod decelerates the advance of spring phenology of six deciduous tree species under climate warming. Global Change Biology, 27 (12), 2914–2927. https://doi.org/10.1111/gcb.15575 Michalak, J. L., Stralberg, D., Cartwright, J. M., & Lawler, J. J. (2020). Combining physical and species-based approaches improves refugia identification. Frontiers in Ecology and the Environment, 18 (5), 254–260. https://doi.org/10.1002/fee.2207 Millar, C. I., & Westfall, R. D. (2019). Geographic, hydrological, and climatic significance of rock glaciers in the Great Basin, USA. Arctic, Antarctic, and Alpine Research, 51 (1), 232–249. https://doi.org/10.1080/15230430.2019.1618666 Morelli, T. L., Daly, C., Dobrowski, S. Z., Dulen, D. M., Ebersole, J. L., Jackson, S. T., Lundquist, J. D., Millar, C. I., Maher, S. P., Monahan, W. B., Nydick, K. R., Redmond, K. T., Sawyer, S. C., Stock, S., & Beissinger, S. R. (2016). Managing Climate Change Refugia for Climate Adaptation. PLOS ONE, 11 (8), e0159909. https://doi.org/10.1371/journal.pone.0159909 Morelli, T. L., Barrows, C. W., Ramirez, A. R., Cartwright, J. M., Ackerly, D. D., Eaves, T. D., Ebersole, J. L., Krawchuk, M. A., Letcher, B. H., Mahalovich, M. F., Meigs, G. W., Michalak, J. L., Millar, C. I., Quiñones, R. M., Stralberg, D., & Thorne, J. H. (2020). Climate-change refugia: Biodiversity in the slow lane. Frontiers in Ecology and the Environment, 18 (5), 228–234. https://doi.org/10.1002/fee.2189 Moss, W. E., Crausbay, S. D., Rangwala, I., Wason, J. W., Trauernicht, C., Stevens-Rumann, C. S., Sala, A., Rottler, C. M., Pederson, G. T., Miller, B. W., Magness, D. R., Littell, J. S., Frelich, L. E., Frazier, A. G., Davis, K. T., Coop, J. D., Cartwright, J. M., & Booth, R. K. (2024). Drought as an emergent driver of ecological transformation in the twenty-first century. BioScience, biae050. https://doi.org/10.1093/biosci/biae050 National Operational Hydrologic Remote Sensing Center (NOHRSC). (2004). Snow Data Assimilation System (SNODAS) Data Products at NSIDC, Version 1 . https://doi.org/10.7265/N5TB14TC Nolin, A. W. (2010). Recent advances in remote sensing of seasonal snow. Journal of Glaciology, 56 (200), 1141–1150. https://doi.org/10.3189/002214311796406077 Overpeck, J. T., & Udall, B. (2020). Climate change and the aridification of North America. Proceedings of the National Academy of Sciences, 117 (22), 11856–11858. https://doi.org/10.1073/pnas.2006323117 Pastore, M. A., Classen, A. T., D’Amato, A. W., Foster, J. R., & Adair, E. C. (2022). Cold-air pools as microrefugia for ecosystem functions in the face of climate change. Ecology, 103 (8), e3717. https://doi.org/10.1002/ecy.3717 Peven, G., Engels, M., Eitel, J. U. H., & Andrus, R. A. (2024). Montane springs provide regeneration refugia after high-severity wildfire. Ecosphere, 15 (9), e70009. https://doi.org/10.1002/ecs2.70009 Piao, S., Liu, Q., Chen, A., Janssens, I. A., Fu, Y., Dai, J., Liu, L., Lian, X., Shen, M., & Zhu, X. (2019). Plant phenology and global climate change: Current progresses and challenges. Global Change Biology, 25 (6), 1922–1940. https://doi.org/10.1111/gcb.14619 Planet Labs PBC (2024). Planet Application Program Interface: In Space for Life on Earth. San Francisco, CA. https://api.planet.com Potter, C. (2020). Snowmelt timing impacts on growing season phenology in the northern range of Yellowstone National Park estimated from MODIS satellite data. Landscape Ecology, 35 (2), 373–388. https://doi.org/10.1007/s10980-019-00951-3 PRISM Climate Group, Oregon State University, https://prism.oregonstate.edu, data created 4 Feb 2014, accessed 10 Nov 2024. R Core Team. 2024. R: A Language and Environment for Statistical Computing. Vienna: R Foundation for Statistical Computing. Renner, S. S., & Zohner, C. M. (2018). Climate Change and Phenological Mismatch in Trophic Interactions Among Plants, Insects, and Vertebrates. Annual Review of Ecology, Evolution, and Systematics, 49 (Volume 49, 2018), 165–182. https://doi.org/10.1146/annurev-ecolsys-110617-062535 Siirila-Woodburn, E. R., Rhoades, A. M., Hatchett, B. J., Huning, L. S., Szinai, J., Tague, C., Nico, P. S., Feldman, D. R., Jones, A. D., Collins, W. D., & Kaatz, L. (2021). A low-to-no snow future and its impacts on water resources in the western United States. Nature Reviews Earth & Environment, 2 (11), 800–819. https://doi.org/10.1038/s43017-021-00219-y Slatyer, R. A., Umbers, K. D. L., & Arnold, P. A. (2022). Ecological responses to variation in seasonal snow cover. Conservation Biology, 36 (1), e13727. https://doi.org/10.1111/cobi.13727 Springer, A. E., Stevens, L. E., Ledbetter, J. D., Schaller, E. M., Gill, K. M., & Rood, S. B. (2015). Ecohydrology and stewardship of Alberta springs ecosystems. Ecohydrology, 8 (5), 896–910. https://doi.org/10.1002/eco.1596 Stevens, L. E., Schenk, E. R., & Springer, A. E. (2021). Springs ecosystem classification. Ecological Applications, 31 (1). https://doi.org/10.1002/eap.2218 Stewart, D. E., R. S. Lewis, E. D. Stewart, and P. K. Link (2013). Geologic Map of the Central and Lower Big Creek Drainage, Central Idaho: Idaho Geological Survey Digital Web Map, Scale 1:75,000. https://www.idahogeology.org/product/dwm-161. Tang, J., Körner, C., Muraoka, H., Piao, S., Shen, M., Thackeray, S. J., & Yang, X. (2016). Emerging opportunities and challenges in phenology: A review. Ecosphere, 7 (8), e01436. https://doi.org/10.1002/ecs2.1436 Tercek, M. T., Gross, J. E., & Thoma, D. P. (2023). Robust projections and consequences of an expanding bimodal growing season in the western United States. Ecosphere, 14 (5), e4530. https://doi.org/10.1002/ecs2.4530 Thoma, D. P., Munson, S. M., & Witwicki, D. L. (2019). Landscape pivot points and responses to water balance in national parks of the southwest US. Journal of Applied Ecology, 56 (1), 157–167. https://doi.org/10.1111/1365-2664.13250 Thornthwaite, C. W. (1948). An Approach toward a Rational Classification of Climate. Geographical Review, 38, 55-94. https://doi.org/10.2307/210739 Tsinnajinnie, L. M., Frisbee, M. D., & Wilson, J. L. (2021). Groundwater from perennial springs provide refuge from wildfire impacts in mountainous semiarid watershed. Journal of Hydrology, 596, 125701. https://doi.org/10.1016/j.jhydrol.2020.125701 U.S. Environmental Protection Agency, 2013, Level III and IV ecoregions of the continental United States: Corvallis, Oregon, U.S. EPA, National Health and Environmental Effects Research Laboratory, map scale 1:3,000,000, https://www.epa.gov/eco-research/level-iii-and-iv-ecoregions-continental-united-states. U.S. Geological Survey, 2024, USGS 3D Elevation Program Digital Elevation Model, accessed November 10, 2024 at URL https://elevation.nationalmap.gov/arcgis/rest/services/3DEPElevation/ImageServer. Weissinger, R., Philippi, T. E., & Thoma, D. (2016). Linking climate to changing discharge at springs in Arches National Park, Utah, USA. Ecosphere, 7 (10), e01491. https://doi.org/10.1002/ecs2.1491 White, J. G., Sparrius, J., Robinson, T., Hale, S., Lupone, L., Healey, T., Cooke, R., & Rendall, A. R. (2022). Can NDVI identify drought refugia for mammals and birds in mesic landscapes? Science of The Total Environment, 851, 158318. https://doi.org/10.1016/j.scitotenv.2022.158318 White, D. C., & Lewis, M. M. (2011). A new approach to monitoring spatial distribution and dynamics of wetlands and associated flows of Australian Great Artesian Basin springs using QuickBird satellite imagery. Journal of Hydrology, 408 (1–2), 140–152. https://doi.org/10.1016/j.jhydrol.2011.07.032 Whiting, J. A., & Godsey, S. E. (2016). Discontinuous headwater stream networks with stable flowheads, Salmon River basin, Idaho. Hydrological Processes, 30 (13), 2305–2316. https://doi.org/10.1002/hyp.10790 Winkler, D. E., Butz, R. J., Germino, M. J., Reinhardt, K., & Kueppers, L. M. (2018). Snowmelt Timing Regulates Community Composition, Phenology, and Physiological Performance of Alpine Plants. Frontiers in Plant Science, 9 . https://doi.org/10.3389/fpls.2018.01140 Wood, S.N. (2011) Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society (B) 73(1):3-36 Wu, C., Peng, J., Ciais, P., Peñuelas, J., Wang, H., Beguería, S., Andrew Black, T., Jassal, R. S., Zhang, X., Yuan, W., Liang, E., Wang, X., Hua, H., Liu, R., Ju, W., Fu, Y. H., & Ge, Q. (2022). Increased drought effects on the phenology of autumn leaf senescence. Nature Climate Change, 12 (10), 943–949. https://doi.org/10.1038/s41558-022-01464-9 Xie, Y., Wang, X., Wilson, A. M., & Silander, J. A. (2018). Predicting autumn phenology: How deciduous tree species respond to weather stressors. Agricultural and Forest Meteorology, 250–251, 127–137. https://doi.org/10.1016/j.agrformet.2017.12.259 Zani, D., Crowther, T. W., Mo, L., Renner, S. S., & Zohner, C. M. (2020). Increased growing-season productivity drives earlier autumn leaf senescence in temperate trees. Science (New York, N.Y.), 370 (6520), 1066–1071. https://doi.org/10.1126/science.abd8911 Zeng, Y., Hao, D., Huete, A., Dechant, B., Berry, J., Chen, J. M., Joiner, J., Frankenberg, C., Bond-Lamberty, B., Ryu, Y., Xiao, J., Asrar, G. R., & Chen, M. (2022). Optical vegetation indices for monitoring terrestrial ecosystems globally. Nature Reviews Earth & Environment, 3 (7), 477–493. https://doi.org/10.1038/s43017-022-00298-5 Zeng, L., Wardlow, B. D., Xiang, D., Hu, S., & Li, D. (2020). A review of vegetation phenological metrics extraction using time-series, multispectral satellite data. Remote Sensing of Environment, 237, 111511. https://doi.org/10.1016/j.rse.2019.111511 Zhang, Y., Parazoo, N. C., Williams, A. P., Zhou, S., & Gentine, P. (2020). Large and projected strengthening moisture limitation on end-of-season photosynthesis. Proceedings of the National Academy of Sciences, 117 (17), 9216–9222. https://doi.org/10.1073/pnas.1914436117 Zhou, H., Min, X., Chen, J., Lu, C., Huang, Y., Zhang, Z., & Liu, H. (2023). Climate warming interacts with other global change drivers to influence plant phenology: A meta-analysis of experimental studies. Ecology Letters, 26 (8), 1370–1381. https://doi.org/10.1111/ele.14259 Information & Authors Information Version history Peer review timeline Published Ecohydrology Version of Record3 Jul 2025Published Copyright This work is licensed under a Non Exclusive No Reuse License. Collection

Authors Metrics & Citations Metrics Article Usage 308views 217downloads Citations Download citation Grace Peven, Jan U.H. Eitel, Timothy Link, et al. The role of spring ecosystems as climate refugia in a semi-arid environment. Authorea. 18 April 2025. DOI: https://doi.org/10.22541/au.174495449.90011594/v1 DOI: https://doi.org/10.22541/au.174495449.90011594/v1 If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download. For more information or tips please see 'Downloading to a citation manager' in the Help menu.