High Resolution SnowModel Simulations Reveal Future Elevation‐Dependent Snow Loss and Earlier, Flashier Surface Water Input for the Upper Colorado River Basin

Continued climate warming is reducing seasonal snowpacks in the western United States, where >50% of historical water supplies were snowmelt‐derived. In the Upper Colorado River Basin, declining snow water equivalent (SWE) and altered surface water input (SWI, rainfall and snowmelt available to enter the soil) timing and magnitude affect streamflow generation and water availability. To adapt effectively to future conditions, we need to understand current spatiotemporal distributions of SWE and SWI and how they may change in future decades. We developed 100‐m SnowModel simulations for water years 2001–2013 and two scenarios: control (CTL) and pseudo‐global‐warming (PGW). The PGW fraction of precipitation falling as snow was lower relative to CTL, except for November–April at high elevations. PGW peak SWE was lower for low (−45%) and mid elevations (−14%), while the date of peak SWE was uniformly earlier in the year for all elevations (17–23 days). Currently unmonitored high elevation snow represented a greater fraction of total PGW SWE. PGW peak daily SWI was higher for all elevations (30%–42%), while the dates of SWI peaks and centroids were earlier in the year for all elevations under PGW. PGW displayed elevated winter SWI, lower summer SWI, and changes in spring SWI timing were elevation‐dependent. Although PGW peak SWI was elevated and earlier compared to CTL, SWI was more evenly distributed throughout the year for PGW. These simulated shifts in the timing and magnitude of SWE and SWI have broad implications for water management in dry, snow‐dominated regions.

• Projections show lower peak snow water equivalent (SWE) below 3,000 m and earlier peak SWE, peak surface water input (SWI) at all elevations • Greater future peak SWI and reduced annual snow-derived SWI for all elevations, with a more even SWI distribution throughout the year • A greater fraction of future SWE will be in high elevations that are currently unmonitored

