This study assesses the ability of a newly developed high-resolution coupled model from the Geophysical Fluid Dynamics Laboratory to simulate the cold-season hydroclimate in the present climate and examines its response to climate change forcing. Output is assessed from a 280-yr control simulation that is based on 1990 atmospheric composition and an idealized 140-yr future simulation in which atmospheric carbon dioxide increases at 1% yr−1 until doubling in year 70 and then remains constant. When compared with a low-resolution model, the high-resolution model is found to better represent the geographic distribution of snow variables in the present climate. In response to idealized radiative forcing changes, both models produce similar global-scale responses in which global-mean temperature and total precipitation increase while snowfall decreases. Zonally, snowfall tends to decrease in the low to midlatitudes and increase in the mid- to high latitudes. At the regional scale, the high- and low-resolution models sometimes diverge in the sign of projected snowfall changes; the high-resolution model exhibits future increases in a few select high-altitude regions, notably the northwestern Himalaya region and small regions in the Andes and southwestern Yukon, Canada. Despite such local signals, there is an almost universal reduction in snowfall as a percent of total precipitation in both models. By using a simple multivariate model, temperature is shown to drive these trends by decreasing snowfall almost everywhere while precipitation increases snowfall in the high altitudes and mid- to high latitudes. Mountainous regions of snowfall increases in the high-resolution model exhibit a unique dominance of the positive contribution from precipitation over temperature.
Snow serves as a natural water reservoir when it collects on the land surface as snowpack. Warm-season snowmelt releases water into the soil and produces runoff with the potential for far-reaching negative effects if the timing or amount of snowmelt is altered. Snowpack replenishes soil moisture when snowmelt begins; lower snowfall values can thus result in earlier onset of soil moisture drying in the springtime (Mahanama et al. 2012; Manabe et al. 1981). Because of this mechanism, earlier snowmelt timing has been linked to droughts (Mahanama et al. 2012; Mishra et al. 2010) and increased fire risk (Westerling et al. 2006). Reductions in snowpack due to lower snowfall values can also alter the amount of water availability for hydroelectric power generation in spring and summer months (Vicuna et al. 2008; Christensen et al. 2004). These implications drive a societal need for sound projections of future snow metrics in a warming climate.
The sign of the response of snow variables is a nuanced combination of temperature and precipitation changes. Absent any changes in total precipitation, uniform increases in temperature will lead to more precipitation falling as rain instead of snow, reducing total snowfall. Adding in a theoretical uniform increase in precipitation globally will result in increased snowfall only in regions where temperature remains below freezing. Actual changes in future precipitation amount and distribution with respect to temperature can therefore either enhance or negate snowfall changes that result from warming temperature alone. Adding complexity to this thought experiment, mountains serve as special isolated points where surface temperature is naturally lower (and often below freezing at the highest peaks) than surrounding low-lying areas because of the mean vertical temperature profile of the troposphere. Resolving orography in model simulations is therefore a necessary component to fully capture how snowfall will change spatially in a warming climate. The amount of snow left on the ground during the snow season is even more complicated than snowfall as it is determined by the influx of snowfall minus the outflow of snowmelt-driven runoff and sublimation.
To approach the problem of determining changes in the cold-season hydroclimate over complex terrain with model simulations, most work has been focused in the regional climate field. Global climate model resolutions until recently have been too coarse, at 200 km or greater, to resolve individual mountain ranges. As a result, studies focused on capturing regional snowfall and snowpack have centered on downscaling global climate model data through statistical and dynamic techniques. Statistical methods like bias correction [some of the most popular methods for hydrologic studies include Wood et al. (2004, 2002)] and constructed analogs (Pierce et al. 2008; Hidalgo et al. 2008) can produce output at finer scales but are limited in regions where there are insufficient observations for comparison or where present climate relationships may rapidly change.
Dynamic downscaling methods (e.g., Ma et al. 2010; Qian et al. 2010; Duffy et al. 2006; Leung and Ghan 1999) using regional climate models can produce output at 1–50-km resolution without observations but are computationally more expensive than statistical methods because of the need to create several model nests to downscale data from global model output to regional scales of interest (Maurer 2007). As a result of computational constraints, time periods of coverage or resolution must be limited. For example, the North American Regional Climate Change Assessment Program, which explores regional climate change over North America using several global climate models as boundary conditions for multiple regional climate models, is limited to 30 yr in the past and future and 50-km resolution (Mearns et al. 2009). Moreover, regional climate models require boundary information from a global climate model, resulting in artificial boundary conditions and uncoupled ocean conditions that introduce biases and complicate connecting regional variability back to the large scale. Ideally, higher-resolution global climate model output would provide self-consistency of a closed global system to analyze large-scale mechanism of snow variability over decades to centuries.
This study presents results for a new high-resolution global climate model in three main parts. 1) For current climate conditions, it provides a comparison of hydroclimate simulations in a high-resolution model with those of a coarser-resolution model and reanalysis. 2) It assesses whether increasing resolution has improved the distributions of snow variables. 3) For idealized changes in radiative forcing, it explores future snow metrics and the controlling mechanisms of cold-season hydrologic change.
We successively work through these problems. In section 2, we describe the two Geophysical Fluid Dynamics Laboratory (GFDL) global climate models and two separate simulations. In section 3, we describe the reanalysis products used to assess the present climate output. Section 4 describes how a specific snow variable is calculated and documents the multivariate regression used to decompose snowfall trends into trends in precipitation and temperature. Sections 5 and 6 describe the snowfall-related climatologies in the present and future climate. Section 7 provides a summary and discussion of our findings.
2. Model simulations
We utilize data from two GFDL global climate models (GCMs) originally described in Delworth et al. (2012). This paper furthers this analysis by focusing on the aspects of these models relating to snow to explore its global variability and provide the foundation for further hydroclimate and water resource studies.
a. CM2.1 model description
We utilize model output from the GFDL Climate Model, version 2.1 (CM2.1), originally developed and widely used as part of the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report (AR4). The model has an approximate horizontal resolution of 200 km in the atmosphere and 100 km in the ocean. Full model documentation is available in Delworth et al. (2006), and output from several experiments was available online at the time of writing from GFDL (http://nomads.gfdl.noaa.gov/CM2.X/).
b. CM2.5 model description
We also use the GFDL Climate Model, version 2.5v1 (CM2.5), the highest-resolution global climate model with similar atmospheric physics, for comparison with CM2.1. The model has an approximate horizontal resolution of 50 km in the atmosphere. The parameterized atmospheric physics are almost identical to those of CM2.1, with some additional cloud-parameter tuning to achieve a radiative balance at finer spatial resolution and more atmospheric vertical levels (24 in CM2.1 vs 32 in CM2.5). The ocean model has a horizontal resolution that varies from 28 km at the equator to 8–11 km at high latitudes. The land model in CM2.5 has been significantly enhanced, including a multilayer model of snowpack and new soil moisture and runoff schemes. A more complete description of the differences in the CM2.1 and CM2.5 model setups is found in the documenting paper (Delworth et al. 2012) and is omitted here for brevity. We will highlight in subsequent sections where these model characteristics are most likely to cause non-resolution-dependent differences in simulated results.
c. Description of simulations
For both CM2.1 and CM2.5, two separate simulations were conducted, referred to as 1) the 1990 control, which is a 280-yr simulation with atmospheric composition (greenhouse gases and aerosols) and external forcing (solar irradiance) fixed at 1990 levels, and 2) the doubled-CO2 climate change experiment, which is a 140-yr simulation that starts with initial and boundary conditions from year 101 of the 1990 control simulation but in which atmospheric carbon dioxide (CO2) increases at a rate of 1% per year until reaching double its initial value after 70 years and is then held constant at the doubled value for the remaining 70 years. The annual-mean characteristics of these simulations were examined in Delworth et al. (2012). Both model versions were found to exhibit some positive model drift in temperature over the duration of their control simulations; care is taken in subsequent sections to express how this is accounted for where appropriate.
3. Analysis products
The lack of continuous global measurements for snowfall and its land surface–dependent variable, snow water equivalent (SWE—the amount of snow on the ground measured by its depth in water if melted), poses a problem for finding a well-vetted product to assess GCM cold-season hydroclimates. This problem leaves two options for assessing snow variables in a GCM: 1) develop a product from available observations of total precipitation or 2) use available reanalysis products. Previous studies have developed proxy snowfall datasets by using temperature data to reclassify total precipitation products into rain or snow (Krasting et al. 2013, hereinafter KBDL; Roesch 2006). We have decided not to use these existing products because they are limited to the Northern Hemisphere, they do not include land surface variables, and they have an assumed relationship among temperature, precipitation, and snowfall. Here, we have chosen to use the European Centre for Medium-Range Weather Forecasts (ECMWF) interim global forecast reanalysis product known as “ERA-Interim” (and subsequently referred to in this paper as ERAI) because of its global coverage and collocated availability of other hydroclimate variables. ERAI is one of the few global reanalysis products with snowfall output and has recently been used by Screen and Simmonds (2012) to assess snowfall variability in the Arctic. We have also chosen a gridded observation-model enhanced SWE product that was specifically developed for GCM assessment in the Northern Hemisphere (Brown and Mote 2009) and a monthly snow-covered-area-extent (SCA) product to analyze seasonality. These products are among the best available off-the-shelf options for assessing global snowfall variability.
ERAI is the latest global atmospheric reanalysis product from ECMWF. It was released in 2009 and covers 1979–present with a horizontal resolution of 75 km. ERAI is an interim product in preparation for replacing the previous 40-yr version of the ECMWF reanalysis product, known as ERA-40 (Dee et al. 2011). The model is improved with new model physics, higher spectral resolution, new data assimilation methods, and—important to our study—a better representation of the hydrologic cycle (Dee and Uppala 2009; Dee et al. 2011). Early papers analyzing the product have shown significant improvements in the global hydrologic cycle when compared with previous reanalysis versions (Dee et al. 2011; Simmons et al. 2010). Given the lack of spatial and temporal coverage for snowfall observations, this product serves as one of the few available proxies for snowfall climatologies, especially in polar regions, where observations are scarce (Screen and Simmonds 2012).
For our comparison of the control model simulations of hydrologic fluxes with reanalysis, we use ERAI snowfall, precipitation, runoff, and SWE. These variables are all archived for daily or monthly means at 0000 and 1200 UTC for 3-, 6-, 9-, or 12-h forecasts. For precipitation, snowfall, and runoff, we created monthly total values by summing the monthly 12-h forecasts for 0000 and 1200 UTC. This method was similarly used in Screen and Simmonds (2012) to explore snowfall variability over the Arctic. For SWE, we created monthly-mean values by averaging the monthly 12-h forecasts for 0000 and 1200. We have documented known issues with ERAI snow variables in the appendix and will reference it throughout this study.
b. Canadian Meteorological Centre SWE analysis data
To complement the ERAI SWE product, we also obtained Northern Hemisphere (excluding Greenland) gridded monthly-mean SWE climatologies from the Canadian Meteorological Centre (CMC; Brown and Brasnett 2010; this database is updated annually). The CMC data were downloaded directly from the Internet (http://nsidc.org/data/docs/daac/nsidc0447_CMC_snow_depth/index.html). They consist of polar-stereographic gridded SWE data at 24-km resolution covering October–June from 1998 to 2010, with 1998 values starting in April. SWE climatologies are estimated from monthly-average snow-depth gridded data in combination with a lookup table of mean monthly snow-density values. The lookup table was originally derived from Canadian snow-course observations corresponding to snow climate classes according to the Sturm et al. (1995) classification. Further explanations of the data product derivation can be found in Brown et al. (2003) and Brasnett (1999). The SWE climatologies have previously been used to analyze spatial patterns in mean annual maximum monthly SWE from 2001 to 2006 from 13 GCMs from phase 3 of the Coupled Model Intercomparison Project (CMIP3; Brown and Mote 2009). We refer to this product as a reanalysis because of the assumptions that are necessary to convert direct snow-depth point measurements into gridded SWE data.
c. Rutgers global snow laboratory product
For understanding the seasonal cycle of Northern Hemisphere snow, we also use the Northern Hemisphere SCA product created by the National Oceanic and Atmospheric Administration (NOAA) and archived at the Global Snow Laboratory of Rutgers, The State University of New Jersey (hereinafter referred to as the Rutgers product). Monthly SCA from November 1966 to the present is available online (http://climate.rutgers.edu/snowcover). This product is produced by integrating available snow stations with interpretations of visible satellite imagery by trained meteorologists (Robinson et al. 1993). These data are somewhat subjective and give pause for using the data for trend analysis; however, statistically significant trends have been detected in previous work (Choi et al. 2010; Liston and Hiemstra 2011). Here we chose simply to compare the Rutgers SCA climatology with that of the control and future-climate simulations.
a. Annual SWE metric
We use SWE as an additional variable for understanding snow biases in the GFDL models. We created our metric by first converting the monthly-mean SWE products into annual product by finding the maximum monthly-mean SWE in a given year at each grid cell. We chose to take the annual maximum to minimize the influences of runoff and sublimation on snowpack loss. This provides the best surface-variable representation of snowfall for comparisons with snowfall analyses and reduces errors caused by early snowmelt onset found in the ERAI product (see the appendix). After converting the SWE record from a monthly product into an annual product, we took the mean of the annual metric (over defined years provided in the analyses) to compare the climatologies of the models and reanalysis products. This metric is referred to as the annual maximum monthly SWE and is used in Fig. 5, which is described in more detail below. A similar calculation was used in Brown and Mote (2009) to compare the CMC product with CMIP3 SWE. This is not a common measurement to be compared with snowfall, but we believe it is appropriate given the availability of additional SWE data products and usefulness of maximum SWE comparisons for hydrological applications.
b. Multivariate regression
Precipitation and temperature are the two main components affecting snowfall. Previous work using station observations has employed simple regression techniques to determine the influence of temperature and precipitation on SWE (Kapnick and Hall 2012, 2010; Mote 2006; Mote et al. 2005). Here, we chose to use a simple multivariate regression to decompose simulated future snowfall S at individual grid cells into components of precipitation P, temperature T, and residuals E,
We employ the simple multivariate regression equation over each grid cell using detrended annual-mean values for each year to retrieve the regression coefficients a and b corresponding to precipitation and temperature, respectively. This process allows us to calculate the contribution of precipitation and temperature to the simulated snowfall trends. We can then calculate the regression-model-derived snowfall using simulated values of precipitation and temperature to compare it with trends in simulated snowfall to quantify the skill of the regression model.
Whereas annual means are used in the control and climate change plots, we use seasonal-mean temperature, precipitation, and snowfall values for the regression calculations. The time frame for calculating seasonal means was determined on a gridcell by gridcell basis by calculating the consecutive 3 months with the greatest snowfall in years 101–200 in the control simulation. This method has been used by Mahlstein et al. (2012) for analyzing future climate changes in precipitation. Observation studies of influences on temperature and precipitation have shown skill in using subannual time scales for attribution calculations (Kapnick and Hall 2012, 2010; Mote 2006; Mote et al. 2005). Subannual regression calculations remove variability in temperature and precipitation during months without snowfall. This is especially important in regions where there can be strong summertime rainfall or temperature fluctuations that are independent of wintertime snowfall.
5. Control climate
a. Topographic characteristics
Topography is plotted in Fig. 1 to highlight the regions of significant orographic variability that are dependent on model resolution. Resolution ranges from approximately 50 km in CM2.5 (Figs. 1a,b) to 75 km in ERAI (Figs. 1e,f) and 200 km in CM2.1 (Figs. 1c,d) in the horizontal plane in the atmosphere and land surface. Of note to snowfall, moving to lower resolutions changes topographic variability in key mountainous regions, leading to the loss of the coastal Scandinavian mountains, a reduction in the European Alps from 2000+ to less than 1000 m, a loss of mountains from the Middle East through Southern Asia, a reduction of western North American mountains to a single sub-3000-m elevation plateau over the Rockies, and a severe truncation of the Andes in the north–south direction. Fifty kilometers is the horizontal resolution necessary to bring these mountain ranges above 1000 m. Changes in topography in the positive tail are the most pronounced: 3.13% of the land surface is above 3500 m in CM2.5 versus only 2.52% in CM2.1. We will show that the topographic differences between models can influence the spatial distribution of elevation-dependent climate characteristics.
b. Time-mean precipitation characteristics
Mean annual total snowfall is compared across the two GFDL models and the ERAI in Fig. 2. Snowfall is given in centimeters of water equivalent, which is the unit tracked in the models and reanalyses. This is not directly comparable to centimeters of depth equivalent (the physical depth of fresh snow when it is on the ground), which requires a snow-density value for conversion into its water equivalent. Grid cells with less than 5 cm of snowfall are shown in white; this threshold was chosen to eliminate locations of light snow dusting of the surface. In the Southern Hemisphere, snowfall values are latitudinally dependent over the ocean between 45° and 60°S, with relative maxima occurring along the coasts. CM2.5 (Fig. 2b) and CM2.1 (Fig. 2d) have relatively greater maxima along coastal Antarctica than in the reanalysis (Fig. 2f) by roughly 20 cm and higher. Station observations have shown a steep gradient in snow accumulation in this region, with more than 110 cm along the coasts, similar to that in CM2.5; past versions of the ECMWF reanalysis have also been observed to have negative snow biases in this region that were attributed to smoothed model topography (Turner et al. 1999). Over the continent of Antarctica, there is a relative maximum of snowfall occurring in the western longitudes (60°–150°W) over the Transantarctic Mountains, with the models again exhibiting greater values of snowfall. This regional maximum is also supported by observed climatologies (Turner et al. 1999; Bromwich 1988). In nonpolar regions in the Southern Hemisphere, the main region with maximum snowfall across the models is over the Andes. Because of differences in topography, the region of snowfall in CM2.1 (Fig. 2d) is smaller and has a lower magnitude (following the topographic mountains seen in Fig. 1d) than in the higher-resolution CM2.5 and reanalysis. Snowfall should occur in a continuous line across the Andes from around the equator down to the tip of South America (Foster et al. 2009). CM2.5 exhibits snowfall across all of the Andes in a continuous line, whereas even ERAI has a break between 15° and 30°S, following the topography.
In the Northern Hemisphere, similar snowfall patterns also emerge, with notable lower values in the reanalysis and CM2.1. Regional maxima are found along coastal regions in the mid- to high latitudes, with the eastern coasts and surrounding oceans exhibiting relatively more snowfall than the western coasts (excluding mountains). Over land, snowfall exceeds 40 cm in high-elevation regions, with the greatest values over southern Greenland, western North America, the Himalaya region, and Scandinavia. Larger values in all mountainous regions are found in CM2.5 when compared with the reanalysis and CM2.1. Maximum snowfall in western North America is pushed inland in CM2.1 (Fig. 2c) when compared with the other climatologies because of the loss of coastal mountains at low resolution (Fig. 1c). Snowfall over the European Alps and the mountains in the Middle East is also notably minimized in CM2.1 because of resolution-dependent topographic smoothing of narrow mountain ranges.
To further explore the sources of differences in global snowfall found in these climatologies, we first look at a map that shows mean annual total precipitation (Fig. 3) only for those regions where annual snowfall exceeds 5 cm; regions with annual snowfall of less than 5 cm have no color shading. Overall, the global precipitation patterns appear to be very similar across the models, with a negative gradient poleward. In South America, similar maxima around 150 cm are found across the Andes in the different climatologies. East–west gradients are seen over North America and Eurasia, with relative maxima on either coast and minima in the continental interior. In mountainous regions with significant snowfall differences in CM2.5 and ERAI (Scandinavia, the European Alps, the Himalaya, and western North America), the corresponding precipitation differences are not as pronounced (Figs. 3a,e). In contrast, coastal and mountainous regional maxima are lower in CM2.1 than in CM2.5 and ERAI. Smoothing in precipitation gradients that is similar to what was shown in snowfall (Fig. 2c) is also seen in CM2.1. Past comparisons between CM2.1 and CM2.5 have attributed this to differences in topography and improvements in the oceanic dynamics in the high-resolution model that lead to reductions in precipitation biases caused by inaccurate coastal sea surface temperatures (Delworth et al. 2012).
In Fig. 4, the ratio of snowfall to total precipitation is shown. Gradients of greater snow as a percentage of total precipitation are seen toward the poles and regionally for high elevations. There appears to be a strong positive difference in the CM2.5 and CM2.1 simulations in the Northern Hemisphere when compared with ERAI. However, the version of the ECMWF model used to produce the ERAI data is known to have an underestimate of snowfall because of the way that the precipitation is partitioned between liquid and solid. The regions with the greatest error caused by this underestimation are found over the western United States, Scandinavia, the Alps, the Middle East, and the Himalaya (see the appendix for a detailed description of this issue). These are the precise locations where there is a negative difference in ERAI when compared with CM2.5 in precipitation partitioning (Fig. 4) and snowfall (Fig. 2). This difference is also seen in the Southern Hemisphere and is manifested by a southward shift in the snowfall region. CM2.5 likely has a positive bias in the proportion of precipitation falling as snow (see next section), but it is not as large as the difference shown in Fig. 4a. This finding is a word of caution when assessing simulated snowfall; snowfall-scheme differences between models can create regional snowfall differences, the sources of which are not always easily identified.
c. Time-mean land surface characteristics
In Fig. 5 we compare the annual maximum monthly SWE climatologies across the two GFDL models, ERAI, and the CMC reanalysis. CM2.1 mainly has a negative difference in SWE (Fig. 5c) when compared with the two reanalysis products. Notably absent are maximum monthly SWE values in the most narrow mountain ranges (the Middle Eastern mountains, the European Alps, Scandinavia, and the west coast of North America), which similarly did not exhibit local maximum snowfall values (Fig. 2c). SWE values above 40 cm are present in the Himalaya region in five grid cells in CM2.1, but without the full geographic span found in CMC. ERAI has a known negative bias in Himalayan snow (Dutra et al. 2011; Tseng et al. 2011; see also the appendix) and only resolves SWE in the highest northwestern mountains (Fig. 2b). Despite these negative differences, however, the spatial patterns of relative maxima in CM2.1 match up with the reanalysis products.
The SWE climatology in CM2.5 more closely resembles that of the reanalyses. In the Himalaya region, CM2.5 has the closest agreement with CMC. In the interior of the continents, CM2.5 and ERAI have greater SWE values than CMC. Maximum values in the Scandinavian mountains, the European Alps, the Middle East, and eastern Canada are also slightly higher in CM2.5. In the western coastal mountains of North America, CM2.5 exhibits a positive difference in SWE when compared with both reanalysis products, with larger SWE magnitudes in the western United States. In some regions, there could be negative differences in SWE in CMC and ERAI as a result of assimilating station observations because of their remoteness as discussed in the appendix. However, there is a positive precipitation bias over the western United States in CM2.5 [shown with careful inspection of Fig. 3, and regionally explored more closely with a separate precipitation product in Delworth et al. (2012)], which likely contributes to the SWE difference. This difference could also arise from land model differences where winter snow does not melt as frequently during the cold season in CM2.5, resulting in greater maximum monthly SWE values than in observations. For the purpose of resolving mountain snowfall for geographic comparisons in the two GFDL models, however, CM2.5 far outperforms CM2.1 because of topographic differences leading to mountain snow improvements in CM2.5.
In addition to looking at annual snow variables, we also calculated the monthly climatologies of Northern Hemisphere SCA to explore seasonality in Fig. 6. SCA appears to have a positive bias from February through October in CM2.5 when compared with the Rutgers product, with the greatest absolute difference in April. This difference is caused by a later peak SCA date in CM2.5 (February and March) than in the Rutgers product (January), which may contribute to this difference as grid cells are still accumulating snow in CM2.5 while melt begins in the Rutgers product. This also mirrors a positive bias in CM2.5 total annual snowfall (Fig. 2) and greater annual maximum monthly SWE (Fig. 5). CM2.5 and the Rutgers product are almost identical during peak Rutgers SCA months (November–January). CM2.1, however, exhibits a greater negative SCA bias, with the maximum absolute bias happening in December and January because CM2.1 does not accumulate snow at a sufficient rate to reach the same magnitude as CM2.5. All of the products reach minimum SCA in August, signaling that the snow season begins across all products at the same time.
We have also calculated the mean annual snowfall to mean annual total runoff ratio in Fig. 7 to compare the modeled influence of snowfall on surface hydrology. This metric was previously used in Barnett et al. (2005) using the Variable Infiltration Capacity model forced at approximately 50-km resolution. The pattern found in Barnett et al. (2005) closely resembles that of CM2.5 (Figs. 7a,b), where maximum values are found in the interior of the northern continents with relative nonzero minima along the coasts, western Europe, and the Midwest into the East Coast of the United States. Where snowfall as a fraction of total precipitation largely diverges between the models and ERAI in the Northern Hemisphere (Fig. 4), the ratio of snowfall to runoff also diverges, supporting the assertion given in the previous section and in the appendix that there is too much liquid precipitation in the ERAI product.
6. Future climate response to 2 × CO2
This section assesses the difference in the climate response found in CM2.5 and CM2.1 under an idealized experiment of increasing greenhouse gas concentrations. As explained in section 2c, atmospheric CO2 in the climate change experiment is increased at a rate of 1% per year until reaching double its initial value after year 70 and is held constant thereafter. Here we examine the transient and postdoubling responses in key snowfall-related variables.
a. Global mean transient response
The global-mean annual-mean near-surface temperature, total precipitation, and snowfall are shown in Fig. 8 for CM2.5 and CM2.1. Values are calculated as area-weighted means over grid cells with 5 cm or more of snowfall in the control climate; these regions correspond to those depicted in Figs. 2a–d. The rate of warming and total warming over snowfall regions is greater in CM2.5 than in CM2.1 (Fig. 8a). As discussed in section 2 and Delworth et al. (2012), there are somewhat different atmosphere–land components in CM2.5 than in CM2.1, and these differences likely contribute to the different rates of warming. Total precipitation in snowfall regions follows a path that is similar to that of temperature, with CM2.5 having a faster rate of increase in precipitation and larger overall precipitation value after CO2 doubling (Fig. 8b). At the end of the simulation, mean global precipitation in CM2.1 rapidly increases; it is unclear if this is due to decadal variability or to a climate system that is continuing to approach equilibrium. More simulation years are needed to further explore these differences.
Global total annual snowfall is shown to decrease in both CM2.5 and CM2.1 and appears to have more rapidly approached equilibrium than did precipitation in years 80–120 but again shows decadal variability at the end of the run (Fig. 8c). In contrast to precipitation gains, CM2.5 more rapidly loses global snowfall and has a greater overall snowfall loss at the end of the simulation. For both models, the magnitude of annual snowfall and precipitation variability is less than the overall climate change response. Future work will examine seasonal-cycle differences in snowfall and precipitation to determine if this signal is robust across all months and seasons and will explore subglobal variability. For example, changes in Northern Hemisphere SCA described in the next section show some differences in regionally averaged snowfall.
b. Postdoubling changes in snowfall characteristics
Here, we examine changes in snowfall characteristics as a result of doubled CO2. These maps differ from the transient response graphs (Fig. 8) in that they depict the change in mean climatologies in years 91–140 in the future simulation versus the corresponding section of the control simulation.
The absolute change in mean annual snowfall in response to doubling CO2 is shown in Fig. 9. Broadly, both CM2.5 and CM2.1 show reductions in snowfall between 30° and 60° latitude, with the greatest relative loses on the western coasts of Europe and North America. The sign of snowfall change flips above approximately 60° latitude, with the largest gains in snowfall along coastal Greenland and Antarctica. The meridional gradient in snowfall change is similar to that of precipitation (shown in Delworth et al. 2012) and corresponds to the often discussed “wet get wetter/dry get drier” pattern that is robustly exhibited across IPCC AR4 climate models (Held and Soden 2006). This snowfall pattern is also shown in the climate models of phase 5 of the Coupled Model Intercomparison Project (CMIP5) by KBDL.
The unique difference in hydrologic change in this study is found over high-elevation regions. In CM2.1, the only increase in mountain snowfall below approximately 60° latitude occurs over the Himalaya region and is weakly positive to neutral (Fig. 9c). In CM2.5, however, there are gains in several locations: gains of 10–25 cm over the northwestern Himalaya region, gains of approximately 5 cm in a small localized region in the southwestern corner of the Yukon province in Canada (at approximately 61°N, 140°W, corresponding to the highest mountain in the region shown in Fig. 1a), and gains of 14 cm in the Chilean Andes (at approximately 27°S, 69°W). Relative to CM2.1, larger snowfall increases are also found over Greenland, eastern Siberia, and coastal Antarctica and greater snowfall losses are shown in the tight western coastal mountain ranges of North and South America, Scandinavia, and the Alps.
Whereas absolute snowfall changes show a change in sign with latitude, changes in precipitation partitioning between snow and rain are more uniform globally (Fig. 10). The only weakly positive gains in snowfall as a percent of total precipitation in CM2.5 occur in two locations with absolute snowfall gains (Figs. 9a,b): eastern Siberia and a few grid cells in the Andes (Figs. 10a,b). There are no positive gains in CM2.1 (Figs. 10c,d). Most regions exhibit a shift from snow to rain, with CM2.5 exhibiting a greater shift in the Arctic and the Southern Ocean. Most of the heating that causes the difference in global temperature response between the two models that is shown in Fig. 8a is confined to these regions (not shown for space; see Fig. 21 of Delworth et al. 2012). It is possible that non-resolution-dependent differences in the two models cause more heating in CM2.5, resulting in a greater shift from snow to rain in these similar regions. We will explore the apparent influence of temperature on snowfall in section 6c. It is striking that, despite the differences in heating between the two models, both have almost universal reductions in snowfall as a percent of total precipitation, even when regional snowfall is increasing.
We calculated the SCA climatology in the double-CO2 simulation to complement the annual changes with a seasonal metric (Fig. 6). We find that the loss in mean global snowfall (Fig. 8c) is indicative of Northern Hemisphere SCA changes: there is a reduction in SCA across all months in both CM2.5 and CM2.1. In CM2.5, Northern Hemisphere SCA loss (signaling snowmelt in more grid cells than accumulation) begins earlier (from February to March) and at a greater rate in the future simulation. The absolute change between the control and future simulation in SCA in both models in the winter (December–February), signaling reductions in snowfall spatially during the cold season, is comparable.
c. A simple decomposition of controlling mechanisms
We have employed the simple regression model described in section 4b to decompose seasonal snowfall trends into additive trends in precipitation and temperature. Figures 11 and 12 present results for CM2.1, and Figs. 13 and 14 present results for CM2.5. Note that, unlike the absolute change in snowfall found in Fig. 9, these regression-model figures do not have a final adjustment for model drift by subtracting mean data from the control and are also for the wettest three consecutive snowfall months. We trained the regression model on simple linearly detrended variables to create a model for seasonal variability. Our results are shown for trends in seasonal snowfall over years 1–80 in the future-climate simulation and are thus not a perfect match to changes shown in Figs. 9 and 10, which are based on the mean annual climatologies of years 91–140 and 191–204 in the future and control climate simulations, respectively.
In CM2.1, precipitation is shown to contribute to increases in snowfall in the high northern latitudes (Fig. 11a) and over the Southern Ocean along the coast of Antarctica (Fig. 12a). Interestingly, despite negative changes in absolute annual snowfall seen over the Himalaya region and Southern Ocean in Figs. 9c and 9d, the regression analysis shows that precipitation in this region contributes to weak positive gains in seasonal snowfall. Temperature almost universally contributes to snowfall losses (Figs. 11b, 12b). Losses are the greatest over the eastern Himalaya region, the North Pacific, western Europe, the Andes, and a portion of the Southern Ocean. The general patterns of snowfall trends found by the regression model (Figs. 11c, 12c) are very similar to that of the actual snowfall trends derived from the simulation (Figs. 11d, 12d). Figures 15c and 15d present the variance explained by the CM2.1 regression model to further quantify the skill of the decomposition in representing snowfall variability. More than 40% of variance is explained by the simple regression model poleward of 30°N and 45°S and with the greatest values over continental interiors, the Himalaya, Greenland, and Antarctica.
In CM2.5, precipitation is similarly shown to contribute to increases in snowfall in the high northern latitudes and over the highest-elevation regions (Figs. 13a, 14a). The magnitudes of trends in snowfall caused by precipitation are greater in CM2.5 than in CM2.1. Remembering the difference in globally averaged precipitation changes in Fig. 8b, CM2.5 exhibits greater increases in global precipitation, which contributes to the magnitude difference seen between Figs. 13a and 11a. Temperature also contributes almost universally to reductions in snowfall, with the greatest trends found in the midlatitudes. Temperature was also shown to increase at a greater rate in CM2.5 than in CM2.1 (Fig. 8a), which contributes to the strength in snowfall trend caused by temperature in Figs. 13 and 14b. The regression model captures magnitudes and spatial distributions of snowfall trends (Figs. 13c, 14c) that are similar to those of the simulation (Figs. 13d, 14d), with an even greater degree of skill than in CM2.1 (Fig. 15; particularly over high elevation regions). In these regions, seasonal-mean temperature and precipitation are good indicators of snowfall variability.
Overall, the regression models appear to perform well in both the higher latitudes and higher altitudes, indicating that mean temperature may be a good indicator of regression-model skill. This is likely due to snowfall being more sensitive to precipitation and temperature fluctuations if climatological-mean seasonal temperatures are at or below freezing than when seasonal temperatures are above freezing. This, in turn, may be due to the following mechanisms: 1) When climatological temperatures are significantly below freezing, temperature fluctuations have little impact and precipitation dominates snowfall variability. 2) When climatological temperatures are close to freezing, temperature and precipitation fluctuations contribute to snowfall variability. 3) When climatological temperatures increase significantly above freezing, precipitation fluctuations will not have an impact on snowfall at all and minor temperature fluctuations will also not influence snowfall as temperatures must depart significantly negatively for snowfall to exist. This phenomenon has been observed over western North America in previous work with monthly observations (Kapnick and Hall 2012). We quantified this phenomenon by regressing variance explained with elevation, mean temperature, and mean precipitation. We found that mean seasonal temperature is the best indicator for grid cells for determining regression-model skill (correlation coefficient r = −0.77 for CM2.5 and r = −0.79 for CM2.1; p < 0.01, Student's t test) whereas elevation (r = 0.44 for CM2.5 and r = 0.48 for CM2.1; p < 0.01, Student's t test) and mean seasonal precipitation do not perform as well (r = −0.43 for CM2.5 and = −0.37 for CM2.1; p < 0.01, Student's t test).
7. Summary and discussion
We have analyzed a low-resolution global climate model and a high-resolution one to assess their ability to simulate snowfall in the present climate and to examine the response of snowfall to idealized radiative forcing changes. We find that global snowfall patterns are very similar in both models, with magnitude differences in snowfall and precipitation found in regions of high elevation. When compared with reanalyses, we find that the high-resolution model has a positive difference in SWE over western North America and central Eurasia, signaling a similar positive difference in snowfall. There are significant improvements in snowfall magnitude and distribution in the high-resolution model. For example, the low-resolution model fails to resolve many mountain ranges in Europe, individual western U.S. mountain ranges, and the full Andes. This new high-resolution model allows us to assess snowfall variability in topographic regions that were previously unresolved in coarser global climate models and represents a significant advance in our ability to simulate regional-scale snowfall and its response to climate change. This is an important step toward providing better water-resource applications.
Changes in snowfall climatologies are explored in both models using an idealized double-CO2 experiment. Figure 10 reflects simulated temperature trends: as temperatures increase, less precipitation tends to fall as snow (Figs. 11–14 all confirm this). The strongest snowfall losses are exhibited in regions that are characterized both by reduced precipitation and increased temperature in the low to midlatitudes in both models (Figs. 11–14). Interestingly, there are sign differences in changes in absolute snowfall in the high- and low-resolution models over the highest-elevation peaks (Fig. 9). The northwestern Himalaya region, a small location in the Andes, and the highest peak in the southwestern corner of the Yukon all exhibit gains in absolute snowfall in CM2.5 and losses in CM2.1. Trends in total precipitation dominate this difference, with greater precipitation enhancement in CM2.5 than in CM2.1. The mountainous regions experiencing enhanced snowfall are unique to CM2.5 as they are high enough to capture moisture and be less susceptible to rising temperatures because of naturally colder climatologies. Even as nearby mountains at lower elevations experience gains in total precipitation, snowfall is reduced because of rising temperature.
The cryosphere observational record shows positive present climate trends in similar isolated regions of increasing snowfall that were found in the idealized double-CO2 experiment. Over the Himalaya, the overall region shows little positive or no trend in recorded precipitation (Shrestha et al. 2000), whereas the Karakoram region in the northwest has exhibited gains in winter precipitation and snowfall (Shekhar et al. 2010; Fowler and Archer 2006). Similarly, observations since the mid-1800s indicate that most Himalayan glaciers are experiencing mass loss, whereas the Karakoram region shows stability or expansion (Bolch et al. 2012). Snowfall is one input variable that contributes to changes in glacier mass balance; future work should further explore how snowfall will change in these regions and its contribution to variability in glacial mass. Positive and accelerating snow accumulation trends have also been observed since the mid-1800s at Mount Logan, the highest peak in northwestern Canada (Moore et al. 2002), which is roughly resolved in CM2.5. Observations have also shown increases in mean snow depths (an indication of snowfall gains) in eastern Siberia (Bulygina et al. 2009). Together, these records provide physical evidence that, in the present warming climate, snowfall is increasing at some of the same unique locations found in CM2.5, implying that the model has captured a continuation of present climate trends. Further simulations at high resolutions are needed to study past variations in snowfall to confirm this result.
Resolutions on the order of CM2.5 (50-km atmospheric resolution) or greater are required to begin to properly resolve mountains to simulate positive trends in snowfall in complex topographic regions in the future climate and resembles observed snowfall trends in recent history. In studies using climate projections produced by GCMs that are similar to CM2.1 (atmospheric resolution of 200 km or more), gains in annual snowfall in these unique regions are not exhibited (KBDL; Brown and Mote 2009). The change in sign of future snowfall between the two models is an important finding. Studies relying on coarse-resolution global models (such as the CMIP5 archive) are presently unable to capture such changes in annual snowfall in key high mountain regions.
The present modeling study does not project positive trends in snowfall in the narrowest western U.S. mountain ranges known for SWE gains. The southern California Sierra Nevada and the southern Rocky Mountains have exhibited increases in SWE in April (Mote 2006; Mote et al. 2005) in the last several decades that have been mainly attributed to increases in snowfall in February (Kapnick and Hall 2012). This regional signal has not been found in past GCM (Pierce et al. 2008) or hydrologic model work (Mote et al. 2005) either. This fact may be due to natural internal climate variability of the models, requiring multiple ensembles to capture such phenomena (it may be a rare event), or may require further resolution enhancement or dynamical upgrades to properly capture the mechanisms relating to western U.S. snowfall. In this study, the jump from 200-km resolution in CM2.1 to 50-km resolution in CM2.5 noticeably increases the model's ability to represent the variability and height of key mountain ranges, resulting in improved spatial distribution of snowfall. As resolution increases in global models beyond 50 km, future downscaling experiments may be able to reproduce observed snowfall trends in extremely narrow mountain ranges because of improved boundary conditions. The higher-resolution GCMs may even be able to reproduce more numerous regions of increased snowfall than are presently produced in this study's idealized climate change experiments. Subseasonal analysis of simulated snowfall trends may also reveal that, over small regions, cold winter months are experiencing snowfall gains while warming springs have led to reduced late-season snowfall as seen in observations (Kapnick and Hall 2012). On a hemispheric-mean basis, however, this is not shown in our calculation of Northern Hemisphere SCA (Fig. 6). This should be noted when interpreting GCM snowfall signals directly; resolution, averaging periods, and physics enhancements may lead to a reversal in projected signals, particularly in high-elevation regions.
This work also highlights the need for further examination of the model schemes used to represent snowfall processes. Without proper dynamics, the distribution of snowfall can be greatly altered. Using the ERAI data as an example, when improvements were made to the mixed-phase cloud scheme in the ECMWF atmospheric model to make it more realistic, snowfall climatologies were improved (discussed in the appendix and seen in Figs. 2–4). The differing treatment of snowfall in global atmospheric models can result in subtle differences in precipitation partitioning, resulting in noticeable differences in snowfall and SWE, the source of which is not easily identified without examining snowfall schemes.
Seasonal snowfall plays an important role in high-elevation and mid- to high-latitude hydrologic budgets. Barnett et al. (2005) calculated the ratio of snowfall to runoff globally using a hydrologic model to identify the regions at risk for runoff timing changes and resulting water supply shocks in response to shifts from snow to rain. CM2.5 and, to a lesser extent, CM2.1 are able to reproduce the results seen in Barnett et al. (2005) in the present climate (Fig. 7). The double-CO2 experiment shows a very different distribution of this ratio (Fig. 16). The shift in precipitation type from snow to rain results in a shrinking of the regions where snowfall dominates runoff with Europe, northwestern Russia, and northern Canada being particularly affected (overlay with circles and a purple border shows the present climate 0.9+ region). CM2.5 shows smaller regions with significant reductions in this ratio. Both also show subtle reductions in the ratio along the coasts in regions where the present climate ratio was lower. In the Southern Hemisphere, the plots focus on South America to show that there are significant changes in the hydrologic budget in narrow coastal mountain bands that are very different between the two models. There is also a uniform reduction in SCA across all months (Fig. 6), implying a significant shift in the Northern Hemisphere seasonal cycle. This may also affect regional warming trends through albedo feedback processes relating to the loss of spring snow cover (Hall 2004). More work is needed to explore the seasonality of SWE, snow cover, and runoff changes to understand the extent of the impact of reduced snowfall.
This work has two model-related limitations. The first limitation is a result of our using only two idealized simulations in two separate global climate models to explore snowfall variability. While the latitudinal gains and losses in total precipitation are robust across IPCC AR4 climate models (Held and Soden 2006), there are significant regional variations, especially over midlatitude mountain ranges (Kharin et al. 2007). Given that the sign of snowfall changes depends on precipitation changes, a larger pool of models and ensembles at high resolution may result in a wider range in snowfall changes. For example, one study by KBDL using the CMIP5 archive has shown that models had significant spread in projected future snowfall over northeastern Greenland. The second limitation is caused by the idealized future-climate experiment being truncated after only 70 yr of CO2 reaching equilibrium. Given that the deep ocean will continue to equilibrate, it is unclear whether the future simulation has a sufficient time scale to fully capture the response of snowfall to a double-CO2 environment.
The main objective of this study was to explore cold-season hydroclimate in two similar GFDL models with different spatial resolutions and to assess future changes; future work will expand on these findings at different scales. Future studies will focus on the unique regions where resolution improvements increased model skill in snow variables more closely by examining the seasonal variability of snowfall and the role of atmospheric circulation. Except for validation purposes in the present climate and a brief discussion of runoff and SCA implications, this study has also not explored the surface characteristics of snowpack. Given the improvements in the land surface model in CM2.5 and various impacts of snowfall variability, future work will also explore the role of the land surface in driving regional snowpack dynamics and water budgets.
The authors are grateful to Sarah Keeley for help in acquiring the ERAI data and in discussing its use. The authors also thank Richard Forbes and Patricia Rosnay for helpful discussions of the ECMWF model components and snowfall schemes. The writers also acknowledge Alexis Berg and Sergey Malyshev for helpful comments and discussions of the manuscript and acknowledge the improvements made by three anonymous reviewers.
Known ERAI Issues Affecting Snow Variables
We have compiled known sources of ERAI snow variable biases relative to observations. These issues help to explain some of the difference between the snow climatologies of the reanalysis products and control simulations.
a. Precipitation partitioning between liquid and solid states
The version of the ECMWF model used to produce the ERA-Interim data is known to underestimate snowfall as a result of the way the precipitation is partitioned between liquid and solid in the mixed-phase temperature zone between 0° and −23°C (Dutra et al. 2010, 2011). In the presented model version, the phase of cloud condensate is diagnosed as a function of temperature, making a transition from all liquid at 0°C to all ice at −23°C. In this temperature range, precipitation generation is also partitioned between liquid and solid in the same way so that an increasing proportion of the precipitation falls as rain as the temperature increases from −23° to 0°C. This assumption is not similarly valid for stratiform clouds; later model versions were modified so that all precipitation generated at temperatures colder than 0°C was in the form of snow. This change leads to an increase in snowfall wherever there is precipitation forming from clouds in the mixed-phase transition zone. The largest snowfall increases due to the mixed-phase precipitation changes are found over the western United States, Scandinavia, the Alps, the Middle East, and the Himalaya (Dutra et al. 2011).
b. Snow water equivalent
The ERAI SWE product is a combination of ECMWF output corrected with available observations. For example, since July 2003, the NOAA/National Environmental Satellite, Data, and Information Service (NESDIS) daily Interactive Multisensor Snow and Ice Mapping System snow cover product has been used to correct grid cells in ERAI (Dee et al. 2011). Grid cells with missing snow cover in ERAI when compared with the NOAA/NESDIS product have 10 cm of SWE added to the background field. At locations that are snow covered according to the ERAI but that are snow free in the NOAA/NESDIS product, a pseudo observation of 0 cm SWE is presented to the reanalysis. This method was adopted after analysis showed that the total area of grid boxes with snow cover was approximately 10% larger in the earlier ECMWF model output than in the NOAA/NESDIS snow-extent product (Drusch et al. 2004).
While the addition of the NOAA/NESDIS product helped to improve SWE in ERAI, it also introduced a small error. Dutra et al. (2010) and Dee et al. (2011) found that SWE has spurious snow-free zones along the coastal grid boxes in some Northern Hemisphere regions. This was created by a geolocation error in which snow was inaccurately mapped from the NOAA/NESDIS product over mismarked ocean instead of land grid cells. As a result, there is a slight negative bias in SWE in some coastal regions where snow was removed.
Point surface snow-depth data are also assimilated into the SWE product using a Cressman spatial interpolation as described in Drusch et al. (2004). It is important to note that point-source snow stations are not evenly distributed and are often more prevalent in regions easily accessible to humans, resulting in a bias of stations located in mid- to low elevations and midlatitudes, resulting in snow-free or low-snow values in the assimilated remote locations (Pan et al. 2003). Potentially because of the lack of observations, snow cover has been shown to be systematically underrepresented in the ERAI product over Scandinavia, western Russia, and central/eastern Canada by up to 1 month during the spring (i.e., snow melts too early in ERAI; Dutra et al. 2010). As a result of these errors, we chose to analyze the annual maximum monthly SWE at each grid cell to minimize the error of missing snow cover during the late-spring melt season seen in Dutra et al. (2010). This calculation is described in section 4a.