Illuminating Snow Droughts: The Future of Western United States Snowpack in the SPEAR Large Ensemble

Seasonal snowpack in the Western United States (WUS) is vital for meeting summer hydrological demands, reducing the intensity and frequency of wildfires, and supporting snow‐tourism economies. While the frequency and severity of snow droughts (SD), that is, anomalously low snowpacks, are expected to increase under continued global warming, the uncertainty from internal climate variability remains challenging to quantify with observations alone. Using a 30‐member large ensemble from a state‐of‐the‐art global climate model, the Seamless System for Prediction and EArth System Research (SPEAR), and an observations‐based data set, we find WUS SD changes are already significant. By 2100, SPEAR projects SDs to be nearly 9 times more frequent under shared socioeconomic pathway 5‐8.5 (SSP5‐8.5) and 5 times more frequent under SSP2‐4.5, compared to a 1921–2011 average. By investigating the influence of the two primary drivers of SD, temperature and precipitation amount, we find the average WUS SD will become warmer and wetter. To assess how these changes affect future summer water availability, we track late winter and spring snowpack across WUS watersheds, finding differences in the onset time of a “no‐snow” threshold between regions and large internal variability within the ensemble that are both on the order of decades. We attribute the inter‐regional variability to differences in the regions' mean winter temperature and the intra‐regional variability to irreducible internal climate variability which is not well‐explained by temperature variations alone. Despite strong scenario forcing, internal climate variability will continue to drive variations in SD and no‐snow conditions through 2100.


Introduction
Mountains play an indispensable role in Western United States (WUS) water supply, as their low temperatures and high precipitation capture significant water reserves in the form of snowpack.Often referred to as the "water towers" of the West, mountains store enormous amounts of winter precipitation which is measured as snow-water equivalent (SWE), or the depth of water if all snow melted instantaneously.During the dry spring and summer, the SWE is released as meltwater and supplies human populations whose water needs continue to rise (Bonsal et al., 2020).A reliable snowpack provides security to human populations across the WUS by providing water for increasing agricultural demands (Barnett et al., 2005), reducing the severity and intensity of wildfires (Trujillo et al., 2012;Gergel et al., 2017), and improving snow tourism economics (Wobus et al., 2017).According to Wobus et al. (2017), ski resorts are expected to lose 50% of ski season length by 2050 and 80% by 2090.Despite large seasonal variability, climate change has already been found to have significantly decreased SWE globally and across the WUS, particularly in late winter (Barnett et al., 2005;Kapnick & Hall, 2012;Fontrodona Bach et al., 2018;Huning & AghaKouchak, 2020).When SWE is abnormally low, the region is said to experience a snow drought (SD).SDs are driven by either warming, as a phase change from frozen to liquid, or reduced precipitation amounts.They affect the WUS's economy and human activity, even in areas far from mountain snowpack that rely on spring and summer melt waters for crop production and human consumption.
The adverse effects of SDs on a region's hydrology vary depending on the type of SD.Dry SDs, characterized by low precipitation and near-or below-normal temperatures, result in low streamflow throughout the melt season.In contrast, warm SDs occur under near-or above-normal precipitation and warm temperatures and often lead to early season snowmelt, increased spring flood risk, and summer hydrological drought (Harpold et al., 2017).While deviations from normal temperature and precipitation dictate SD occurrence, their absolute conditions impact how SDs are expected to respond to climate change.Shrestha et al. (2021) demonstrate that additional warming above a critical average winter temperature threshold of -6 to -5°C decreases snowpack.As all WUS large hydrologic unit code (HUC2) regions have historical average winter temperatures at or above -5°C, we expect their snowpack to be vulnerable to any level of warming.
To study SD across the WUS, we focus on comparing changes in SWE.Large observational uncertainty in WUS SWE measurements implies high biases are likely between any two datasets or models (Wrzesien et al., 2019).Observational model bias is driven by low sampling rates and terrain complexity, present in mountain regions, and is further magnified by assumptions in models used to generate SWE estimates (Wrzesien et al., 2019).Coupled global climate models (GCMs) are expected to produce snowpack estimates that are biased compared to observations because they have a lower spatial resolution and have temperature and precipitation biases (McCrary et al., 2017;Wrzesien et al., 2019;Kim et al., 2021;McCrary et al., 2022).Despite these biases, Matiu and Hanzer (2022) show that many models exhibit uniformity in simulating robust decreases in WUS SWE.Huning and AghaKouchak (2020), for example, have shown that SD total duration, average duration, and intensity in the WUS have increased by 28% between 198028% between and 201828% between , and Shrestha et al. (2021) ) adds that these conditions are expected to continue to worsen because of the WUS's low latitude.These previous results imply that although GCMs are typically biased in their SWE base state, changes relative to their base states are still informative.As a result, we will primarily focus on comparing changes in SWE across data sets.
To investigate historical and future changes in SD frequency and intensity we use 30-member initial condition large ensembles from a state-of-the-art coupled global climate model, called the Seamless System for Prediction and EArth System Research (hereafter SPEAR) (Delworth et al., 2020).To assess SD intensity relative to the historical period, we focus on SPEAR's simulation of severe to exceptional snow droughts (D2+ SD) and follow the classification framework used by the US Drought Monitor (Svoboda et al., 2002).We first show that SPEAR accurately simulates changes in WUS SD by comparing it to an observationally based dataset and with previous studies across the historical period   (Livneh et al., 2013;Huning & AghaKouchak, 2020).The classifications in SPEAR show both an increase in D2+ SD occurrence across the historical period and a continued increase under future warming scenarios.To understand the conditions driving these SDs, we examine the average temperature and precipitation conditions for the study period, finding that temperature and not lack of precipitation is the main driver of the D2+ SD increase at monthly time resolution.We then provide a region-level assessment of the transition to a "no-snow" environment by the end of the 21 st Century that accounts for scenario uncertainty and internal climate variability.
By separating the uncertainty into the portion attributable to internal climate variability and emissions uncertainty, we can determine the distribution of D2+ SD changes until 2100, the variability in the conditions that generate drought/non-drought conditions, and the probability distribution of the transition timing to a no-snow regime.We assess these changes under two scenarios in the SPEAR projections (2014-2100): a middleof-the-road scenario (Shared Socioeconomic Pathway 2-4.5, hereafter SSP2-4.5), and a high-emissions scenario (SSP5-8.5)(Delworth et al., 2020).While the two emissions scenarios allow us to explore the effects of emissions uncertainty, the 30-member ensembles enable the estimation of internal climate variability.

