Model Spread and Bias in Snowpack in the American Mountain WestTravis Aerenson,a Daniel McCoy,aGregory Elsaesser,bca University of Wyoming, Laramie, WYb Columbia University, New York, NYC NASA-GISS, New York, NYCorresponding author : Travis Aerenson, taerenso@uwyo.eduABSTRACTAs the climate is rapidly changing, we are faced with questions about which regions will become drier and which wetter. In many regions such as the American Mountain West the amount of moisture available for agriculture, drinking water, and industry is dependent on the depth of the wintertime snowpack that melts throughout the spring and summer, filling the downstream rivers and reservoirs. We provide an overview on how snowpack in the Mountain West relates to interannual climate variability in observations and reanalysis and how snowpack in reanalysis has changed during the past century of warming. We then compare these changes with the historical simulations from an ensemble of CMIP6 Global Climate Models (GCMs) to determine how skillful they are at predicting snowpack in the Mountain West, which models are the most skillful, and how we can constrain future predictions by selecting only the models with the best simulation of historical snowpack. We perform a decomposition of Snow Water Equivalent (SWE) trends and variability in reanalysis and CMIP6 models to determine if the model deficiencies in their SWE simulation are due to issues relating to their simulation of precipitation, fraction of precipitation occurring as snow, or the melt-rate of the snowpack. Finally, we construct machine learning emulator models trained on reanalysis that are useful diagnostic tools for the snowpack biases and show that most of the spread in the GCMs’ prediction of historical snowpack change is due to spread in their simulation of regional surface temperature and moisture convergence.1. IntroductionAs the Earth is warming due to anthropogenic emissions of heat-trapping gasses, we anticipate numerous changes to global circulations and local water cycle that will alter regional moisture availability and may cause regional water stress (e.g., Held & Soden, 2006; Pendergrass & Hartmann, 2014; Reed et al., 2023). Hence, it is important to determine if a region is expected to become more moist or arid, and if any adaptive moisture storage methods may be necessary to provide ample water for drinking, agriculture, and industry in a warmer future. In this study we examine the trends and variability in moisture availability in the American Mountain West, defined by the region bounded by 37° to 48° latitude and 247° to 256° longitude, which roughly encompasses the mountainous regions of Colorado, Utah, Wyoming, and Montana, with 3 primary questions in mind: 1) how does large scale variability in the atmospheric and oceanic circulations impact regional moisture availability in the Mountain West, 2) how has moisture in this region changed during the historical global warming during the 20th century, and 3) what aspects of GCMs could be improved to most effectively increase the fidelity of their predictions of regional water stress?The American Mountain West is a continental climate with complex terrain, where most of the run-off throughout the year originally precipitates as wintertime snow (Li et al., 2017). The seasonal snowpack acts as a reservoir, allowing the wintertime precipitation to run-off throughout the spring and summer. Understanding how precipitation and snow-storage is changing is important for future water management plans (Middelkoop et al., 2001), so we focus our analysis on changes in snow water equivalent (hereafter SWE). As the climate warms, one might anticipate a decrease in SWE as it becomes too warm for the surface to support seasonal snowpack for a greater portion of the year. The influence of changing precipitation patterns on SWE also needs to be considered, where there may be increase in precipitation (and snowfall) that can increase the SWE in a warming climate. Due to the complexity of these drivers of future SWE, it is important to thoroughly evaluate the causes of previous trends in SWE and model predictions of future SWE changes with global warming. Long-term snow-monitoring stations have reported widespread snow declines throughout the western United States, including maritime climates such as the Pacific Northwest and California (Chavarria & Gutzler, 2018; Clow, 2010; Fritze et al., 2011; Harpold et al., 2012; Knowles et al., 2006; Mote, 2003), however the observed changes in the intercontinental regions such as the Mountain West are less consistent compared with maritime regions (Harpold et al., 2012; Regonda et al., 2005). In this study we use a combination of in situ and satellite observations, reanalysis, and state-of-the-art global climate model (hereafter GCM) simulations to determine if simulations of historical SWE in the Mountain West region are consistent with observations. We then use a machine learning model trained on historical reanalysis to determine which aspects of the GCMs contribute most to their bias and spread in SWE predictions.Throughout this paper we use the relationship between SWE in the Mountain West and the phase of the El Niño Southern Oscillation (ENSO) to benchmark if models are correctly simulating the relationship between SWE in the Mountain West and large-scale climate variability. ENSO is an interannual oscillation of sea-surface temperature in the equatorial eastern pacific and is the largest mode of interannual variability in the historical sea-surface temperature record (Battisti et al., 2019; Chen & Wallace, 2015). It has numerous impacts on atmospheric circulation (Bayr et al., 2014; Bjerknes, 1969; Burgers & van Oldenborgh, 2003; Oort & Yienger, 1996; Wang & Su, 2015) and some studies have linked weather patterns and hydrology over the continental United States to the phase of ENSO (Hidalgo & Dracup, 2003; Singh et al., 2021). There have also been a handful of studies that have shown that the climate’s response to ENSO may be a useful proxy for understanding the climate’s response to global warming (Klein & Hall, 2015; Norris et al., 2022). We expect that models which can simulate the response in the Mountain West to ENSO with high fidelity are likely most skillful at simulating future climate change within that region, so we evaluate the efficacy of our models’ predictions of SWE in the Mountain West based on if they correctly predict the relationship with ENSO. After which, we use the models to evaluate the response to past warming, and what we might expect from future warming. We note that the analysis used in this study could easily be applied to other regions too, so we hope this paper lays the framework for future assessments of models’ ability to forecast regional SWE changes.This paper is organized as follows, in section 2 we describe the data used in this study, including observational datasets, reanalysis, and GCM output. Then section 3 describes a decomposition to determine how three mechanisms contribute to the SWE change and the machine learning framework that is described and validated. The main results of this study are discussed in section 4, followed by concluding remarks in section 5.2. Dataa. Observational DatasetsIn this paper, we use a combination of satellite estimates of precipitation and ground-based measurements of SWE to evaluate if the reanalysis and GCM response to climate variability is consistent with observations.Precipitation estimates from the Integrated Multi-satellitE Retrievals for Global Precipitation Measurement (IMERG; Huffman et al., 2020) product serve as our precipitation observations over the Mountain West. IMERG incorporates, merges, and inter-calibrates various spaceborne infrared (IR) and microwave retrievals, which are then calibrated by ground-based rain gauges. Most of the spaceborne observations were either used in the Tropical Rainfall Measuring Mission Multi-Scale Precipitation Analysis (TMPA; Huffman et al., 2007), or were launched in the follow-on mission, the Global Precipitation Measurement (GPM) constellation, which was meant to improve the accuracy and resolution of TMPA and expand its measurement region from 35° N-S to 65° N-S (Hou et al., 2014). We use the sixth version (V06) of the IMERG dataset, which includes substantial updates to the morphing system that relies on a Lagrangian approach to propagate precipitation vectors into regions without instantaneous precipitation observations (Tan et al., 2019). The newest V06 dataset allows for precipitation to be estimated globally (year 2000 – present) at 0.1°x0.1° resolution, and when compared to surface observations, V06 performs considerably better than previous versions (Pradhan et al., 2022). Though IMERG continues to advance, it still exhibits biases including an underestimation of precipitation over mountainous regions (Anjum et al., 2018; Asong et al., 2017; Dezfuli et al., 2017; Huang et al., 2018; Kim et al., 2017; Sharifi et al., 2016). Due to this known bias, we do not compare the precipitation averages in GCMs and reanalysis directly against the precipitation average measured by IMERG. Instead, we evaluate if the GCMs and reanalysis correctly simulate the relationship between ENSO and precipitation variability in the Mountain West relative to IMERG. In their review of IMERG performance Pradhan et al. (2022) showed that over mountainous terrain there is considerable bias in the mean-precipitation, but with sufficient aggregation (monthly or greater) it can capture ground-based observed variability.In this study we are interested in seasonal snowpack, and so to evaluate the reanalysis and GCMs predictions of SWE we use the SNOw TELemetry (SNOTEL) network of snowpack related climate sensors operated by the Natural Resources Conservation Service of the United States Department of Agriculture that are situated throughout the mountainous terrain of the Western United States. Some of their earliest uses in research was by McGinnis (1997), who estimated the impact of climate change on Colorado snowpack with a machine-learning model trained on the SNOTEL network. The SNOTEL network has since been used throughout the scientific literature as the measure of snowpack in North America (e.g., Fassnacht et al., 2003; Molotch & Bales, 2006; Schneider & Molotch, 2016; Serreze et al., 1999). Using ground-based observations to estimate changes throughout an entire region has some inherent problems; namely that there is typically not sufficient sampling over the entire region to accurately measure the average SWE level in the entire Mountain West. For this reason, like our usage of IMERG, we focus on the SNOTEL variability as opposed to its average value during any given period. SNOTEL observations are also subject to biases related to uneven sampling. Specifically, the SNOTEL locations are not evenly spaced, and do not sample all aspects, slope angles, or regions evenly, so an average over all SNOTEL sites would over-weight the regions that are well sampled and under-weight those with poor sampling. Hence, we do not average the SNOTEL sites over the Mountain West, and instead look at them on a station-by-station basis. Additionally, since the start of the SNOTEL project numerous sites have been added and some decommissioned, so we limit our scope to the SNOTEL sites within the Mountain West that have been consistently providing data from 1980 to 2015.b. ReanalysisOften, observations cannot be directly compared with GCM output due to a combination of the aforementioned observational biases and their relatively short observational periods in comparison to the timescale of climate change (decades/centuries). Hence, we use the European Centre for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5 (ERA5; Hersbach et al., 2020), as a comparison point to determine the efficacy of GCM predictions of the previous century. The ERA5 dataset is a global climate model that assimilates roughly 24 million observations per day from a variety of satellite instruments, and in situ observations of 10 m winds over sea, 2m humidity over land, surface pressure, measurement of upper-air temperature, humidity and wind from dropsondes, radiosondes, and aircraft measurement, and rain rate from ground-based radar-gauge composite observations (Hersbach et al., 2020). All variables are assimilated together to yield a simulated atmosphere that is as consistent as possible with that of the real world. Later in this manuscript we show that ERA5 precipitation and SWE in the Mountain West are rather consistent with observations, enabling us to perform direct comparisons between GCM output and reanalysis to determine the deficiencies in GCMs and their potential causes.c. Global Climate ModelsTo evaluate our current ability to simulate SWE in the Mountain West we also use GCMs that are not constrained by observations (as in ERA5). We primarily use output from two types of simulations performed as a part of the sixth phase of the Coupled Model Intercomparison Project (CMIP6; Eyring et al., 2016). We use the historical simulations, which are fully coupled model simulations forced by the historical levels of atmospheric constituents from the year 1850 to 2014 and were performed as a part of the CMIP6 DECK experiments (Eyring et al., 2016). We also use output from the Scenario Model Intercomparison Project (ScenarioMIP), which provides multi-model projections of future climate based on possible future emissions and land-use changes produced with integrated assessment models (O’Neill et al., 2016). The full list of CMIP6 models that we use is provided in the Supplementary Materials .In addition to the CMIP6 models, we also include the latest model developed by the NASA (National Aeronautics and Space Administration) Goddard Institute for Space Science (NASA-GISS), Model E version 3 (hereafter GISS-E3), which has not yet been published on the CMIP database. This model provides unique insight into the drivers of model spread because NASA-GISS has produced three equally viable versions of their climate model, each characterized by different parameter tunings, but with emergent mean climatological states that broadly match satellite observations. Later in this paper we highlight the results of these three versions of the GISS-E3 model to indicate how much model spread could be a result of parametric uncertainty as opposed to structural errors.3. Methodsa. SWE DecompositionThroughout this study we decompose changes in SWE (either due to variability or climate change) into the portion of the anomaly due to changes in precipitation rate, the fraction of precipitation which occurs as snow, and the fraction of the previous snowfall that remains in the snowpack (and has not melted, evaporated, or blown away). We apply this decomposition to changes in SWE driven by internal variability (ENSO) or climate change. The decomposition method used here was initially introduced by Räisänen (2008), where they show that SWE of seasonal snowpacks can be expressed as a product of the fraction of accumulated snowfall that remains on the ground (\(G\)), and the amount of accumulated snowfall, which is the product of the precipitation (\(P\)) and fraction of precipitation that occurs as snow (\(F\)) integrated from sometime before the start of the snowy season (in our case we use August) to a selected time during the snowy season, as expressed in Equation 1.\begin{equation} \begin{matrix}SWE=G\int_{\text{Aug}}^{t}{\text{FPd}t^{{}^{\prime}}}\ \#\left(1\right)\\ \end{matrix}\nonumber \\ \end{equation}Then in Equation 2 the SWE anomaly is decomposed into the portion of the anomaly due to precipitation anomaly (first term on RHS of Equation 2), the fraction of the precipitation to occur as snow anomaly (second term on RHS of Equation 2), the fraction of the accumulated snowfall to remain on the ground at time t (third term on RHS of Equation 2), and the non-linear combinations of \(G\), \(F\), and \(P\) which is typically at least one order of magnitude less than the other terms (Räisänen, 2008, 2023).\begin{equation} \begin{matrix}\Delta SWE=\overline{G}\int{\overline{F}\Delta}Pdt^{\prime}+\overline{G}\int\Delta F\overline{P}dt^{\prime}+\Delta G\int{\overline{F}\ \overline{P}dt^{{}^{\prime}}}+\frac{1}{4}\Delta G\int{\Delta F\Delta\text{Pd}t^{{}^{\prime}}\ }\ \#(2)\\ \end{matrix}\nonumber \\ \end{equation}In Equation 2 the overbar denotes an average over the entire time domain and \(\Delta\) denotes the anomaly in some subset of the time domain (either a specific set of years or composites of El Niño and La Niña years). Throughout this paper, Equation 2 is applied to each month of the year individually, meaning that \(\Delta\) denotes the anomaly of a specific month in a specific year (or average over some subset of years) from the monthly average over all years within the domain. This allows us to show how global warming and climate variability impact each term in Equation 2 at months of the year. Unfortunately, we cannot perform this calculation with observational datasets because we do not have any observational retrievals of \(F\).b. Gaussian Process EmulationWe use machine learning to determine if deficiencies in GCMs’ SWE and their inter-model spread are due to their simulation of surface temperature and moisture convergence, or if they are due to the numerous small-scale processes that can impact SWE, such as aerosol snowpack effects, evaporation and sublimation rates, and runoff. We use Gaussian Process Emulation (hereafter GPE), a machine-learning algorithm that has robust uncertainty estimates and performs particularly well in cases with relatively small amounts of training data (Watson-Parris et al., 2021a). To train and use our GPE models we use the Earth System Emulator (ESEm) software package, which is designed for use with geospatial data, and allows for complex calibration and emulation to be done in a simple analysis framework (Watson-Parris et al., 2021b). Using ESEm, we train a GPE model to predict March (the typical start of the melting season) SWE and its components from Equation 2 using only surface temperature and moisture convergence from October to March as the predictor variables. We train the model on ERA5 output from the year 1940 to 2015, where a subset of 20 randomly selected years is removed from the training data to test the efficacy of the GPE model.We argue that \(P-E\) is an acceptable approximation for moisture convergence. Moisture convergence is not explicitly archived for GCMs participating in CMIP6, so we approximate moisture convergence as\(P-E\), which does neglect changes in column water vapor; however, changes in column water vapor are typically small compared to \(P-E\). We validate our use of \(P-E\) with ERA5 (which does directly save output of the true moisture convergence). In the Supplementary Materials we show that the seasonal cycle and ENSO variability of\(P-E\) and moisture convergence in ERA5 are quite similar, and in fact, we find that from 1940 to 2015 the mean values of ERA5 moisture convergence and \(P-E\) have a correlation coefficient of 0.95 and during Autumn, Winter and Spring are statistically indistinguishable at\(2\sigma\) confidence.GPE requires selection of a kernel, which can be thought of as a “best guess” of the functional form of the relationship between the input predictors and the output predictand. We select our kernel by looping through all possible additive combinations of the kernels available in ESEm and choosing that for each component of (and total) SWE where the relationship between the predicted SWE and the ERA5 SWE for the testing dataset explains more than 50% of the SWE variance (r2 greater than 0.5), where the slope of the linear relationship is closest to 1, and where a sufficient number of the GPE predictions for the testing data cross the 1:1 line.