Supporting Information:
Supporting Information may be found in the online version of this article.
supply, irrigation, mining, hydropower, agriculture, and recreation (e.g., Huss et al., 2017). Continued global warming is reducing seasonal snowpacks in the western United States (Mote et al., 2018), where a majority of the region's water supply has historically been generated by snowmelt, especially in semi-arid basins (Li et al., 2017). These changes are modifying the timing and magnitude of streamflow, with (a) mid-winter melt occurring more frequently, (b) declines in peak snow water equivalent (SWE) (Musselman et al., 2021), and (c) earlier streamflow timing in snowmelt driven areas (Dudley et al., 2017;Stewart et al., 2005). Across the western United States, SWE declines of ∼25% are expected by the middle of the 21st century, and seasonal snowpacks are likely to be replaced by low-to-no snow conditions at low and mid elevations (Fyfe et al., 2017;Siirila-Woodburn et al., 2021).
Snowmelt in the Upper Colorado River Basin (UCRB) generates as much as 92% of the streamflow for the larger Colorado River Basin, which supports over 40 million people and irrigates 5.5 million acres across seven western U.S. states (Butler et al., 2015;Lukas & Payton, 2020). Thus, better understanding the present spatiotemporal distribution of SWE and the associated hydrologic response to global warming in this region is critical to prepare for future adjustments to water resource management. Existing allocations between upper and lower basin states may be unrealizable with continued snow loss and heightened interannual climatic variability, especially for persistent drought periods (Adler, 2008;Udall & Overpeck, 2017;Woodhouse et al., 2006). Prior work in the Colorado River Basin has shown annual streamflow to be most sensitive to warming during the warm season as compared to other seasons for the last several decades (Ban & Lettenmaier, 2022), with annual discharge decreasing 9.3% per °C of warming over the period 1913-2017 (Milly & Dunne, 2020). Other recent studies report that increases of temperature of 1-2°C in the Colorado River Basin without changes in precipitation could result in streamflow reductions of 10%-20% (Christensen et al., 2004;Christensen & Lettenmaier, 2007;Wolock, 2007, 2011). Projections using regional climate models indicate reduced, earlier peak SWE in the UCRB (Rasmussen, Baker et al., 2011) and increases in streamflow for the 2030s followed by decreases for the 2080s compared to present, with potential streamflow increases insufficient to meet growing water demand (Miller et al., 2021).
Accurately characterizing mountain precipitation and SWE remains a key task for responding to projected streamflow reductions and intensified droughts (Kim et al., 2021). This characterization has been a challenge because most global climate models (GCMs) or mesoscale climate models produce simulations at coarse spatial resolutions and are unable to represent complex orographic precipitation patterns and small-scale snow processes that are important in mountainous regions such as the UCRB, resulting in substantial biases (Salzmann & Mearns, 2012). Fine spatial scale is crucial for regional predictions of future hydroclimate (Ikeda et al., 2021;Wrzesien & Pavelsky, 2020) to accurately represent the influence of topographic complexity on orographic precipitation patterns and snow process dynamics (Liston, 2004;López-Moreno et al., 2013;Musselman, Clark, et al., 2017, Musselmann, Molotch, & Margulis et al., 2017Rasmussen, Baker et al., 2011). In this study, we extend prior mesoscale snow modeling work for areas of the western United States (e.g., Ikeda et al., 2021;Li et al., 2017;Musselman, Clark, et al., 2017) to a finer resolution by using process-based snow modeling (Snow-Model; Liston & Elder, 2006a) at 100-m (m) spatial resolution driven with fine-resolution convective-permitting and orography-resolving Weather Research and Forecasting Model (WRF) simulations for a control (CTL) and a pseudo-global-warming (PGW) scenario. The PGW method Rasmussen, Baker et al., 2011, Rasmussen, Liu et al., 2011Sato et al., 2007;Schär et al., 1996) is best used to ask the question, "what will today's weather look like in a future climate state." The difference between two 30-year monthly mean conditions in a future (2071-2100) and current  are applied to reanalysis data used to drive the mesoscale model that allows for historical events and patterns (i.e., extreme storm events, dry vs. wet seasons) to be directly studied in contrast to other methods that dynamically downscale from GCMs. The approach well characterizes the thermodynamic aspects of climate change, but does not represent large-scale changes in atmospheric circulation. However, sub-monthly mesoscale variability is well captured and allows for significant improvements in the representation of precipitating systems in regions of complex terrain (Rasmussen, Baker et al., 2011. This study optimally combines two cutting-edge models at a fine spatial scale that better captures the effects of complex topography and elevation gradients on snow variability and snow processes than prior studies at 4-km or coarser resolutions, and it allows us to investigate how snowpack processes change with climate warming without introducing uncertainty inherent to other methods that include potential changes in storm tracks. Additionally, in this study we analyze surface water input (SWI), which is the sum of rainfall and snowmelt 10.1029/2022EF003092 3 of 23 available to enter the soil, to characterize the effects of a changing snowpack on the hydrologic cycle. Prior studies at plot and catchment scale have demonstrated linkages between SWI and observed and modeled streamflow (Hammond et al., 2019;Kiewiet et al., 2022;Kormos et al., 2014), and the use of SWI in this study enables linking snow changes to hydrologic effects not fully captured in prior studies using metrics including snowfall fraction (Berghuijs et al., 2014) and snowmelt rate (Barnhart et al., 2016).
With the aim to reduce uncertainty in present and future spatiotemporal distributions of SWE and SWI, we ask our primary study question (a) how will SWE change under a future warmer climate by elevation, and how might this affect future water availability via changes to the timing and magnitude of SWI? (Sections 3.3, 3.4, 4.1).
Using PGW as an indicator of potential future conditions associated with the thermodynamic aspects of climate change, we hypothesize future SWI will occur earlier and at lower magnitude for mid to low elevations coincident with decreases in SWE (Barnhart et al., 2016;Musselmann, Molotch, & Margulis et al., 2017), but peak SWE may increase at times for high elevations where precipitation is expected to increase; temperatures during late fall, winter and early spring will still be conducive for snowfall (e.g., Musselman, Clark, et al., 2017). Despite similar or elevated SWE at high elevations, we hypothesize the timing of future SWI will be earlier in the year than for historical estimates (Musselmann, Molotch, & Margulis et al., 2017). To place these results in context of current snow monitoring locations following in depth analysis of the present and future distributions of SWE and SWI, we ask a secondary study question (b) how well do current snow observation stations represent the spatial distribution of SWE and SWI, and how might these networks need to change to continue to provide predictive information for future water availability? (Section 4.2). We hypothesize future snowpacks will be more limited to higher elevations than at present (Ikeda et al., 2021), and that low to mid elevations will provide a smaller snowmelt contribution to streamflow despite substantial present-day snowmelt contribution given their large spatial footprint (Hammond et al., 2018;Harrison et al., 2021).

Methods
To answer our study questions, we selected the Colorado portion of the UCRB (Section 2.1), implemented CTL and PGW SnowModel simulations (Section 2.2), and performed a multi-source model evaluation using station and remotely sensed data (Section 2.3). We then aggregated model inputs and outputs to elevation bands and daily, monthly, and annual timesteps, and calculated several annual metrics to describe differences in SWE and SWI between CTL and PGW scenarios (Section 2.4).

Study Area
The model domain was a 311 × 300 km area in Colorado, United States, encompassing the headwaters of Colorado and Gunnison River Basin portions of the UCRB (Figure 1). Elevations within the domain ranged from 1,313 to 4,385 m ( Figure 1). The model domain is a snow-dominated mountainous region with deep persistent snowpacks and high runoff generation ( Figure 1; Hammond et al., 2018;Miller et al., 2016); seasonal snowpacks in the region serve as the primary source of annual streamflow with snowmelt contributing to an estimated 71% of total UCRB runoff (Li et al., 2017).

Model Description and Simulations
To simulate detailed snow evolution processes for both current and future climate conditions, we utilized Snow-Model, a spatially distributed physics-based snow evolution modeling system (Liston & Elder, 2006a;Liston et al., 2020), forced with WRF convection-permitting and orography-resolving regional climate simulations on a 4-km grid resolution . WRF climate simulations provided the atmospheric forcing conditions to drive SnowModel in both a current and future climate scenario. A pair of continuous 13-water-year (2001-2013) WRF model simulations was used: (a) a current climate simulation CTL forced using ERA-Interim reanalysis, and (b) a future climate simulation using the PGW method. The PGW simulation uses the ERA-Interim reanalysis for the same period as (a) and adds an ensemble mean climate delta (increase in radiative forcing) from 100 years in the future for the most extreme Representative Concentration Pathway (RCP) 8.5 scenario. The two-dimensional 4-km grid spacing hourly WRF simulation output  for both CTL and PGW scenarios were subset to our study area model domain and used as atmospheric forcing for running SnowModel.
10.1029/2022EF003092 4 of 23 SnowModel simulations were run for both CTL and PGW scenarios with a three-hour time step and 100-m spatial grid resolution (9.32 million grid cells across the model domain; Figure 1; Sexstone et al., 2022). Elevation and land cover data (30-m spatial resolution) were provided by the U.S. Geological Survey (USGS) National Elevation Data set (U.S. Geological Survey, 2020) and North American Land Change Monitoring System (Canada Centre for Remote Sensing, 2020) and were resampled to the 100-m spatial grid resolution using average resampling and mode resampling (GDAL/OGR, 2022), respectively. Hourly WRF precipitation, air temperature, relative humidity, and wind speed and direction data were aggregated to three-hour values to correspond with the model simulation time step. Three-hour WRF forcing data along with the mean elevation of each 4-km WRF grid cell were used by MicroMet (Liston & Elder, 2006b), a high-resolution meteorological distribution submodel of SnowModel, to downscale and create the 100-m spatial resolution meteorological forcing data required to run SnowModel. The Dai (2008) rain-snow precipitation phase fraction parameterization  was used to calculate a rain-snow fraction of precipitation inputs for this study. Other submodels that make up the SnowModel modeling system were used to simulate spatially distributed snowpack evolution across the study area including: EnBal (Liston, 1995), which computes surface energy exchanges between the snow and atmosphere; SnowPack (Liston & Hall, 1995), which simulates the seasonal evolution of snow depth and SWE; and SnowTran-3D (Liston & Sturm, 1998), a three-dimension model that simulates snow redistribution by wind over topographically variable terrain. CTL and PGW SnowModel simulations include snow sublimation from the surface, canopy, and blowing snow (Sexstone et al., 2018). SnowModel simulation outputs (100-m spatial resolution) of SWE, SWI, snowmelt, precipitation, fraction of precipitation falling as snow, and air temperature were aggregated from a three-hour to a daily time step and used in this study.

Model Evaluation
SnowModel simulations for the CTL scenario were evaluated using two independent data streams, snow telemetry (SNOTEL) station observations and Moderate Resolution Imaging Spectroradiometer (MODIS) satellite observations. Simulated precipitation and SWE values were evaluated using SNOTEL station observations of precipitation and SWE at 71 sites ranging in elevation from 2,591 to 3,523 m that included at least three complete water years of observations during the study period (Table S1 in Supporting Information S1). Simulated snow-covered Black and orange contours in (a, b) and the vertical lines of the same color in the insets of (a and b) separate the domain into three elevational categories that are used for synthesizing patterns throughout our analyses: "high," >3,300 m; "mid," 2,300 m < elevation ≤ 3,300 m; "low," 2,300 m. area and snow persistence (SP) were evaluated using daily cloud gap filled MODIS observations, and derived SP grids. A grid cell and day in the model output was considered to be snow covered if SWE ≥10 mm; and SP was calculated as the fraction of time snow was present from 1 January to 3 July. Evaluation using MODIS data was completed for all areas of the domain in the seasonal snow zone, where SP is greater than 50% ( Figure S3 in Supporting Information S1). Prior work has evaluated the WRF forcing ) used in our model simulations, and the model estimates reliably reproduced the annual, seasonal, and sub-seasonal precipitation and surface temperature climatology in mountain regions of the western United States Rasmussen, Baker et al., 2011). Likewise, SnowModel has been rigorously evaluated and shown to perform well in seasonally snow-covered environments similar to the study area (e.g., Greene et al., 1999;Hiemstra et al., 2006;Liston & Elder, 2006a, 2006bListon & Hiemstra, 2008;Prasad et al., 2001;Sproles et al., 2013;Sexstone et al., 2018;Sexstone, Penn et al., 2020). For details on the methods used to evaluate the SnowModel simulations used in this study, please refer to Text S1 in Supporting Information S1.