SPEAR Large Ensemble Global Climate Model
To assess changes in the probable distribution of historical and future SD, we analyzed WUS SWE in multiple 30-member SPEAR large ensembles (Delworth et al., 2020).
SPEAR is a coupled global climate model recently developed at the NOAA Geophysical Fluid Dynamics Laboratory (GFDL) that is designed for improved prediction and projection on seasonal-to-multidecadal timescales.SPEAR is composed of GFDL's AM4 atmosphere, LM4 land, MOM6 ocean, and SIS2 sea-ice models.These component models are the same as GFDL's Global Climate Model version 4 (CM4) (Held et al., 2019), which is a contributor to the Coupled Model Intercomparison Project phase 6 (CMIP6).SPEAR's configuration differs from CM4 as its physical parameterization choices are optimized for climate prediction on seasonal to centennial timescales.SPEAR has a moderately high atmospheric and land-surface resolution (approximately 50 km) and a coarser ocean and sea-ice horizontal resolution of about 1°, which has meridional refinement to 0.33°at the equator.For this study, we use SPEAR's monthly SWE, temperature, and precipitation across the historical period and projections from 2014-2100 under both SSP2-4.5 and SSP5-8.5 emissions scenarios.

Observational Data
To evaluate SPEAR's historical simulation of SWE, temperature, and precipitation, we use an observations-based dataset (Livneh et al., 2013), available from 1915 to 2011, hereafter the Livneh dataset.Livneh uses statistically gridded in situ daily precipitation and temperature observations on a 1/16°grid to generate SWE estimates (among other land surface variables) using the Variable Infiltration Capacity (VIC) land model (Liang et al., 1994).To compare the Livneh dataset with the SPEAR ensemble members, we re-gridded Livneh to SPEAR's 1/2°grid and re-sampled it to SPEAR's monthly timescale.Despite incorporating observational data, gridded datasets, like Livneh, retain large uncertainties across variables including temperature, precipitation and SWE (Walton & Hall, 2018;Wrzesien et al., 2019).Many recent papers have found SWE estimates to vary widely, by upwards of a factor of 3 in some cases (Wrzesien et al., 2019), leading us to expect significant absolute biases between SWE estimates (McCrary et al., 2017(McCrary et al., , 2022)).To overcome this issue, we focus our analysis on proportional changes, comparing SWE values to their own historical distributions within each dataset, and then comparing these relative changes across datasets.
We chose 1921-2011 as our historical period as it is the overlapping period of the Livneh and historical SPEAR datasets.We use the 90 complete winters to validate SPEAR and develop a baseline against which to compare the modeled future climatology.We chose to consider data at monthly resolution intervals for the following three reasons: (1) data availability, as SPEAR only recorded SWE at monthly intervals; (2) consistency with previous studies (Huning & AghaKouchak, 2020); and (3) because the monthly resolution is an appropriate timescale for monitoring snow drought.As ENSO and PDO drive inter-annual variability across the region, assessing SPEAR's representation of these teleconnections is important for understanding how accurately the model may reproduce other extremes across the region, like SDs.Delworth et al. (2020) shows SPEAR accurately captures the relationship between PDO and North American precipitation, while Maher et al. (2022) (Deser et al., 2020) and is essential for modeling extremes.When evaluating model bias, however, it means that, short of a taking a long-term average as shown in Figure 1, we do not expect biases between observations and either a single SPEAR ensemble member or the SPEAR ensemble mean to be reflective of SPEAR's accuracy in simulating the climate.While we do not expect a single SPEAR ensemble member or the ensemble mean to reproduce Livneh exactly, we do expect SPEAR to simulate a realization of the climate at least as extreme as the observed historical climate over most regions.However, with only 30 ensemble members it is still reasonable to expect an occasional observation to fall outside of the SPEAR spread.Thus, if the change in SD frequency observed in Livneh falls within the SPEAR ensemble spread, we can assume SPEAR produces a realistic historical climate.Our analysis reveals that the majority of the Livneh SWE statistics fall near clusters of SPEAR ensemble members, further strengthening the conclusion that SPEAR accurately represents the WUS climate as demonstrated in Figures 2 and S3.

