Projections of twenty-first-century Northern Hemisphere (NH) spring snow cover extent (SCE) from two climate model ensembles are analyzed to characterize their uncertainty. Phase 5 of the Coupled Model Intercomparison Project (CMIP5) multimodel ensemble exhibits variability resulting from both model differences and internal climate variability, whereas spread generated from a Canadian Earth System Model–Large Ensemble (CanESM-LE) experiment is solely a result of internal variability. The analysis shows that simulated 1981–2010 spring SCE trends are slightly weaker than observed (using an ensemble of snow products). Spring SCE is projected to decrease by −3.7% ± 1.1% decade−1 within the CMIP5 ensemble over the twenty-first century. SCE loss is projected to accelerate for all spring months over the twenty-first century, with the exception of June (because most snow in this month has melted by the latter half of the twenty-first century). For 30-yr spring SCE trends over the twenty-first century, internal variability estimated from CanESM-LE is substantial, but smaller than intermodel spread from CMIP5. Additionally, internal variability in NH extratropical land warming trends can affect SCE trends in the near future (R2 = 0.45), while variability in winter precipitation can also have a significant (but lesser) impact on SCE trends. On the other hand, a majority of the intermodel spread is driven by differences in simulated warming (dominant in March–May) and snow cover available for melt (dominant in June). The strong temperature–SCE linkage suggests that model uncertainty in projections of SCE could be potentially reduced through improved simulation of spring season warming over land.
Seasonal snow cover is a crucial component of the climate system, with major impacts on the surface energy budget and water balance. At its peak, snow covers approximately 47 × 106 km2 of Northern Hemisphere (NH) land (about 40% of the land area) each year (Hall 1988; Robinson and Frei 2000). The reflective properties of snow mean that it has a very strong influence on land surface albedo, controlling its seasonal evolution. This high albedo has a cooling influence on climate, with the contribution from terrestrial snow to cryospheric cooling being roughly equal to that of sea ice (Flanner et al. 2011). Furthermore, snow cover can indirectly impact atmospheric circulation (Fletcher et al. 2009; Cohen et al. 2012). Water resources across most NH extratropical (poleward of 30°N) land areas rely on natural water storage provided by snowpack (Diffenbaugh et al. 2013), with approximately one-sixth of Earth’s population dependent on snowmelt for a portion of their water supply (Barnett et al. 2005; Mankin et al. 2015). Earlier spring snowmelt across the western United States has been linked to increased summer heat extremes (Diffenbaugh et al. 2005; Hall et al. 2008) and wildfires (Westerling et al. 2006). Snow cover also has a low thermal conductivity, meaning that it can have an insulating effect on soil temperatures, with important impacts on permafrost extent (Zhang 2005; Zhang et al. 2008; Lawrence and Slater 2010). Variability in snow conditions also has implications for travel and tourism (e.g., Scott et al. 2008). It is crucial, therefore, that we understand how projected warming will affect snow cover.
Extensive climatological snow cover, and relatively high insolation, make the climate system most sensitive to snow and albedo changes during spring (changes during fall and winter are less important because of decreasing insolation across the NH) (Ingram et al. 1989; Hall 2004). Snow albedo feedback (SAF) peaks during March–May (Qu and Hall 2014). Numerous observational studies have shown that Northern Hemisphere spring snow cover extent (SCE) has been decreasing rapidly over recent decades (Groisman et al. 1994; Déry and Brown 2007; Flanner et al. 2009; Brown et al. 2010; Brown and Robinson 2011; Hernández-Henríquez et al. 2015). Most of the loss in snow occurs over “temperature-sensitive regions,” where changes in SCE are closely linked to temperature variability (Karl et al. 1993; Déry and Brown 2007). March–April SCE is decreasing at a rate of −3.4% ± 1.1% decade−1 (1979–2005) (Brown and Robinson 2011), and June SCE has decreased by approximately −18% decade−1 from 1979 to 2014 (Derksen et al. 2015). This rate of decline in SCE exceeds the well-publicized declining trend in September Arctic sea ice (−13% decade−1). It should be noted that the absolute areal SCE loss in March–April is comparable to that observed during June because there is much less snow-covered area in June (section 3a).
By contrast, the suite of climate models contributing to phase 5 of the Coupled Model Intercomparison Project (CMIP5) simulate March–April SCE trends roughly one-third as large as observed for the same time period (−1.0% ± 0.3% decade−1) (Brutel-Vuilmet et al. 2013). This underestimation, also found for the CMIP3 models (Roesch 2006), is associated with underestimated extratropical spring warming (Brutel-Vuilmet et al. 2013). Derksen and Brown (2012) and Mudryk et al. (2014) illustrate other aspects of the observed trends that are not well captured by general circulation models (GCMs).
Northern Hemisphere SCE is expected to continue decreasing under future warming. Several studies over the past 30 years have utilized GCMs to show that SCE decreases substantially in a doubled-CO2 climate (Manabe and Wetherald 1987; Boer 1993; Essery 1997). More recently, the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report (AR5) stated that early spring (March–April) SCE is likely to decrease by 7%–25% (from RCP2.6 to RCP8.5 scenario) by 2100 (Collins et al. 2013; Brutel-Vuilmet et al. 2013). However, this projection was only assigned a medium confidence level because of a large intermodel spread, and a lack of sophistication in the representation of snow processes in many models (e.g., single-layer snowpack, and snow schemes that assume an equal distribution of snow mass across a grid cell; Collins et al. 2013; Slater et al. 2001). Furthermore, no projections of SCE were provided for late spring (May–June), when snow cover is largely restricted to the Arctic, but still represents a significant area (mean 1982–2002 Arctic SCE for May and June was 11.0 and 3.9 × 106 km2, respectively; Brown et al. 2010). The Arctic has experienced the greatest warming in recent decades (Bekryaev et al. 2010), and that trend is expected to continue due to positive feedback mechanisms such as the lapse rate and albedo feedbacks (Pithan and Mauritsen 2014) with implications for spring snow cover.
CMIP (Meehl et al. 2007; Flato et al. 2013) and Snow Models Intercomparison Project (SnowMIP; Etchevers et al. 2004; Rutter et al. 2009) studies have demonstrated that a large intermodel spread exists when simulating snow properties (extent and mass). This limits our confidence in future projections, and is likely to be caused by a combination of internal climate variability and model uncertainty (intermodel variability among GCMs in response to the same forcing) from physical processes that are either missing or oversimplified (Hawkins and Sutton 2011). For example, simulated snow mass (and similarly SCE) is sensitive to different parameterizations for snowfall, albedo, snow–vegetation masking, and snow cover fraction (see Slater et al. 2001; Bartlett et al. 2006; Dai 2008; Rutter et al. 2009; Essery et al. 2013).
Internal climate variability arises primarily from naturally occurring nonlinear atmospheric and oceanic processes, and their interactions (Deser et al. 2012b; Kay et al. 2015). The processes can be nearly instantaneous, or take several years (Hegerl et al. 2007). Internal variability influences projected regional trends in temperature and precipitation, even in the presence of a background trend in CO2 forcing (Hawkins and Sutton 2009, 2011; Deser et al. 2012b), both of which are crucial factors for future winter snow accumulation patterns (Räisänen 2008; Krasting et al. 2013; Mankin and Diffenbaugh 2015; Shi and Wang 2015). Therefore, it is likely that internal climate variability could also influence projected trends for spring SCE.
The primary goal of this work is to investigate the spread in twenty-first-century changes to spring snow cover as projected by the CMIP5 suite of climate models. We use the recent past to help understand the spread in trend strength between simulations and observations. We also seek to determine the influences of both internal variability and model uncertainty in these simulations, to answer the following research questions: (i) What impact does the representation of snow–climate processes [e.g., the sensitivity of snow cover to warming (henceforth, snowmelt sensitivity)] have on simulations of SCE? (ii) What are the respective roles of temperature and precipitation changes in governing SCE trends? The data and methods are described in section 2. In section 3, we present historical and projected SCE trends, and diagnose the respective roles of model uncertainty and internal variability. Last, section 4 highlights the key findings of this research and provides a discussion of how our findings relate to previous research.
2. Data and methods
a. Climate model data
We use monthly mean output from the suite of historical (1850–2005) and future (2006–2100) simulations from the CMIP5 archive (Taylor et al. 2012; http://cmip-pcmdi.llnl.gov/cmip5/) to evaluate snow cover in 15 models (Table 1). Only models that archived the variable snc [snow cover fraction (SCF)] are included in this analysis so as to avoid introducing uncertainty through the estimate of SCE from snow water equivalent (snow mass, the variable snw) as done in previous studies (e.g., Roesch 2006). Future snow cover projections are forced using the RCP8.5 scenario because it most closely resembles the observed emissions pathway over the past decade (Peters et al. 2013). We use all available realizations (n = 1, …, 10) from each of the 15 CMIP5 models for a total of 61 historical runs and 41 runs for the representative concentration pathway 8.5 (RCP8.5) scenario. We compute individual trends and averages for each realization and then take the interrealization average across each model to calculate ensemble means. These values are then used to determine the CMIP5 multimodel mean values.
We also use output (variables snc, snw, tas, and psl) from a large (50 realization) ensemble of the Second-Generation Canadian Earth System Model (CanESM2; Arora et al. 2011). To produce the 50-member Canadian Earth System Model–Large Ensemble (CanESM-LE) 10 realizations are initiated in 1950 from each of the 5 original realizations of CanESM2 submitted to CMIP5. Each new ensemble member is perturbed by changing the seed of a random number generator used in the parameterization of radiative transfer through clouds. Following the CMIP5 design, historical forcing is applied from 1950 to 2005, followed by RCP8.5 from 2006 to 2100. This methodology is appropriate for sampling the statistics of climate variability within CanESM2, because the initial ocean conditions down to 300-m depth have very little influence on the simulation after 5–10 years (Branstator and Teng 2012). This implies that, by the time our analysis period begins in 1981, CanESM-LE represents 50 statistically independent climate states.
As a result, CanESM-LE is used to quantify the component of internal variability within future projections from a single GCM, while the CMIP5 ensemble includes a combination of model uncertainty and internal variability. Similar large ensembles have previously been used to separate the components of future climate patterns related to forced and internal variability (Deser et al. 2012a,b; Wettstein and Deser 2014; Swart et al. 2015). Consistent with previous research (Deser et al. 2012b, 2014), we estimate the “forced response” to greenhouse gas (GHG) forcing as the ensemble mean response from all 50 realizations of CanESM-LE. As in Hawkins and Sutton (2011) we estimate the contribution of “internal variability” to each realization by subtracting the ensemble mean of a particular quantity from the values in that realization. This approach is effective when there are enough simulations so that climate noise can be sufficiently diluted (Deser et al. 2012b, 2016).
b. Observational data
Seven observation-based estimates (Table 2) are used to evaluate the CMIP5 models during recent decades (1981–2010). The use of an observational ensemble reduces the chance of incorrectly identifying a model bias due to errors in a single observational analysis. The seven observation-based estimates are derived from reanalyses, satellite data, and in situ measurements. Each dataset must have complete Northern Hemisphere coverage and a data record spanning the 1981–2010 period. Those datasets that meet this criteria are 1) the NOAA snow chart climate data record (NOAA CDR) (Brown and Robinson 2011; ftp://eclipse.ncdc.noaa.gov/cdr/snowcover/), 2) the Brown snow cover product derived from a combination of climate station data and a simple snow model (Brown et al. 2003), 3) the combined in situ and satellite passive microwave derived Global Snow Monitoring for Climate Research (GlobSnow) SWE product (Takala et al. 2011; www.globsnow.info), 4) SWE from the Modern-Era Retrospective Analysis for Research and Applications (MERRA; Rienecker et al. 2011), 5) SWE from the European Centre for Medium-Range Forecasts (ECMWF) interim reanalysis global land surface dataset (ERA-I-Land; Balsamo et al. 2015), 6) SWE from the Global Land Data Assimilation System, version 2 (GLDAS-2), product (Rodell et al. 2004), and 7) SWE output from the snowpack model Crocus driven by ERA-Interim meteorology (Brun et al. 2013). Three of the snow products use the same atmospheric forcing data from ERA-Interim (Brown dataset, Crocus, and ERA-I-Land). However, despite this similarity, they exhibit very different SCE trends due to differences in the snow parameterizations between Crocus and HTESSEL (ERA-I-Land), while the Brown dataset only uses temperature and precipitation to force a simple snow model (Table 2).
The NOAA CDR is derived primarily from optical satellite data (Brown and Robinson 2011). This dataset provides monthly fractional snow cover, which is calculated as the percent of days per month that a grid cell is at least 50% snow covered. The Brown dataset (Brown et al. 2003) uses ground-based snow measurements and a simple snowpack model to produce SCF from daily SWE thresholds exceeding 4 mm. The five remaining datasets (MERRA, ERA-I-Land, GLDAS-2, GlobSnow, and Crocus) were used in the SWE product intercomparison described in Mudryk et al. (2015). For these products, SWE is initially interpolated to a 1° × 1° grid and SCF is then derived from daily SWE thresholds exceeding 4 mm. The 1981–2010 time period is a shorter record than available from individual datasets (e.g., the NOAA CDR starts in 1967) but the compromise in time series length is mitigated by the advantages of a multidataset perspective, which has typically not been used in previous snow–climate studies.
An observational ensemble of temperature is used to determine spring snow extent sensitivity. We select five datasets for temperature: the Climatic Research Unit (CRU) land station temperature database (Jones et al. 1999; 2012), the Goddard Institute for Space Studies (GISS) surface temperature analysis (Hansen et al. 2010), the National Climatic Data Center (NCDC) temperature product (Smith et al. 2008), and the National Centers for Environmental Prediction (NCEP) surface and 2-m (NCEP2m) temperature datasets (Kalnay et al. 1996).
c. Analysis methods
The CMIP5 models output data at a variety of resolutions (1°–3° grid boxes), and to account for this we must create a consistent land–sea mask such that land area biases are reduced [particularly in the Canadian Arctic Archipelago (CAA), where many narrow channels may not be resolved at coarser resolution; Laliberté et al. 2016]. For each CMIP5 model a land–sea mask extracted from the MERRA product is remapped to the native model resolution to isolate simulated NH land-only snow cover and temperature data. This ensures that we reduce discrepancies in the amount of land area between models (mainly in the CAA). Using this mask along with the model-specific land mask reduces climatological SCE in the models but has minimal impact on their trends (not shown). Furthermore, the study area for this analysis is spatially restricted to the Northern Hemisphere extratropics (>30°N), with Greenland excluded, and temporally restricted to March–June (MAMJ).
SCE is calculated from model output by multiplying grid cell snow cover fraction (%) by the area of each grid cell (m2) and then taking the hemispheric sum for each month. Similarly, snow water mass (SWM) is computed by multiplying grid cell level snow water equivalent (SWE) by the area of each grid cell and summing over the NH. This is applied to SWE data from the CanESM-LE to allow for an illustration of the influence that changes in winter precipitation have on SCE trends. The premelt SWM is a useful measure of snowfall totals over the winter months, particularly across the Arctic, because wind-driven snow processes are not represented in current models (Turner et al. 2006; Lawrence et al. 2012). We find that 1981–2010 winter (October–March) snowfall trends are strongly correlated with March SWM trends within CanESM-LE (r = 0.92; not shown). Along with measures of correlation, we also use the coefficient of determination R2 to recognize the contribution from precipitation and temperature to SCE variability.
Time series of SCE and SWM data are used to calculate climatologies and linear trends. We calculate these values for each of the four 30-yr climatological time frames during the study period: historical (1981–2010), near-future (2011–40), midcentury (2041–70), and long-term (2071–2100). In some cases, a twenty-first-century trend (2011–2100) is used to simplify the discussion of results. Since the historical CMIP5 data end in 2005, we use RCP8.5 data to complete the 1981–2010 period so that a comparison with recent observations can be made [similar to Derksen and Brown (2012)]. SCE trends are computed as both absolute area (106 km2 decade−1) and percent changes (% decade−1). Absolute area calculations are useful in the context of comparing different months, whereas percent changes may be more suitable for interensemble comparisons because they account for potential differences in snow cover climatology. Throughout, all trend values reported are accompanied by one standard deviation (1σ) to represent uncertainty. Note that trends are calculated at each model’s native resolution, and regridding to a 1° × 1° grid is only used for spatial mapping of the snow cover from the CMIP5 models. Last, the term “bias” will be used solely for comparing models in relation to observation-based estimates.
a. Historical spring SCE trends
Considering first the entire spring season for the Northern Hemisphere, SCE has decreased at a mean rate of −0.55 ± 0.21 × 106 km2 decade−1 from 1981 to 2010, according to the seven observation-based estimates evaluated here (Table 2). Dividing this rate by the climatological spring SCE produces a percent change of −3.3% decade−1. The strongest trend in terms of absolute area occurs in March (mean: −0.66 ± 0.26 × 106 km2 decade−1) and the weakest in June (−0.41 ± 0.30 × 106 km2 decade−1) (Fig. 1). However, direct measures of trend magnitude do not account for the much greater total snow area in March (32.6 ± 2.5 × 106 km2) than June (2.6 ± 1.9 × 106 km2; not shown). When viewed as a percent change relative to the monthly climatology, March SCE is decreasing at a rate of −2.0% ± 0.8% decade−1, whereas the rate of June SCE loss is −16% ± 11% decade−1 (not shown). Both early spring and June trends found here are weaker than those previously reported that were based only on the NOAA CDR; Brutel-Vuilmet et al. (2013) reported trends of −3.4% decade−1 for March–April over the 1979–2005 period, and Derksen and Brown (2012) reported trends of −18% decade−1 in June over the 1979–2011 period.
These reported differences result because there is a substantial spread among the observation-based estimates of SCE trends, and of the seven products evaluated here the NOAA CDR trends are the largest in magnitude (Table 2). Mudryk et al. (2015) have recently shown an analogously large spread in SWM trends from various snow analysis products. Such spread likely occurs as a result of differences in methodology and the type of data used to construct each dataset (e.g., in situ, reanalysis, satellite-derived). Three of the four reanalysis products (GLDAS-2, ERA-I-Land, and MERRA) exhibit the weakest spring SCE trends over recent decades, with GLDAS-2 losing the least SCE in each spring month. On the other hand, two of the three products that utilize either satellite-derived or in situ information (NOAA CDR and Brown dataset) exhibit the strongest trends in spring SCE over the recent past.
Simulated spring trends from the CMIP5 models are approximately 22% weaker than observed on average [multimodel ensemble mean (MM): −0.43 ± 0.17 × 106 km2 decade−1 or −2.5% ± 1.0% decade−1]. This is also weaker than that of the CanESM-LE (mean = −0.62 ± 0.18 × 106 km2 decade−1), which demonstrates much greater late-spring snow loss. From a monthly perspective, the CMIP5 mean SCE trend is largest during April and weakest in June (Fig. 1). The agreement between CMIP5 models and observations is very good during April and May, but less so during March and June, when the models have weaker trends (MM: −1.4% ± 0.8% decade−1 and −8.4% ± 5.4% decade−1, respectively). However, March is the only month with a statistically significant difference between the observed and simulated mean trends (p < 0.05; using a two-sided Student’s t test). Furthermore, nearly all models exhibit negative SCE trends throughout the spring, except for INM-CM4.0, which has a slight increasing trend during March. The CMIP5 models range from losing very little snow in MAMJ at 13% of the observed baseline (−0.07 × 106 km2 decade−1) to 142% of observed (−0.78 × 106 km2 decade−1) (Table 1). As a demonstration that single model contributions to the CMIP5 archive may underestimate internal variability, CanESM2 is the model with the greatest spring snow loss of any CMIP5 model from 1981 to 2010, yet we find that CanESM-LE (the large ensemble produced using the same model) is much closer to the CMIP5 average (Fig. 2).
The variability in 30-yr trends from the CMIP5 ensemble is equally large for March, April, and May (standard deviation σ = 0.24 × 106 km2 decade−1) and slightly lower in June (σ = 0.17 × 106 km2 decade−1). However, when we examine 10-yr trends the spread widens dramatically for all spring months. For example, in May the simulated range (maximum–minimum) for decadal trends is more than 7 times that for 30-yr trends (from −3.7 to 2.7 × 106 km2 decade−1; see Fig. S1 in the supplemental material). This highlights the larger contribution to the SCE trends from internal variability, compared to the forced response to GHG increases, over shorter time periods. This result is consistent with similar findings for temperature and precipitation trends (Hawkins and Sutton 2009, 2011). The following section will investigate possible factors contributing to the large intermodel spread within historical simulations from the current generation of climate models.
b. Sources of model uncertainty: Historical trends
1) Sensitivity of SCE to warming trends
First we evaluate how differences in warming could be affecting the simulated intermodel spread in SCE trends. Although a very important contributor (R2 = 0.74), differences in simulated warming do not explain all of the intermodel spread in spring SCE trends for the 1981–2010 period. We use mean extratropical land warming rather than local warming because unnaturally high sensitivities can occur for some models in areas where the warming trend is close to zero. The observed spring NH extratropical land warming trend over the 1981–2010 period is 0.34 ± 0.04 K decade−1. On average, the CMIP5 models accurately capture spring warming over recent decades, with an ensemble mean of 0.36 ± 0.11 K decade−1 (Fig. 2). In contrast, the CanESM-LE has a greater mean warming trend (0.52 ± 0.08 K decade−1), which overlaps with the CMIP5 warming trend but falls outside the uncertainty range of the observed warming trend. The mean trend seen in the CanESM-LE ensemble reflects the majority of CMIP5 models (10 out of 15) that simulate recent (1981–2010) spring warming that is greater than or equal to the warming found in observations.
Despite realistically reproducing observed temperature trends, only two models produce more snow loss than the observation-based estimates. This suggests that snow in some models tends to be less sensitive to temperature variations than in observations. To quantify this property of the models, we compute a snowmelt sensitivity λsmelt = 〈ΔSCE〉/〈ΔTs〉, which measures how much SCE is reduced per degree of warming, averaged over the NH land area (averaging is denoted by the angle brackets). Observed spring λsmelt during the spring months (MAMJ mean) is −1.62 ± 0.61 × 106 km2 K−1, with the large uncertainty driven mainly by disagreement among the observed SCE trends (illustrated by the shaded rectangle in Fig. 2), which creates a large spread among the 35 possible combinations of observed temperature (5 products) and SCE (7 products). This exceeds the estimate of λsmelt that is computed for both the CanESM-LE (−1.18 ± 0.15 × 106 km2 K−1), and the CMIP5 ensemble (MM = −1.19 ± 0.31 × 106 km2 K−1). The 68% confidence interval for the CMIP5 estimate of λsmelt overlaps that from observations, indicating that the two estimates of mean snowmelt sensitivity are not significantly different. However, similar to Brutel-Vuilmet et al. (2013), we find that the weaker-than-observed SCE trends from the CMIP5 ensemble are likely due in part to a weaker-than-observed snow response to warming.
2) Climatological mean snow cover
A secondary cause of the weaker-than-observed SCE trends in CMIP5 is biases in the simulated climatological (1981–2010 mean) snow cover () for a given month. The ability to accurately represent present-day snow cover is important for simulated SCE trends because of a positive correlation between snow extent and SAF strength (Levis et al. 2007). This relationship indicates that models with greater SCE produce stronger SAF for a given rise in temperature, because larger SCE implies a greater potential area over which albedo can be reduced from its snow-covered to its snow-free value. In the CMIP5 multimodel mean, nearly all land poleward of 45°N is at least 50% snow covered in March (Fig. 3a; MM of ~30.5 ± 3.5 × 106 km2). However, there is a significant spread in March within the CMIP5 ensemble: the model with the highest has 18% more snow-covered area than the mean (red line in Fig. 3a), whereas the model with the lowest has 28% less snow cover than the mean (green line in Fig. 3a). Much of the disparity between these models is found across western North America, western Eurasia, and the Tibetan Plateau (similarly for April; Fig. 3b). Comparatively, the observation-based estimates show ranging from 30 to 36 × 106 km2 in March, with a mean of approximately 32 × 106 km2. Disparity between the minimum and maximum observation-based snow cover products is greatest over eastern Eurasia and western North America (Fig. 3a). Of note is the good agreement over western Eurasia, where both of the extreme observation-based estimates exceed or closely resemble the maximum model extent. This implies that the CMIP5 models may be systematically underestimating early spring snow cover in this region.
Late spring (May–June) snow cover resides primarily across the Arctic (>60°N), with much of the high latitudes still more than 50% snow covered during May (Fig. 3c). On average, the CMIP5 models simulate May of 11.8 ± 3.9 × 106 km2, while the observation-based products range from 6.9 to 14.7 × 106 km2, with a mean of 10.6 ± 2.6 × 106 km2. In June any remaining snow cover is restricted to Siberia, Arctic Canada, and Alaska and is characterized by local snow cover fractions lower than 50% (CMIP5 mean is 3.2 ± 2.0 × 106 km2; Fig. 3d). Similar to the models, the observation-based estimates have a mean of 2.6 ± 1.9 × 106 km2, and a large spread spanning from 0.5 (GLDAS-2) to 5.9 × 106 km2 (NOAA CDR). Note that differences between models in their values is likely the result of model uncertainty, rather than internal variability, as demonstrated by a very small range within the CanESM-LE (<1 × 106 km2 for all spring months; not shown).
Biases in June have the most significant impact on SCE trends of any month: those models with minimal SCE in June tend to show very weak SCE trends because in the future there is so little snow left to melt. The models with low June (e.g., BCC_CSM1.1, CSIRO Mk3.6, INM-CM4.0, and MPI-ESM-LR; see Fig. S2 in the supplemental material) exhibit a mean SCE trend of only −0.06 × 106 km2 decade−1 (not shown), a factor of 6 weaker than the other CMIP5 models (−0.35 × 106 km2 decade−1; not shown). These same models with low June have previously been shown to have mean late spring near-surface air temperatures that are substantially warmer than the other CMIP5 models (Thackeray et al. 2015). Therefore, biases in can affect the SCE trend in seasons when SCE becomes limited (e.g., late spring).
There is a slightly weaker correlation (r = 0.43) between trends in 1981–2010 March SWM (used as a proxy for variability in winter snowfall) and MAMJ SCE trends within the CMIP5 ensemble (not shown). A moderately strong correlation also exists between March SWM trends and λsmelt (r = 0.70). However, because SWM does not have the same subseasonal importance as (there is a weak correlation with SCE trends for all months other than March) it is not investigated further. Furthermore, although there are many other potential sources of model uncertainty in simulated SCE (e.g., model resolution, land surface scheme complexity, and climatological temperature biases), we do not find any clear linkages between these parameters and spring SCE trends, so they are not discussed further.
c. Projected trends in spring SCE
To evaluate projections of future spring SCE, we use two model ensembles: the multimodel CMIP5 ensemble and the CanESM-LE. As previously noted, the latter only contains (land–atmosphere–ocean induced) spread due to internal variability, so it provides a useful benchmark to compare with the estimate of intermodel spread from CMIP5. It should be noted that CanESM-LE provides one model’s estimate of internal variability, which could vary for other CMIP5 models (e.g., Kay et al. 2015). First, we discuss the median projected change for the spring as a whole, followed by early (March–April) and late spring (May–June). On average, the CMIP5 models project that spring SCE trends will strengthen during the twenty-first century relative to the recent past. The mean rate of spring snow loss over the twenty-first century (2011–2100 trend) is approximately −3.7% ± 1.1% decade−1, 33% greater than in the period 1981–2010 (Fig. S3 in the supplemental material). Similarly, the CanESM-LE exhibits a strengthening of 41% (more negative) compared to its 1981–2010 trend. However, CanESM2 exhibits the strongest 1981–2010 trend of the CMIP5 models (Table 1) so the median rate of twenty-first-century spring snow loss from CanESM-LE is also larger (−4.3% decade−1).
The CMIP5 models project that early spring (March–April) SCE trends will strengthen in the twenty-first century relative to the recent past (Figs. 4a,b). The mean SCE responses over the twenty-first century (2011–2100 trend) from CMIP5 (−0.80 ± 0.23 × 106 km2 decade−1) and CanESM-LE (−1.02 ± 0.22 × 106 km2 decade−1) are more than 65% stronger than their respective simulated rates for the period of 1981–2010. This is consistent with greater simulated rates of warming during the twenty-first century (not shown). Since these two ensembles have different mean (Table 3), we also calculate the ensemble mean percentage of snow loss over the twenty-first century to account for differences in the amount of snow cover available for melt: the CMIP5 models lose −3.0% ± 0.9% decade−1, while CanESM-LE loses 3.3% ± 0.7% decade−1.
Within the CMIP5 ensemble, mean twenty-first-century May SCE loss (2011–2100 trend) is projected to strengthen slightly (by ~25%) compared to the 1981–2010 trend. Trends in May also exhibit the greatest discrepancy between CMIP5 (−5.5% ± 2.0% decade−1) and CanESM-LE (−6.5% ± 1.3% decade−1) (Fig. 4c). Unlike the other months, June SCE trends are projected to weaken over the course of the twenty-first century (Fig. 4d). A gradual weakening within CanESM-LE is tied to a significant reduction in the amount of snow area remaining for melt (mean SCE < 0.5 × 106 km2 by 2071–2100). This same reasoning explains why the simulated June trends (even under the most aggressive GHG forcing scenario) of −8.1% ± 4.9% decade−1 from CMIP5 and −8.7% ± 1.9% decade−1 from CanESM-LE are weaker than observed in recent decades (−16% decade−1). In summary, rates of projected snow cover loss are expected to increase for all spring months, with the exception of June, where a relatively small SCE remains in the latter half of the twenty-first century.
Last, we consider the intermodel spread of land surface warming trends as a possible source of uncertainty for projected spring SCE trends within the CMIP5 ensemble. We estimate the uncertainty in the CMIP5 projections using the multimodel standard deviation (Hawkins and Sutton 2011), and we average σ over three different 30-yr periods (2011–40, 2041–70, and 2071–2100). This procedure yields uncertainties in SCE trends for March (σ = 0.28 × 106 km2 decade−1), April (σ = 0.36 × 106 km2 decade−1), May (σ = 0.34 × 106 km2 decade−1), and June (σ = 0.18 × 106 km2 decade−1). Taking the example of the relatively large uncertainty in April, we find a wide range of projected SCE trends for the 2011–40 period, from a small gain in one model to a loss of −1.6 × 106 km2 decade−1 in another model (Fig. S4 in the supplemental material). In this case, the model with the strongest (weakest) SCE loss also warms the most (least) over this period. Warming trends explain much of the intermodel spread in early spring SCE trends (R2 = 0.79; Fig. 5), whereas June SCE trends are heavily influenced by biases in (R2 = 0.93; Fig. 5). For the spring as a whole, variability in twenty-first-century NH land warming explains approximately 80% of the intermodel spread in SCE trends. Therefore, reducing variability in simulated future warming should in turn reduce uncertainty in SCE trends (further discussed in section 4).
d. The contribution of internal variability to projected trends in SCE
Many of the CMIP5 models that project extremely strong or weak spring SCE trends contributed only a single realization in the CMIP5 archive, while considerably better agreement in projected SCE trends exists among the set of four models that contributed n ≥ 5 realizations (particularly in early spring when warming trends dominate; not shown). This motivates an important question as to the role of internal variability in SCE trends; however, the majority of CMIP5 models completed fewer than five realizations, which is likely insufficient for estimating internal variability (Kay et al. 2015).
The 50-realization CanESM-LE exhibits a smaller spread in MAMJ SCE trends throughout the twenty-first century than the CMIP5 ensemble (σ = 0.18 and 0.29 × 106 km2 decade−1, respectively, averaged over three epochs 2011–40, 2041–70, and 2071–2100). Internal variability, as indicated by the shading for CanESM-LE in Fig. 6, is likely a very important contributor to the intermodel spread in the near term (2011–40), but the fraction of total variance within the CMIP5 ensemble attributable to internal variability decreases on longer time scales as a relatively larger fraction is explained by model uncertainty (Fig. 6). This same finding has also been shown for precipitation and temperature trends, where internal variability has a greater influence in the near future than at the end of the century (Hawkins and Sutton 2009, 2011).
Last, to demonstrate the interplay between internal variability and model uncertainty, we can compare the results from CanESM-LE with the intraensemble spread for all CMIP5 models with n ≥ 5 realizations available for the RCP8.5 experiment (CanESM2, CCSM4, CNRM-CM5, and CSIRO Mk3.6.0). This subset of models, which we assume provides an improved estimate of model uncertainty isolated from internal variability, contains substantial variation, both in the intermodel SCE trends (n = 4, σ = 0.22 × 106 km2 decade−1) and the interrealization variability (n = 26, minimum and maximum σ are 0.12 and 0.23 × 106 km2 decade−1). The interrealization spread of trends for each model is thus of similar magnitude to the intermodel spread, making it plausible that a significant fraction of the intermodel spread is caused by internal variability. We find considerable similarities between this analysis and the work on trends in September Arctic sea ice extent (SIE) by Swart et al. (2015). For example, there is a remarkable similarity in the contributions from internal variability and model uncertainty to projected trends of SIE and SCE (not shown). The conclusion for near-term projections is that the large contribution from internal variability presents a challenge to determining the physical cause of 30-yr SCE trends.
Case study: May SCE trends in CanESM-LE
In the CanESM-LE, near-future (2011–40) springtime (MAMJ) SCE trends range from −0.26 to −1.08 × 106 km2 decade−1, with the largest monthly spread occurring in May (from −0.42 to −1.49 × 106 km2 decade−1). For May, this represents a more than doubling of the range exhibited by the five CanESM2 runs contributed to CMIP5 (from −0.69 to −1.11 × 106 km2 decade−1). We will therefore use May as a case study for better understanding the primary physical factors contributing to the spread.
First, we examine the contribution from variations in trends of near-surface air temperature. Only a relatively small fraction (17%) of the interrealization variability in CanESM-LE projected near-future May SCE trends is explained by annual mean global surface warming (land plus ocean), with r = −0.41 (Table 4) and the negative sign implying that enhanced global warming is associated with greater snow loss. However, nearly half of the interrealization variability (45%) is explained when we restrict the analysis to include only contemporaneous (e.g., May) and local (e.g., NH extratropical land averaged) temperatures (Fig. 7a; r = −0.67). The majority of May snow cover resides across the Arctic (Fig. 3c), so one would expect an even stronger correlation with temperature there if local warming were the only contributor to differences in SCE trends. Yet, when we restrict the temperatures to the Arctic (>60°N) region, the relationship becomes only slightly stronger (r = −0.69). This demonstrates that differences in simulated warming cannot fully account for variability within CanESM-LE SCE trends, and so next we examine the roles for changes in precipitation and atmospheric circulation.
Whereas in CMIP5 we find that June SCE trends are highly correlated with the intermodel spread in June climatological SCE (Fig. 5), the spread in climatological SCE in CanESM-LE is minimal (<1 × 106 km2) for all spring months and there is no correlation with SCE trends (not shown). However, we do find a relationship in CanESM-LE between spring SCE trends and snow accumulation during the previous winter. We use March SWM as a proxy for simulated snowfall totals over the winter months, and find a weak positive correlation (r = 0.44 for 2011–40; see Table 4; note that the R2 values for other climatological periods are shown in Fig. 7b) between twenty-first-century trends in SWM and May SCE. To illustrate the importance of SWM for interrealization differences in SCE trends, in Fig. 7 we compare May land warming, May SCE loss and March SWM loss over the NH extratropics for the period 2011–40. Across the 50 realizations May land warming varies from 0.35 to 0.78 K decade−1, and the realization that warms the most (run 40) also produces the greatest SCE loss (Fig. 7a). However, the realization with the weakest SCE loss (run 1) is not the realization with the least warming (run 49). The reason is that a weaker-than-average decreasing trend in March SWM in run 1 contributes to a weaker-than-average SCE trend (Fig. 7b). Positive SWM anomalies extend snow cover duration because greater melt energy is required to remove deeper snow.
Last, we explore the role of local temperature and atmospheric circulation changes in explaining the interrealization spread in CanESM-LE SCE trends. We compute local correlations between near-future May SCF trends and contemporaneous temperature and sea level pressure (SLP) changes. Trends in temperature and snow cover have a very strong negative association over most NH areas with substantial May snow cover (Fig. 8a). For SLP, the relationship with snow cover is of hemispheric spatial scale, with moderate correlations of either sign that project onto the North Atlantic Oscillation (NAO) pattern (Barnston and Livezey 1987) (Fig. 8b). Greater Eurasian snow loss is associated with an increased meridional pressure gradient across the North Atlantic typical of a positive NAO phase, bringing enhanced warm advection into Eurasia. In contrast, contributions to changes in western North American SCE largely stem from North Pacific SLP patterns. Figure 8b consequently represents a combination of several unique circulation patterns. This demonstrates that atmospheric circulation responses associated with internal variability exert an influence on near-term SCE trends over much of NH land.
4. Discussion and conclusions
This study uses seven observation-based estimates of snow cover, five surface temperature datasets, and two climate model ensembles to characterize the uncertainty in simulations of NH spring snow cover extent. We find that weaker than observed historical (1981–2010) SCE trends from the CMIP5 ensemble can be partially explained by biases in climatological spring snow extent within these models. However, biases in simulated SCE trends during recent decades are much smaller than previously shown from studies that relied on a single observation-based reference dataset (Derksen and Brown 2012; Brutel-Vuilmet et al. 2013). These studies used the NOAA CDR because of its long time series (1967–present), which we find to have the strongest spring SCE trend of the seven observation-based estimates. SCE in some models appears to lack sensitivity to warming, but the ensemble means are not significantly different.
Spring snow cover is projected to decrease by −3.7% ± 1.1% decade−1 within the CMIP5 ensemble over the twenty-first century. This represents a strengthening of 33% relative to the rate simulated over recent decades (1981–2010). Projected snow cover loss is expected to increase for all spring months over the twenty-first century, with the exception of June (when nearly all remaining snow has melted by the latter half of the twenty-first century). For 30-yr spring SCE trends over three time periods in the twenty-first century (2011–40, 2041–70, and 2071–2100), we find that internal variability, as estimated from the CanESM initial condition ensemble (CanESM-LE; σ = 0.18 × 106 km2 decade−1), is substantial but smaller than the intermodel spread from CMIP5 (σ = 0.29 × 106 km2 decade−1). In contrast, the spread in SCE trends from CanESM-LE and CMIP5 are very similar for the historical period (Fig. 1). The main physical drivers of intermodel differences in projected spring SCE trends are differences in simulated warming trends (R2 = 0.80) and biases in mean SCE, with the latter more important in late spring. In theory, a reduction in the variability of projected warming should lead to a decrease in the spread of spring SCE trends. Internal variability is a major contributor to intermodel spread (total variance) in the near-term, but the fraction of total variance attributable to internal variability decreases on longer time scales because of greater model uncertainty. We find large internal variability in near-term (2011–40) warming trends over NH extratropical land, which explains almost 50% of the variability in projected SCE trends, even in the presence of a strong forced trend from GHGs. Furthermore, internal variability in winter snowfall trends has a significant (but lesser) impact on SCE trends (R2 = 0.20).
There are a number of ways to potentially reduce the uncertainty in projections of NH SCE. The first involves increasing the number of realizations from each model as a part of future modeling efforts. Following the approach of Deser et al. (2012b), we calculate the minimum number of realizations Nmin required to detect the near-future forced May SCF trend at the 5% significance level, given by Nmin = 8/(X/σ)2, where X is the ensemble mean trend, and σ is the standard deviation of the 50 trends. Regions with stronger snow responses generally need between 3 and 10 realizations to detect a significant trend, whereas areas with weaker responses (eastern Siberia and Canadian Arctic; Fig. 9a) require upward of 50 realizations (Fig. 9b). The implication here is not that hundreds of realizations are necessary, but that over some regions the near-term forced response is so weak that it cannot be captured. However, a majority of the models contributing to CMIP5 provided fewer than three realizations for RCP8.5.
Second, there is a very strong relationship between projected spring SCE trends and warming trends (R2 = 0.80). However, under RCP8.5 the CMIP5 models exhibit a rather large spread in twenty-first-century spring warming (from 0.39 to 0.95 K decade−1). Therefore, a reduction in the uncertainty of the forced component of projected warming could lead to a decrease in the spread for spring SCE trends (the component due to internal variability is essentially random, and therefore unconstrained). Previous research has shown that 40%–50% of the spread in CMIP5 twenty-first-century spring warming over NH extratropical land can be explained by variability in simulated snow albedo feedback (SAF) (Qu and Hall 2014). Furthermore, Thackeray and Fletcher (2016) demonstrated that selecting only models with SAF closest to observed estimates reduces the spread in CMIP5 twenty-first-century NH land warming by about 40%. Therefore, model development focused on alleviating process-level biases—particularly those related to SAF— could help to reduce model uncertainty in future projections of warming and snow cover.
We acknowledge funding from the Natural Sciences and Engineering Research Council of Canada’s Climate Change and Atmospheric Research Initiative via the Canadian Sea Ice and Snow Evolution (CanSISE) Network (Grant 210832161). We also acknowledge Environment and Climate Change Canada’s Canadian Centre for Climate Modelling and Analysis for executing and making available the CanESM2 Large Ensemble simulations used in this study. We thank Eric Brun for providing data from the Crocus snowpack model, and Ross Brown for providing the Brown dataset. We also thank three anonymous reviewers for their helpful comments.
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-16-0341.s1.