Projections of North American snow from NA-CORDEX and their uncertainties, with a focus on model resolution

Snow is important for many physical, social, and economic sectors in North America. In a warming climate, the characteristics of snow will likely change in fundamental ways, therefore compelling societal need for future projections of snow. However, many stakeholders require climate change information at finer resolutions that global climate models (GCMs) can provide. The North American Coordinated Regional Downscaling Experiment (NA-CORDEX) provides an ensemble of regional climate model (RCMs) simulations at two resolutions (~ 0.5° and ~ 0.25°) designed to help serve the climate impacts and adaptation communities. This is the first study to examine the differences in end of twenty-first-century projections of snow from the NA-CORDEX RCMs and their driving GCMs. We find that the broad patterns of change are similar across RCMs and GCMs: snow cover retreats, snow mass decreases everywhere except at high latitudes, and the duration of the snow covered season decreases. Regionally, the spatial details, magnitude, percent, and uncertainty of future changes vary between the GCM and RCM ensemble but are similar between the two resolutions of the RCM ensembles. An increase in winter snow amounts at high latitudes is a robust response across all ensembles. Percent snow losses are found to be more substantial in the GCMs than the RCMs over most of North America, especially in regions with high-elevation topography. Specifically, percent snow losses decrease with increasing elevation as the model resolution becomes finer.


Introduction
Terrestrial snow plays a key role in the climate, ecology, hydrology, and economy of North America (NA). Snow's high albedo alters the surface energy budget consequently influencing both long-term climate and short-term weather (e.g., Vavrus 2007). It also provides an important habitat for wildlife that are adapted to living in snow conditions (Campbell et al. 2005; Barsugli et al. 2020). Seasonal snow accumulation is a natural reservoir for water storage, and the timing and amount of snowmelt is critical for water supply (Barnett et al. 2005), agriculture (Qin et al. 2020), and hydropower production (Markhoff and Cullen 2008). The timing and amount of snowmelt is linked to droughts (Harpold 2016) and wildfires (Westerling et al. 2006). Snow is crucial for winter transportation (Palko and Lemmen 2017) and tourism (e.g., skiing, snowmobiling, snowshoeing; Chin et al. 2018;Wobus et al. 2017) which drive regional economies (Burakowski and Magnusson, 2012). In addition to the many benefits of snow, it also contributes to a wide range of hazards including damages to roads and buildings (Palko and Lemmen 2017; Jeong and Sushama 2018), avalanches (Campbell et al. 2007), and spring flooding (Berghuijs et al. 2016).
Future changes in snow conditions that are expected to be associated with climate change will have important implications for all of these sectors. This drives a strong societal need for regional projections of snow and their uncertainties. Researchers and stakeholders alike need such information to determine the regional and local impacts of future changes in snow as well as to inform decision-makers regarding how to adapt to future changes.
Changes in snow result from the combined interactions between increasing temperatures and changing precipitation patterns. Future projections of snow for NA have been investigated using modeling techniques across multiple time and space scales including global climate models (GCMs; Räisänen 2008;Brown and Mote 2009;Mudryk et al. 2020;Krasting et al. 2013), regional climate models (RCMs; McCrary and Mearns 2019; Rhoades et al. 2018a;Rasmussen et al. 2011), variable resolution climate models (Rhoades et al. 2018b), and statistical downscaling applied to hydrologic models (Christensen and Lettenmaier 2007;Notaro et al. 2014). These studies show that increasing temperatures will dominate the climate change signal over most of NA resulting in widespread decreases in snowfall, snow cover extent and duration, and snow water equivalent (SWE). However, mid-winter snowfall and SWE may increase over the cold high latitudes and high elevations (Räisänen 2008;McCrary and Mearns 2019;Rasmussen et al. 2011).
The North American Coordinated Regional Downscaling Experiment (NA-CORDEX; Mearns et al. 2017) consists of an ensemble of regional climate projections for NA where multiple RCMs were driven with boundary conditions from multiple GCMs to produce downscaled climate projections at two resolutions (~ 0.5° and ~ 0.25°). NA-CORDEX fills a need for scientists and stakeholders who desire spatially uniform and consistent climate change data at higher resolutions than GCMs can provide and with enough models to explore uncertainty. RCM ensembles like NA-CORDEX are heavily used across multiple disciplines in order to study climate change and its impacts (Mearns et al. 2015;McGinnis and Mearns 2021). While there exists an abundance of papers examining temperature and precipitation projections in RCM ensembles, far fewer studies have looked at snow, even over Europe where CORDEX simulations have been available for longer.
Our goal here is to evaluate snow and examine future changes and their uncertainties in the NA-CORDEX RCMs and their driving GCM simulations. Since NA-COR-DEX provides downscaled simulations at two resolutions, we focus on the differences between the RCM and GCM ensembles to identify what, if any, additional information is gained by increasing resolution from GCM scales (ranging from 1.25° to 2.8°) down to 0.5° and 0.25°. To do this, we perform a side-by-side comparison of the GCMs and RCMs used in NA-CORDEX, parsing NA-CORDEX by resolution. We focus on how the spatial distribution, magnitude, and percent change of future projections and their uncertainties differ across the ensembles. Our analysis starts broadly over all of NA but then narrows down to three unique regions (Fig. 1a) to further explore regional Fig. 1 Representation of topography from the 5-min ETOPO5 (1988) dataset (a), the ensemble-average topography from the CMIP5-Driver ensemble (b), the NA-CORDEX-0.5° ensemble (c), and the NA-COR-DEX (0.25°) ensemble (d). The three sub-regions examined are outlined in (a) where (1) is the US Intermountain West, (2) is North-Central Canada, and (3) is the Northeast USA and Southeast Canada differences in model fidelity and future change. In this work, we define uncertainty as the spread across either the observations or the individual climate model ensembles.