Drought Classification
Before we can assess changes in SD, we first introduce our SD classification method.
To ensure that only regions which typically have snow are eligible for classification, we restrict our region of study to the "historically snowy" region, areas that historically have average seasonal SWE maxima above 20 mm, based on the SPEAR ensemble mean.We then assign a classification based on how extreme each month is compared to the historical distribution of SWE across all grid cells and months.
Our methodology assigns standardized indices to each location by month and uses the US Drought Monitor's (USDM) drought classification method for hydrological drought to categorize observations into six descriptive bins: near normal (NN), abnormally dry (D0), and moderate (D1), severe (D2), extreme (D3), and exceptional (D4) drought.Wet conditions are classified analogously, with labels W0-W4 for increasingly wet months; see Figure S2 (Svoboda et al., 2002;Huning & AghaKouchak, 2020).We use a non-parametric empirical model to classify SWE, temperature, and precipitation values for each month.
Without assuming the underlying distributions, a non-parametric model allows us to efficiently capture the variability without imposing subjective constraints on the data.
We begin by assigning each extended winter month of the year (Oct-April) a score based on the historical conditions at that location.Our time indices are by year (y) and month (m), e.g.t 1931,1 for January 1931, and spatial indices are in degrees latitude (i) and longitude (j).For example, s t1931,1 40.5,250 corresponds to a SWE value at latitude-longitude pair (40.5, 250) during January 1931.We now compute an empirical distribution over , representing the historical SWE values during month m at location (i, j).We then assign a value in (0, 1) to each SWE measurement using the empirical cumulative distribution function, F m i,j , based on the proportion of the observed data in S m i,j that fall below it.In equation 1, I(•) takes the value 1 if SWE measurement x is larger than the historical SWE measurement, S ty,m i,j , and 0 otherwise.We sum over the historical period which ranges from 1921 to 2011, which is 91 complete years.
For each observed or simulated SWE value, s ty,m i,j , we can then compute the z-score by plugging the SWE value into the corresponding F and then into the inverse normal distribution, Φ.We refer to these z-scores as ZSWE, which are indexed by location, month, and year.We can now classify snow droughts from the SWE value, s ty,m i,j , using Each month is then assigned a classification (W4-W0, NN, D0-D4) which can now be compared across regions.While we primarily use this framework to classify SDs, we extend the classification scheme to temperature and precipitation as needed.
A similar empirical methodology is used by Huning and AghaKouchak (2020) to classify snow droughts across the Alps, Himalayas, and WUS.Their framework is inspired by the USDM which uses the same D0-D4 classification.However, the USDM approach is not purely statistical, relying on experts to incorporate regional sensitivity into the published drought classification.Without experts, our model attempts to match the frequency of meteorological droughts in the US Drought Monitor (USDM) with snow drought frequency because the USDM is the widely accepted standard, despite its subjectivity (Svoboda et al., 2002).While our method may result in a mismatch of SWE values and impact in some locations, it provides a statistically-rigorous way to quickly capture extremes without gathering detailed human and environmental data for each pixel.

Computing Changes in Snow Drought
We can now apply our drought classification scheme to evaluate how well SPEAR reconstructs historical changes.We define two 41-year windows containing 40 complete winters to assess change, and after applying our drought classification scheme to snowpack data aggregated to the HUC2-level, we count the number of D2+ SD occurrences across the early and late historical periods, given by a ZSWE of less than −1.3, e.g.I(Z t R < −1.3) for HUC2 region R at time t.The percent change for a given region, ∆ R , is de- For example, in the Upper Colorado region, 27 months of Livneh-derived D2+ SD occur in the early historical period and 28 in the late historical period, translating to an increase of 3.7%.Next, we leverage the SPEAR ensemble spread to determine whether the overall trend is significant.

Snow Transition Threshold
In addition to evaluating drought climatology, we are also motivated to determine how a changing SWE will affect water resources.We seek to discern when a shifting climate will begin to severely and persistently impact snow as a water resource.Long-term droughts are particularly damaging, as one or two years of low snow-pack can be buffered by groundwater, above-ground reservoirs, or stored in live biomass, but these buffers dwindle with extended exposure to drought conditions.Thus, we are particularly interested in determining when no-snow conditions are expected to become systemic (Siirila-Woodburn et al., 2021;Harpold et al., 2017).
To determine this transition, we focus on April SWE because April typically corresponds to peak SWE.By first calculating the fraction of April 15 th SWE remaining in the historically snowy portion across each of the 5 HUC2 regions: Upper Colorado, Lower Colorado, Great Basin, Pacific Northwest, and California (abbreviated UC, LC, GB, PNW, and CA), we can classify an April (m = 4) grid cell s t,4 i,j as no-snow for that year if there is at most 10% of the historical snowfall average remaining at the location (Siirila-Woodburn et al., 2021).We then calculate the regional no-snow area proportion as the fraction of the historically snowy region which experiences those conditions.Formally, we let N Y R denote this no-snow area proportion, where R represents the region, for our application a WUS HUC2, and Y the year.As before, St Y ,4 i,j is the average historical SWE value for the grid cell and s t Y ,4 i,j the SWE value for the specific year.Using 10% as our no-snow threshold, T = 0.1, then N Y R can be written as: Thus we have a fraction of the historically snowy region that is snow free in a given year in April.To assess when no-snow conditions become endemic, we apply a 10-year movingwindow mean and then define the no-snow transition time as the year when the movingwindow mean last crosses the area threshold, A, before 2100.Applying this procedure to all ensemble members, we compute a distribution for when these conditions are likely to become endemic.Formally, the no-snow transition time for an ensemble member, T , is given by: where Ñ t ′ R gives the moving-window mean fraction of region R that experiences no-snow conditions at time t ′ .By requiring the moving-window average to be above A for all subsequent years (until 2100), T is uniquely determined.For a graphical explanation of this method, please refer to Figure S5.historical period -in fact, a 95% confidence interval for the SPEAR ensemble mean across four of the five regions contains the 28% benchmark for drought intensity increases found in Huning and AghaKouchak (2020), with only the UC interval exceeding the benchmark with a 30% lower bound on historical D2+ SD increases.While we could not use the same historical period due to data constraints, the agreement helps to further validate the SPEAR ensemble.See supplemental Text S1 and Figure S3 for an analysis of changes in precipitation and temperature across the historical period.