Analysis of Model Inputs and Outputs
We examined the precipitation, temperature, and fraction of precipitation falling as snow on monthly, annual, and mean annual timescales across 13 years of simulation. We further examined SWE and SWI at daily resolution. We focus many of our analyses on variability in mean annual SWE and SWI by elevation, as prior studies have used coarser spatial resolutions less able to resolve mountain topography and elevational snow variability (Ikeda et al., 2021;Musselmann, Molotch, & Margulis et al., 2017). Because the 13-year window is too short for evaluating temporal variability and trends, we evaluated 13-year average differences between CTL and PGW scenarios. We aggregated daily SWE and SWI to 100-m elevation bands ( Figure S1 in Supporting Information S1) using the exactextractr (Baston, 2020) R package to obtain daily area-weighted averages for each unit. Spatial units across the domain were classified into 31,100-m elevation bands ranging from 1,300 to 4,300 m. For mean annual scale analyses, the mean of annual totals (e.g., precipitation) or the mean of annual means (e.g., temperature) were computed. For analyses plotting results by elevation, we computed the median value of a given metric for all spatial units falling within each 100-m elevation band. For example, when displaying the CTL peak SWE for 3,000 m in Figure 3, we first calculated the peak SWE for each simulation year (2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) and each spatial unit, then took the mean of all simulation years for that unit, and finally took the median value of peak SWE for all spatial units with the 100-m elevation classification of 3,000 m.
We used a threshold in rainfall/SWE to differentiate between snowmelt dominated input and mixed input because without a threshold, snowmelt events with an insignificant amount of rain would be categorized as mixed. A 10% threshold in rainfall/SWE, where input with rainfall/SWE >= 10% is mixed input and rainfall/SWE < 10% is snowmelt input, was adopted here and is functionally similar to those used in other rain on snow studies (Cohen et al., 2015;Guan et al., 2016;Lundquist et al., 2008). We calculated six annual SWE and SWI metrics for each spatial unit ( Figure S1 in Supporting Information S1) using median daily values of model outputs as described above. Our focus is on SWI, which includes rainfall and snowmelt, because it is more indicative of hydrologic response in snow dominated areas than precipitation alone. Annual metrics calculated were: (a) peak SWI, which represents peak daily SWI, (b) peak SWE, (c) peak SWI DOWY, the day of the water year on which the peak SWI occurs, (d) peak SWE DOWY, the day of the water year when peak SWE occurs, (e) SWI50 DOWY, the day of the water year on which 50% of the annual SWI has occurred, and (f) input seasonality (IS), which is similar to the precipitation seasonality metric (Dingman, 2002;Markham, 1970) but calculated with SWI time series rather than precipitation . IS is an index of how much mean annual input occurs within a concentrated period or is equally distributed throughout the year and is calculated using monthly SWI totals. where: A is the annual sum of SWI.
M i is the monthly sum of SWI.
Finally, we computed the fraction of SWI sourced from snowmelt (f snow ) and mixed input (f mix ) for each day using summed total values of SWI for each spatial unit (rather than mean values) separated into three categories: snowmelt input (rainfall/SWE ≤ 0.1 and SWE > 0), mixed input (rainfall/SWE > 0.1 and SWE > 0), and rainfall input (rainfall > 0 and SWE = 0), and these values were then totaled for each spatial unit to calculate f snow and f mix as: In addition to assessing SWI changes by elevation between CTL and PGW scenarios using annual scale metrics, we implemented two subsequent analyses to examine domain-wide impacts to total SWI. To better understand patterns and changes in total study area mean annual SWI by elevation and time of the year, we calculated the fraction of total study area mean annual SWI generated by each elevation band for each month for CTL and PGW separately. To investigate snowmelt contributions to SWI, and how they differ between PGW and CTL, we mapped the mean annual fraction of SWI contributed from snowmelt and calculated the mean monthly snowmelt contributed SWI for both scenarios, and how this changed between scenarios. Finally, we investigated differences between mean annual change and annual change within the same scenario and between scenarios for two simulation years (2011, greatest peak SWE for CTL, and 2012, lowest peak SWE for CTL) that highlight high interannual variability of snow and SWI.
In reporting the results of our analyses, we use three elevation categories to synthesize patterns in SWE and SWI change by elevation: "high," >3,300 m; "mid," 2,300 m < elevation ≤ 3,300 m; "low," ≤2,300 m. These elevation categories were chosen because they broadly reflect the satellite-derived boundaries of the present intermittent and seasonal snow zones in the study region, zones that reflect similar snow patterns (Moore et al., 2015) and similar patterns of soil moisture and streamflow response to snowmelt (Hammond et al., 2018;Harrison et al., 2021). These elevation categories are also consistent with regional regression equation elevation ranges for streamflow prediction in ungaged areas (Capesius & Stephens, 2009;Eurich et al., 2021).

Results
To examine differences in SWE and SWI between CTL and PGW scenarios, we first evaluated performance of SnowModel in the Colorado portion of the UCRB for the CTL scenario (Section 3.1) and assessed differences in the meteorological forcing between CTL and PGW (Section 3.2). We then quantified differences in SWE (Section 3.3) and SWI (Section 3.4) between CTL and PGW using daily model output and annual-scale synthesis metrics for elevation bands across the study area.

Model Evaluation
We evaluated simulated SnowModel SWE for the historic (water years 2001-2013) CTL scenario using SNOTEL stations observations ( Figure S2, Table S1 in Supporting Information S1) and simulated SnowModel snow covered area using MODIS observed snow covered area for water years 2001-2005. Comparison of simulated SWE with all available observations from 71 SNOTEL stations within the study domain showed a root mean square error (RMSE) of 98 mm and Nash-Sutcliffe efficiency (NSE) of 0.72 for periods when both observed and simulated SWE were not equal to zero. NSE performance statistics computed individually for simulation results at each SNOTEL station (Table S1 in Supporting Information S1) show 30% of stations displayed a very good performance rating (0.75 < NSE < 1.00), 25% of stations had good performance (0.65 < NSE < 0.75), 17% of stations were satisfactory (0.50 < NSE < 0.65), and 28% of stations were unsatisfactory (NSE < 0.50) according to criteria established in Moriasi et al. (2007). Simulated SWE at SNOTEL stations was biased low overall (PBias −6.7%, mean error −14 mm) compared to observations, following a low bias in simulated precipitation (i.e., WRF precipitation downscaled to a 100 m spatial resolution by SnowModel) relative to observations (PBias −6.5%, mean error −27 mm). Simulated precipitation exhibited an RMSE of 97 mm and NSE of 0.88 when compared to SNOTEL observations. The SnowModel snow covered area daily time series by elevation band showed agreement with remotely sensed values with a median study area NSE of 0.33 and PBias of 1.9% for areas with seasonal snowpack (SP > 50%). SnowModel differed most from remotely sensed estimates during abrupt accumulation and melt events, but showed close alignment for accumulation and melt periods ( Figure S3a-d in Supporting Information S1). The NSE of modeled and remotely sensed snow-covered area values for 100-m elevation band units across the study area, performance was generally best for mid and high elevation areas in the seasonal snow zone ( Figure S3e in Supporting Information S1), with poorer performance in low elevation, intermittently snowy areas. Finally, in comparison to existing SWE products at 4-km and 1-km spatial resolutions, the 100-m simulations for current climate conditions used in this paper better matched the spatial variability of near peak SWE for the UCRB (Figures S17, S18 in Supporting Information S1). These results suggest the forcing data and level of performance of SnowModel is adequate for our analyses tracking snow accumulation and melt.