Models
In NA-CORDEX, multiple RCMs were driven with boundary conditions from multiple GCMs that were part of phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012). Many of the RCM simulations in NA-CORDEX were performed at two resolutions (0.44° or 50 km and 0.22° or 25 km, depending on the model configuration). All RCM simulations cover at least 1951-2098, and future projections (2006-2098) follow the Representative Concentration Pathway 8.5 (RCP 8.5). In this work, the historical time period spans 1976-2005 and the end of twenty-first-century future time period spans 2070-2098.
We analyze the subset of NA-CORDEX simulations that have SWE output (Table 1). Although SWE is available from the RegCM4 NA-CORDEX simulations (na-cordex. org), unbounded snow accumulation was found to occur in many mountainous regions, so this RCM was excluded from the analysis (Supplemental Information (SI) Section S1). We split the NA-CORDEX models into two climate ensembles based on resolution. As discussed more in Section 2.3.2, we regrid the 0.44°/50 km simulations to a common 0.5° grid and the 0.22°/25 km simulations to a common 0.25° grid. Throughout the paper, we refer to these ensembles as either NA-CORDEX-0.5° (8 members) or NA-CORDEX-0.25° (11 members).
Of the 7 driving GCMs used in NA-CORDEX, daily SWE was available from 6 ( Table 1). Throughout the paper, we refer to this set of 6 GCMs as the CMIP5-Driver ensemble. We also look at broad changes in SWE from all of the models in CMIP5 with daily SWE output which we refer to as the CMIP5-All ensemble (SI Table S1, 18 members). In our analysis, the CMIP5-All ensemble is included on all time series plots, but not on the spatial maps.

Snow datasets
A major challenge for evaluating SWE in climate models is a lack of long-term, highresolution (spatial and temporal), well-vetted gridded observations (e.g., McCrary et al. 2017). The insufficiency of snow observations has led many to create gridded SWE datasets that are observationally constrained and informed by models, which we call Modeled-Observations, or MObservations (MObs). These include atmospheric and land surface reanalysis products and statistical and physical models that are constrained by in situ snow observations. Following McCrary et al. (2017), we use a multi-dataset approach to capture the uncertainty in observed snow by creating an ensemble of MObs datasets ( Table 2). All of the MObs datasets included are gridded products with 0.25° or finer resolution, at least 5 years of data between 1981 and 2010, and cover CONUS or North America. These datasets have considerable uncertainties related to sparse observational networks and surface meteorological forcing, satellite retrieval algorithms, and the use of models that must parameterize snow processes. While this ensemble will not capture the full uncertainty in snow observations, it serves as a reference dataset in which to assess the climate models used in this study.
Optical satellite products can be used to identify the presence of snow. To evaluate snow cover metrics, we include the Interactive Multisensor Snow and Ice Mapping System (IMS) 24 km daily snow cover dataset (U.S. National Ice Center 2008) in combination with snow cover estimated from the SWE MObs.