Analyzing SWE into the 21 st Century: Accelerating Loss
We next shift our attention to projected changes in 21 st century D2+ SD, focusing first on changes in droughts classified with our ZSWE metric.We construct our empirical CDF F m i,j distributions from the historical period  and calculate corresponding ZSWE scores for each winter month across the historically snowy west (2014-2100) for all 30 ensemble members.Projected changes in SWE are dramatic, with rapid increases in D2+ SD occurring at mid-century (Figure 3).Under SSP5-8.5, we find that towards the end of the century, all regions are projected to experience severe, extreme, or exceptional SD during most months.Under SSP2-4.5, SD increases are less severe, with conditions by the end of the century resembling conditions under SSP5-8.5 by midcentury.As expected, the higher forcing scenario corresponds with an accelerated timeline for increases in snow drought frequency.SD frequencies for all 18 study decades are shown in Figure S4.
Examining the spatial distribution of D2+ SDs in Figure 3, a pattern of regional "hot spots" emerges through time.D2+ SD frequency is consistently higher in certain regions beginning in 2030 in SSP5-8.5 and SSP2-4.5.For example, the Washington Cascades and Colorado Rockies are projected to experience more frequent D2+ SD across all decades than regions in south-central Idaho and the California Sierra Nevada.We expected to see more dramatic D2+ SD increases in the southern basins, including the Cal- ifornia and Lower Colorado regions, as Shrestha et al. (2021) found that even low amounts of warming at southern latitudes result in strong SWE loss signals.We assume the hot spot pattern emerges because we are looking over a narrow enough range of latitudes that the latitude signal is overshadowed by regional variation, perhaps coming from elevation variability.Shrestha et al. (2021) examined basins ranging from the Yukon to Columbia River basins that have average winter temperatures of -8°C to +4°C, finding that below -5°C to -6°C warming temperatures did not reduce SWE.Our HUC2 regions had mean winter temperatures in historically snowy regions ranging from -5.1°C (UC) to 0.3°C (California).Therefore, we expect any amount of warming will decrease SWE and correspondingly increase D2+ SD.
While Figure 3

Temperature and Precipitation Controls on SWE
As changes in SWE are primarily driven by changes in temperature and precipitation climatology (McCrary et al., 2017;Harpold et al., 2017), we next examine changes in SWE in the phase space spanned by temperature and precipitation.By aggregating over the entire historically snowy WUS, we can determine how temperature and precipitation anomalies are driving the dramatic increase in SD.In Figure 5, each dot represents the average temperature and precipitation anomaly by decade and is colored according to the average ZSWE score.By definition, the average all-month historical (1921-2011) temperature and precipitation mean is (0, 0).However, by breaking the century down by decade we can see variation within the 20 th century.
As expected, all-month decadal averages in the historical period cluster around a zero temperature and precipitation deviation.We observe small changes in anomalies before 2000, a finding consistent with our understanding of changing D2+ SD frequency.
Beginning in the 2000s, the all-month decadal-average rapidly shifts towards warmer and wetter conditions.By 2050 under SSP5-8.5, the average temperature and precipitation are 1.50 and 0.25 standard deviations higher than the 20 th century average, respectively.This corresponds to a dramatic warming and slight wetting across the WUS and indicates the average month in 2050 to be warmer than 93% of months in the historical period for a given location.For SSP2-4.5, the values are 1.18 and 0.20, respectively, reflecting a moderate increase in temperature and precipitation by mid-century, with the average month in 2050 being warmer than 88% of historical months.
To understand changes in SDs, we also track the underlying climatology of months that experience D2+ SD.Outlined in grey in Figure 5, we find historical D2+ SD averages are both dry and warm with an average temperature and precipitation anomaly of 0.6 to 0.8 and -0.6 to -0.8, respectively, indicating historical snow droughts are primarily driven by a near equal combination of both warm and dry conditions.These conditions suggest that an average historical D2+ SD month is both warmer and drier than 75% of months.However, when examining SPEAR's future climate, we find the average drought is both warmer and wetter.By 2050 under SSP5-8.5, the temperature deviation is 1.84 while the precipitation deviation is -0.015, indicating that future D2+ SDs are significantly warmer than the historical ones and that dry conditions are no longer needed to produce a SD.We conclude future D2+ SD conditions are driven by the increasingly high-temperature average, which is warmer than 97% of historical conditions.
By 2090, the average drought month has a temperature deviation of 2.18 and a precipitation deviation of 0.27, close to the all-month anomalies of 2.10 and 0.36 for temperature and precipitation, respectively.Average monthly temperature for both D2+ and all-month averages are in the 98th percentile of historical conditions, indicating that future winter conditions will, on average, be extremely warm and that the difference between average conditions for all months and SD months has decreased.Examining the ZSWE scores for 2090 under SSP5-8.5 confirms that the convergence is also reflected in SWE changes, with the average all-month ZSWE being -1.79 and the average D2+ month having a ZSWE of -2.10.Thus, the 2090 all-month average is expected to be a D3 SD, while the average month classified as a SD is D4.Under SSP2-4.5, conditions do not reach such an extreme, with average all-month conditions by 2090 reaching 1.48 for temperature, 0.27 for precipitation, and -1.10 ZSWE.The temperature, precipitation, and ZSWE deviations for the months that experience D2+ SD are 1.75, 0.064, and -1.91, respectively.
Although the gap between drought months and all-months shrinks, the difference is far less extreme than under SSP5-8.5; the average month under SSP2-4.5 is only given a D1 snow drought classification.The convergence of the all-month and drought-month temperature and precipitation anomalies, particularly under SSP5-8.5 emphasize that D2+ SDs will require increasingly smaller deviations from normal conditions to produce.This underscores that SDs will become a "new normal" for the WUS by the end of the 21 st century.

