Future snowfall and snowpack changes over the mountains of Southern California are projected using a new hybrid dynamical–statistical framework. Output from all general circulation models (GCMs) in phase 5 of the Coupled Model Intercomparison Project archive is downscaled to 2-km resolution over the region. Variables pertaining to snow are analyzed for the middle (2041–60) and end (2081–2100) of the twenty-first century under two representative concentration pathway (RCP) scenarios: RCP8.5 (business as usual) and RCP2.6 (mitigation). These four sets of projections are compared with a baseline reconstruction of climate from 1981 to 2000. For both future time slices and scenarios, ensemble-mean total winter snowfall loss is widespread. By the mid-twenty-first century under RCP8.5, ensemble-mean winter snowfall is about 70% of baseline, whereas the corresponding value for RCP2.6 is somewhat higher (about 80% of baseline). By the end of the century, however, the two scenarios diverge significantly. Under RCP8.5, snowfall sees a dramatic further decline; 2081–2100 totals are only about half of baseline totals. Under RCP2.6, only a negligible further reduction from midcentury snowfall totals is seen. Because of the spread in the GCM climate projections, these figures are all associated with large intermodel uncertainty. Snowpack on the ground, as represented by 1 April snow water equivalent is also assessed. Because of enhanced snowmelt, the loss seen in snowpack is generally 50% greater than that seen in winter snowfall. By midcentury under RCP8.5, warming-accelerated spring snowmelt leads to snow-free dates that are about 1–3 weeks earlier than in the baseline period.
Streamflow from snowfall and snowpack in mountain regions is a key freshwater source for California. The natural reservoir that snow provides is a critical element of water supply management and flood control. Climate change poses challenges to California’s water infrastructure through its large effects on snow. Higher surface temperatures due to increased greenhouse gas emissions will likely cause more precipitation to fall as rain instead of snow and lead to faster and earlier snowmelt. The resulting snowpack loss is likely to mean less streamflow in late spring and early summer and possibly more streamflow in winter. In addition to the implications for water supply management and flood control, this could have negative impacts on agriculture and recreational activities, as well as plant and wildlife ecology. Meanwhile, water demand is expected to increase because of a hotter climate and growing population. It is therefore critical to assess regional snowpack changes at space and time scales relevant for regional and local resource management decisions.
Many studies have documented observed changes in snow measures over the past several decades and assessed impacts of global and regional warming on snow (e.g., Barnett et al. 2005; Knowles et al. 2006; Barnett et al. 2008; Pederson et al. 2011; Pierce and Cayan 2013). Several observational and numerical modeling studies have also investigated potential effects of warming on Sierra Nevada snowpack (e.g., Howat and Tulaczyk 2005; Mote et al. 2005; Mote 2006; Maurer 2007; Cayan et al. 2008; Kapnick and Hall 2010; Pavelsky et al. 2011) and other mountainous regions in the western United States (Kim et al. 2002; Knowles and Cayan 2002; Leung and Qian 2003; Mote 2003; Leung et al. 2004; Snyder et al. 2004; Bales et al. 2006; Salathé et al. 2008; Minder 2010; Abatzoglou 2011; Kapnick and Hall 2012; Ashfaq et al. 2013; Klos et al. 2014). However, few studies have focused on snowfall and snowpack variability and change in Southern California mountains.
Climate change information is available from projections by general circulation models (GCMs), such as those archived in phase 5 of the Coupled Model Intercomparison Project (CMIP5; Taylor et al. 2012). However, the spatial scale of GCM output is far too coarse to provide accurate estimates of snow variability and projections for regional or local studies. GCMs poorly represent topography, even in the largest mountain ranges, minimizing topographic influences on circulation and hence rainfall and snowfall in mountain regions (Luce et al. 2013). This presents an issue in the study of climate change effects on snowpack, because snowpack changes may vary by region and elevation. Whereas snowpack may decline in some areas, it could be less sensitive to warming in others, with no change or even increases where total precipitation is projected to increase. In addition to changes in mean snow quantities, changes in seasonality and timing of snowfall and snowmelt are expected. Understanding which areas of mountain ranges are most vulnerable to climate change is critical for regional and local climate assessments and water management planning. To address the limitations of coarse-resolution GCMs, dynamical and statistical techniques are commonly used to downscale GCM projections.
Dynamical downscaling employs a regional climate model with much higher spatial resolution than GCMs, while relying on GCMs for boundary information. Like GCMs, a regional climate model is based on fundamental physical laws governing the climate system. A high-resolution regional climate model can simulate regional- and local-scale climate variations, including orographic precipitation, rain shadows, and snow processes in mountainous regions. Dynamical downscaling has been widely applied over many regions (e.g., Leung and Ghan 1999; Giorgi et al. 2001; Wang et al. 2004; Chin 2008; Gao et al. 2012; Mallard et al. 2015) and is valuable in providing information on regional climate change, including in California (Cayan 1996; Cayan et al. 2001; Leung et al. 2003; Caldwell et al. 2009; Qian et al. 2009; Pan et al. 2011). Although dynamical downscaling probably provides the most physically realistic approach to downscaling low-resolution climate data and provides a comprehensive suite of climate variables, it is enormously computationally expensive.
Alternatively, statistical downscaling is computationally cheap. Statistical models are based on empirical relationships between known climate predictors and climate variables of interest at the regional scale. These relationships are presumed to hold true for the future and are used to project regional climate change given the change in the climate predictors (von Storch et al. 1993; Wilby et al. 2004). Statistical downscaling approaches have been applied to various regions of interest. Temperature and precipitation are the climate variables most commonly used in statistical downscaling studies (e.g., Hayhoe et al. 2004; Pierce et al. 2013). Snowfall and snowpack are statistically downscaled less often, in spite of their importance for hydroclimate.
To obtain reliable climate change information at the regional scale, this study develops and applies a novel hybrid dynamical–statistical framework. Unlike previous downscaling studies, which use either a dynamical or a statistical technique, this study combines both, developing statistical relationships directly from dynamically downscaled output. This study undertakes dynamical downscaling for a reanalysis-driven baseline simulation (1981–2000) and a small representative sample of GCMs forced by a “business as usual” greenhouse gas emissions scenario [i.e., representative concentration pathway 8.5 (RCP8.5; Meinshausen et al. 2011)] for a mid-twenty-first-century time slice (2041–60). A statistical model is then developed to reproduce the snowfall and snowpack variations in the baseline period using other climate variables (i.e., surface temperature and precipitation) as predictors. We confirm the accuracy of the statistical model by comparing its future projections to those of the dynamically downscaled mid-twenty-first-century simulations. Finally, we apply the validated statistical model, using temperature and precipitation projections from previous studies (Walton et al. 2015; Sun et al. 2015; Berg et al. 2015) to project regional snowfall and snowpack changes for all the available GCMs.
Combining dynamical and statistical downscaling techniques in this way allows us to incorporate the most important dynamical processes shaping regional snowfall and snowpack variations and change without dynamically downscaling each GCM. The hybrid technique provides ensemble-mean estimates and quantifies uncertainties associated with differences across various GCM projections. This study marks the first time a full ensemble of GCM output has been downscaled to high-resolution regional scales to create both snowfall and snowpack projections. The hybrid framework also allows us to assess snowfall and snowpack changes associated with different emissions scenarios and time periods: We downscale the available GCMs for the mid-twenty-first-century and end-of-twenty-first-century (2081–2100) periods under both the “business as usual” RCP8.5 scenario and the “mitigation” RCP2.6 scenario, in which greenhouse gas emissions are aggressively reduced in the coming decades.
The paper is organized as follows: section 2 describes the dynamical downscaling experiment design and evaluation. Section 3 describes the statistical downscaling framework and its performance in baseline prediction. Section 4 presents the evaluation of performance of statistical downscaling under future climate change. Section 5 presents the future snowfall and 1 April snow water equivalent (SWE) changes, followed by comparison and discussion of sensitivities to choice of greenhouse gas emissions scenario. Major findings are summarized and discussed in section 6.
2. Dynamical downscaling
a. Experimental design
Dynamical downscaling simulations are performed using the Weather Research and Forecasting (WRF) Model (Skamarock et al. 2008), version 3.2. WRF is a community mesoscale model developed by the National Center for Atmospheric Research (NCAR). It is designed for use on regional grids for a range of applications, including weather forecasts and climate simulations. WRF consists of a fully compressible nonhydrostatic dynamical core with high-order, conserving numerical techniques. WRF provides a number of options for its various physics parameterization schemes. The following parameterization schemes are chosen: Kain–Fritsch (new Eta model) cumulus (Kain 2004), Yonsei University boundary layer (Hong et al. 2006), Purdue Lin microphysics (Lin et al. 1983), Rapid Radiative Transfer Model longwave radiation (Mlawer et al. 1997), and Dudhia shortwave radiation (Dudhia 1989) schemes. We use the Noah land surface model (Chen and Dudhia 2001) to simulate land surface processes, including vegetation, soil, snowpack, and exchange of energy, momentum, and moisture between the land and atmosphere.
We use three nested domains (18, 6, and 2 km) to reach a high enough resolution to represent the region’s topography and coastlines. As Fig. 1a shows, the outermost domain, at 18-km resolution, covers the entire state of California and the adjacent ocean. The middle domain (6 km) covers roughly the southern half of California, including the southern Sierra Nevada. The innermost domain (2 km) encompasses the greater Los Angeles region and its surrounding mountains. In each domain, all variables within five grid cells from the horizontal lateral boundary are relaxed toward the corresponding values at the boundaries. This procedure ensures smooth transitions across these boundaries. Each domain has 43 sigma levels in the vertical with model top level at 50 hPa. To provide a better representation of surface and boundary layer processes, the model’s vertical resolution is enhanced near the surface, with 30 sigma levels below 3 km. Figure 1b shows the topography for the innermost domain at its native 2-km resolution. The main features of both the topography and coastlines are represented well at this resolution.
Using this model configuration, we first perform a reanalysis-driven baseline simulation (1981–2000), the purpose of which is threefold: 1) to evaluate the model’s ability to simulate regional climate, 2) to provide a baseline climate state against which a future climate simulation can be compared, and 3) to build the statistical relationships in the statistical downscaling framework (section 3a). For the baseline simulation, WRF is forced along the boundaries of the outermost domain by the National Centers for Environmental Prediction’s 3-hourly North American Regional Reanalysis (NARR; Mesinger et al. 2006) for the 1981–2000 period. Using the same model configuration, we then perform a series of dynamical downscaling experiments to simulate future climate associated with five CMIP5 GCMs for the midcentury period under the RCP8.5 forcing scenario. The selected GCMs are the NCAR Community Climate System Model, version 4 (CCSM4; Gent et al. 2011); the Centre National de Recherches Météorologiques Coupled Global Climate Model, version 5 (CNRM-CM5; Voldoire et al. 2012); the NOAA/Geophysical Fluid Dynamics Laboratory Climate Model, version 3 (GFDL CM3; Donner et al. 2011); the Atmosphere and Ocean Research Institute (AORI) at the University of Tokyo, National Institute for Environmental Studies (NIES), and JAMSTEC Model for Interdisciplinary Research on Climate, Earth System Model, Chemistry Coupled (MIROC-ESM-CHEM; Watanabe et al. 2011); and the Max Planck Institute for Meteorology Earth System Model, low resolution (MPI-ESM-LR; Brovkin et al. 2013).
To produce future climate boundary conditions for the WRF simulations, we first quantify the differences in GCM monthly climatology between the midcentury (2041–60) and baseline (1981–2000) periods. Monthly differences are calculated for three-dimensional variables, including temperature, relative humidity, zonal and meridional winds, and geopotential height, and two-dimensional surface variables, including surface temperature, relative humidity, winds, and pressure. On a monthly varying basis, we add these climate change signals to the NARR reanalysis data corresponding to the baseline period. Thus, we perturb the NARR baseline boundary conditions (1981–2000) with monthly averaged climate change signals between the midcentury and baseline (2041–60 minus 1981–2000) provided by each GCM. These perturbed NARR data are then used as boundary conditions for the outermost domain of the regional model. This technique has been previously used to downscale GCM signals to fine spatial scales for other regions (e.g., Schär et al. 1996; Sato et al. 2007; Hara et al. 2008; Knutson et al. 2008; Kawase et al. 2009; Lauer et al. 2010; Seo and Xie 2011; Rasmussen et al. 2011). Sun et al. (2015) addressed this perturbation approach’s caveats and limitations, including the unchanged interannual variability and weather and transient signals on the model domain’s boundaries. CO2 levels are also increased in WRF to match the changes in CO2-equivalent radiative forcing in the RCP8.5 scenario averaged over the midcentury period compared to the baseline.
Computational limitations do not allow us to perform full 20-yr dynamical downscaling simulations for each of the five GCMs for the mid-twenty-first-century period. So we first perform a 20-yr (2041–60) dynamically downscaling of climate change signals in CCSM4. Because we perturb each year in the future 20-yr period with the same monthly varying climate change signals from CCSM4, we expect the climate change patterns for each year to be similar. In fact, any 4-yr period within the full 20-yr period captures the snow change and other climate change signals found in the full 20-yr period of CCSM4 downscaling very well (not shown). Therefore, we dynamically downscale the other four GCMs for only four years. In each simulation, the perturbed boundary conditions are created by adding the 20-yr GCM climate change signal (2041–60 values minus 1981–2000 values) to the 1997–2000 NARR data.
b. Model evaluation
Figure 2 presents the dynamically downscaled spatial distributions of the baseline (1981–2000) snowfall water equivalent (SFE; units in millimeters) climatology over the Los Angeles region from November through April. Winter months December, January, February, and March (DJFM) have the greatest snowfall for most areas and have the largest spatial extent of snowfall. DJFM snowfall accounts for more than 80% of annual accumulated snowfall across the region. Snowfall is mainly found in mountain regions at elevations of 1500 m and higher. Snowfall generally increases with elevation, with larger amounts on the coastward-facing side of the ranges. Climatological snowfall is negligible (less than 1 mm month−1) at high-elevation desert regions (e.g., the Mojave Desert). Such a tiny amount of snowfall would probably not survive long enough on the ground to lead to any substantial accumulation, especially when the surface temperature is not cold enough. Snowfall greater than 200 mm month−1 is seen in high-elevation mountain regions (2000 m and higher), including the southern rim of the Sierra Nevada and the Tehachapi, San Emigdio, San Gabriel, San Bernardino, and San Jacinto Mountains (refer to Fig. 1b for mountain locations). At the highest elevations, monthly accumulated snowfall can reach 300 mm month−1 in the peak season.
The dynamical model’s ability to reproduce snowfall climatology and its temporal and spatial variations are assessed by comparing output from the 2-km baseline climate simulation to available observational measurements. Quality-controlled daily snowfall, precipitation, and maximum temperature data are obtained from the Western Regional Climate Center (WRCC), which collects daily data from the National Weather Service (NWS) Cooperative Observer Program (COOP).
Data from each NWS COOP station within the 2-km WRF domain are evaluated, and stations with no recorded snowfall are excluded. From the remaining stations, we select those with daily snowfall and temperature measurements for at least 75% of the baseline period, which allows for the assessment of both climatology and interannual variability. Four NWS COOP stations meet these criteria: Big Bear Lake and Lake Arrowhead in the San Bernardino Mountains, Idyllwild in the San Jacinto Mountains, and Tehachapi in the Tehachapi Mountains. Table 1 summarizes the identifying information associated with each observational station, including COOP station number, location, elevation, and observation period of available data. All four stations are in mountain regions and represent a variety of elevations ranging from 1224 to 2070 m. At the four stations, individual months with more than five missing days of data are not used for monthly statistics and are not included in annual totals for that year.
There are several barriers to comparing observed and simulated snowfall data in a straightforward way (Kunkel et al. 2007, 2009). First, the model grid cells in the vicinity of a measurement station may not be at the exact elevation of the measurement station. Because of the strong dependence of snowfall on elevation, this can lead to a slight mismatch between observed and simulated data. To minimize this issue, we consider the four model grid cells nearest each station and select the one with an elevation that is in closest agreement with that of the station.
Second, a direct comparison of SFE simulated by WRF to observed snowfall is not possible because of an absence of in situ liquid water equivalent of snowfall measurements in the study area for the baseline period. To compare WRF-simulated SFE to fresh snowfall observations at NWS COOP stations, we convert the observed depth of fresh snowfall to its water equivalent using the following relationship:
The density of freshly fallen snow is not directly measured. However, total daily precipitation is available. We use the above equation to estimate density of freshly fallen snowfall at each NWS COOP station. SFE is taken to be the observed daily precipitation value on days with nonzero snowfall and maximum temperatures less than or equal to 0°C. We exclude days with nonzero snowfall and maximum temperatures above 0°C, as some of the measured precipitation on those days is likely to have been in the form of rain.
The threshold 0°C is chosen to be consistent with the Noah land surface model, which uses a step function with a temperature threshold of 0°C for the partitioning of precipitation into snowfall and rainfall. The model evaluation also uses the 0°C temperature threshold step function when estimating the density of freshly fallen snowfall at each NWS COOP station in our study domain. However, we do acknowledge the sensitivity of SFE to snow versus rain partitioning schemes (Nolin and Daly 2006). To address the uncertainty in our model evaluation associated with using a step function with a temperature threshold of 0°C, we performed two additional tests using temperature thresholds of 2.5° and 5°C. In these cases, the estimated observed snow density only changes about 5%–10% for a range of possible rain-versus-snow temperature threshold values, which would not qualitatively affect our results. Additionally, because the estimated snow density is not significantly changed when considering a threshold of 0° or 5°C for each of the four NWS COOP stations in our domain, it is likely that a ramp function would give similar results, considering that commonly used ramp thresholds (e.g., Jordan 1991) are functions of temperatures colder than 5°C.
Across the four stations estimated snow densities range from 63 to 129 kg m−3, well within the range of previous studies. These studies found that snow density varies from 10 to 500 kg m−3, depending on location, meteorological condition, crystal size, crystal shape, degree of riming, and other snow metamorphosis processes (Pomeroy and Brun 2001; Roebber et al. 2003; Baxter et al. 2005; Kay 2006). Our estimated snow densities are also comparable to the most commonly observed values, between 60 and 100 kg m−3, suggested by Judson and Doesken (2000). We acknowledge there are more issues with the observational snowfall as well as the snowfall density (Kunkel et al. 2007, 2009). Nevertheless, using our calculation of observed snowfall density, we are able to estimate the observed SFE using Eq. (1), and then compare the estimate to simulated snowfall at the selected model grid cells.
We now present model evaluation results, focusing on the climatological seasonal cycle and interannual variability. Figure 3a compares monthly climatological simulated snowfall values (i.e., SFE) to observations. WRF’s seasonal cycle of snowfall is consistent with observations for each of the four stations. The model also accurately simulates spatial variations in climatological snowfall, as the simulated and observed values are very well correlated, with an average correlation coefficient of 0.95 across the four stations. There is very little systematic bias across the four stations, with a mean bias of 0.46 mm and a normalized mean bias of 5.3%. The root-mean-square error of all data points in Fig. 3a is 4.4 mm. A large fraction of this error has to be because of the unavoidable assumption of constant snow density. It is noteworthy that only four NWS COOP stations are available for the validation, and for each station only a monthly climatology is presented so that the sample size is rather small. Figure 3b compares annual accumulated snowfall (from September to August of the following year) between each NWS COOP station and the corresponding WRF grid cell for all years in the baseline period with available data. For each station, the simulated and observed values are significantly correlated, with an average correlation of 0.59. The overall correlation of the data points in Fig. 3b is 0.76, providing a combined evaluation of the spatial and interannual snowfall variability. The fact that this number is higher than the correlation associated with any individual station is an indication that the model captures spatial variability somewhat better than temporal variability. There is no systematic bias in WRF-simulated snowfall, with a mean bias of 6.2 mm, a normalized mean bias of 5.8%, and a root-mean-square error across all points of 66 mm. Again, the assumption of constant snow density probably contributes significantly to the model–observation discrepancies in Fig. 3b.
Overall, Fig. 3 shows that our WRF framework simulates the temporal and spatial variations of snowfall during the baseline period with reasonable accuracy at specific mountain locations where reliable observational data are available. Based on this evidence, it is very likely that the model is able to reproduce the temporal and spatial snowfall variations across the whole domain, even at very high elevations where there is certainly substantial snowfall but observations are sparse or unavailable.
3. Statistical downscaling
In this section, we present the statistical framework to reproduce snowfall and snowpack variations. The framework is based on the multiple linear regression analysis between interannual variations in snow properties and climate variables from the reanalysis-driven, dynamically downscaled baseline simulation.
a. Snowfall model
Precipitation and temperature are key factors affecting snowfall. The relationship between precipitation and snowfall estimates is fairly straightforward. In the Noah land surface model, precipitation phase is determined by a simple partitioning scheme based on estimated surface air temperatures. (Note that both surface air temperature and surface skin temperature are used in the analysis, and they produce similar results.) Temperatures below the freezing point of water (0°C) are assumed to result in precipitation that is 100% snowfall, and those above are assumed to result in 100% rainfall.
Snowfall is highly sensitive to elevation and season. To create the statistical framework, we first bin all elevations in 100-m increments. Then, for all grid cells within each bin, we average variables for each month. Figure 4a shows a vector representation of interannual correlations between snowfall and precipitation (x direction) and snowfall and temperature (y direction) in each elevation bin and winter month (DJFM) in the baseline period. Nearly all elevation bins (except for very low elevations below 1500 m) show significantly positive correlations (rightward direction) between precipitation and snowfall in each of the winter months. This is expected, as precipitation places a fundamental limit on how much snowfall is possible and is completely determinative of snowfall at elevations above the freezing line. Significantly negative correlations (downward direction) between snowfall with temperature are seen in low- to midelevation bins and high-elevation bins only in the two shoulder months. Especially at low elevations, temperature plays the stronger role. For instance, the correlation between snowfall and temperature is as large as 0.8 around 1700 m in February. At these lower elevations, regional warming in the absence of precipitation changes would result in snowfall declines. In contrast, at high elevations, because temperatures remain below the freezing point, temperature fluctuations have less impact on precipitation phase. The highest elevations might be also susceptible to warming, but most are so cold that the warmth of a particular winter has only a small effect on snowfall and snowpack, particularly in January and February.
We next explore how well WRF-simulated interannual variability in monthly accumulated snowfall is represented using monthly averaged temperature and accumulated precipitation. To do this, we use the following multiple linear regression equation:
where α and β are the regression coefficients for monthly averaged temperature T and monthly accumulated precipitation P, respectively, and γ is the residual term. The monthly accumulated snowfall values associated with temperature and precipitation are represented by , which is not allowed to be negative. As shown in Fig. 4a, the sensitivity of snowfall to temperature and precipitation varies by elevation and month. We construct the best-fit multiple linear regression model for each winter month and elevation bin in the region and apply this multiple linear regression model to the WRF-simulated baseline snowfall, temperature, and precipitation.
Using this regression model, a snowfall value for each winter month and each elevation bin can be predicted from WRF’s temperature and precipitation values; this value can then be compared with WRF-simulated snowfall. Figure 4b shows that the regression-derived snowfall is significantly correlated with WRF-simulated snowfall in almost all elevation bins in each winter month. Less significant correlation at low elevations (below 1500 m) in January may be the result of the fact that neither temperature nor precipitation is a very good predictor of snowfall at this particular time (Fig. 4a). Figure 5 further compares the regression-derived snowfall and WRF-simulated snowfall for four sample elevation bins in each of the winter months. It suggests temperature and precipitation can predict the interannual variability of mountain snowfall with a high degree of success. At elevations above 3000 m (Fig. 5, top row), the WRF-simulated snowfall in each winter month can be nearly perfectly reproduced by the regression model, especially in January–March (JFM). For a high-elevation bin (2500–2600 m; Fig. 5, second row from top), the regression-derived snowfall and WRF-simulated snowfall is highly correlated (r = 0.95 averaged over the winter months). At these elevations, snowfall is determined almost entirely by total precipitation (see Fig. 4a). For a midelevation bin (2000–2100 m; Fig. 5, third row from top), the correlation between the regression-derived snowfall and WRF-simulated snowfall is 0.90. For a low-elevation bin (1500–1600 m; Fig. 5, bottom row), the correlation drops to 0.72. This indicates the statistical model’s quality declines somewhat as elevation decreases. Overall, Figs. 4b and 5 suggest snowfall derived from the multiple linear regression model represents the interannual variation at all elevations quite well.
We then evaluate how well the regression model reproduces the seasonal cycle of snowfall in the baseline period. Figure 6 first presents the regression-derived (dashed black) baseline seasonal cycles of snowfall for elevations above 1500 m, compared to the dynamically downscaled results (solid black). It shows that the regression-derived snowfall overestimates the dynamically downscaled snowfall in and before February, while underestimating thereafter. However, the differences between the two are small, with a normalized mean bias for the annual accumulated snowfall of approximately −3.5%. This underestimation is partially attributable to the fact that the regression model uses monthly mean temperature and precipitation as predictors. The monthly mean climates likely smooth out the day-to-day or event-to-event variability of snowfall. We further discuss this underestimation as an additional uncertainty when interpreting the future changes derived from the statistically downscaled projections in section 6.
b. SWE model
As with snowfall, we build an analogous statistical framework for the snowpack on the ground (represented by SWE) using multiple linear regression analysis. SWE is the depth of water that would result if the entire snowpack were to melt instantaneously. The 1 April SWE is commonly used to assess snowpack and its variability. It is the most frequent observation date and is extensively used for spring streamflow forecasting and analysis (Howat and Tulaczyk 2005; Mote et al. 2005), as it is an indicator of the interannual variation in snowpack.
The relationship between 1 April SWE and winter mean temperature is more complex than that between 1 April SWE and winter accumulated precipitation. Additionally, temperature affects SWE in a more complex way than it does snowfall, as it directly impacts snowmelt through convective heat transfer from the air to the snowpack. Temperature is also indirectly consequential for humidity and water vapor pressure, which both contribute to snow sublimation and snowmelt (Hamlet et al. 2005). Previous studies have assessed the contributions of temperature and precipitation to observed snowpack variations and trends (Serreze et al. 1999; Howat and Tulaczyk 2005; Mote et al. 2005; Mote 2006). Informed by these studies, we predict 1 April SWE from mean temperature and accumulated precipitation during the preceding winter months (December through March) using the following multiple linear regression equation:
As with snowfall, SWE sensitivity to temperature and precipitation varies with elevation. So we construct the regression model for each elevation bin. A plot similar to Fig. 5 (not shown) reveals that the multiple linear regression model reproduces the interannual variations of WRF-simulated baseline 1 April SWE at all elevations as successfully as the model predicting snowfall.
The success of this statistical framework indicates that the interannual variability in regional snowfall and snowpack are well explained by summaries of monthly regional climate. This implies that day-to-day or event-to-event weather fluctuations might be of secondary importance.
4. Statistical model performance under climate change
To project future snowfall and SWE using the statistical framework, it is necessary to verify that the relationships between baseline snow properties and baseline temperature and precipitation hold in the future climate. In this section, we validate the statistical downscaling framework by comparing its predictions of future snowfall and SWE against output from the multiple dynamical downscaling future simulations described in section 2a.
Figure 6 shows the statistically downscaled mid-twenty-first-century seasonal cycles of snowfall under the RCP8.5 scenario (green), compared to the corresponding dynamically downscaled results for each future simulation (red) and the baseline WRF simulation (solid black) and statistical downscaling results (dashed black). Data are shown for snowfall averaged over elevations above 1500 m within the domain. There are significant snowfall declines at midcentury in all months for four of the five dynamical downscaling experiments. The only exception is the CNRM-CM5 simulation, which shows no change or even a slight increase in snowfall. This is the result of a projected precipitation increase for CNRM-CM5 (Berg et al. 2015), which apparently cancels out any snowfall reduction due to warming. The statistically downscaled results reproduce dynamically downscaled snowfall changes in each wet month with an error of less than 10%. The dynamical snowfall declines for CCSM4, GFDL CM3, MIROC-ESM-CHEM, and MPI-ESM-LR are all generally well captured by the statistical framework. The statistical model also reproduces snowfall’s insensitivity to climate change in the CNRM-CM5 projection.
Figure 7 compares statistically and dynamically downscaled results for winter (DJFM) accumulated snowfall as a function of elevation for all five experiments. The dynamically downscaled snowfall (solid black) and statistically downscaled snowfall (dashed black) for baseline are shown as a background reference. Dynamically downscaled results show snowfall reductions at all elevations for all experiments, with the exception of the very high elevations in CNRM-CM5. For each experiment, the statistical model tracks the dynamically downscaled results well at most elevations. For CCSM4, the statistical result reproduces the WRF-simulated snowfall almost perfectly in every elevation bin. For GFDL CM3, CNRM-CM5, and MPI-ESM-LR, the statistical results overestimate overall snowfall, but the bias is generally less than 10%. It is noteworthy that, above 2800 m, CNRM-CM5 simulates a snowfall increase, suggesting the stronger control is exerted by increased precipitation. This increase is captured in the statistically downscaled results as well. In general, the statistical model captures dynamically projected snowfall changes at all elevations in all five experiments with a reasonable degree of accuracy.
Next, we evaluate the statistical model’s ability to estimate 1 April SWE. Figure 8 compares statistically and dynamically downscaled results for 1 April SWE as a function of elevation for all five experiments. For each GCM, the statistical SWE estimate generally matches dynamically downscaled SWE, with less than 10% error. The dwindling of snowpack on or even before 1 April at low and moderate elevations in GFDL CM3 and MIROC-ESM-CHEM underscores the dominant effect of regional warming at these elevations in these models. The statistical model captures this effect well. The statistical projections also capture the increased SWE at very high elevations for CNRM-CM5. Overall, Figs. 6, 7, and 8 demonstrate that the statistical framework based on baseline relationships can be used to project future snowfall and SWE.
5. Projections for all GCMs
In this section, we estimate future winter snowfall and 1 April SWE using the aforementioned statistical framework, which efficiently approximates snow changes that would have been produced had dynamical downscaling been performed on all available GCMs (Table 2). We then assess the sensitivity of snow outcomes to emissions scenario.
a. Snowfall projection
Two future time slices under the RCP8.5 emissions scenario are chosen for comparison: mid-twenty-first century (2041–60) and end-of-twenty-first century (2081–2100). Statistical model predictors—temperature and precipitation—are taken from Walton et al. (2015), Sun et al. (2015), and Berg et al. (2015), who developed hybrid dynamical–statistical downscaling approaches to project surface air temperature and precipitation changes for the two future time slices. These studies downscaled available CMIP5 GCMs and estimated the ensemble mean as well as the associated intermodel range of future surface warming and precipitation changes in the greater Los Angeles region. Their main findings include warming at midcentury and continued warming at end of century, although the warming amplitude varies significantly across the region and GCMs (Walton et al. 2015; Sun et al. 2015). Winter precipitation projections vary in both sign and amplitude across models. Some GCMs project moistening, and others project drying in the region. But overall precipitation change signals are weak, yielding no significant ensemble-mean precipitation change (Berg et al. 2015). With these studies’ temperature and precipitation projections for each GCM as inputs, we use the snowfall statistical model to downscale and quantify snowfall changes in the region.
Figures 9a and 9b present CMIP5 ensemble statistically downscaled DJFM accumulated snowfall as a function of elevation for the mid-twenty-first century and end-of-twenty-first century, respectively, under RCP8.5. Projected snowfall is shown as a percentage of the baseline snowfall. Under RCP8.5, the ensemble mean shows a snowfall loss everywhere. Low elevations have the greatest reductions in snowfall, with less than 50% of baseline snowfall remaining on average. At lower elevations, surface air temperatures during precipitation events are more likely to breach the freezing point of water as the climate warms. Hence, a snow event is more likely to be converted to a rain event. Midelevation (2000–2500 m) snowfall is somewhat less sensitive to climate change, retaining about 70% of baseline snowfall in the ensemble mean. High-elevation snowfall is projected to be relatively resilient, with roughly 90% of baseline snowfall remaining. Below 2400 m, every GCM projects a snowfall decline compared to the baseline under RCP8.5. At the highest elevations (above 3100 m), about two-thirds of the GCMs predict snowfall loss. High-elevation snowfall is relatively insensitive to warming because of the insensitivity of snowfall to temperature fluctuations (see Fig. 4) and is instead dominated by the precipitation change. The GCMs showing increased snowfall above 3100 m are those with a projected increase in total precipitation. Figure 9a also shows the spread across GCM projections for each elevation bin. The spread is substantial across GCMs at each elevation, roughly 50%–60% points. For example, at the highest elevations, projected snowfall percentages range from about 70%–120%. This range, which can be taken as a measure of uncertainty, is nearly half the ensemble-mean projection (about 95%).
End-of-twenty-first-century snowfall under RCP8.5 (Fig. 9b) shows a further reduction from midcentury values at every elevation. Only about a quarter of GCM projections show a snowfall increase at the highest elevations. Ensemble-mean end-of-century snowfall is less than 20% of baseline snowfall at elevations below 1800 m, about 50% at moderate elevations, and about 80% at high elevations. Model spread becomes somewhat smaller at end of century (40%) than at midcentury (50%). The greater consistency among GCMs at end of century may be because of the increasing magnitude of the temperature signal in all models and its increasingly powerful effect on snowfall reduction.
b. SWE projection
We next apply the statistical framework for 1 April SWE. Figures 9c and 9d present the CMIP5 statistically downscaled 1 April SWE as a function of elevation for the mid-twenty-first century and end-of-twenty-first century, respectively, under RCP8.5. Ensemble-mean 1 April SWE decreases at all elevations for both time periods. In general, higher elevations have more remaining snowpack, in accordance with elevation-dependent snowfall projections. As shown in Fig. 9c, every GCM projects snowpack reduction below 2400 m at midcentury, while above 2400 m, about one-tenth to one-third of the downscaled GCMs project increased snowpack. At very high elevations, ensemble-mean midcentury remaining snowpack is about 90% of baseline snowpack, similar to the ensemble-mean snowfall projections. Regional warming does not cause the temperature to breach the freezing point of water at very high elevations, so its impact on snow ablation and snowmelt is minimal at these high elevations.
In contrast, at low and moderate elevations, the percentage of SWE lost is substantially larger than that of snowfall. Figure 10 shows the ensemble-mean percentages of 1 April SWE compared to those of winter snowfall for the three sampled elevations. The remaining SWE percentage is as low as 26% at low elevations and about 54% at moderate elevations, while the corresponding snowfall percentage is nearly 48% and 71% at low and moderate elevations, respectively. Thus, 1 April SWE is reduced by an additional 20% points from the already-reduced winter snowfall. This suggests that, in addition to its impacts on snowfall loss, warming plays a further role in enhancing ablation and melting of snow at low-to-moderate elevations. In all models, warming is large enough either to exaggerate the snow decline seen in models with reduced precipitation or to overcome any snow accumulation increase in models with increased precipitation. Figure 9c also shows that there is a significant spread in projections of midcentury SWE across the GCMs, particularly at moderate and high elevations. For example, at 2700 m, the SWE percentage ranges from about 35% of baseline to about 120% of baseline. Model spread is larger than that seen in snowfall projections shown in Fig. 9a, indicating that variations among GCMs in snow ablation and melting add to their variations in snow deposition.
At end of century under RCP8.5, the 1 April SWE is further reduced from midcentury values at all elevations, including the very high elevations (Fig. 9d). This implies that further warming more than compensates for any precipitation increases. Moderate elevations see the largest further reduction of snowpack, and the uncertainty range across GCMs is generally smaller than at midcentury. At elevations lower than 1700 m, all GCMs project a complete disappearance of snowpack on or before 1 April by end of century. End-of-century 1 April SWE reduction is larger than that of snowfall, particularly at moderate elevations. Figure 10 shows that, whereas about 52% of baseline winter snowfall remains at the end of century for moderate elevations, only about 31% of the snowpack remains on 1 April, yielding an additional SWE reduction of 20% points. This further demonstrates that warming at end of century significantly enhances snow ablation and melting processes. A rule of thumb is that roughly two-thirds of the 1 April SWE loss is due to snowfall reduction, while about one-third is due to enhanced melting.
c. Sensitivity to choice of emissions scenario
To account for uncertainty associated with choice of future emissions scenario, we project midcentury and end-of-century winter snowfall and 1 April SWE for the CMIP5 ensemble under RCP2.6, which assumes greenhouse gas emissions peak around 2030 then decline substantially thereafter. Figure 9 and Fig. 10 assess the sensitivities of snow changes to emissions scenario. The cross markers in Fig. 9 denote the ensemble-mean snowfall and 1 April SWE projections under RCP2.6 for the corresponding time period. The ensemble-mean projection of snowfall under RCP2.6 is greater than that under RCP8.5 at all elevations for midcentury and for nearly all elevations for end of century. However, at midcentury, RCP2.6 shows only about 10% points more snowfall in low and middle elevations than RCP8.5, and the difference between the two scenarios in higher elevations is minimal. End-of-century projections show a greater contrast between the two scenarios, especially in low and moderate elevations. For instance, at very low elevations, the difference is as large as about 40% points. In contrast to RCP8.5, the RCP2.6 scenario sees a negligible snowfall change at end of century compared with midcentury. Similarly, the remaining snowpack on 1 April in RCP2.6 represents minimal change from midcentury to end of century at all elevations, leading to significant contrast at end of century between these two scenarios.
6. Summary and discussion
In this study, we develop a hybrid dynamical–statistical downscaling technique to produce 2-km-resolution projections of future snowfall and snowpack changes in the mountains of Southern California at the middle and end of the twenty-first century. This new hybrid technique combines both dynamical and statistical downscaling and develops the statistical relationships directly from dynamically downscaled output. The first step is to perform a dynamical downscaling for a baseline simulation (1981–2000) and a representative sample of GCMs forced by the RCP8.5 emissions scenario for a midcentury time slice (2041–60). The choice of GCMs and discussions on the caveats of this dynamical downscaling approach were addressed in Sun et al. (2015) and Berg et al. (2015). A statistical model is then developed to reproduce the snowfall and snowpack variations in the baseline period using surface temperature and total precipitation as predictors. The accuracy of the statistical model is evaluated by comparing its predictions to those of the dynamically downscaled baseline and midcentury simulations. Using surface temperature and precipitation projections from previous studies (Walton et al. 2015; Sun et al. 2015; Berg et al. 2015), we apply the validated statistical model to downscale regional snowfall and snowpack changes corresponding to all available GCMs. We further downscale GCM output for the end-of-century time slice (2081–2100) under both the “business as usual” RCP8.5 scenario and the “mitigation” RCP2.6 scenario to assess snow changes associated with different emissions scenarios and time periods.
We project that, in the future, the Southern California mountains are likely to receive substantially less snowfall and have less snowpack on the ground than in the baseline period. Under RCP8.5, midcentury area-mean snowfall is just 70% of the corresponding baseline value. Under RCP2.6, the amount is somewhat higher (80% of baseline snowfall). After midcentury, however, the two scenarios diverge significantly. By end of century under RCP8.5, snowfall sees a dramatic further reduction from midcentury levels; area-mean snowfall is only about half the baseline value. On the other hand, under RCP2.6 snowfall sees only a minimal further reduction from midcentury values. Because of the spread in the GCM climate projections, these values are all associated with large intermodel uncertainty, in the range of 50%–60% points. For both time slices, the snowfall loss is consistently greatest at low and moderate elevations. At higher elevations, snowfall totals are similar to those in the baseline, and about one-third of GCMs project a snowfall increase. High-elevation snowfall is insensitive to warming because temperatures are well below the freezing point of water. Instead, its changes are dominated by total precipitation change.
In addition to intermodel and emissions scenario uncertainties, the future snowfall loss projection is subject to uncertainty resulting from our statistical downscaling technique. We have shown that the statistically downscaled snowfall generally underestimates the dynamically downscaled results. However, the underestimation is small, with a normalized mean bias of a few percentage points. The bias is also relatively minor compared to the percentage loss because of climate change, especially in low and middle elevations. This same bias applies to the snowpack projections. If this underestimate did not exist, the corresponding snowfall and snowpack percentage loss could be a few percentage points less than what we present. However, this does not qualitatively change our main conclusion.
We then project that the percentage reduction of snowpack, represented by 1 April SWE, is larger than that of snowfall, especially at low and moderate elevations. The difference between winter snowfall reductions and 1 April SWE reductions is about 15%–20% points at low and moderate elevations for both periods and both scenarios. In addition to its impacts on winter snow accumulation, warming further enhances snowmelt at these elevations. However, the further reduction is only a few percentage points at high elevations. For low and moderate elevations, about two-thirds of the 1 April SWE loss is due to snowfall reduction, while about one-third is due to enhanced melting. The greater percentage snowfall and snowpack loss at low and moderate elevations in all future climate states probably accounts for the variation in snowfall and snowpack loss across the region’s mountain complexes, which vary in their average elevations.
The effect of snowfall decline on streamflow from mountain snow will be magnified by warming-accelerated melting. A comprehensive assessment of the effect of snowmelt changes on streamflow in the region is beyond the scope of this study. However, it is possible to make meaningful inferences based on simulated snowpack from the dynamically downscaled baseline and midcentury climate under RCP8.5. Figure 11a presents the date when the ground becomes snow-free in the baseline. This snow-free date is defined as the day when SWE reaches a critically low value. A subjective value of 2 mm is used here, though the results are not sensitive to this threshold. The spatial distribution of the snow-free date matches the snowfall distribution. On mountain peaks, the seasonal snow cover disappears from the landscape after 1 June, while at lower elevations, snow cover disappears as early as February. Figures 11b–f show that in all five midcentury dynamical simulations the dates on which snow completely disappears generally occurs earlier than during the baseline period. A spread is evident among the GCMs in how much earlier the snow-free dates occur. On average, the snow-free date occurs 16 days earlier. For each GCM, snowmelt timing is sensitive to winter and spring temperature, with the greatest changes apparent at low elevations, where winter and spring temperatures are warmer. In contrast, significantly earlier snow-free dates are not seen at high elevations, where warming likely has limited impact on snowfall or snowmelt.
Our projections reveal how the stark contrast between the global warming outcomes of the two emissions scenarios by century’s end corresponds to a dramatic difference in snowfall and snowpack outcomes in the mountains of Southern California. From our projections, it is clear that roughly one-third of snowfall and a somewhat greater amount of snowpack are likely to be lost by midcentury, no matter how aggressively greenhouse gas emissions are reduced. By end of century, however, the choice of emissions scenario does make a difference. The amount of snowfall likely to be lost at end of century (roughly half of baseline snowfall), and the corresponding further reduction of the snowpack, can be substantially mitigated by aggressively reducing greenhouse gas emissions in the coming decades.
Support for this work was provided by the City of Los Angeles and the U.S. Department of Energy as part of the American Recovery and Reinvestment Act of 2009. Additional funding was provided by the National Science Foundation (Grant EF-1065863: “Collaborative Research: Do Microenvironments Govern Macroecology?”) and the Southwest Climate Science Center. The authors thank Dr. Sarah Kapnick for reviewing an early draft of this work and the two anonymous reviewers for their invaluable comments.