Calculation of snow cover
Similar to McCrary and Mearns (2019), we calculate snow cover from SWE by applying a 5 mm threshold to daily SWE fields from the MObs and the simulations to produce a binary yes-no snow cover field. Snow cover extent (SCE) is calculated by first averaging this binary field over each month to produce monthly snow cover fraction (SCF) at each gridbox and then summing over NA. The daily binary snow cover field is also used to calculate snow cover duration (SCD), defined as the number of days with snow on the ground. SCE, SCF, and SCD are also calculated from the satellite IMS dataset which provides estimates of yes-no snow cover.  1998Brown and Brasnett (2010 Snow-on-ice (sea ice or land ice) is a complex process that is not well simulated by many climate models. We remove any points which may be ice covered or strongly influenced by land/sea ice (see SI Section S3 for a description of how these points were removed).

Ensemble analysis
As regridding snow onto different grids can greatly impact mass budgets, regionally averaged time series plots are calculated on the native grids of the climate models. However, since each model uses a different grid, for ensemble mean spatial calculations, we regridded the MObs and models to common grids using conservative remapping. The CMIP5 models have been regridded to a 1.5° grid, and the NA-CORDEX ensembles have been regridded to 0.5° and 0.25°.

Evaluation of climate models
In the following section, we compare the historical simulations from the CMIP5 and NA-CORDEX ensembles with the MObs ensemble.

Snow cover extent
The annual cycle of SCE from the MObs and climate models is shown in Fig. 2. Only IMS and ERA5-land have daily data covering all of NA. The timing of the annual cycle of SCE in the MObs follows each other closely. SCE is near zero in July-August, starts to increase in September, reaches a maximum in January, and declines throughout the spring and early summer. Between December and April, IMS has higher values than ERA5-land; otherwise, the two datasets are closely matched the rest of the year. Spatially, the largest differences in January snow cover fraction (SCF) occur in the Central Plains and Great Basin (SI Fig. S4.) Ensemble mean SCE for the CMIP5-Driver and CMIP5-All ensembles is similar to each other and is lower than both MObs ( Fig. 2; individual models results SI Fig. S8). The spread in historic SCE in the both CMIP5 ensembles is much larger than the MObs spread. In both NA-CORDEX ensembles, the ensemble mean annual cycle of SCE is on the higher end of the MObs, following the IMS dataset almost exactly, with half the individual RCMs slightly overestimating SCE compared to IMS (SI Fig. S8). The spread in the NA-COR-DEX-0.5° ensemble is larger than the NA-CORDEX-0.25° ensemble, primarily due to the HIRHAM5 EC-EARTH simulation being a low outlier which is only present in the 0.5° ensemble (SI Fig S8). The spatial distribution of January snow cover fraction (SCF) highlights where differences in simulated SCE arise (SI Fig. S5-S7).

Climatological SWE
Maps of ensemble mean annual maximum monthly SWE (AM-SWE) from the MObs and climate simulations are shown in Fig. 3. For the MObs, ensemble mean AM-SWE is calculated across all available datasets at each gridbox ( Fig. 3a-c). Although the spatial patterns of AM-SWE are similar across the three resolutions, distinct topographic features such as the Sierra Nevada's in California, the Cascade Range in the Pacific Northwest, and the Rocky Mountains deteriorate with decreasing spatial resolution. Results from the individual MObs and MObs ensemble statistics highlight the uncertainty across the MObs ensemble (SI Figs. S9-S12). In this study, we use them as a reference to qualitatively evaluate the models.
The spatial patterns of AM-SWE in the CMIP5-Driver and NA-CORDEX ensembles are broadly similar to the MObs ensemble ( Fig. 3; individual model results SI Figs. S13-15). Compared to the MObs ensemble mean, the CMIP5-Driver GCMs underestimate SWE in the mountains and overestimate SWE in the lower elevation regions of western NA (Fig. 3d, g, j; SI Fig. S16). As the mountains in these coarse GCMs are relatively smooth and low (Fig. 1b), there is limited orographic enhancement of precipitation in the mountains resulting in too much moisture penetrating inland, likely contributing to positive biases east of the mountain ranges in the west (e.g., Rasmussen et al. 2011). The smoothed topography in the GCMs also Snow Cover Extent (10 6 km 2 )

NA-CORDEX CMIP5 % Decrease in Cover Extent
The annual cycle of monthly mean NA SCE from the MObs (dashed lines on a and b) and the historical climate simulations from the two CMIP5 ensembles (a) and two NA-CORDEX ensembles (b) examined in this study. Also shown is the annual cycle of the percent decrease in NA SCE projected for the end-of-century from the CMIP5 ensembles (c) and the NA-CORDEX ensembles (d). The spread of each ensemble is displayed with colored shading. The average of each ensemble is plotted with a corresponding solid line. Percent decreases in SCE for July-September have been masked as they are skewed by small number division. NA SCE has been calculated using the native grid of all models and datasets

Fig. 3
Maps of the average annual monthly maximum SWE (AM-SWE) from the MObs ensemble mean which has been regridded to the common 1.5°, 0.5°, and 0.25° resolution (a-c). Ensemble mean AM-SWE from the historical time period for the three model ensembles (d-f). Also shown are the magnitude (g-i) and percent (j-l) of the simulated bias (model -MObs) of AM-SWE. The MObs ensemble mean is calculated independently at each gridbox using datasets with available data (see SI Figs. S9-S11) results in artificially high terrain in the interior valleys of western NA, which likely contributes to cold biases and excess SWE accumulation and suppressed ablation. In the CMIP5-Driver ensemble, negative biases also occur over most of Canada. AM-SWE in the CMIP5-Driver ensemble falls within the range of the MObs over much of CONUS; however, biases fall outside the MOBs range near the western mountains and Northeast Canada (SI Fig. S19). AM-SWE in the two NA-CORDEX ensembles differs from the MObs ensemble in similar ways (Fig. 3h, i, k, l; individual model results SI Figs. S14-S15, S17-S18). The only real difference is that the spatial details of SWE patterns are finer with increasing resolution. In both ensembles, positive biases dominate over the domain, with negative biases near some mountains and across central Canada. Over the western half of the domain, AM-SWE is greatly overestimated on the western side of the mountains and underestimated just east of the highest peaks of the mountains (SI Fig. S20). Relative to the MObs ensemble mean, the magnitude of positive SWE biases are larger in the mountains than the lower elevation regions; however, percent SWE biases are much larger at lower altitudes. The percent bias figures ( Fig. 3j-l) highlight regions where the climate models simulate snow, but the MObs do not. When compared to the range of the MObs ensemble, AM-SWE values in both RCM ensembles are greater than the MObs on the western side of Mountain ranges, in the Central USA, and Northern Canada/Alaska (SI Fig. 19).
The similarities in the AM-SWE bias patterns in the RCMs suggest that biases in largescale forcing and RCM configuration/parameterizations play a similar role in the simulation of SWE at both resolutions. As the mountains in the NA-CORDEX simulations are higher than the GCMs, orographic precipitation is larger in the RCMs due to enhanced lifting (see Mahoney et al. 2021) resulting in higher SWE values. However, winter precipitation in NA-CORDEX far exceeds observations in the RCMs (Mahoney et al. 2021). This is possibly because even at 0.25°, convection is insufficiently resolved and convective parameterizations play a large role in precipitation biases (Hughes et al. 2021). The bias in precipitation likely translates to a positive bias in SWE; however, SWE biases will also be linked to temperatures the evolution of snowpack in the RCMs (McCrary et al. 2017).

Snow cover duration
The length of the snow covered season or SCD may also change in the future. In the MObs, SCD increases with latitude and elevation ( Fig. 4a-f; individual MObs in SI Figs. S21-S23). The broad spatial patterns of SCD are similar across the different resolutions, but the details are lost as resolution coarsens. Much like winter SWE values, SCD is reduced in the mountains when aggregated to coarse scales (Fig. 4a).
The CMIP5-Driver GCMs underestimate SCD in the mountains and most of the eastern half of the domain, although SCD is positively biased at low-elevation regions in the western half of the domain ( Fig. 4g; individual models results SI Fig S24). Both RCM ensembles overestimate SCD over most of the domain, with larger positive biases at lower elevation in the Great Basin and east of the Middle and Southern Rockies and negative biases over northcentral Canada (Fig. 4h, i; individual model results SI Figs. S25-S26).

Snow cover extent
SCE is projected to decrease in all months of the year in all of the models examined ( Fig. 2c-d; individual model results SI Fig. S8). The largest percent losses are projected to occur in October, May, and June when snow cover is marginal in the historic climate period. Average SCE losses are larger in both CMIP5 ensembles than both NA-CORDEX ensembles. Although ensemble mean changes in SCE are similar for the CMIP5-Driver and CMIP5-All model ensembles, the uncertainty (measured here as the multi-model spread) is considerably larger in the CMIP5-All ensemble. This may indicate that our subset of CMIP5 models does not capture the full potential of the combination of temperature and precipitation changes that drive changes in snow. The NA-CORDEX 0.5° and 0.25° ensembles have nearly identical projections for SCE loss, both of which are smaller than the CMIP5-Driver ensemble mean. The uncertainty in future SCE changes is slightly larger in the NA-CORDEX ensembles than the CMIP5-Driver ensemble during October-February but slightly smaller in March-July.  Räisänen 2008). Absolute losses are larger in the mountains over the western and eastern portions of the domain, while percent losses are higher at low latitudes and lower elevations. Total snow losses are projected in all of the ensembles along the southern edge of the snow boundary ( Fig. 5d-f). The individual models in all of the ensembles also show total losses along the southern snow boundary (SI Figs. S28, S30, S32).