Timeline for Snow-Free Conditions
In addition to changes in D2+ SD frequency, we also examine how total SWE availability is expected to change, by assessing the timing of Western regions' transition to a no-snow regime.A no-snow regime, characterized by a 10-year moving average of April SWE consistently below 10% of the historical April average, indicates severely limited summer water supply from SWE.To understand when a no-snow regime is likely to affect a HUC2 region, we examine the distribution of transition times to no-snow across SPEAR's ensemble members.By varying the area threshold, A, we can assess how quickly conditions are expected to deteriorate.Figure 6 shows the distribution of the transition to no-snow regimes for 3 different area thresholds, A: 50%, 75%, and 90%, for the historically snowy HUC2 regions.Note that by construction, an individual ensemble member's transition year always occurs later for higher A. However, the ensemble distributions can overlap, which indicates large variability in the severity of conditions, especially later this century.When aggregated to the entire historically snowy WUS ("West-Wide"), the average transition time for A = 50% is 2071 for SSP2-4.5 and 2048 for SSP5-8.5.However, when considered as separate regions, transition times for A = 50% varied from as early as 2025 (CA) to 2088 (UC) under SSP2-4.5 and 2018 (CA) to 2056 (UC) for SSP5-8.5.
The snow-free transition distribution center occurs later for all regions under SSP2-4.5 scenario than SSP5-8.5.However, the difference is less pronounced in regions that experience a no-snow transition earlier, such as California.We conclude that while following a lower emissions trajectory improves the probability that transitioning to a no-snow regime will occur later, large irreducible internal climate variability could result in a transition to no-snow much sooner than the ensemble mean projects.
Another notable feature of Figure 6 is the large range of transition times within each region of the 30-ensemble member transition times.We find that in some ensemble members, the earliest transition occurs over 15 years earlier than the mean transition for many regions.For example, under the SSP5-8.5 and 90% area threshold, the first ensemble member in the Lower Colorado region transitions to no-snow in 2069 while the mean transition time of the ensemble members is not until 2086.The shape of the transition time distribution under SSP2-4.5 is also more spread out than the high emissions scenario indicating larger uncertainty in the onset of no-snow conditions.The compressed timeline is a byproduct of the rapid warming accelerating the transition to no-snow because the forcing of temperature and precipitation changes happens more quickly.Thus, internal climate variability is particularly influential in SSP2-4.5 when determining nosnow transition times, while in SSP5-8.5, the accelerated radiative forcing is the dominant effect.Furthermore, while emissions reductions improve the probability that the no-snow transition will occur later in the 21 st century, they do not guarantee a later arrival.For example, in the PNW, a quarter of the SSP2-4.5 SPEAR members transition to no-snow before the median ensemble member under SSP5-8.5.This is particularly true for regions where the transition is projected to occur earlier in the 21 st century, likely because scenario forcing is much more similar.
To assess the probability that a region becomes snow free over the next century, we examine the fraction of ensemble members that transition to no-snow before 2100.
We model the likelihood of the transition by the maximum likelihood estimator (MLE), or fraction of ensemble members that hit the transition threshold by 2100, and display these values in Table 1.By further splitting across the low and high emissions scenarios, we can model how the likelihood also changes as a function of the radiative forcing scenario.In Table 1, we see that under SSP5-8.5,A = 75% is guaranteed by 2100 across all regions.The highest threshold (A = 90%) is guaranteed only for California, while uncertainty remains for the other 4 HUC2s.Conditions by 2100 are much less severe under SSP2-4.5, with only A = 50% likely or certain for all regions, while for A = 75%, only California is very likely to transition to a low-snow regime; the other regions have low probability of doing so.For A = 90% it is unlikely that any region will have transitioned by 2100 under SSP2-4.5.
Furthermore, when we compare the likelihood of transition to no-snow conditions with the historical regionally averaged winter temperature, we find the coldest regions are least likely to transition while the warmest are most likely.For example, under SSP5-8.5 with A = 90%, the order of regions by cold to warm average winter temperature and lowest to highest transition probability is the same: UC (-5.1°C, 30%), PNW (-3.9°C, 53%), GB (-2.4°C, 70%), LC (-0.7°C, 83%), and CA (0.3°C, 100%).Like Shrestha et al.
(2021), we find that warming any region with a winter average temperature to greater than −5°C negatively impacts SWE.We also find that warmer regions are expected to experience a greater increase in no-snow conditions, emphasizing the role historical temperature has in determining not only whether a region will see decreased SWE but also the magnitude of the change.
Table 1 indicates that under either SSP2-4.5 or SSP5-8.5 we expect at least half of the historically snowy WUS to have less than 10% of its historical April SWE by 2100.
Both columns where A = 50% show greater than 80% probability for all regions, with the threshold guaranteed under SSP5-8.5.We also find that under SSP5-8.5, 4 of the 5 Western watersheds are more likely than not to cross the A = 90% no-snow threshold by 2100.Upper Colorado is the exception with only a 30% chance, likely driven by lower average winter temperatures.While severe, it is important to consider how snow-covered area and total snow volume differ.As SWE declines are dominated by losses at lower elevations that are closer to the freezing point (Mote et al., 2005;Minder, 2010), we expect the topological smoothing of SPEAR may result in an overestimate of the total amount of SWE storage lost.Therefore we expect the area-based no-snow transition to over-predict the hydrological impact of warming.