Differences in Meteorological Forcings Between Scenarios
The PGW scenario was warmer and wetter than the CTL scenario and exhibited reductions to the fraction of precipitation falling as snow. Mean monthly PGW temperatures for all months ranged from 4.9 to 5.6°C warmer than for CTL for elevation bands across the domain. We observed the greatest positive changes in temperature in the eastern, high elevation region of the study domain, with the greatest changes in fall (Figure 2b), followed by spring. In these seasons, high elevations exhibited the greatest positive changes, with changes in October up to 6.3°C for high elevations. The smallest long-term changes occurred during winter on the order of 4°C.
Following the thermodynamic relation between temperature and specific humidity, mean monthly PGW precipitation tended to be higher, by approximately 146 mm/y relative to CTL across all elevations. The greatest precipitation changes occurred for high elevations in the central and northeast parts of the study domain. Precipitation was higher compared with CTL in all seasons except spring (Figure 2a). The greatest positive precipitation changes occurred during winter at high elevations, resulting in increased snowfall, and the greatest negative changes occurred during spring months at high elevations.
The fraction of precipitation delivered as snow was reduced by as much as 40% between CTL and PGW, with every 100-m elevation band indicating lower snowfall fraction during at least one month (Figure 2c). In cold season months between November and April, high elevations indicated little change in the fraction of precipitation delivered as snow, whereas the low elevations indicated losses of up to 40%, especially in December through February. High elevations indicated the greatest losses in proportion of precipitation delivered as snow during May and October, likely owing to positive changes in temperature projected during those months (e.g., Klos et al., 2014;Scalzitti et al., 2016aScalzitti et al., , 2016b.   Figure S9 in Supporting Information S1 contains annual SWI metrics calculated using only snowmelt-derived SWI. Transparent orange, yellow and blue coloring in the background of each panel reflects three elevation categories used to synthesize patterns in snow and SWI changes by elevation: high, >3,300 m; mid, 2,300 m < elevation ≤ 3,300 m; low, ≤2,300 m.

High Elevations
We evaluated changes in the timing and magnitude of SWE between the CTL and PGW scenarios using annual metrics that synthesize relevant shifts in daily time series of SWE. The greatest mean annual peak SWE for both scenarios occurred at the highest elevations ( Figure 3), with the PGW scenario exhibiting, on average, peak SWE values that were only 3% lower (13 mm) relative to CTL. Despite higher mean annual precipitation for PGW at all elevations during most of the snow accumulation season (Figure 2), peak SWE was only higher over a small, high elevation area in the PGW scenario relative to CTL (Figure 3, Figure S13 in Supporting Information S1). The date of peak SWE was latest at the high elevations for CTL and earlier at high elevations in the PGW scenario (−17 days, Table S4 in Supporting Information S1). Additional comparisons using mean daily values ( Figure 4) and mean monthly values (Table S2 in Supporting Information S1) confirmed SWE magnitude and timing changes observed using annual metrics. The months where high-elevation SWE differed most between the scenarios were May-July (−44 to −191 mm, Table S2 in Supporting Information S1) and October-December (−21 to −34 mm). On average, during these months, the PGW scenario exhibited lower SWE than the CTL scenario. High elevations showed elevated February to April SWE from CTL to PGW (+1 to +13 mm, Table  S2 in Supporting Information S1), with patterns in average monthly snowmelt matching the patterns observed for SWE change between CTL and PGW. Notably, while snow did not melt at high elevations during the CTL scenario for the months December to February, small amounts of snowmelt were generated for these months in the PGW scenario (Table S2 in Supporting Information S1).

Mid Elevations
When PGW was compared with CTL, pronounced reductions in peak SWE occurred at mid elevations (−26 mm, −14%, Table S4 in Supporting Information S1, Figure 3). The date of peak SWE was earlier during the water year (−23 days, Table S4 in Supporting Information S1), a larger temporal shift than observed for high or low elevations. Changes in mean monthly values revealed that mid elevations displayed the largest reductions in SWE for April (−53 mm, Figure S15 in Supporting Information S1) and May (−51 mm, Table S2 in Supporting Information S1), and exhibited no change for the months July-September when snow was absent in both scenarios. Monthly mean snowmelt was higher for the mid elevations during winter months, as earlier peak SWE corresponded with greater winter melt and reduced spring and summer snowmelt (Table S2, Figure S8 in Supporting Information S1).

Low Elevations
Low elevations displayed the greatest percentage reductions in peak SWE for the PGW scenario compared to CTL (−30 mm, −45%, Table S4 in Supporting Information S1) and when comparing low elevation reductions to those at mid and high elevations. The date of peak SWE occurred 21 days earlier, a similar advancement to the earlier peak SWE observed for mid elevations. Low elevations displayed reductions in mean monthly SWE (−60% to −100%, Table S2 in Supporting Information S1) for all months except for July and August which lacked SWE in both scenarios. These reductions were further highlighted by patterns displayed using mean daily SWE values for the two scenarios ( Figure 4). Although all elevations experienced earlier peak SWE with greater snowmelt earlier in the year, low elevations were unique in displaying a nearly complete loss of SWE and subsequently snowmelt in the PGW scenario ( Figure 4, Figure S8, Table S2 in Supporting Information S1).

High Elevations
Because we observed differences in SWE and snowmelt between the two scenarios, we then analyzed differences in SWI between CTL and PGW scenarios by elevation. Mean annual peak daily SWI was highest in the high elevation band for CTL (Figure 3), with higher peak daily SWI in the PGW scenario compared to CTL (+30%, +12 mm/d, Table S4 in Supporting Information S1). The date of peak daily SWI for high elevations occurred 14 days earlier in the PGW scenario. The IS metric, which indicates whether input is focused primarily within a few months or is evenly distributed throughout the year, showed that the most temporally concentrated input was at high elevations during CTL. That input became less concentrated for PGW at the same elevations (−13%, Table S4 in Supporting Information S1). For the CTL scenario, the mean annual centroid date of SWI was earlier in the western part of the domain, and later at high elevations and in the eastern part of the study area. The centroid date of mean annual SWI was earlier for PGW than for CTL for all elevations (Table S4, Figure S14 in Supporting Information S1).
Mean monthly SWI was higher in the PGW scenario than the CTL scenario for all months except for June and July (Table S3 in Supporting Information S1). The magnitude of the difference was greatest in April and May, with the negative differences in June and July nearly perfectly offsetting the positive differences in April and May from a mass balance perspective. Changes observed in mean monthly SWI were further evidenced by patterns in mean daily values. For CTL, high elevations experienced long and continuous periods of elevated snowmelt-derived . Mean annual daily snow water equivalent (SWE, top) and mean annual daily surface water input (SWI, bottom) by scenario (control -CTL, pseudo-globalwarming -PGW) and elevation band. A 7-day moving mean was applied for better visualization between bands, but this compresses peak daily SWI. SWE and SWI patterns are separated into three elevational categories: high, >3,300 m; mid, 2,300 m < elevation ≤ 3,300 m; low, ≤2,300 m. Figure S7 in Supporting Information S1 displays greater details in daily SWE and SWI for the low elevation category.
SWI in May and June (Figures 4 and 5), whereas PGW SWI at high elevations displayed higher, concentrated peaks earlier in the year. From here forward, we use the term "flashier" to represent more rapidly occurring and higher peak SWI observed for PGW as compared to CTL. At high elevations, snowmelt contributions to SWI for CTL were greatest in April and May ( Figure 5), with PGW snowmelt contributions to SWI generally higher earlier in the spring and lower later in spring and summer.

Mid Elevations
Peak daily SWI was similarly greater for PGW at mid elevations relative to CTL as compared to high elevations (+34%, +10 mm/d, Table S4 in Supporting Information S1). However, the date of peak daily SWI was 20 days earlier in the year at mid elevations, a greater advance than observed for high elevations. IS declined the most at mid and low elevations (−18%, Table S4 in Supporting Information S1). Mean monthly total SWI values for mid elevations revealed higher values in monthly SWI for October-April, similar to what was observed for high elevations, but with reductions observed for May, June, July, and September. Elevated SWI for March and April in mid elevations was offset by reductions to May and June SWI. Assessing differences in CTL versus PGW April SWI, mid elevations displayed the greatest reductions ( Figure S15 in Supporting Information S1). Similar to the PGW-CTL differences observed in mean daily values for high elevations, mid elevations displayed elevated daily SWI ( Figure 4) and earlier daily snowmelt contributions to SWI ( Figure 5). This period was followed by reductions in SWI and snowmelt contributions, which occurred earlier in the year than at high elevations. While peak daily SWI for high elevations consistently displayed higher values for PGW, some mid elevation bands displayed lower and longer elevated periods potentially due to higher SWI earlier in the season (Figures 4 and 5).

Low Elevations
At low elevations, peak daily SWI was higher for PGW, and with a greater percent difference than for mid and high elevations (+42%, +10 mm/d, Table S4 in Supporting Information S1). The date of peak daily SWI occurred 23 days earlier in the year and the centroid date of SWI occurred 30 days earlier in the year. Both metrics indicate greater differences between scenarios than at mid or high elevations. At low elevations, the difference in mean monthly total SWI values between the two scenarios was larger than observed at the mid and high elevations, with Figure 5. Time series of daily surface water input (SWI) by elevation colored by the fraction of SWI sourced from snowmelt (f snow ) for control (CTL) and pseudo-global-warming (PGW) scenarios. A 7-day moving mean was applied for better visualization between bands, but this compresses peak daily SWI. SWI patterns are separated into three elevational categories: high, >3,300 m; mid, 2,300 m < elevation ≤ 3,300 m; low, ≤2,300 m.
SWI lower in the PGW scenario compared with the CTL scenario. We observed lower SWI in the PGW scenario relative to the CTL scenario in October, March, April, May, July, and August (Table S3 in Supporting Information S1), with the greatest reductions for April and May. Changes in mean daily SWI reinforced the patterns observed at monthly scale. While mid and high elevations showed (a) larger magnitude shifts in late winter and spring SWI ( Figure 4) and (b) elevated snowmelt contributions to SWI earlier in the year and reduced contributions later in the year ( Figure 5), low elevations displayed more muted changes both to the change in daily SWI as well as the change in daily snowmelt contribution to SWI.

Differences Between Control and Pseudo-Global-Warming Scenarios for High and Low Snow Years
Peak SWE and IS in both scenarios differed greatly between 2011 and 2012 ( Figure 3), with substantially higher IS, peak SWE and peak SWI during the high snow year (2011) than the low snow year (2012) for CTL and PGW. The timing of these metrics also varied greatly, with peak SWE timing occurring considerably earlier for the low snow year, especially at high elevations. Peak daily SWI for the high snow year for PGW showed greater differences between PGW and CTL values whereas peak SWI was more similar between mean annual and low snow year values for the two scenarios. Peak daily SWI appeared to change the most for PGW high snow years as compared to CTL values. The centroid of SWI (SWI50) for the low snow year for PGW showed much earlier timing as compared to CTL values for the same year.

Combined High, Mid and Low Surface Water Input Changes and Basin-Wide Response
Total monthly contributions to study area SWI were lowest for both CTL and PGW during the months of December, January, and February, with spring and summer months exhibiting the highest contributions ( Figure 6). However, June contributions for PGW were substantially lower than other spring and summer months, and substantially lower than for CTL, corresponding with lower precipitation for PGW in the preceding month (Figure 2a). Subtracting CTL monthly SWI contributions from PGW monthly contributions revealed higher October to March contributions to total study area SWI for PGW, and lower monthly contributions for June-September for all elevations as compared to CTL. For the months March-May, the direction of monthly contribution changes between PGW and CTL depended on elevation. In March, low elevations showed lower PGW monthly contributions to total annual study area SWI than for CTL, whereas mid and high elevations showed higher values. In April at low and mid elevations, the PGW scenario contributed less to total mean annual study area SWI than the CTL scenario. However, at high elevations, the PGW scenario contributed more to total mean annual study area SWI than the CTL scenario. In May, mid elevations showed substantially lower PGW contributions than during CTL, while high elevations had greater PGW monthly contributions to total annual study area SWI than CTL.
Though snowmelt-derived SWI showed similar spatial distributions across the study area for both scenarios (Figures 7a and 7b), less of the total SWI in the PGW scenario came from snowmelt ( Figure 7c). Elevation was a major CTL on the percent differences in mean monthly snowmelt-derived SWI between PGW and CTL ( Figure 7d, Table S3 in Supporting Information S1), with high elevations showing higher snowmelt-derived SWI Figure 6. The fraction of total mean annual study area surface water input (SWI) generated by each elevation category for each month for control (CTL) and pseudo-global-warming (PGW) scenarios. Elevation categories as follows: high, >3,300 m; mid, 2,300 m < elevation ≤ 3,300 m; low, ≤2,300 m.
values for PGW relative to CTL in fall and spring months, and reductions in snowmelt-derived SWI relative to CTL for summer months. Mid elevations showed higher contributions in fall, winter, and early spring months for PGW as compared to CTL, and reductions during late spring and summer months. Low elevations showed higher values for December and January, declines for October, November, February, March, April, and May, with no change for June to September.

Discussion
The results of high spatial resolution SnowModel simulations used in this study revealed novel elevational patterns of changes to SWE and SWI under continued global warming. Simulations indicate that high elevation SWE is retained under PGW (Figures 3 and 4, Section 4.1.1). Peak daily SWI was earlier and elevated under PGW, but SWI was also more evenly distributed throughout the year (Figures 3 and 4, Section 4.1.2). PGW led to altered timing, but not complete loss, of snowmelt contributions to SWI (Figures 5-7), with similar annual elevation Figure 7. Snowmelt contributions to total water year surface water input (SWI) for control (CTL) (a) pseudo-global-warming (PGW) (b) and PGW minus CTL (c), and the mean monthly PGW minus mean monthly CTL snowmelt contributions to total SWI by elevation range (d). Elevation categories in (d) are as follows: high, >3,300 m; mid, 2,300 m < elevation ≤ 3,300 m; low, ≤2,300 m. The gray contour in (a,b, and c) is for 2,300 m and separates the low and mid elevation category, while the black contour is for 3,300 m and separates the mid and high elevation category.
contributions of SWI to total study area SWI for CTL and PGW scenarios despite seasonal changes (Figures 6 and 7, Figure S12 in Supporting Information S1, Section 4.1.3). These findings suggest that high resolution snow models that (a) better capture mountain topography (e.g., Figures S17, S18 in Supporting Information S1) and (b) model forcings that better capture convective and orographic processes will lead to different predictions of snow patterns than models without these processes. This is an important consideration for future regional to national scale hydrologic modeling applications in snow dominated areas (Garousi-Nejad & Tarboton, 2021;Sexstone, Driscoll et al., 2020). Additionally, the IS metric  associated with SWI is an important advance in analyzing snowmelt and rainfall inputs to the soil surface. IS captures both timing and intensity of snowmelt and rainfall combined, which previously used metrics such as snowfall fraction (Berghuijs et al., 2014) and snowmelt rate (Barnhart et al., 2016) fail to fully capture. Finally, we show that existing snow monitoring stations become less representative of future SWE (Figure 8, Section 4.2), with implications for the predictive capacity of future streamflow forecasting.

Elevation-Dependent SWE and SWI Changes
PGW precipitation was greater for all elevations and all months (0.4-0.9 mm/d), consistent with observed trends for the Northern Front Range of Colorado (Fassnacht et al., , 2020, except for spring, when precipitation declined as compared to CTL (Figure 2). Temperature was higher for all elevations and months (4.9-5.6°C), with the greatest positive changes for late spring and fall months. Precipitation and temperature changes from SnowModel simulations were generally consistent with regional changes from CMIP6 scenarios (Almazroui et al., 2021). The fraction of precipitation falling as snow declined for all elevations and months (up to 40%), except for November-April at high elevations. This decline in snowfall fractions is consistent with results across the western United States for the mid-21st century (Klos et al., 2014), though our work demonstrates high elevations retain high snowfall fractions in PGW simulations in contrast to prior work showing ubiquitous declines (Rhoades, Jones et al., 2018;Wi et al., 2012).

Retention of High Elevation SWE Under Pseudo-Global-Warming
In response to changes in precipitation and temperature between the two scenarios, simulations revealed peak SWE declined between CTL and PGW for low and mid elevations, on the order of −14% to −45%, while the date of peak SWE uniformly shifted earlier in the year for all elevations. The magnitude of this shift was greatest for mid elevations (∼23 days, Figure 3). Declines in peak SWE for PGW in the study area are consistent with (a) Ikeda et al. (2021) who assessed SWE changes between 2070-2099 and 1976-2005 at 4-km resolution and found 30%-60% decreases in April 1 SWE in the region with peak SWE 20-30 days earlier and (b) Siirila-Woodburn et al. (2021) who reported a 25% decline in peak SWE for the western United States for the period 2025-2049, though in our study declines range from −45% for low elevations to −3% for high elevations in the study domain for our period of analysis. Our finding that low elevations show the greatest declines in SWE between CTL and PGW in the study region (Figures 3 and 4) is consistent with historical declines in SWE for similar elevations from 1984 to 2017 in the neighboring Rio Grande Basin (Sexstone, Penn et al., 2020). Limited high elevation snow change for PGW is consistent with Barsugli et al. (2020) who reported limited changes in high elevation April and May SWE using a 250-m resolution model forced with CMIP5 projections for 2041-2070. SnowModel simulations for the Northern Rocky Mountains of Canada similarly project snowpack loss that is strongly elevation and season dependent for the periods 2045-2059and 2085-2099as compared to 1979-1994(Mortezapour et al., 2022. Taken together, these comparisons emphasize that inland, high elevation areas will likely remain conducive to substantial snow accumulation for decades to come, especially when assessing high resolution model output that better resolves mountain topography and snow dynamics. In maritime, snowy regions of North America including the Sierra Nevada and Cascades, high resolution simulations that resolve mountain topography are needed to further understand the degree to which snow loss will occur by elevation.

Earlier and Elevated Peak Daily Surface Water Input, yet More Evenly Distributed Input Under Pseudo-Global-Warming
IS, which reflects the concentration of input in time, declined between CTL and PGW across all elevations, with similar reductions at all elevations due to higher winter input and lower spring input that result in a more uniform distribution of input throughout the year (Figure 3). Despite these reductions in IS, elevated peak daily SWI for PGW was observed for all elevations. The dates of peak SWI and SWI centroid moved to earlier in the year for all elevations and with the biggest advances for low and mid elevations, consistent with Xu et al. (2021). While prior work has shown earlier streamflow timing in the UCRB for future decades (Wrzesien & Pavelsky, 2020) our work highlights the potential for higher peak flows and a shorter period of elevated snowmelt SWI during snowmelt ( Figure 5). Earlier SWI timing for PGW as compared to CTL follow the results from Musselman et al. (2021), who show that snowmelt is already occurring prior to peak SWE in recent years as compared to early SNOTEL records of the 1980s and suggest that these observed patterns will continue in future decades. For low and mid elevations, our results generally align with Kampf and Lefsky (2016) who observed reduced snowmelt contributions to peak flows in areas where peak flows were historically sourced both from rain and snow contributions. Mixed input was a small component of total SWI for both CTL and PGW, consistent with rain on snow analysis by McCabe, Betancourt, and Hidalgo (2007), with mixed input only exceeding 10% of total input for approximately 1 week per year during CTL at the highest elevations and otherwise remaining at 0% of total input.

Altered Timing, but Not Complete Loss, of Snowmelt Contributions to Surface Water Input
In our study, the fraction of total annual SWI generated by snowmelt was lower for all elevations between CTL and PGW ( Figures 5 and 7), similar to the results of Li et al. (2017) who showed declines in snowmelt-derived streamflow across the UCRB for RCP4.5 and RCP8.5 out to 2,100. Seasonal snowmelt contributions to SWI were higher for fall, winter, and spring months at mid elevations ( Figure 5, Figure 7d), but lower for summer months in PGW, differing from the results of Li et al. (2017) who showed reductions for snowmelt generated streamflow for all months. Shifts in the timing of total SWI during spring months was dependent on elevation ( Figure 6), with low and mid elevations displaying these shifts earlier in the year than high elevations. Half of the SWI for in the study area CTL and for PGW was generated by elevations >2,700 m (only ∼37% of study area, Figure 1), generally consistent with prior research in the UCRB and the neighboring Cache la Poudre watershed where half of the streamflow is generated from elevations >3,000 m (Hammond et al., 2018;Harrison et al., 2021). Cumulative contributions to total study area SWI at water year scale retain a nearly identical relation with elevation for PGW as compared to CTL ( Figure S12 in Supporting Information S1) despite the shifts in timing of sub-annual SWI.

Implications for the Future of Timing and Magnitude of Water Availability
The combined, and potentially countervailing, effects of greater precipitation and earlier, flashier SWI complicate predictions for the future timing and magnitude of streamflow. The potential for greater precipitation in future decades, similar to what is shown by the PGW scenario, may offset streamflow reductions caused by the loss of efficiently generated snowmelt from seasonal snowpacks at mid elevations (Reclamation, 2016). However, future precipitation changes are uncertain (Almazroui et al., 2021), and while simulations suggest future droughts in western North America will be more spatially pervasive and extreme, uncertainties in vegetation processes and precipitation representation in models complicate the interpretation of these results (Cook et al., 2020). Additionally, warming during spring and summer months could partition summer SWI to soil and vegetation evaporative losses rather than streamflow or groundwater recharge Gordon et al., 2022;Hammond et al., 2019;Harpold & Molotch, 2015;Tague & Peng, 2013), as earlier snow disappearance can lead to earlier ET fluxes.
Prior work has shown that snowmelt rate and timing are linked (Alonso-González et al., 2022;Serreze et al., 1999), with slower snowmelt associated with decreased streamflow generation (Barnhart et al., 2016;Musselman, Molotch, & Molotch, 2017). However, recent research also highlights that earlier snowmelt timing and associated higher runoff efficiency (runoff generated per unit of input) may compensate for lower future snowmelt rates Gordon et al., 2022;Jeton et al., 1996). The competition of these two mechanisms may dictate how runoff from snow-dominated sites might change in the future as snowmelt rate and timing change under future climate conditions. Barnhart et al. (2020) report that snowmelt season runoff was governed by the competing influence of snowmelt timing and snowmelt rate via a modeling study where each factor was considered separately due to the linkage between snowmelt rate and timing (Trujillo et al., 2014). The simulated changes in subsurface storage and melt season runoff were more sensitive to changes in snowmelt timing than snowmelt rate for a site in Colorado, suggesting greater runoff from earlier snowmelt may counteract runoff losses due to slower snowmelt and that this process is mediated by plant available water storage and energy availability. In the Lower Colorado River Basin, the effects of temperature-induced snow loss on annual streamflow appear to have been moderated by enhanced winter hydrological inputs and streamflow production (Robles et al., 2021).
In the current hydroclimate, streamflow generation for many western U.S. mountain watersheds tends to peak in spring whereas runoff efficiency peaks in winter. As watersheds become warm enough to experience winter melt and rain on snow, this enhanced winter runoff efficiency may act as a temporary buffer against reductions to streamflow generated by concentrated spring and summer snowmelt . Results of the present study indicate earlier SWI with higher peak SWI values at all elevations ( Figure 3). This result, which integrates changes in snowmelt and precipitation phase through SWI, contrasts with Musselman, Clark, Liu et al. (2017), Musselman, Molotch andMargulis (2017), who, using the same meteorological forcing data as this study, reported declines in future snowmelt input to the hydrologic system and possible concomitant declines in streamflow. Earlier SWI may lead to efficient streamflow generation if soil moisture is conducive to infiltration excess overland flow, saturation excess overland flow, or shallow subsurface flow. However, for warmer low and mid elevations where snowmelt is a smaller component of SWI, earlier SWI was simulated to be associated with reduced annual total SWI . Integration of high resolution SnowModel output with a hydrological model would be useful to assess study-area-wide hydrologic partitioning and streamflow generation for future periods. This integration is needed to examine the combined effects of (a) changes occurring to deep snowpacks in relatively small areas of the western United States and (b) the near loss of shallow snowpacks across large areas on basin-wide water availability.
The combined effects of shifts in timing and magnitude in SWE and SWI also have broad implications for reservoir management and wildfire dynamics. Continued warming is diminishing role of snowy mountains as natural reservoirs, and thus water resource management will likely need to adjust to shifting reservoir inflow timing and magnitude . With drought severity and extent expected to increase in coming decades (Cook et al., 2020;Gutzler & Robbins, 2011), interannual variability in annual SWE and SWI metrics (highlighted by years 2011 and 2012, Figure 3) could test coordinated reservoir management for multiple objectives. Even if annual streamflow remains the same or increases slightly, it could be hard for coordinated reservoir management to adjust to changes in the timing of inflows (Cohen et al., 2020;Li et al., 2010;Sterle et al., 2020), given a greater offset between water supply and demand timing. These changes will be most impactful and relevant for water management decision making when hydroclimatic shifts for individual years compound and result in changes to decadal and multi-decadal variability (Garrick et al., 2008;Hurkmans et al., 2009;McCabe, Clark, & Hay, 2007;Nowak et al., 2012;NRC, 2007).
Hedging in reservoirs, the practice of storing more water during the spring to meet the widening gap between inflows and outflow demand during the summer low flow period, may be needed to respond to shifting seasonal reservoir inflows (Jones & Hammond, 2020). Other potential adaptations to shifting SWI timing and magnitude for downstream water availability include the construction of new surface water reservoirs or the enhancement of existing infrastructure (Perry et al., 2017;Reclamation, 2021), the adjustment of current reservoir management strategies (Delaney et al., 2020;Reclamation, 2021;Sumargo et al., 2020), and recharging groundwater during periods of excess (Scanlon et al., 2016;Zhang, 2015), and according to Kellner and Brunner (2021) immediate action is required to plan for balanced upstream-downstream water availability in water infrastructure and management decision making. Over the past several decades, earlier snowmelt occurrence and reduced peak SWE have corresponded with increased wildfire activity and wildfire severity (Westerling, 2016;Westerling et al., 2006) as well as with increasing solar forcing on snow in the western U.S (Gleason et al., 2019). Wildfires have increasingly occurred in seasonally snow covered areas from 1984 to 2020 throughout the western United States (Kampf et al., 2022). Continued snow loss and earlier melt may lead to larger and more intense fires in these areas in future decades, though potential increases in spring and summer precipitation could have a countervailing influence.

Snow Measurement Network Coverage and the Distributions of SWE and SWI
Current snowpack monitoring does not represent the range of land covers, elevations and terrain classes found in mountainous regions (Fassnacht et al., 2012;Kampf et al., 2020;Sexstone & Fassnacht, 2014), partially due to the difficulties of measuring snow in alpine areas. Stations do not always reflect their nearby surroundings (Meromy et al., 2013), and conditions in mountainous regions experience substantial inter-annual variability (Fassnacht et al., 2003(Fassnacht et al., , 2012 especially in the alpine (Sexstone et al., 2018). The 75 existing SNOTEL sites in the study region range in elevation from 2,605 to 3,532 m, with a mean elevation of 3,094 and median of 3,127 m. Although the mean and median elevation of the SNOTEL sites is greater than the median and median elevations of the modeling domain, the SNOTEL site elevation range does not cover much of the persistent snow zone above 3,500 m in the high elevation category (Figure 8). Prior work has suggested these networks will have reduced predictive capacity for streamflow forecasting as snowpacks continue to decline (Livneh & Badger, 2020). A lack of observations at both the high and low extents of where snow accumulates could inhibit the ability to predict snow and streamflow conditions.
For both CTL and PGW, the existing network of SNOTEL stations does not monitor areas with the highest mean annual peak SWE values (Figure 8a), though it does sample the elevations with the greatest volumetric SWE in both scenarios (Figure 8b). With the greatest declines in SWE for low and mid elevations for PGW, greater future volumetric SWE will occur at high elevations not sampled by the existing network. Adding sites at higher elevations not presently monitored could allow for better tracking the remaining water content within snow to provide information for downstream streamflow forecasting and coordination. However, in future decades, treeline may advance in elevation, shrinking alpine areas (Kummel et al., 2021), potentially making existing stations more reflective of study area snowpacks in the future sub-alpine. The study area examined here is generally higher in elevation and snowier than other areas of the UCRB (Figures 1a and 1b), so additional snow monitoring in this area may provide more useful information on the distribution of SWE and snowmelt-derived SWI that drives water availability than instrumentation in other areas of the basin.
Considering the volumetric distribution of mean annual SWI (Figure 8c), the existing SNOTEL network generally samples elevations above those with greatest SWI. Adding additional sites to slightly lower elevations could better capture the timing and rate of rainfall and snowmelt input both in the near-term and long-term. Low and mid elevation sites would also better capture mid-winter melt not currently measured by exiting monitoring stations in low and mid elevations that measure precipitation but not SWE or soil moisture. Given that a greater fraction of precipitation in future years will fall as rain rather than snow (Figure 2, Klos et al., 2014), precipitation measurements at low and mid elevation sites may suffer less from under catch since snow is more susceptible to this phenomenon than rain (Groisman & Easterling, 1994). USGS gaged headwater watersheds are generally more uniformly distributed across elevations than SNOTEL sites, and although there are a number of gaged headwater watersheds draining only the high elevation persistent snow zone, there are fewer gaged watersheds isolating the transitional or low snow zones; rather, gages at lower elevations tend to integrate the hydrologic response from multiple snow zones . Thus, additional gaging better isolating the hydrologic response to different snow zones could also provide information for downstream flow forecasting and coordination. Additionally, soil moisture is sparsely monitored in the basin, but additional soil moisture monitoring along elevation gradients in the basin may help to mitigate lessened predictive capacity due to declining snowpacks at existing snow monitoring stations (Livneh & Badger, 2020).

Limitations, Uncertainties and Areas for Continued Research
Though the results of this study were coherent between analyses implemented at multiple time scales (daily, monthly, annually) and broadly consistent with prior research, several factors limit the scope of inference of our study. First, our model evaluation revealed poorer performance in low snow areas, areas with limited existing snow monitoring to use for model calibration and validation. Further work is needed to evaluate whether modeling snow accurately in these low snow areas is important for assessing water availability at management relevant spatial and temporal scales. Second, though SNOTEL sites use Alter wind shields to increase precipitation catch efficiency (Serreze et al., 1999), under catch has been approximated at ∼11% (Yang et al., 1998), with up to 30% in an intense snowstorm with high winds in Colorado (Rasmussen, Liu et al., 2011). Our use of SNOTEL data for model evaluation is thus susceptible to this source of error (Scalzitti et al., 2016a(Scalzitti et al., , 2016b. Additionally, several of our analyses were based on aggregation of model output to elevation bands, with local patterns lost in the process. Further analysis of model output for sub-regions may identify nuanced patterns compared to those identified using study-area-wide elevational aggregation. There are also limitations associated with the regional climate model data used as input to this study's SnowModel simulations. Liu et al. (2017) broadly characterized the biases inherent in the WRF PGW data set and noted a warm and dry bias in the central U.S. in the late summer months. However, given that this issue is related to convective storms in the late summer, it is unlikely that this would substantially affect the main results of our analyses. There are also limitations to using a single model representation of future climate precipitation given that it is not possible to characterize uncertainty with this method. As ensemble convection-permitting simulations become available on longer timescales, future research could analyze changes in snowpack from an ensemble framework to assess uncertainty. Finally, although the PGW methodology does not represent large-scale changes in storm tracks associated with climate change, the methodology allows assessment of specific time periods and see how actual weather events in that time period might have been different if they experienced thermodynamic climate warming and enhanced moisture.
Subsequent snow modeling efforts could explore including dynamical vegetation to account for increased forest cover at high elevations (e.g., Barnhart et al., 2021), to explore the effects of significantly changed canopy structure due to more prevalent wildfire in snowy areas (Kampf et al., 2022), and to address the effects of forest disturbance on snowpack and streamflow (Goeking & Tarboton, 2020;Varhola et al., 2010). The existing simulations might also be used to examine differences in SP by aspect and landcover (e.g., Barsugli et al., 2020;Sexstone, Penn et al., 2020). While SWI as used in this analysis provides an initial estimate of runoff resulting from rainfall and snowmelt input available at the soil surface, this estimate does not account for the various pathways that SWI could take into the soil, across the subsurface, or into the atmosphere. Thus, connecting the surface water inputs from this study to a hydrological model would enable detailed exploration of changes to evapotranspiration, soil moisture, groundwater recharge, and streamflow generation due to snow loss, changes in SWI timing and magnitude, and changes to the synchrony of energy and water availability (e.g., Gordon et al., 2022).

Conclusion
The results of high spatial resolution PGW SnowModel simulations for the Colorado portion of the UCRB revealed widespread SWE loss, with reductions in peak SWE at low (−45%, −30 mm) and mid (−14%,-26 mm) elevations. However, high elevation, continental areas in the bounds of this study, and elsewhere in the UCRB, are projected to experience a more limited peak SWE change (−3%) and will likely remain conducive to seasonal snow accumulation and melt for decades to come. Peak SWE and peak SWI occurred earlier in the year (peak SWE, 17-23 days; peak SWI, 14-23 days) for all elevations under PGW, with higher peak daily SWI magnitude (30%-42%, 10-12 mm/d) also noted at all elevations. Despite earlier and elevated peak daily SWI, reduced IS revealed more evenly distributed SWI under PGW. Snowmelt-derived surface water inputs were lower at the annual scale, but experienced seasonal growth and shifts in timing dependent on elevation for PGW as compared to CTL, with low elevations experiencing shifts earlier in the year. Finally, existing snow monitoring stations became less representative of future SWE, because there are few monitoring sites at high elevations where future SWE will contribute a larger fraction of total SWI compared to current climatic conditions. Water stored in snowpacks has historically served as a reliable extension of manmade reservoir storage, providing sustained streamflow generation throughout spring and summer months. The combined effects of shifts in timing and magnitude in SWE and SWI have broad implications for (a) reservoir management given a greater offset between water supply and demand timing under PGW, and (b) future wildfire dynamics because large fires are already occurring in snowy mountain areas where SWE has declined, and snowmelt has shifted earlier in the year. These results and added insights suggest high resolution snow models that better capture mountain topography and model forcings that better capture convective and orographic processes may improve representation of seasonal snow accumulation and melt dynamics, elevational variability in snow patterns, and changes in precipitation phase under future climate conditions. This is an important consideration for future regional to national scale hydrologic modeling applications in snow dominated areas. Additional work that connects high resolution SnowModel outputs to a hydrological model could benefit downstream coordinated water management by to evaluating the combined effects of widespread snow loss and earlier, flashier SWI.

Data Availability Statement
The data that support the findings of this study, including model input and output, are available in USGS Science-Base at https://doi.org/10.5066/P9K0QUCN. Availability and Use Science Program as part of the Water Resources Mission Area Snow Hydrology Research Project. Additional funding support was provided by the Colorado Water Center at Colorado State University. Computing resources were provided by USGS Core Science Systems Advanced Research Computing Center. Thanks to Gregory McCabe (USGS, Water Mission Area) and two anonymous reviewers for their suggestions and edits that improved this article. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.