Annual maximum SWE
While the three ensembles tell broadly the same story, details emerge in the RCMs that are not found in the GCMs. Focusing on percent change, as it reduces the influence of simulated difference in baseline historical AM-SWE amounts, it is apparent that percent losses are generally larger in the CMIP5-Driver ensemble than in both of the NA-CORDEX RCM ensembles, especially in regions of complex topography. By the end of the twenty-first century, the CMIP5-Driver ensemble projects that 50.5% of NA will experience AM-SWE losses of greater than 50% and that 15.9% of NA will experience losses of greater than 90%. While the NA-CORDEX-0.5° ensemble projects greater than 50% losses over 38.9% of NA and 90% losses over 9.16% of the domain, the NA-CORDEX-0.25° ensemble projects greater than 50% losses over 36.4% of NA and 90% losses over 8.8%. The reduced losses in the RCMs is partially due to the fact that they have more higher-elevation mountain points ( Fig. 1 and Section 3.1.1) where temperatures can remain below freezing during winter. But also, percent change is affected by the historical snow amount where for the same magnitude loss, smaller percentage loss will be found if there is more SWE in the historical baseline climate.

Snow cover duration
Along with losses in SCE and SWE, the duration of the snow covered season is also projected to decrease (Fig. 4j-l; individual model results SI Figs S33-S35). The largest decreases in SCD are found over the mountains in the western half of the domain, over Southwestern Alaska, and the eastern half of Newfoundland/Labrador and New England. While the broad spatial patterns of changes in SCD are similar across the ensembles, again the spatial details are lost in the GCMs. For example, in the western mountains, the RCMs demonstrate that SCD will decrease more at higher elevations than lower elevations. While AM-SWE is projected to increase at high latitudes (Fig. 5) and at high elevations in a few of the models (SI Figs. S28, S30, S32), SCD is found to decrease everywhere, indicating while AM-SWE may increase in some locations, the snow covered season will still contract.

Regional changes
In the previous section, we explored the continental-scale patterns of changes in snow conditions over NA. While important from a large-scale climate perspective, most researchers and stakeholders often want to know what will happen over smaller-scale regions. Here, we zoom in on three unique climate regions, to explore more deeply how resolution influences future projections of SWE. These regions (Fig. 1a)