Summary
In this study, we analyze large ensembles from a coupled global climate model, SPEAR, to understand changes in SWE across the 20 th and 21 st centuries.According to SPEAR, the frequency of D2+ SD has already increased dramatically across the historical period, with an average increase across all regions of 51%.While higher than the estimate of 28% in observational data found by Huning and AghaKouchak (2020), the large amount of internal climate variability of WUS SWE within the SPEAR large ensemble indicates that chaotic climate variability could account for some of the difference.SPEAR projects even more dramatic changes to come by 2100, classifying over 35% of winter months as snow droughts under RCP2-4.5 and 60% under RCP5-8.5 compared with a normalized 9.6% across the historical period.End-of-the-century projections suggest the average monthly temperature will exceed the 93 rd and 97 th percentiles of historical conditions under RCP2-4.5 and RCP5-8.5, respectively, and were found to be the primary driver of increased D2+ SD.To understand when future conditions will deviate significantly from 'normal,' we applied the no-snow classification defined in Siirila-Woodburn et al. (2021) to each grid cell and across years for all SPEAR ensemble members, and aggregated on the HUC2 level.We found that for the most severe threshold, A = 90%, a no-snow transition was more likely than not in four out of the five WUS HUC2s, the UC region being the exception.Under RCP2-4.5,only A = 50% was likely for all regions.Furthermore, our Table 1.Probability of a snow free transition occurring before 2100 at the 3 thresholds A based on the fraction of ensemble members who transition to a no-snow regime by 2100.We show the probabilities by area threshold, 50%, 75%, and 90%, across SSP2-4.5 and SSP5-8.5 for the historically snowy portions of each of the 5 Western HUC2 regions.
finding that California is expected to transition to no snow earlier than most regions, and Upper Colorado later, is consistent with Siirila-Woodburn et al. ( 2021) who use different climate models in their analysis.These conclusions emphasizes the role of future emissions in determining the no-snow transition timing.
We found regions with higher average winter temperatures were more likely to experience a transition to no-snow.The Lower Colorado and California regions, which have the highest average winter temperatures, also had the highest probability of reaching nosnow conditions across both emissions scenarios and all area thresholds.The Pacific Northwest and Upper Colorado, the regions with the coldest average temperatures, had the smallest transition probabilities.This finding parallels Shrestha et al. (2021), who found a strong correlation between average basin temperatures and the sensitivity of the region's snow to warming.

Remarks
By using initial condition large ensembles from a state-of-the-art GCM to study SD, we can conduct a region-wide study that accounts for both radiatively forced changes and the uncertainty attributable to internal climate variability.However, while SPEAR has higher atmospheric and land resolution than most current GCMs, its 1/2°horizontal resolution is low when compared with many mountain snowpack models (Minder, 2010), which makes it unable to resolve complex mountain topography.This limitation can result in significant warm biases and less snow (Matiu & Hanzer, 2022).We expect this may make SPEAR snowpack estimates particularly sensitive to warming, and therefore likely to overestimate increases in SD.Furthermore, Hoylman et al. (2022) asserts that using timescales longer than 30 years for drought baseline climatology, as has been done here and in the vast majority of previous literature (Svoboda et al., 2002), can result in over-estimating the drought threat in a climate that is shifting towards (in this case) a less snowy state -although they argue that the reference period should take into consideration the adaptive capability of the system in question.Further work should investigate both the sensitivity of SD estimates to GCM resolution and the effect of reference climatology choice on drought severity estimation.
Here, we have assessed changes in SD across the WUS in a GCM, focusing on validating historical changes, assessing changes to the underlying climatology, and determining when WUS regions may essentially become snow-free.For this latter objective, we developed a metric, the no-snow transition time, to track both how soon a region is expected to change and the uncertainty of this timing attributable to internal climate variability.One promising avenue for future research is to examine SD changes over smaller regions, such as HUC4s, to determine the most vulnerable locations on a sub-region scale.
This would also allow further exploration of SWE's sensitivity to latitude and elevation, although at smaller watershed scales the GCM's horizontal resolution will become more problematic.Also, estimating total SWE losses and melt timing across each region would allow us to better estimate the impacts of snow droughts on the West's hydrological system.The impacts of future SDs will be felt across the entire country, both directly from the hydrological and tourism resources that consistent snowpack provides and indirectly through loss of agricultural output from summer water shortages or drifting wildfire smoke.
Understanding the probable severity and timing of when these conditions are projected to become most damaging, alongside uncertainty from emissions and internal climate variability, will allow policymakers and infrastructure planners to best prepare the West for a future with less snow.

Introduction
We include additional text and figures to complement the discussion in the main text.
August 1, 2023, 11:34pm To assess SPEAR's ability to represent historical temperature and precipitation extremes, we applied the methodology used in section 2.4 of the main text for snow drought classification to temperature and precipitation.We first aggregated SPEAR and Livneh to monthly time scales, taking the maximum high and minimum low daily temperature for each month.We do so because we assume these metrics are more likely to capture heat waves and cold extremes than an average of daily high and low temperatures.For example, a severe multi-day heat wave in January has the potential to melt snowpack quickly but might fail to show up on a 30-day average.To measure changes in meteorological drought conditions, we utilize the D2+ severity notation for dry extremes, and introduce W2+ for wet extremes.For temperatures, we introduce H2+ to indicate heat extremes and C2+ to indicate cold extremes; refer to Figure S2 for notation.As our analysis focuses on both monthly maximum and minimum temperatures, we append the subscripts "max" and "min" to distinguish between these conditions, respectively.Therefore, H2+ min represents the change in frequency of warm extremes for monthly minimum temperatures.
To assess whether changes are significant across the historical time period, we evaluate a 95% confidence interval for the SPEAR ensemble mean, assuming the underlying changes were distributed normally.If the interval does not contain zero change, then the forced component is significant in SPEAR.We present these calculations in Figure S3 which assesses changes to monthly meteorological drought (D2+), warm temperature extremes (H2+ max ), and cold temperature extremes (C2+ min and H2+ min ).Together, these panels provide a method to validate the SPEAR ensemble against changes in extreme temperature and precipitation observations.
August 1, 2023, 11:34pm : X -3 Across the historical period , we found that while SPEAR's changes in D2+ meteorological drought were not significant, several measures of temperature extremes were.Across the five HUC2 regions, we find that H2+ max extreme heat increased on average between 59% and 73%, C2+ min extreme cold decreased between 18-21%, and H2+ min extreme heat increased between 41% and 60% (Figure S3(b, c, d)).These changes indicate significant warming in both maximum and minimum temperatures and both of these trends are expected to negatively impact WUS snowpack.In both SPEAR and Livneh, extreme heat events have increased in frequency while extreme cold events have decreased on average.When we assess agreement between SPEAR and Livneh, we find that all but one Livneh observation falls within the SPEAR ensemble range, suggesting SPEAR is able to accurately reproduce changes in precipitation and temperature extremes across the historical period.Examining which trends are significant in SPEAR, we find all temperature trends to be significant while only the increase in D2+ meteorological drought in the LC region is significant.The LC saw an average increase in meteorological drought of 16% in SPEAR, while in Livneh the increase in the LC was 48%.Amongst the five HUC2 regions, the LC region increase was also the most extreme increase among Livneh meteorological drought observations (Figure S3(a)).Together, these observations may indicate the LC is drying more rapidly than other regions.When examining temperature trends, the PNW stood out as it experienced the smallest changes in extreme temperatures and was the only region to observe a decrease in meteorological drought.We assume these underlying colder, wetter conditions across the latter half of the historical period explains the decrease in D2+ SD frequency over the historical period in the PNW seen in Figure States.The wet bias is consistent with Delworth et al. (2020).By examining the maximum and minimum temperature biases we see that SPEAR has a significant cold bias for maximum temperatures across the entire WUS, while it has a systematic cold bias for minimum temperatures over mountainous regions and slight warm bias over the rest of the WUS.We expect that some of the bias can be explained by the differences in model resolution: SPEAR is on a 1/2°grid while Livneh is on a 1/16°grid.We also note that Livneh has a particularly high lapse rate of 6.5°C/1000m which may contribute some additional bias (Walton & Hall, 2018).
August 1, 2023, 11:34pm each, the corresponding Z-Score, and the probability of an event being at least as extreme.This probability captures how likely a random historical month is to be classified in that category or one that is more extreme.This table uses identical Z-Score ranges to Huning and AghaKouchak (2020) and attempts to mimic frequencies of hydrological drought given by the US Drought Monitor (Svoboda et al., 2002).
August 1, 2023, 11:34pm August 1, 2023, 11:34pm Delworth et al. (2020)  andMaher et al. (2022) demonstrate that SPEAR accurately reproduces temperature and precipitation patterns across the US and outperforms many other state-of-the-art large ensemble climate models.Delworth et al. (2020) finds that SPEAR has negligible temperature bias and a slight positive precipitation bias across the WUS.As temperature and precipitation inform snowfall, Delworth et al. (2020) lends confidence that the underlying conditions at SPEAR's approximately 50 km resolution are well-simulated.Delworth et al. (2020), Johnson et al. (2022), and Maher et al. (2022) assess SPEAR's accuracy in representing teleconnections of the El Niño-Southern Oscillation (ENSO) and the Pacific Decadal Oscillation (PDO) to North American climate.