US Intermountain West (IMW)
The IMW region is large and contains portions of the Middle and Southern Rocky Mountains and the Great Basin (Fig. 1a). We chose this region because of its complex topography including high-elevation mountains where seasonal snowpacks and spring snowmelt are critical for water supply, ecosystem health, forest fire risk, and recreation.
First, we examine the annual cycle of monthly averaged total snow mass (SM) for the region (Fig. 6a, d; Fig. 6 The first two rows show the annual cycle of SM from the CMIP5-All and CMIP5-Driver ensembles (a-c), the NA-CORDEX ensembles (d-f), and available MObs, for the IMW (a, d), NC-Canada (b, e), and Northeast (c, f) regions. The last two rows show the percent decrease SM projected for the end of the century for the CMIP5 (h-j) and NA-CORDEX ensembles (k-m). All regional averages are calculated on each model's native grid (excluding August and September) and the timing of peak SM (showing either a February or March maximum). There is a clear separation between the 4 highest MObs and the 3 lowest MObs over the region. Lundquist et al. (2020) demonstrated that most snow reanalysis datasets underestimate SWE in the mountains, so our judgment here is that datasets with higher SM values are more realistic.
The spread in historical SM in both CMIP5 ensembles is larger than the spread of the MObs with a few GCMs greatly underestimating peak SM. Ensemble mean SM in the GCMs falls within the middle of the MObs range. The NA-CORDEX ensembles have considerably more snow than their driving GCMs and about half the spread during peak months. Between December and April, ensemble mean SM in the RCMs falls within the 4 highest MObs, but spring snowmelt occurs more rapidly in the RCMs than those same datasets.
To examine future changes in regional SM, we again primarily focus on percent changes (Fig. 6h, k); however, the annual cycle of future snow and the magnitude of snow changes are shown in SI Fig. S37. Average IMW SM is projected to decrease for all months of the year in all of the models, except for one CMIP5-All ensemble member (Fig. 6h, k). Regional percent SM losses are smaller in the NA-CORDEX ensembles than the CMIP5 ensembles. For example, in March (around the timing of peak snow amounts), SM is projected to decrease by 84.1% in CMIP5-All, 76.8 in CMIP5-Driver, 58.2% in NA-COR-DEX-0.5°, and 52.4% in NA-CORDEX-0.25°. The uncertainty in future losses is also much larger in the GCM ensembles than the RCM ensembles. For example, in March, the spread in future change is 2.05 times larger in the CMIP5-All ensemble than the NA-COR-DEX-0.25° ensemble. In most of the models, the largest absolute SM losses occur when historical maximum SM occurs, and the timing of peak SM occurs 1 month earlier in the future (SI Fig. S36-S37). Figure 7 examines the spatial distribution of historical and future changes in AM-SWE over the IMW. Terrain plays a large role in determining precipitation, snowfall, and SWE patterns in the IMW. In the coarse CMIP5-Driver GCMs, topography is fairly smooth and the Rocky Mountains are captured as one mountain feature (Fig. 7a) although this varies slightly with GCM (see SI Fig. S38). With increasing resolution, individual mountain ranges begin to appear and become more distinct (Fig. 7b-c). Comparisons of the distribution of elevation over the IMW (Fig. 8a) highlights that all of the ensembles, but in particular the CMIP5-Driver ensemble, have too many grid points at mid-elevations between 1500 and 2500 m and too few points at lower elevation valleys and higher-elevation mountains.
In both the historical and future simulations, AM-SWE generally increases with elevation and latitude (Fig. 7d-i; individual models SI Figs. S39-S41). As with topography, distinct spatial patterns of SWE become more refined with increasing resolution. Comparison to the MObs AM-SWE (SI Fig. 20) reveals unrealistic patterns in the CMIP5-Driver ensemble related to the GCMs' unrealistic underlying terrain, as the Middle and Southern Rockies should have distinct peaks in SWE.
By the end of the century, ensemble mean AM-SWE is projected to decrease at all grid points over the IMW region in all of the ensembles (Fig. 7g- (Future-Historical)

Middle Rockies
SouthernRockies G r e a t B a s in Wyoming Basin Fig. 7 Maps over the IMW region of the ensemble mean topography (a-c), ensemble mean historic AM-SWE (d-f), ensemble mean future AM-SWE (g-i), and the magnitude (j-l), and percent change (futurehistoric) in AM-SWE from the three model ensembles higher topography. Although temperatures are projected to increase everywhere over the IMW, we might expect higher-elevation SWE to be partially preserved as temperatures can still remain below freezing during the heart of winter.
To further explore the relationship between elevation and snow over the IMW, we bin annual maximum snow mass (AM-SM) over the IMW by elevation (Fig. 8). To calculate AM-SM, we take the climatological AM-SWE (e.g., Fig. 7) from each model and the four highest MObs in Fig. 6a and calculate the mass of snow stored in each elevation bin for each dataset.
None of the datasets have grid points below 500 m, the CMIP5-Driver models have no grid points above 3000 m, and only the MObs and the NA-CORDEX-0.25° ensemble have grid points above 3500 m (Fig. 8a). Most of the observed AM-SM occurs between 2000 and 3000 m (Fig. 8b), with lower values at higher and lower elevations. In the climate models, most of the AM-SM occurs at lower elevations, between 1500 and 2500 m. On average, the RCMs overestimate AM-SM at mid elevations (1500-2500 m) and underestimate AM-SM above 2500 m, with larger biases in the NA-CORDEX-0.5° ensemble. The GCMs skew toward negative AM-SM biases, except between 1500 and 2000 m, where ensemble mean values are slightly higher than observed. There is also large uncertainty in historical AM-SM in the CMIP5-Driver ensemble between 1500 and 2500 m, likely linked with the large spread horizontal resolution and topography (SI Fig S38 and Fig. 8a).
In terms of magnitude, the largest losses in AM-SM for all the models occur between 1500 and 2500 m, corresponding with the elevation bins with the largest historical AM-SM. Percent losses in the CMIP5-Driver GCMs are also highest between 1500 and 2500 m, but in the RCMs, percent losses are highest at lower elevations (1000-2000 m).
In all the models, percent losses steadily decrease above 2000 m, where temperatures will remain below freezing more frequently than the lower elevation bins. Since vast majority of grid points in the GCMs occur between 1500 and 2500 m where the largest AM-SM losses (magnitude and percent) occur and the GCMs have no grid points at the higherelevation bins, this supports the idea that reduced relative snow losses occur in the RCMs The spread for each ensemble shows the minimum and maximum values calculated for each bin. In (b-d), the spread of each ensemble is shown, and the white space between the bars represents the ensemble mean change because they have higher mountain elevations, which help to buffer snow losses. However, percent AM-SM losses are higher in the GCMs at all elevation bins, which suggests differences are not solely due to elevation, but also related to baseline SM amounts (Fig. 8b) and the magnitude of total SM loss that occurs in each elevation bin (Fig. 8c).

North-Central Canada (NC-Canada)
We next examine changes over North-Central Canada (NC-Canada, Fig. 9) as this is the one area of NA where the climate model ensembles project potential increases in winter snowpacks. In this region, increases in SWE could have important implications for winter transportation, wildlife, and indigenous populations. Here, we examine how robust these projected changes are.
Only three members of MObs ensemble have SWE data for NC-Canada, and there are large difference between them (Fig. 6b, e; individual model results SI Fig. S48). Also, given the very few in situ observations (snow or surface meteorology) over the region (e.g., Mekis et al. 2018), observed snow amounts are highly uncertain. Over NC-Canada, much of the year is snow covered, with only a short snow-free period in summer. Snow accumulates between October and March/April and declines quickly between March/April and July, with the timing of peak SM varying across the datasets. Over the region, ERA5-land has the highest SM values compared to CMC and GlobSnow. Values are likely underestimated in GlobSnow, as the mountains in the southwest part of the domain are masked (see SI Fig. S11). SM in all four of the model ensembles lie on the upper end of the MObs estimates (following ERA5-land). The spread is higher in the GCM simulations than the RCM simulations, and the RCMs tend to have more snow than the GCMs.
In the future, losses are projected in all of the models examined between May and November. However, from December to April, most of the RCMs and many of the GCMs project up to 20% increases in SM for the region. Snow increases in this region are likely associated with increases in the amount of moisture in the atmosphere associated with warming temperatures which result in increases in precipitation and snowfall, as winter temperatures remain well below freezing (Räisänen 2008). In the future, the timing of peak SM occurs 1 month earlier in most of the models, and the largest magnitude decreases in SM occur in May for all ensembles (SI Figs. S48-S49).
Spatially, most of the increases in AM-SWE occur over the northern and eastern portions of the domain (Fig. 9) although this is model dependent (SI Figs. S50-S58). At individual points, ensemble mean AM-SWE increases by 1-50 mm or 1-20%. The areal extent over which AM-SWE is projected to increase is largest in the NA-CORDEX-0.5° ensemble and smallest in the CMIP5-DRIVER ensemble. All of the models in all of the ensembles project increase in SWE along the northern edge of continental Nanavut (Fig. 9p-r); however, model agreement regarding where increases in SWE may occur is lower for the west and south of the domain.

Northeast USA and Southeast Canada (Northeast)
While the vast majority of previous studies have examined future projections for snow over western NA, snow is also important over the Northeast (Fig. 1a). In this region, heavy snowfall and snow loads are hazards for transportation and building infrastructures and snowmelt plays a key role in spring flooding. The Northeast region is also home to over 180 ski resorts (https:// www. skice ntral. com/). Lake effect snow may play a role in driving  Fig. 9 Maps over the NC-Canada region of the ensemble mean topography (a-c), ensemble mean historic AM-SWE (d-f), ensemble mean future AM-SWE (g-i), the magnitude (j-l) and percent (m-o) change (future-historic) in AM-SWE, and the percent of models in each ensemble that agree AM-SWE will increase in the future (p-r) SWE amounts in this region, and the ability of the models to capture lake effect snow will depend on if the models use a lake model for the Great Lakes or interpolate lake temperatures from sea surface temperatures (see RCM Characteristics at https:// na-cordex. org/).
As in the previously examined regions, the spread across the MObs ensemble is substantial over the Northeast (Fig. 6c). There is also disagreement on whether SM peaks in February or March over the region. Compared to the MObs, the spread in both the CMIP5 ensembles is larger than the estimated observed spread, with the ensemble mean values falling in the middle of the MObs (Fig. 6c, SI Figs. S59-S60). The spread across the NA-CORDEX RCMs is smaller than the MObs, with the RCMs following the upper end of the MObs (Fig. 6f). It is possible that all of the MObs underestimate SWE in the region, as observations are limited.
Substantial losses are projected for the region in all of the models (Fig. 6j-m; individual model results SI Figs. S59-S60). Percent losses are largest October-November and April-May and smaller in the middle of winter. However, even in winter, most models project that the region will lose more than 40% of its total SM. The uncertainty range for SM losses is twice as large in the CMIP5-All ensemble than the CMIP5-Driver ensemble, again suggesting CMIP5-Driver models may not capture the full range of climate possibilities. Regional scale losses are nearly identical in the two NA-CORDEX ensembles. The uncertainty of the change in the RCMs is smaller than the CMIP5-All ensemble but larger than the CMIP5-Driver ensemble.
At first glance, the spatial representation of winter SWE in the CMIP5-Driver and NA-CORDEX RCM ensembles appears to be very similar, with higher values in the Northeast portion of the region and lower values to the Southwest (Fig. 10d-f; SI Figs. 61-69). However, an examination of the spatial details and comparison with topography highlights that even though topographic variations are less extreme in this region than in the IMW, mountains still play a role in driving SWE patterns (e.g., Adirondack Mountains in New York).
The broad spatial patterns of SWE changes for the end of the century are similar across the ensembles, with larger total losses to the northeast and smaller total losses to the south and larger percent losses over the southern portion of the domain and smaller relative losses over the north. However, at closer inspection, we see that percent snow losses are smaller with increasing elevation and higher resolution, indicating topography also dampens snow losses here. While beyond the scope of this study past RCM studies suggest warming of the Great Lakes may result in increased lake-effect snowfall possibly mitigating some snow losses in the future (e.g., Notaro et al. 2015).

Summary and Discussion
RCM ensembles like NA-CORDEX are widely used by scientists and stakeholders across multiple fields. While a plethora of studies have examined temperature and precipitation changes, far fewer have examined critical variables such as snow. In this study, we performed a side-by-side comparison of historical and future snow over NA between the NA-CORDEX dynamically downscaled 0.5° and 0.25° RCM simulations and their driving GCMs (1.25°-2.8°). The primary goals of this study were to evaluate model performance and examine how end-of-century projections for snow differ between the different resolution ensembles. To evaluate model performance, we used an ensemble of observationally constrained SWE datasets. We demonstrate that the uncertainty in gridded snow datasets is large, even across datasets with high resolutions. This uncertainty is associated with difficulties in measuring snow and is a major challenge for the snow science community.
In their historical climate simulations, the CMIP5-Driver and CMIP5-All ensembles underestimate NA, while both NA-CORDEX ensembles tend to follow the higher end of the MObs. Simulated biases in SCE can have an impact on the radiation budget, thereby influencing surface temperatures and weather patterns (Vavrus 2007). On average, the CMIP5-Driver ensemble underestimates AM-SWE and SCD over eastern Canada and at high elevations but overestimate them everywhere else. In both RCM ensembles, AM-SWE and SCD biases are positive everywhere except the highest elevations and across tracts of central Canada. While the 0.25° simulations have greater spatial details and higher mountains than the 0.5° models, the spatial pattern of SWE and the magnitude of the biases are similar between the ensembles. Over the three regions examined, the spread in historical SM is greater across the GCM ensembles than the RCM ensembles, especially in the IMW and Northeast regions.
End-of-century projections for snow over NA are broadly similar across the ensembles considered. In all ensembles, SCE, AM-SWE, and SCD are projected to decrease over most of NA, with the exception of increases in AM-SWE at high latitudes. In terms of magnitude, the largest losses in AM-SWE and SCD occur over the mountains in the western half of the domain and over coastal-eastern Canada. However, percent losses in AM-SWE are largest at low-elevation and low-latitude regions. Comparison of these ensembles shows that in terms of percent change, which can be more useful in the application of climate model data to climate change impacts, the CMIP5 GCMs tend to project a more severe picture of total snow loss for NA. For example, the CMIP5-Driver ensemble projects that just over half of NA will experience greater than 50% losses in AM-SWE while the RCMs project that only 36-38% of NA will experience greater than 50% losses in AM-SWE. These differences are likely largely related to the poor representation of topography in the GCMs but are also related to differences in baseline snow amounts (see more below).
While the large-scale picture of future changes in snow are similar between the ensembles, zooming in on individual regions helps highlight where differences between the ensembles occur. Over the IMW and Northeast regions, percent SM losses are larger in the GCMs than the RCMs, while over NC-Canada, percent increases in winter SM are smaller in the GCMs. The uncertainty in these future changes is larger in the CMIP5-Driver GCMs than the RCM ensembles over the IMW and NC-Canada but smaller in the Northeast region. While SM is projected to decrease during all months of the year in all of the models over the IMW and Northeast region, 70% of the models examined show winter SM increasing over NC-Canada. The largest differences across the ensembles are found over the IMW region where topography plays a large role driving snowfall and SWE amounts through orographic enhancement of precipitation and lower temperatures. As the GCMs oversample low-to-middle elevations and undersample the higher elevations, historical SM is underrepresented at most elevations, percent snow losses are larger in the GCMs than the RCMs at most elevations, and the GCMs have no information about snow at the higher elevations. Our results suggest that a more accurate representation of snow (especially at high elevations) allows for the buffering of snow losses, which we do not see in coarse models. In contrast to the IMW, over NC-Canada, the GCM and RCM ensembles also show the greatest agreement in historical SM and percent SM changes in this region, which we suspect could be due to the lack of significant topographic features in the region.
Overall, while we find interesting differences in the specific regional details between the CMIP5-Driver GCMs and the NA-CORDEX RCMs, we do not see significant differences between the two NA-CORDEX ensembles. The largest differences between the GCMs and RCMs are found in regions of complex topography, but even over the IMW, the two RCM ensembles have very similar climate change responses. So while the spatial details of snow are more refined in the 0.25° ensemble, the overall impact to regional snow is small. However, as shown in Walton et al. (2021), fine-scale details associated with snow in the mountainous may be important for end-users who statistically downscale temperature from climate models, as incorrectly capturing snow to no-snow transitions in the future can result in the incorrect amplification of surface temperatures associated with the snow albedo feedback.
While these results are similar to other studies which have examined snow over NA (e.g., McCrary and Mearns 2019; Rasmussen et al. 2011), our study highlights that the severity of future changes in snow, their uncertainties, and regional details are a function of the size and configuration of the model simulations/ensembles examined and the resolution of the simulations. The result that the GCMs tend to project a more severe picture of relative snow losses, in part because they have fewer high-elevation points, is important to remember when considering studies such as Diffenbaugh et al. (2012), which used the CMIP5 models to assess hydrologic extremes and water availability over regions of the Northern Hemisphere, including the western USA.
There are also a few notable limitations to this study. First, while our focus has been on how increases in resolution impacts the representation of historical and future snow over NA and their uncertainties, the NA-CORDEX-0.5° and NA-CORDEX-0.25° ensembles consist of different combinations of RCM/GCM pairs. These differences in ensemble configuration may also contribute to the differences we found in this study. As discussed in McGinnis and Mearns (2021), funding was extremely limited for NA-CORDEX, and the choice of simulations included in the experiment was opportunistic and required leveraging other modeling activities (McGinnis and Mearns 2021). While the archive consists of simulations with available SWE from 5 RCMs driven with boundary conditions from 7 GCMs with simulations performed at 2 resolutions (Table 2), the simulation matrix itself is both sparse and unbalanced, limiting our ability to dive deeply into the different roles the choice of RCM, GCM, or resolution have on climate change uncertainty. In this work, we chose to include all available ensemble members to highlight the differences in the available datasets that end users may consider. In SI Section S17, we compare results from using the full ensemble with the 7 simulations that have matching RCM/GCM pairs and both resolutions. We find the uncertainty in the historical simulations to be lower in the smaller subset of models, but that the projections of future change and their uncertainties for all snow variables are similar between the full and subset ensembles.
Second, while the NA-CORDEX simulations are higher resolution than their GCM counterparts, they are still relatively coarse for capturing precipitation, snowfall, and SWE especially in topographically complex areas. Many of the sectors discussed in the "Introduction" require very high-resolution data to study the impacts of future changes in snow. Statistical downscaling techniques are often employed to get at very high-resolution climate information, but this can break the coupling between the atmosphere and the land surface leading to inconsistent results especially in snow dominated regions (Walton et al. 2021).
Past RCM studies have found resolutions of 4-6 km to be necessary to match in situ point observations of precipitation, snowfall, and SWE (Garvert et al. 2007; Rasmussen et al. 2011). The argument for this is that terrain-induced convection and local air circulation patterns associated with smaller ridges and valleys that are important for snowfall patterns are better resolved, and surface temperatures are better represented. Wind redistribution of snow is also important, which is also not captured in many models (Musselman et al. 2015). Many high-resolution modeling studies have examined changes in snow over regions of NA (e.g., Sun et al. 2019;Rasmussen et al. 2011;Musselman et al. 2017). While these studies have been able to look at detailed process level changes that are important, they have been limited in either domain size or by the use of only one RCM or one GCM, limiting the examination of uncertainty. As computing power, storage, and analysis of big data continues to advance, we expect to see the creation of larger ensemble convection permitting simulations (CPSs) over larger domains. The coordination of regional CPSs is ongoing over Europe in one CORDEX Flagship Pilot Study (Coppola et al. 2021), but such studies are not being coordinated over NA yet. Until that time, NA-CORDEX fills a need in the community as it provides spatially uniform, higher-resolution simulations with enough model diversity to explore uncertainty while covering a large enough domain to be useful for many regional interests and are adequate for efforts such as the US National Climate Assessment and the IPCC.