Figure 1 .
Figure 1.Winter average SPEAR SWE deviation from Livneh (%).Red indicates regions where SPEAR has a negative SWE bias while blue indicates regions with a positive bias.The five HUC2 regions are outlined in black.

Figure 2 .
Figure 2. Comparison of SPEAR-estimated D2+ SD increases across the 1921-2011 historical period to Livneh observed increases.The SPEAR distribution is given by the box and whisker plot.The lower and upper bounds of the box correspond to the 25 th and 75 th percentiles, respectively, and points more than 1 interquartile range away from the box are denoted with a "+".The observed change in D2+ SD frequency in the Livneh dataset is marked with a red circle.

Figure 3 .
Figure 3. SPEAR D2+ SD frequencies between 1960-2100 under (a) low (SSP2-4.5)and (b) high (SSP5-8.5)emissions scenarios.The plots are masked to historically snowy regions and shaded by the percentage of winter months that the grid-cell experiences D2+ SD across a two decade period.Historically snowy regions are characterized by having an average peak SWE of at least 20mm.All 18 study decades are shown individually in Figure S4.
reveals the expected changes in D2+ SD frequency under different emissions scenarios, it does not show the impact of internal climate variability.The analysis of a large ensemble allows us to examine this effect by looking at the distribution of SD frequency across the ensemble.To study internal climate variability at the level of the entire WUS, we consider D2+ SD across the WUS in each ensemble member separately.The individual trajectories, shown in Figure4, reveal large tail probabilities that emphasize the region may experience worse drought conditions much earlier than the ensemble mean.For example, under both future warming scenarios, the ensemble mean D2+ SD frequency is reached in some ensemble members a decade or two earlier.This emphasizes that the WUS must be prepared for D2+ SD conditions well before the ensemble mean expects them.

Figure 4
Figure 4 also reveals just how dramatic the increases in D2+ SD frequency may be.SPEAR ensemble members experience an average of 5-12% D2+ SD frequency during the historical period and an average of 6.5% before 2000.However, the probability of D2+ SD is projected to be over 35% by 2050 under SSP5-8.5, while under SSP2-4.5, the 35% D2+ SD probability is projected for 2070.Examining the shape of the two curves, we see an inflection point in 2000.Before 2000, both curves do not show a noticeable increase in D2+ SD frequency while after 2000 the increase is dramatic and sustained.Under SSP2-4.5, the increase in D2+ SD has a second inflection point in 2070, where the increase in snow droughts flattens.We assume the slowdown parallels the changes in the underlying climatology discussed in 3.3.Contrary to the simulations, Livneh does not show the same uptick in drought frequency in 2000.When examining the observed changes in Figure2, we find a 53% decrease in D2+ SD frequency in the PNW.While within the SPEAR ensemble range, this decrease is far from the SPEAR ensemble mean and perhaps explains the deviation.

Figure 4 .
Figure 4.Each thin curve represents the percentage of historically snowy months classified as D2+ SDs and averaged by decade in Livneh (pink) and for each member of the three SPEAR ensembles; historical (blue), SSP2-4.5 (green), and SSP5-8.5 (yellow).The dark curves and surrounding shaded regions represent the ensemble mean and 95% confidence interval for the historical (blue), SSP2-4.5 (purple), and SSP5-8.5 (red) scenarios.

Figure 5 .
Figure 5. Temporal evolution of average temperature and precipitation anomalies with respect to the historical conditions (1921-2011).Each dot represents the average temperature and precipitation condition for historically snowy locations during winter (Oct-April) for a given decade either for all months and locations (outlined in green) or only for months classified as D2+ (outlined in gray).Each point is shaded by its average ZSWE score; thus because D2+ SD months are restricted to have a ZSWE of less than −1.3, these points average snow drought conditions are less than −1.3.Both all-month and D2+ SD-month points are surrounded by a contour which captures 95% of ensemble members.Panel (a) depicts these changes under SSP2-4.5 while (b) depicts changes under SSP5-8.5.

Figure 6 .
Figure 6.Distribution of SPEAR-simulated transition times to no-snow regimes, or T , byWestern HUC2 region, split between SSP5-8.5 and SSP2-4.5 scenarios.The 3 subplots represent the different thresholds A = 50%, 75% and 90%.Meeting a higher threshold corresponds with an increased proportion of the region experiencing perennial no-snow conditions, and implies more severe conditions.The vertical lines in the distributions represent the quantiles of the ensemble members that transition.We also include a transition time for the entire historically snowy WUS, labeling it "West-Wide".
Changes in Temperature and Precipitation Extremes

Figure S2 .
Figure S2.List of drought and temperature classification abbreviations, a text description for