Abstract

Statistically downscaled climate projections from nine global climate models (GCMs) are used to force a snow accumulation and ablation model (SNOW-17) across the central-eastern North American Landscape Conservation Cooperatives (LCCs) to develop high-resolution projections of snowfall, snow depth, and winter severity index (WSI) by the middle and late twenty-first century. Here, projections of a cumulative WSI (CWSI) known to influence autumn–winter waterfowl migration are used to demonstrate the utility of SNOW-17 results. The application of statistically downscaled climate data and a snow model leads to a better representation of lake processes in the Great Lakes basin, topographic effects in the Appalachian Mountains, and spatial patterns of climatological snowfall, compared to the original GCMs. Annual mean snowfall is simulated to decline across the region, particularly in early winter (December–January), leading to a delay in the mean onset of the snow season. Because of a warming-induced acceleration of snowmelt, the percentage loss in snow depth exceeds that of snowfall. Across the Plains and Prairie Potholes LCC and the Upper Midwest and Great Lakes LCC, daily snowfall events are projected to become less common but more intense. The greatest reductions in the number of days per year with a present snowpack are expected close to the historical position of the −5°C isotherm in December–March, around 44°N. The CWSI is projected to decline substantially during December–January, leading to increased likelihood of delays in timing and intensity of autumn–winter waterfowl migrations.

1. Introduction

Snow is an important element of the hydrologic cycle, replenishing soil moisture and contributing runoff to river basins through spring melt (Lettenmaier and Gan 1990; Groisman et al. 2001). Snow cover enhances longwave emissivity and induces an insulating effect on soil temperatures (Zhang et al. 2008). Continued increasing temperatures may lead to a diminished snowpack, which would shift snowmelt from spring to winter (Manabe et al. 1981; Mahanama et al. 2012), reduce soil moisture in spring–summer during low-flow season (Mastin et al. 2011), favor drought development (Mishra et al. 2010; Mahanama et al. 2012), and increase fire risk (Westerling et al. 2006). Furthermore, snow influences many facets of society. Spring snowpack is an essential water-supply reservoir for irrigation, hydroelectric power generation, and municipal demands (Norton and Bolsenga 1993; Kunkel et al. 2000; Christensen et al. 2007; Barnett et al. 2005; Vicuna et al. 2008). Snowfall supports multimillion-dollar winter recreational industries and tourism (Scott et al. 2008). Heavy snowstorms can be a natural hazard to life and property, leading to automobile accidents (Eisenberg and Warner 2005), heart attacks (Rogot and Padgett 1976), disrupted air and surface traffic (Robinson 1989; Rasmussen et al. 1993), and large economic costs (Rooney 1967). The survival and distribution of many animals are strongly impacted by the presence and depth of snow (Severinghaus 1947; Edwards 1956; Verme 1968; Krohn et al. 2006; Hoving et al. 2005; Carroll 2007; Notaro et al. 2011).

North American snow cover has declined since the mid-twentieth century, especially in spring (Brown 2000; Lemke et al. 2007). Snow cover duration across the Northern Hemisphere decreased during 1966–2007 across a zone where the seasonal mean air temperatures were within ±5°C (Brown and Mote 2009). In contrast, a diminishing lake ice cover and more favorable circulation patterns have enhanced lake-effect snowfall across the Great Lakes basin (Norton and Bolsenga 1993; Leathers and Ellis 1996; Burnett et al. 2003). If no change is expected in precipitation, then future warming should favor more rain and less snow. However, given projections of greater cold season precipitation, particularly at mid-to-high latitudes, an increase in snowfall might be expected where temperature continues to remain below freezing (Kapnick and Delworth 2013). The response in snow cover to enhanced greenhouse warming is therefore complicated by projected increases in mid- to high-latitude precipitation, with the response varying by latitude and elevation (Groisman et al. 1993; Räisänen 2008; Brown and Mote 2009). Climate models generally display a transition zone, around the present-day −10°C isotherm, with expected snowfall increases to the north (high latitudes) and decreases to the south (midlatitudes) (Krasting et al. 2013). Projections of snow depth are more complicated than of snowfall, as the former is determined by both snowfall and snowmelt-driven runoff and sublimation (Kapnick and Delworth 2013).

Several modeling studies have generated coarse projections of twenty-first-century snowfall. Long-term planning for water management strategies and reservoir design is driven by the potential for greater winter flows and reduced spring snowpack and spring–summer runoff (Adeloye et al. 1999; Draper and Kundell 2007; Mastin et al. 2011). Simulations by the National Center for Atmospheric Research (NCAR) Community Climate System Model suggest that snow cover will decline across all of North America, except northernmost Canada, during the twenty-first century, particularly during springtime, despite projected increases in cold season precipitation (Peacock 2012). Using the Geophysical Fluid Dynamics Laboratory Climate Model version 2.1 (GFDL CM2.1), with approximate 50-km spatial resolution, Kapnick and Delworth (2013) simulated positive trends in twenty-first century snowfall across the mid-to-high latitudes and negative trends across the mid-to-low latitudes, including reductions across all of North America, except Hudson Bay and northern Quebec. According to an ensemble of coupled atmosphere–ocean global climate models (GCMs) from phase 5 of the the Coupled Model Intercomparison Project (CMIP5), annual snowfall is anticipated to decline during the twenty-first century across much of the Northern Hemisphere, except with high-latitude increases (Brutel-Vuilmet et al. 2013; Krasting et al. 2013). The CMIP3 and CMIP5 GCMs share similar deficiencies in their capacity to simulate mean patterns and trends in Northern Hemispheric snow cover, including an underestimation of the recent negative trend in spring snow cover extent (Roesch 2006; Brutel-Vuilmet et al. 2013).

The purpose of this study is to generate high-resolution projections of various winter snow metrics for the mid and late twenty-first century across the central-eastern North American Landscape Conservation Cooperatives (LCCs) using an operational snow model forced by statistically downscaled CMIP3 climate projections. The LCCs are an applied conservation science partnership among states, federal agencies, universities, nongovernmental organizations, and tribes that provides the expertise needed to support landscape-scale conservation planning. These winter projections will benefit state, provincial, and federal natural resource agencies and nongovernmental agencies in generating survival, distribution, population, and carrying capacity projections for wildlife influenced by winter severity. To date, a lack of quality twenty-first-century snow projections has hindered the capacity to estimate the potential effects of climate change on species that are sensitive to snowpack and temperature.

2. Data and methods

a. Statistical downscaling

Through the Wisconsin Initiative on Climate Change Impacts (WICCI), daily statistically downscaled climate projections were developed to guide climate change impact and adaptation studies (WICCI 2011). Notaro et al. (2011) used this downscaled data to force an operational snow model, SNOW-17, and generate high-resolution snow projections for Wisconsin. The current study uses an improved version of the downscaling covering 10 of the LCCs across the United States and southern Canada, to the east of the Rocky Mountains (Tables 1 and 2), at 0.1° × 0.1° resolution. Downscaled variables include daily maximum temperature, minimum temperature, and precipitation for the late twentieth (1961–2000), mid-twenty-first (2046–65), and late twenty-first centuries (2081–2100), based on three CMIP3 emission scenarios (A2, A1B, and B1) and 13 GCMs. The present study focuses on the nine GCMs (see Table 3, including expansions of the nine model names) for which downscaled projections are available for the late twentieth century (20C3M) and each emission scenario, with results shown for the contrasting high-end (A2) and low-end (B1) scenarios; one ensemble member is considered per model. By the year 2100, atmospheric carbon dioxide levels are projected to reach 840 and 550 ppm, under the respective scenarios.

Table 1.

Summary of results for five LCC regions: South Atlantic, Peninsular Florida, Appalachian, Gulf Coastal Plains and Ozarks, and Eastern Tallgrass Prairie and Big Rivers. Variables include DJFM air temperature, DJFM precipitation, annual snowfall, November–April (NDJFMA) mean snow depth, number of days per year with snowfall of at least 1 cm, average snowfall per event, number of days per year with a snowpack of at least 1 cm, and cumulative WSI. For the first seven variables, Mod represents the late twentieth-century mean value and A2 and B1 represent projected changes (both actual change and percentage change in parentheses) by the end of the century (2081–2100 vs 1981–2000). For cumulative WSI, Mod represents the mean WSI value for 1981–2000, while A2 and B1 represent the mean WSI value for 2081–2100 for two emission scenarios.

Summary of results for five LCC regions: South Atlantic, Peninsular Florida, Appalachian, Gulf Coastal Plains and Ozarks, and Eastern Tallgrass Prairie and Big Rivers. Variables include DJFM air temperature, DJFM precipitation, annual snowfall, November–April (NDJFMA) mean snow depth, number of days per year with snowfall of at least 1 cm, average snowfall per event, number of days per year with a snowpack of at least 1 cm, and cumulative WSI. For the first seven variables, Mod represents the late twentieth-century mean value and A2 and B1 represent projected changes (both actual change and percentage change in parentheses) by the end of the century (2081–2100 vs 1981–2000). For cumulative WSI, Mod represents the mean WSI value for 1981–2000, while A2 and B1 represent the mean WSI value for 2081–2100 for two emission scenarios.
Summary of results for five LCC regions: South Atlantic, Peninsular Florida, Appalachian, Gulf Coastal Plains and Ozarks, and Eastern Tallgrass Prairie and Big Rivers. Variables include DJFM air temperature, DJFM precipitation, annual snowfall, November–April (NDJFMA) mean snow depth, number of days per year with snowfall of at least 1 cm, average snowfall per event, number of days per year with a snowpack of at least 1 cm, and cumulative WSI. For the first seven variables, Mod represents the late twentieth-century mean value and A2 and B1 represent projected changes (both actual change and percentage change in parentheses) by the end of the century (2081–2100 vs 1981–2000). For cumulative WSI, Mod represents the mean WSI value for 1981–2000, while A2 and B1 represent the mean WSI value for 2081–2100 for two emission scenarios.
Table 2.

As in Table 1, but for five other LCC regions: Plains and Prairie Potholes, Great Plains, Gulf Coast Prairie, Upper Midwest and Great Lakes, and North Atlantic.

As in Table 1, but for five other LCC regions: Plains and Prairie Potholes, Great Plains, Gulf Coast Prairie, Upper Midwest and Great Lakes, and North Atlantic.
As in Table 1, but for five other LCC regions: Plains and Prairie Potholes, Great Plains, Gulf Coast Prairie, Upper Midwest and Great Lakes, and North Atlantic.
Table 3.

List of the nine global climate models analyzed in this study.

List of the nine global climate models analyzed in this study.
List of the nine global climate models analyzed in this study.

The station data used to train the statistical models consist of daily maximum and minimum temperature and precipitation from the National Weather Service (NWS)’s Cooperative Observer Program (1950–2009) and Environment Canada’s Canadian Daily Climate Data (1950–2007), for roughly 4000 stations in the United States and Canada. Objective methods are applied to correct the hour of observation and remove stations with large biases in wet day frequency (Daly et al. 2007). Only stations with at least 30 yr of data, and fewer than five missing days per month, are included. The large-scale atmospheric state is obtained from the National Centers for Environmental Prediction (NCEP)–NCAR reanalysis (Kalnay et al. 1996). The downscaling methodology is probabilistic, as a given large-scale state does not determine a single, precise value for a location’s temperature or precipitation, but rather the parameters’ values in a parametric probability density function (PDF). These parameters vary daily with the large-scale fields. A normal distribution is used for maximum and minimum temperature, a generalized gamma distribution (Stacy 1962) is used for precipitation amount on wet days, and the probability of measurable precipitation is modeled with logistic regression. A realistic temperature skewness is simulated, because the final temperature distribution is determined by both the large-scale predictors with normally distributed noise added to the mean. Unlike classic linear regression, the variance is allowed to change, based on large-scale fields.

The reanalysis predictors for daily maximum (minimum) temperature include the 2-m maximum (minimum) temperature itself, 1000-m wind vector, and column-integrated relative humidity. Here, column-integrated relative humidity measures the lower to midtropospheric water holding capacity, computed as the ratio of the sum of specific humidity between sigma levels 0.9 and 0.5 and the sum of the saturated specific humidity between those levels. For the probability of measurable precipitation, the predictors are precipitation rate to the 0.25 power, column-integrated relative humidity, average wind vector between sigma levels 0.9 and 0.5, and a pseudo-K index. The inclusion of wind as a predictor allows for the representation of lake-effect precipitation. The K index is defined as

 
formula

where T is temperature and TD is dewpoint, both in degrees Celsius. Given that this equation contains an imbalance of three temperature values and two dewpoint values, a projected increase in temperature would almost certainly imply increased instability, whether realistic or not. If relative humidity is constant, which is a reasonable first-order approximation for global warming (Held and Soden 2000), then both the dewpoint temperature and actual temperature increase by the same amount. Therefore, the pseudo-K index is computed by eliminating the T850hPa term. To make the index suitable for high elevations, sigma levels are used. The pseudo-K index is defined as

 
formula

The predictors for precipitation amount include precipitation itself, the pseudo-K index, available water, and available water flux. The latter two predictors are motivated by the scaling relationships between moisture and precipitation under climate change (O’Gorman and Schneider 2009). Available water is defined as the amount of water that is condensed out when an air parcel lifts vertically from its current pressure to half its pressure. The available water and available water fluxes (computed using the wind field) are averaged vertically between sigma levels 0.9 and 0.5 before serving as predictors.

The above statistical models determine the spatially and temporally varying PDFs of maximum and minimum temperature and precipitation at each station and day in the observed record. To apply the downscaling to climate models, the cumulative distribution function remapping algorithm of Wood et al. (2004) is used to debias the GCM predictors to match the reanalysis. The statistical models are then applied to the debiased model predictors, creating downscaled PDFs for each station and day in the climate models.

To create a gridded downscaling dataset, the PDF parameters are interpolated to a predefined grid. In the study by Notaro et al. (2011), the PDFs were interpolated to a grid based on distance to surrounding stations. The mean temperature at a grid point g was given by

 
formula

where the summation is across stations, d(g, sj) is the distance between the grid point g and station sj, and D is a parameter that is fit by leaving each station out in turn and using (3) to predict the station’s temperature from the remaining stations. The D that minimizes the error over all stations was used. A different value of D was found for each variable and month. In the new downscaling, (3) is modified so that the weighting also depends on the difference in the mean climate between the stations sj and the grid point g:

 
formula

where ΔT(g, sj) is the mean temperature difference between the grid point and the station. The mean climate is obtained from the Parameter-Elevation Regressions on Independent Slopes Model (PRISM; Daly 2006; Daly et al. 2008) and the Canadian Climate Normals dataset (McKenney et al. 2001). These datasets are interpolated to the downscaling grid and to the stations using bilinear interpolation. Using PRISM and the Canadian Climate Normals dataset adds information because the fitting of the conditional PDF estimates requires stations with extensive data, while the PRISM and Canadian Climate Normals data are intended only to estimate the mean. Therefore, many more stations are used to constrain the mean climate in the PRISM and Canadian Climate Normals datasets. PRISM also estimates the role of topography and proximity to coastlines in its interpolation scheme. In (4), the mean climate effect is intentionally of lower order than the distance term, so that the distance term dominates if one is far removed from the grid point. Stations with a similar climate are given extra weight only if they are relatively close to the grid point. The constants D1 and D2 are found by cross validating over stations.

In regions of complex topography and data-sparse regions of Canada, the downscaling is relatively unconstrained by data, as the applied approach is restricted to stations with long records. Therefore, PRISM and the Canadian Climate Normals dataset are also used to adjust the mean maximum and minimum temperature and precipitation. For temperature, an additive correction is used, which is the difference between the PRISM and/or Canadian Climate Normals dataset and gridded downscaled data. For precipitation, a multiplicative correction is used. This correction is not meant to correct the downscaled mean of the station data. For example, the conditional mean for temperature is determined by ordinary least squares regression, so the mean is exact to numerical precision. Instead, this correction is used to provide additional information in between stations in regions of complex topography and in data-sparse regions, as confirmed by the fact that the corrections are small, except in mountains and parts of Canada.

Spatially and temporally varying PDFs of maximum and minimum temperature and precipitation are now known at each grid point and day. To generate normal data, random numbers are drawn from the PDFs, generating a possible realization of the small-scale state that is consistent with the large-scale fields. The random numbers are correlated in space (and time, in the case of temperature) so that the spatial and temporal correlations of the downscaled variables are similar to observations. One of the randomly generated realizations is used to drive the snow model.

Wintertime variance and covariance of temperature and precipitation from observations and downscaled data are used to gauge the fidelity of the downscaling. The downscaling, by design, reproduces the mean climate, so this analysis does not focus on the mean. All variables are averaged across the nine climate models before plotting the results. For the standard deviation in daily maximum temperature (see Figs. S1a,b in the supplementary material), the downscaling reproduces the strong variability from the Great Plains to the Ohio Valley and weak variability near the Great Lakes and coasts. The largest errors are found in Florida, where the downscaled maximum temperature variability is too large. For the standard deviation in daily minimum temperature (Figs. S1c,d), the downscaling reproduces the strong variability in the north and relatively weak variability near the Great Lakes, Atlantic coast, and southwestern Great Plains.

Given that projected snowfall extremes are examined in this study, precipitation extremes in the observations and downscaling data are compared (Fig. S2 in the supplementary material). The downscaled extremes are calculated directly from the PDFs and thus exhibit less sampling variability (spatial noise) than the observations. The downscaling reproduces the spatial patterns and magnitudes of the observed 99.9th and 99.99th percentiles. The largest differences are found across the southeastern United States for the 99.99th percentile, although given the limited sample size, much of the differences likely results from sampling variability.

The correlation between daily wintertime temperature and precipitation is assessed (Fig. S3 in the supplementary material). Unlike the precipitation and temperature variability, no attempt was made to fit the correlation between temperature and precipitation. The correlations in the downscaling are those that naturally appear as a result of covariability in the large-scale predictors used to individually predict temperature and precipitation. For maximum temperature and precipitation, the downscaling reproduces the general pattern of positive correlations in the Northeast and negative correlations in the western Great Plains. The downscaling also shows decreased correlations to the lee of the Appalachians, although the downscaled correlations are too positive. The downscaled correlations are typically too large across the Northeast. For minimum temperature and precipitation, the downscaling reproduces the general large-scale patterns, although the downscaled correlations are too large across the northwestern Great Plains and around Tennessee.

b. Cumulative winter severity index

Schummer et al. (2010) developed a cumulative winter severity index (CWSI) to explain changes in the relative abundance of dabbling ducks during autumn–winter at midlatitude migration areas of eastern North America. CWSI is defined as

 
formula

where TEMP is the mean daily temperature (degrees Celsius multiplied by −1), TEMPDAYS is the number of consecutive days with a mean temperature less than 0°C, SNOW is the snow depth (converted to inches dividing by 2.54 cm in.−1), and SNOWDAYS is the number of consecutive days with a snow depth of at least 2.54 cm. CWSI is cumulatively summed each day from 1 September to 31 March. For temperate regions, CWSI remains below zero throughout the warm season, but as the air temperature consistently drops below freezing and a snowpack develops from late autumn to early winter, the index accrues and achieves positive values. CWSI reflects the current and cumulative effects of temperature on ducks’ energy expenditure and snow and ice cover on food availability to meet the ducks’ energy needs. Based on the CWSI threshold of 7.2 for mallards, researchers can identify the likely date when duck abundance will begin to decrease at northern latitudes with southward migration (Schummer et al. 2010).

c. Snow model

SNOW-17 is an empirically based, operational snow accumulation and ablation model (Anderson 1973, 2002, 2006; Franz et al. 2008a,b; Raleigh and Lundquist 2012), a component of the NWS River Forecast System. SNOW-17 treats the key physical processes that regulate snow dynamics in a conceptual manner, requiring only commonly available data observations of temperature and precipitation as inputs (Fig. S4 in the supplementary material). In treating a column of snowpack, the model addresses the following principal processes: form of precipitation, snow accumulation, energy exchange at the snow–air interface, heat exchange at the soil–snow interface, heat storage and deficit within the snowpack, and liquid water retention and transmission of water through the snow cover (Anderson 2006). Snow cover is modeled as a single bulk layer with a specified water holding capacity. The density of new snow is treated as a function of air temperature, while the density of existing snow responds to compaction, destructive metamorphism, and the component of melt metamorphism due to the presence of liquid water (Anderson 2006). The simple temperature-based SNOW-17 model performs at least as well as more physically based energy balance models (Anderson 1973, 1976; Ohmura 2001; Zappa et al. 2003; Franz et al. 2008a). There is a precedent to applying SNOW-17 in climate change studies (Notaro et al. 2011). Prior studies have used either SNOW-17 or the Sacramento Soil Moisture Accounting Model (Burnash et al. 1973) coupled to SNOW-17, forced by climate projection data, to investigate potential changes in California water resources and reservoir operations (Miller et al. 2003; Brekke et al. 2009; Maurer et al. 2010; Hanson et al. 2012).

d. Datasets

Hourly air temperature data at 12 weather stations across the United States for 1980–99 are used to compute their mean diurnal cycle for December–February (DJF) and identify the times of maximum and minimum temperature. The station data are extracted from the National Climatic Data Center (NCDC) global and U.S. integrated surface hourly dataset (available at http://www.ncdc.noaa.gov/cdo-web) for Billings, Montana; Havre, Montana; Pueblo, Colorado; and Cheyenne, Wyoming, in the western study region (103°–116°W); for LaCrosse, Wisconsin; Dubuque, Iowa; Madison, Wisconsin; and Oklahoma City, Oklahoma, in the central region (87°–103°W); and for Bangor, Maine; Atlanta, Georgia; Greenville, South Carolina; and Teterboro, New Jersey, in the eastern region (67°–87°W). Based on analysis of these stations, the daily minimum and maximum temperatures in DJF typically occur at 1400 UTC (range of 1400–1500 UTC) and 2200 UTC (2000–2200 UTC) in the western region, 1300 UTC (1100–1400 UTC) and 2100 UTC (2000–2100 UTC) in the central region, and 1200 UTC (1000–1200 UTC) and 2000 UTC (1900–2100 UTC) in the eastern region, respectively. This information aids the interpolation of downscaled daily maximum and minimum temperatures to hourly values, for input to SNOW-17. Based on hourly present weather (rain or snow occurrence) and air temperature for 1980–99 at 135 stations across the study region, from the NCDC global and U.S. integrated surface hourly dataset, the mean snow–rain temperature threshold is computed at different locations and assigned to the SNOW-17 parameter, PXTEMP. Climatologies of snow depth and surface soil temperature for DJF from the North American Regional Reanalysis (NARR; Mesinger et al. 2006) are used to estimate SNOW-17’s antecedent temperature index parameter (TIPM) and the parameter for daily snowmelt at the snow–soil interface (DAYGM), respectively. The SI parameter that quantifies the required amount of liquid snow for 100% cover is determined by applying values used by the Noah land surface model to an International Geosphere–Biosphere Programme land cover map. Detailed descriptions of SNOW-17 parameters are provided in the supplemental material. Annual snowfall and snow depth climatologies for 1981–2000 are computed for a different set of 196 stations across the study region from NCDC local climatological data publications (available at http://www.ncdc.noaa.gov/IPS/lcd/lcd.html) and Environment Canada’s Canadian Climate Normals (available at http://climate.weather.gc.ca/climate_normals) and used to tune and evaluate SNOW-17. Elevation is retrieved from the 2′ Gridded Global Relief Data (ETOP02; National Geophysical Data Center 2013).

e. SNOW-17 simulations and their evaluation

After calibrating SNOW-17 and its parameters across the study region, the model is run for the late twentieth century (1981–2000), mid-twenty-first century (2046–65), and late twenty-first century (2081–2100), according to the A2 and B1 emission scenarios. For each time period, statistically downscaled climate data from nine CMIP3 GCMs are used to force the snow model across the 10 LCC regions of central-eastern North America (Fig. 1). SNOW-17 generally requires hourly climatic inputs. Downscaled daily maximum and minimum temperature is converted to hourly temperature data using a cubic spline under tension, with the minimum and maximum temperature assigned to 1200–1400 and 2000–2200 UTC, respectively (depending on longitude). This interpolation approach is successful in maintaining different projected changes in maximum versus minimum temperature, recognizing that most GCMs simulate a greater increase in minimum temperature. Sensitivity experiments, in which the timing of maximum or minimum temperature is shifted ahead or back by one or two hours, demonstrate that simulated snowfall patterns are robust and not overly sensitive to small adjustments in the timing of the temperature diurnal cycle. Daily precipitation amounts are converted into hourly values by randomly distributing precipitation on wet days into continuous 6-h periods (0000–0600, 0600–1200, 1200–1800, or 1800–2400 UTC), which is the typical duration of a cold-season precipitation event based on studies for the Midwest and northeastern United States (Cox and Armington 1914; Changnon 1969; Knapp et al. 2000; Laird et al. 2009).

Fig. 1.

Map of elevation (m) based on ETOP02 global topographic grid, with red polygons outlining the relevant LCCs and black dots for the 196 stations. Snowfall climatology at these stations is used to assess the performance of SNOW-17 and estimate the snow correction factor (SCF) parameter in SNOW-17. The 10 LCC regions of interest are identified by number, for each LCC region, the largest populated city is identified.

Fig. 1.

Map of elevation (m) based on ETOP02 global topographic grid, with red polygons outlining the relevant LCCs and black dots for the 196 stations. Snowfall climatology at these stations is used to assess the performance of SNOW-17 and estimate the snow correction factor (SCF) parameter in SNOW-17. The 10 LCC regions of interest are identified by number, for each LCC region, the largest populated city is identified.

SNOW-17 requires regionally specific calibration to produce quality snow simulations. The model applies six major parameters (SCF, MFMAX, MFMIN, UADJ, SI, and ADC) and six minor parameters (MBASE, NMF, DAYGM, PLWHC, PXTEMP, and TIPM). The methodology for determining appropriate parameter values is outlined in the supplemental materials. The model is tuned to best capture climatological patterns of annual snowfall and DJF snow depth, compared to observations at 196 stations (Fig. 1).

Using the resulting parameter values, SNOW-17 is forced by statistically downscaled climate data from nine CMIP3 GCMs for 1981–2000. The simulated climatological patterns of annual snowfall and DJF snow depth (averaged among nine GCMs) are evaluated against the set of 196 stations (Fig. 2). The observed and simulated mean annual snowfall closely match, with a correlation of 0.99 (ranging from 0.95 in April to 0.99 in December–January), root-mean-square difference of 5.9 cm (ranging from 2.8 cm in November to 3.4 cm in December), and a mean absolute difference of 3.8 cm. Because of deficiencies in assigned parameter values, the model produces excessive snowfall in December–January (+8.3% in December) and too little in November and February–April (−14.8% in November). Regarding DJF climatological snow depth, the model and observations exhibit slightly less agreement, with a correlation of 0.90 (range from 0.86 in March to 0.95 in November), root-mean-square difference of 4.4 cm (range from 0.4 cm in November to 6.7 cm in February), and a mean absolute difference of 1.3 cm. The simulated snowpack is too deep in spring, because of insufficient snowmelt. Percent biases range from −7% in November to +37% in March. Overall, SNOW-17 performs well in reproducing the modern-day snow climatology across the vast study region when forced with statistically downscaled climate data, with a distinct improvement over the original GCMs (Fig. S7 in the supplementary material).

Fig. 2.

Scatterplots of observed vs simulated (a) annual snowfall (cm) and (b) DJF snow depth for 196 stations based on a climatology for 1981–2000. Simulated snow is based on output from SNOW-17 forced by downscaled climate data. For snowfall (snow depth) the temporal correlation is 0.99 (0.90) and root-mean-square difference is 5.9 cm (4.4 cm). Scatterplots for individual months are shown for snowfall [labeled as (a1)–(a6)] and snow depth [labeled as (b1)–(b6)], with the correlation coefficient provided in each panel (p < 0.001 for all panels).

Fig. 2.

Scatterplots of observed vs simulated (a) annual snowfall (cm) and (b) DJF snow depth for 196 stations based on a climatology for 1981–2000. Simulated snow is based on output from SNOW-17 forced by downscaled climate data. For snowfall (snow depth) the temporal correlation is 0.99 (0.90) and root-mean-square difference is 5.9 cm (4.4 cm). Scatterplots for individual months are shown for snowfall [labeled as (a1)–(a6)] and snow depth [labeled as (b1)–(b6)], with the correlation coefficient provided in each panel (p < 0.001 for all panels).

3. Results

a. Downscaled climate projections

At every location in the study region, and according to both emission scenarios, warming is a robust projection for December–March (DJFM) by the mid and late twenty-first century by all nine downscaled GCMs. According to the B1 and A2 scenarios, the area-average increase in surface air temperature by the mid-twenty-first century is +2.0°C (range from +1.3° to +2.6°C) and +2.9°C (range from +2.0° to +4.3°C) and by the late twenty-first century it is +3.1°C (range from +2.1° to +4.2°C) and +5.3°C (range from +4.1° to +6.2°C), respectively (Fig. 3). The scenarios diverge by the late twenty-first century in terms of the magnitude of projected warming. The magnitude of warming increases with latitude, with a local minimum across the Great Lakes basin. The greatest mean projected warming by the late twenty-first century, from the A2 scenario, is found across the Upper Midwest and Great Lakes LCC (+6.0°C) and Plains and Prairie Potholes LCC (+5.7°C) (Tables 1 and 2; Fig. 3). In this study, at each grid cell, the “low-end projection” and “high-end projection” in temperature refer to the least and most warming, respectively, among the nine downscaled GCMs. The A2-based low-end and high-end projections by the late twenty-first century differ dramatically by roughly 3.5°–4°C across Michigan and southern Canada (Figs. 3a,c,d), indicating consistent projections of warming but much uncertainty regarding its magnitude (particularly in minimum temperatures). This uncertainty may be attributed to 1) greater natural variability with increased latitude, 2) variations in the intensity of snow–ice feedbacks among GCMs, and 3) differing representations of the Great Lakes among the original GCMs (e.g., GCMs excluding the lakes may warm more regionally).

Fig. 3.

Projected change in DJFM mean temperature (°C) by the late twenty-first century, computed as the difference between 2081–2100 and 1981–2000. Results are shown for the (a)–(d) A2 and (e)–(h) B1 emission scenarios and are broken down by (a),(e) low-end projection, (b),(f) mean projection, (c),(g) high-end projection, and (d),(h) spread (difference between high-end and low-end projections) per grid cell. The upper color bar pertains to (a)–(c),(e)–(g) and the lower color bar pertains to (d),(h). Among nine GCMs, the greatest warming at each grid cell is considered high-end and the least warming is considered low-end. All of the differences shown in (a)–(c),(e)–(g) are statistically significant (p < 0.1), except in (e) over central Montana and southern Alberta as identified with hatching.

Fig. 3.

Projected change in DJFM mean temperature (°C) by the late twenty-first century, computed as the difference between 2081–2100 and 1981–2000. Results are shown for the (a)–(d) A2 and (e)–(h) B1 emission scenarios and are broken down by (a),(e) low-end projection, (b),(f) mean projection, (c),(g) high-end projection, and (d),(h) spread (difference between high-end and low-end projections) per grid cell. The upper color bar pertains to (a)–(c),(e)–(g) and the lower color bar pertains to (d),(h). Among nine GCMs, the greatest warming at each grid cell is considered high-end and the least warming is considered low-end. All of the differences shown in (a)–(c),(e)–(g) are statistically significant (p < 0.1), except in (e) over central Montana and southern Alberta as identified with hatching.

There is even greater uncertainty in DJFM projections of precipitation than temperature among GCMs. The area-average projected change in DJFM precipitation by the mid-twenty-first century is +0.09 mm day−1 (range from −0.01 to +0.21 mm day−1), or +5.3%, according to the B1 scenario and +0.13 mm day−1 (range from +0.03 to +0.25 mm day−1), or +7.4%, according to the A2 scenario (Fig. 4); projections for the late twenty-first century are +0.13 mm day−1 (range from −0.01 to +0.31 mm day−1), or +7.4%, for the B1 scenario and +0.20 mm day−1 (range from +0.00 to +0.32 mm day−1), or +11.5%, for the A2 scenario. The model-mean DJFM projections for the late twenty-first century, based on the A2 scenario, range from −0.19 mm day−1 (−10.7%) in the Gulf Coast Prairie LCC to +0.61 mm day−1 (+21.7%) in the North Atlantic LCC, with the largest increases over the northeastern United States (Fig. 4b, Tables 1 and 2). Precipitation projections largely disagree among the GCMs across the Gulf Coastal Plains and Ozarks LCC, ranging from drying of about −1 mm day−1 to moistening of about +1.0 mm day−1 (Figs. 4a,c,d). While a potential decline in precipitation across the southern LCCs would reinforce the effect of warming to reduce snowfall, a robust projection of greater precipitation across the northernmost LCCs should counter the impact of warming on future snowfall.

Fig. 4.

As in Fig. 3, but for projected change in DJFM mean precipitation (mm day−1). Among nine GCMs, the greatest increase in precipitation at each grid cell is considered high-end and the greatest drying (or least increase in precipitation) is considered low-end. Differences that fail to attain statistical significance in (a)–(c),(e)–(g) are dotted (p > 0.1).

Fig. 4.

As in Fig. 3, but for projected change in DJFM mean precipitation (mm day−1). Among nine GCMs, the greatest increase in precipitation at each grid cell is considered high-end and the greatest drying (or least increase in precipitation) is considered low-end. Differences that fail to attain statistical significance in (a)–(c),(e)–(g) are dotted (p > 0.1).

b. Benefits of downscaling for snow simulations

Given the coarse horizontal resolution of the CMIP3 GCMs (mean grid cell size of 7.4 × 104 km2) and insufficient representation of the Great Lakes (either completely absent or represented by a few grid cells), lake-effect processes are largely neglected; unfortunately, this limitation continues to exist in most CMIP5 GCMs. However, by including wind as a precipitation predictor in the statistical downscaling, lake-effect precipitation events are captured in the downscaling product, and therefore SNOW-17 simulations forced by this downscaled data include lake-effect snowstorms. To illustrate this point, two December days are selected from the models GISS-ER and ECHAM5 (see Table 3 for model information), with cold air outbreaks and associated west-northwesterly low-level winds over the lakes at a time of the year when ice cover is usually limited and therefore lake-effect snow would be expected (Figs. 5a,b). However, in both cases, there is no signature of lake-effect precipitation in the GCM (Figs. 5a,b). SNOW-17, when forced by the statistical downscaling data, produces distinct lake-effect snowfall maxima downwind of the lakes on both days (Figs. 5c,d). The downscaling relies on the strength of GCMs, which is the representation of large-scale circulation patterns, and adds accuracy to the projections by relating the GCM’s large-scale circulation patterns to the expected local climatic response, as established through observations.

Fig. 5.

Simulated daily precipitation (mm) and 850-hPa wind vectors from a select day in (a) December 1989 from the 20C3M simulation of GISS-ER and in (b) December 1987 from the 20C3M simulation of ECHAM5, which include northwesterly winds over the Great Lakes that should support lake-effect snowfall. Reference vectors are assigned values of 5 and 10 m s−1 for (a) and (b), respectively. For these two events, simulated snowfall (cm) is shown for (c) GISS-ER and (d) ECHAM5 based on SNOW-17 simulations forced by downscaled climate data. Lake-effect precipitation is absent in (a) and (b) but captured in (c) and (d).

Fig. 5.

Simulated daily precipitation (mm) and 850-hPa wind vectors from a select day in (a) December 1989 from the 20C3M simulation of GISS-ER and in (b) December 1987 from the 20C3M simulation of ECHAM5, which include northwesterly winds over the Great Lakes that should support lake-effect snowfall. Reference vectors are assigned values of 5 and 10 m s−1 for (a) and (b), respectively. For these two events, simulated snowfall (cm) is shown for (c) GISS-ER and (d) ECHAM5 based on SNOW-17 simulations forced by downscaled climate data. Lake-effect precipitation is absent in (a) and (b) but captured in (c) and (d).

Compared to the raw GCMs, the use of statistically downscaled climate data in SNOW-17 results in a more reasonable spatial pattern of annual mean climatological snowfall across the study region. The mean observed snowfall peaks across the northeastern United States, southern Ontario and Quebec, the Great Lakes basin, the Appalachian Mountains, and on the lee side of the Rocky Mountains (Figs. 2 and 6c). The mean snowfall climatology among the raw CMIP3 models crudely captures this spatial pattern but contains significant biases and minimal signature of the Great Lakes or Appalachian Mountains (Fig. 6a). However, when forced by downscaled climate data, SNOW-17 produces minimal biases in annual snowfall and accurately represents the lake-effect snow maxima and elevated snowfall over the Appalachians (Fig. 6b). For the comparisons in Fig. 6 (and Fig. S7), given that the GCMs only saved the liquid equivalent of snowfall, a 10: 1 ratio is applied to crudely convert the liquid-equivalent snowfall to simulated snowfall depth. Krasting et al. (2013) likewise applied a 10: 1 ratio to the CMIP3 snowfall mass flux output but cautioned that this ratio has significant observed regional variations. The lake-effect zones of the Great Lakes basin are typically characterized by snow-to-liquid-equivalent ratios of 14: 1 or 15: 1 (Baxter et al. 2005; Kutikoff 2013).

Fig. 6.

Climatological mean annual snowfall (cm) from the (a) raw CMIP3 GCMs, (b) downscaled SNOW-17, and (c) 196 weather stations for 1981–2000. Since the GCMs only saved the liquid equivalent of snowfall, a 10: 1 ratio was applied to crudely estimate simulated snowfall.

Fig. 6.

Climatological mean annual snowfall (cm) from the (a) raw CMIP3 GCMs, (b) downscaled SNOW-17, and (c) 196 weather stations for 1981–2000. Since the GCMs only saved the liquid equivalent of snowfall, a 10: 1 ratio was applied to crudely estimate simulated snowfall.

c. Snow projections

The model-mean projections, for both emission scenarios and both time periods, indicate a warming-induced decline in annual snowfall across the entire study region (Figs. 7b,e,h,k). The area-average mean projections for annual snowfall, by the mid-twenty-first century, are −17.9 cm (range from −8.4 to −28.6 cm), or −16.6%, according to the B1 scenario and −26.0 cm (range from −19.0 to −37.8 cm), or −24.2%, according to the A2 scenario (Fig. 7); these projections for the late twenty-first century are −28.4 cm (range from −14.6 to −35.9 cm), or −26.4%, for the B1 scenario and −45.0 cm (range from −32.5 to −55.7 cm), or −41.8%, for the A2 scenario. The largest declines, by the late twenty-first century, are projected for the North Atlantic LCC (−58.8 cm, or −30.8%, for B1 versus −92.0 cm, or −48.1%, for A2) and the Upper Midwest and Great Lakes LCC (−53.9 cm, or −24.1%, for B1 versus −89.9 cm, or −40.2%, for A2) (Tables 1 and 2), as precipitation falls more in the form of rain than snow. By the late twenty-first century, the two emission scenarios result in similar patterns of projected snowfall decline, although the high-end projection (defined as the largest decline) for A2 produces about 25% greater snowfall declines across the domain than does B1 (Figs. 7f,l). According to the low-end projection (smallest decline, or a potential increase) for annual snowfall, relatively modest increases are possible across southern Canada and the Plains and Prairies Potholes LCC, and the southern states, as increases in cold season precipitation offset the warming trend (Figs. 7a,d,g,j).

Fig. 7.

Projected change in annual snowfall (cm) by the (a)–(c),(g)–(i) mid- (2046–65) and (d)–(f),(j)–(l) late twenty-first century (2081–2100), computed as the difference from the late twentieth century (1981–2000), according to the (a)–(f) B1 and (g)–(l) A2 emission scenarios. Results are broken down by (top) low-end, (middle) mean, and (bottom) high-end projections per grid cell. Among the nine GCMs, the largest decline in snowfall at each grid cell is considered high end and the greatest increase (or smallest decline) is considered low end.

Fig. 7.

Projected change in annual snowfall (cm) by the (a)–(c),(g)–(i) mid- (2046–65) and (d)–(f),(j)–(l) late twenty-first century (2081–2100), computed as the difference from the late twentieth century (1981–2000), according to the (a)–(f) B1 and (g)–(l) A2 emission scenarios. Results are broken down by (top) low-end, (middle) mean, and (bottom) high-end projections per grid cell. Among the nine GCMs, the largest decline in snowfall at each grid cell is considered high end and the greatest increase (or smallest decline) is considered low end.

The mean frequency of snowfall days, in excess of 1 cm, is projected to decline everywhere by the mid and late twenty-first century, according to model-mean results from both emission scenarios (Figs. 8a–d). During the late twentieth century, SNOW-17 simulates a mean frequency of snowfall days that ranges from 0.01 day in the Peninsular Florida LCC to 39.08 days in the Upper Midwest and Great Lakes LCC. The largest model-mean declines in snowfall days are projected for the Upper Midwest and Great Lakes LCC (for B1: mean = −8.9 days, or −22.9%, range from −2.6 to −13.8 days; and for A2: mean = −15.6 days, or −39.9%, range from −12.2 to −19.6 days) and the North Atlantic LCC (for B1: mean = −8.2 days, or −28.2%, range from −4.4 to −11.8 days; and for A2: mean = −13.7 days, or −46.7%, range from −8.6 to −18.0 days) (Tables 1 and 2, Figs. 8a–d). According to the B1 scenario, a modest increase in snowfall day frequency by the mid-twenty-first century is identified across northern Wisconsin, Minnesota, and Montana for the low-end projection (defined as smallest decline, or greatest increase, in snowfall days among nine GCMs).

Fig. 8.

Mean projected percentage change in the annual (a)–(d) frequency of daily snowfall events of at least 1 cm and (e)–(h) mean snowfall per event by the (a),(c),(e),(g) mid- (2046–65) and (b),(d),(f),(h) late twenty-first century (2081–2100), compared to the late twentieth century (1981–2000). Results are shown both for the (a),(b),(e),(f) B1 and (c),(d),(g),(h) A2 emission scenarios. The median projections are also computed (not shown) and found to be quite consistent with the mean projections shown here.

Fig. 8.

Mean projected percentage change in the annual (a)–(d) frequency of daily snowfall events of at least 1 cm and (e)–(h) mean snowfall per event by the (a),(c),(e),(g) mid- (2046–65) and (b),(d),(f),(h) late twenty-first century (2081–2100), compared to the late twentieth century (1981–2000). Results are shown both for the (a),(b),(e),(f) B1 and (c),(d),(g),(h) A2 emission scenarios. The median projections are also computed (not shown) and found to be quite consistent with the mean projections shown here.

Despite large projected declines in annual-mean snowfall and snowfall frequency across the study region, SNOW-17 produces only slight changes in mean snowfall intensity, or total daily snowfall per event (Figs. 8e–h). Most areas are simulated to experience lighter snowfall events, with declines in mean snowfall intensity simulated across 64% of the study region by the mid-twenty-first century with the B1 scenario, and across 79% of the region by the late twenty-first century with the A2 scenario (Figs. 8e–h). However, an increase in snowfall intensity is possible for southern Canada and the Gulf states (Figs. 8e–h). Southern Canada may experience lower annual snowfall totals, fewer snowfall days, and higher snowfall totals per event (Figs. 7 and 8), which suggests fewer, but more intense, snowfall events. Model-mean projections for snowfall totals per event, by the late twenty-first century, are modest and differ notably across models. Projected intensity changes for the Great Plains LCC are −0.26 cm (range from −1.00 to +0.54 cm), or −4.1%, for the B1 scenario and −0.52 cm (range from −1.29 to +0.24 cm), or −8.4%, for the A2 scenario (Tables 1 and 2). For the Gulf Coastal Plains and Ozarks LCC, these projections include +0.10 cm (range from −0.58 to +1.24 cm), or +2.2%, for the B1 scenario and −0.62 cm (range from −1.37 to +0.03 cm), or −14.3%, for the A2 scenario; in this region, snowfall events are rare, making these conclusions uncertain. A distinct climate change signal in snow event intensity is generally difficult to detect, with prior studies rarely tackling this area of uncertainty.

Projected percentage changes in the frequency of snowfall events of different daily magnitudes are computed both across the study region and individual LCC regions (Fig. 9). For both scenarios and time periods, the projected percentage reduction in the frequency of daily snowfall events is maximized around 15–23 cm, ranging from −21.5% by the mid-twenty-first century under the B1 scenario to −47.6% by the late twenty-first century under the A2 scenario (Fig. 9). Consistent with GCM projections for heavier precipitation events later this century, SNOW-17 simulates more frequent intense snowfall events, particularly by the mid-twenty-first century, under the low-end emission scenario (B1), and across the north-northwest study region. Specifically, daily snowfall events exceeding (under) 70 cm are projected to become more (less) frequent across the study region by the mid-twenty-first century, for the B1 scenario (Fig. 9). Under more extreme warming projections (late twenty-first century or A2), the potential for heavier snowfall events diminishes substantially (Fig. 9). The projected increase in frequency of heavy daily snow events is primarily confined to the Plains and Prairie Potholes LCC and Upper Midwest and Great Lakes LCC (Fig. 9), for events exceeding 56 cm for the mid-twenty-first century and B1 scenario or 76 cm for the late twenty-first century and A2 scenario.

Fig. 9.

Projected percentage change in the frequency of daily snowfall (cm) events, within the entire 10 LCC study region (blue lines) and the region encompassing the Plains and Prairie Potholes LCC and Upper Midwest and Great Lakes LCC (red lines). Projections are shown for the (a),(b) B1 and (c),(d) A2 scenarios and for the (a),(c) mid- and (b),(d) late twenty-first century. The x axis consists of 1-cm daily snowfall bins (0.1–1, 1–2, 2–3, …, 99–100) for the former region and 2-cm bins (0.1–2, 2–4, 4–6, …, 98–100) for the latter region. Thin black curves are also included for each of the following LCC regions: Appalachian, Eastern Tallgrass Prairie and Big Rivers, Plains and Prairie Potholes, Great Plains, Upper Midwest and Great Lakes, and North Atlantic.

Fig. 9.

Projected percentage change in the frequency of daily snowfall (cm) events, within the entire 10 LCC study region (blue lines) and the region encompassing the Plains and Prairie Potholes LCC and Upper Midwest and Great Lakes LCC (red lines). Projections are shown for the (a),(b) B1 and (c),(d) A2 scenarios and for the (a),(c) mid- and (b),(d) late twenty-first century. The x axis consists of 1-cm daily snowfall bins (0.1–1, 1–2, 2–3, …, 99–100) for the former region and 2-cm bins (0.1–2, 2–4, 4–6, …, 98–100) for the latter region. Thin black curves are also included for each of the following LCC regions: Appalachian, Eastern Tallgrass Prairie and Big Rivers, Plains and Prairie Potholes, Great Plains, Upper Midwest and Great Lakes, and North Atlantic.

With larger projected snowfall declines in autumn than spring, SNOW-17 generally simulates a delayed onset of the snow season, more prominently than an earlier termination (Fig. 10). During the cold season (November–April), the study region is projected to warm the most during November–December, with the largest snowfall reductions during December–January (see late twenty-first century and A2 scenario in Fig. 10). The primary exception is found across the northern Great Plains, particularly South Dakota, Nebraska, Colorado, and Wyoming, with the greatest projected decline in snowfall during springtime (March–April); these regions are characterized by high rain–snow temperature thresholds exceeding 1.5°C (Fig. S5a in the supplementary material).

Fig. 10.

Month of maximum projected decline in snowfall, according to the A2 emission scenario by the end of the twenty-first century (2081–2100), compared to the late twentieth century (1981–2000). SNOW-17 is forced by downscaled climate projections from nine GCMs, the mean snowfall projection is computed among these models, and then the month of maximum decline in snowfall is identified.

Fig. 10.

Month of maximum projected decline in snowfall, according to the A2 emission scenario by the end of the twenty-first century (2081–2100), compared to the late twentieth century (1981–2000). SNOW-17 is forced by downscaled climate projections from nine GCMs, the mean snowfall projection is computed among these models, and then the month of maximum decline in snowfall is identified.

For all 10 LCC regions, the projected percentage reduction in mean November–April snow depth exceeds that of snowfall, as projected warming also accelerates snowmelt. For example, the North Atlantic LCC is expected to experience a 48.1% reduction in annual snowfall and a 72.1% reduction in November–April snow depth by the late twenty-first century according to the A2 scenario (Tables 1 and 2, Fig. 11). The largest projected reductions in snow depth are found within the Upper Midwest and Great Lakes LCC, with a decline by the late twenty-first century of −7.63 cm (−36.9%) for the B1 scenario versus −11.94 cm (−57.8%) for the A2 scenario. Across the study region, November–April mean snow depth is projected to decline by −1.79 cm (range from −1.12 to −2.67 cm) according to the B1 scenario and −2.56 cm (range from −1.95 to −3.30 cm) according to the A2 scenario, by the mid-twenty-first century; late twenty-first-century projections include declines of −2.72 cm (range from −1.50 to −4.08 cm) for the B1 scenario and −4.11 cm (range from −3.17 to −4.78 cm) for the A2 scenario. Here, the snow line, or southernmost mean extent of snow cover, is defined by a criterion of at least 0.5 cm of mean snow depth during November–April. The snow line across the study region is anticipated to shift northward by 72 km (128 km) by the mid-twenty-first century and 126 km (229 km) according to the B1 (A2) scenario. The projected northward migration rate of the snow line is less than half of that of the mean isotherms during DJFM (for −5° to 5°C), implying a more gradual response to climate change; the mean isotherms are projected to shift northward at a much more rapid rate than that of the snow line. The low-end projection among the nine GCMs suggests the potential for a slight increase in annual mean snowfall and snow depth across the southern states, with a large percentage increase (Figs. 11a,d,g,j).

Fig. 11.

As in Fig. 7, but for projected percentage change in NDJFMA mean snow depth.

Fig. 11.

As in Fig. 7, but for projected percentage change in NDJFMA mean snow depth.

Projected changes in the annual number of days with existing snowpack (≥1 cm) are assessed (Fig. 12). The regions likely to experience the largest reductions in days with snowpack are the Upper Midwest and Great Lakes LCC, North Atlantic LCC, and Plains and Prairie Potholes LCC, with respective declines of −27.0 days (−19.7%, range from −16.2 to −42.8 days), −28.5 days (−30.3%, range from −18.3 to −36.7 days), and −26.3 days (−22.2%, range from −17.1 to −47.2 days) for the B1 scenario and −47.6 days (−34.8%, range from −35.7 to −62.0 days), −47.3 days (−50.3%, range from −31.8 to −60.0 days), and −42.7 days (−35.9%, range from −29.5 to −55.9 days) for the A2 scenario (Tables 1 and 2). An optimal zone of maximum reduction in days with snowpack is identified near 44°N, closely following the −5°C isotherm of the late twentieth-century mean climatology (Fig. 12). As this isotherm shifts northward, projected trends in declining snowfall accelerate. Ratios between the projected snowfall decline by the late twenty-first century and the projected decline by the mid-twenty-first century are on the order of 1.4–1.6 to the south of this isotherm and 1.6–2.0 to the north (Fig. 13). The northern LCCs can expect an accelerated negative trend in annual snowfall during this century.

Fig. 12.

Mean projected change in the number of days per year with a minimum snowpack of 1 cm by the (a),(c) mid- (2046–65) and (b),(d) late twenty-first century (2081–2100), compared to the late twentieth century (1981–2000). Results are shown for the (a),(b) B1 and (c),(d) A2 emission scenarios. The green curve represents the −5°C mean isotherm for DJFM for 1981–2000. Scatterplots display (y axis) mean climatological DJFM temperature (°C) from 1981–2000 vs (x axis) projected change in number of days per year with a minimum snowpack, with a dot for each grid cell.

Fig. 12.

Mean projected change in the number of days per year with a minimum snowpack of 1 cm by the (a),(c) mid- (2046–65) and (b),(d) late twenty-first century (2081–2100), compared to the late twentieth century (1981–2000). Results are shown for the (a),(b) B1 and (c),(d) A2 emission scenarios. The green curve represents the −5°C mean isotherm for DJFM for 1981–2000. Scatterplots display (y axis) mean climatological DJFM temperature (°C) from 1981–2000 vs (x axis) projected change in number of days per year with a minimum snowpack, with a dot for each grid cell.

Fig. 13.

Ratio of the mean projected change in annual snowfall by the late twenty-first century (2081–2100) to the mean projected change by the mid-twenty-first century (2046–65), compared to the twentieth century (1981–2000). Results are based on the A2 emission scenario. Ratios greater than one indicate accelerated snowfall declines later in the twenty-first century. The black dashed line represents the −5°C mean isotherm for DJFM for 1981–2000.

Fig. 13.

Ratio of the mean projected change in annual snowfall by the late twenty-first century (2081–2100) to the mean projected change by the mid-twenty-first century (2046–65), compared to the twentieth century (1981–2000). Results are based on the A2 emission scenario. Ratios greater than one indicate accelerated snowfall declines later in the twenty-first century. The black dashed line represents the −5°C mean isotherm for DJFM for 1981–2000.

d. Case study application: Projected waterfowl distributions

To illustrate the utility of the SNOW-17 based projections, potential impacts of changing winter conditions on mallard ducks are assessed as a case study. CWSI is computed for the late twentieth, mid-twenty-first, and late twenty-first centuries, according to the A2 and B1 scenarios (Tables 1 and 2, Fig. 14). Analysis focuses on the CWSI = 7.2 isopleth, which is the threshold at which mallards typically migrate southward and locally decrease in abundance (Schummer et al. 2010), thereby providing an index of relative duck distribution. During Januaries in the late twentieth century, this isopleth was typically located across central Kansas, southern Illinois, and southern Ohio (Fig. 14j). In response to warming and declining snowpack, this isopleth is projected to shift northward, reaching northern Nebraska, northern Illinois, and central Michigan by the late twenty-first century under the A2 scenario. The arrival of this isopleth from Canada into the United States could generally be delayed by roughly one month in late autumn, by the late twenty-first century. The largest reductions in mean DJF CWSI are expected for the Upper Midwest and Great Lakes LCC, from 94.9 during the late twentieth century to either 60.5 or 40.3 for the B1 and A2 scenarios (Tables 1 and 2). Within the study region, the probability of WSI > 7.2 (likely southward migration) is projected to decline especially during December across the northern LCCs (e.g., Plains and Prairie Potholes) and during January across the central LCCs (e.g., Eastern Tallgrass Prairie and Big Rivers) (Fig. 14). For example, during Decembers in the late twentieth century, there was a 62.5% probability of WSI > 7.2 across the North Atlantic LCC. By the late twenty-first century, this probability is projected to decrease to 40.4% for the B1 scenario or 19.8% for the A2 scenario (Fig. 14f), suggesting a delayed southward migration during early winter.

Fig. 14.

(a)–(f) Mean probability of the daily CWSI in specific LCC regions exceeding 7.2 from November to March, for the late twentieth century (black), mid-twenty-first century [B1 (blue lines), A2 (green lines)], and late twenty-first century [B1 (orange lines), A2 (red lines)], based on mean projections from nine GCMs. (g)–(l) Mean CWSI for October–March for the late twentieth century. The thick black lines indicate CWSI = 7.2, a critical threshold for mallard ducks. Note that the CWSI is cumulatively summed each day from 1 September to 31 March.

Fig. 14.

(a)–(f) Mean probability of the daily CWSI in specific LCC regions exceeding 7.2 from November to March, for the late twentieth century (black), mid-twenty-first century [B1 (blue lines), A2 (green lines)], and late twenty-first century [B1 (orange lines), A2 (red lines)], based on mean projections from nine GCMs. (g)–(l) Mean CWSI for October–March for the late twentieth century. The thick black lines indicate CWSI = 7.2, a critical threshold for mallard ducks. Note that the CWSI is cumulatively summed each day from 1 September to 31 March.

4. Discussion and conclusions

High-resolution projections of daily snowfall, snow depth, and winter severity for the mid-twenty-first and late twenty-first century are generated for the central-eastern North American LCCs by forcing SNOW-17 with statistically downscaled climate data from nine CMIP3 GCMs and two emission scenarios. These cold season projections will aid natural resource agencies in assessing impacts of changes in snow and winter severity to wildlife. Compared to raw GCM snow output, this approach leads to a better representation of weather extremes, lake-effect processes, and topographic influences, along with smaller biases in snowfall and snow depth.

For both time periods and both emission scenarios, the entire region is expected to experience reduced annual snowfall and a shift toward less snow and more rain, based on model-mean projections. The largest snowfall losses are simulated for the North Atlantic LCC and the Upper Midwest and Great Lakes LCC, with respective declines by the late twenty-first century of −48.1% and −40.2% under the A2 scenario. For much of the region, snowfall is projected to decline particularly during December–January, implying a delayed onset of the snow season, except over the Northern Plains, where the greatest losses are expected in March–April when snowstorms are more common (Changnon et al. 2006). These results are consistent with studies by Kapnick and Delworth (2013) and Krasting et al. (2013), which project a general decline in Northern Hemispheric snowfall, including across North America, outside of the high latitudes. Considering the range of outcomes from nine GCMs for the B1 scenario, there is a potential for greater annual snowfall across southern Canada, north of the historical mean position of the −10°C isotherm in DJFM, by both the mid and late twenty-first century, in response to increased precipitation. Similarly, Krasting et al. (2013) concluded that snowfall would increase to the north of this isotherm and decrease to the south. SNOW-17 simulates a trend toward fewer snowfall days, but with modest changes in mean intensity. A larger percentage decline in mean snow depth is simulated than in annual snowfall, because of enhanced snowmelt rates with warming. The mean snow line is simulated to shift northward, ranging from 72 km by the mid-twenty-first century according to the B1 scenario to 229 km by the late twenty-first century according to the A2 scenario. The projected decline in the number of days with existing snowpack is expected to optimally peak around 44°N, close to the mean −5°C DJFM isotherm. Dramatic decreases are projected in mean CWSI and the probability of CWSI > 7.2, with the arrival of this critical CWSI isopleth from Canada into the United States becoming delayed by a month by the late twenty-first century, indicative of a delayed southward migration of mallards; this shift will impact the $3.4 billion migratory bird hunting industry (U.S. Fish and Wildlife Service 2011; Southwick Associates 2012).

Few studies have addressed the impact of future climate change on the frequency of heavy snow events (Gutowski et al. 2008; López-Moreno et al. 2011), with much uncertainty in the future response in mean snow event intensity. Changnon (2007) found that 88% of the economic losses attributed to United States’ snowstorms during 1949–2003 occurred in the nation’s eastern half. Snowstorms in the United States have become less frequent, but larger and more intense (Changnon 2007). As temperatures continue to rise, the atmospheric moisture-holding capacity will increase, favoring more intense precipitation events, potentially even snowstorms (CMAP 2013). The current study provides support for intensified snowstorms during this century. Heavy daily snow events are projected to increase in frequency, particularly across the Plains and Prairie Potholes LCC and the Upper Midwest and Great Lakes LCC and most notably for low-warming scenarios.

Several limitations are identified in the study. SNOW-17 only requires inputs of ground-level temperature and precipitation, although ground-level temperature is not a perfect indicator of precipitation type (Anderson 2006). SNOW-17 requires an extensive set of spatially explicit parameters, which can be challenging to develop as most prior applications focused on single watersheds. The applied parameter values result in insufficient snowmelt and excessive snow depth in spring. The LCC statistical downscaling is based on output from a small set of nine CMIP3 GCMs, rather than the latest CMIP5 archive. Only one realization for each CMIP3 model is considered; thus, the full range of possible climate outcomes may be insufficiently represented. This is particularly a concern for projections out to the mid-twenty-first century, given the importance of internal natural variability on the decadal-to-multidecadal time scale (Deser et al. 2012a,b, 2014); future studies should consider large ensemble sets to represent uncertainty. The snowfall projections need to be considered with caution for the Great Lakes basin, given that the statistical downscaling does not consider future changes in lake ice cover, which regulates lake-effect snow. Observations have shown a positive trend in the Great Lakes lake-effect snowfall (Norton and Bolsenga 1993; Leathers and Ellis 1996; Burnett et al. 2003), in response to greater lake evaporation, although an observed trend reversal toward less lake-effect snowfall might have recently begun (Bard and Kristovich 2012). To address future changes in lake-effect snow, a combination of dynamical downscaling, using regional climate models coupled to interactive lakes (Notaro et al. 2013), and improved statistical downscaling techniques are encouraged.

Acknowledgments

The research was funded by the Upper Midwest and Great Lakes LCC, Gulf Coastal Plains and Ozarks LCC, Wisconsin Focus on Energy, and the U.S. Environmental Protection Agency’s Great Lakes Restoration Initiative through a contract with the Michigan Department of Natural Resources. CMIP3 data were obtained through the Program for Climate Model Diagnosis and Intercomparison (PCMDI). The authors appreciate the helpful comments from three anonymous reviewers.

REFERENCES

REFERENCES
Adeloye
,
A. J.
,
N. R.
Nawaz
, and
M.
Montaseri
,
1999
:
Climate change water resources planning impacts incorporating reservoir surface net evaporation fluxes: A case study
.
Water Resour. Dev.
,
15
,
561
581
, doi:.
Anderson
,
E. A.
,
1973
: National Weather Service River Forecast System–Snow Accumulation and Ablation Model. NOAA Tech. Memo. NWS Hydro-17, 217 pp.
Anderson
,
E. A.
,
1976
: A point energy and mass balance model of a snow cover. NOAA Tech. Rep. NWS 19, 150 pp.
Anderson
,
E. A.
,
2002
: Calibration of conceptual hydrologic models for use in river forecasting. Office of Hydrologic Development, U.S. National Weather Service, Silver Spring, MD, 372 pp.
Anderson
,
E. A.
,
2006
: Snow Accumulation and Ablation Model—SNOW 17. U.S. National Weather Service, Silver Spring, MD, 61 pp.
Bard
,
L.
, and
D. A.
Kristovich
,
2012
:
Trend reversal in Lake Michigan contribution to snowfall
.
J. Appl. Meteor. Climatol.
,
51
,
2038
2046
, doi:.
Barnett
,
T. P.
,
J. C.
Adam
, and
D. P.
Lettenmaier
,
2005
:
Potential impacts of a warming climate on water availability in snow-dominated regions
.
Nature
,
438
,
303
309
, doi:.
Baxter
,
M. A.
,
C. E.
Graves
, and
J. T.
Moore
,
2005
:
A climatology of snow-to-liquid ratio for the contiguous United States
.
Wea. Forecasting
,
20
,
729
744
, doi:.
Brekke
,
L. D.
,
E. P.
Maurer
,
J. D.
Anderson
,
M. D.
Dettinger
,
E. S.
Townsley
,
A.
Harrison
, and
T.
Pruitt
,
2009
:
Assessing reservoir operations risk under climate change
.
Water Resour. Res.
,
45
,
W04411
, doi:.
Brown
,
R. D.
,
2000
:
Northern Hemisphere snow cover variability and change, 1915–97
.
J. Climate
,
13
,
2339
2355
, doi:.
Brown
,
R. D.
, and
P. W.
Mote
,
2009
:
The response of Northern Hemisphere snow cover to a changing climate
.
J. Climate
,
22
,
2124
2145
, doi:.
Brutel-Vuilmet
,
C.
,
M.
Ménégoz
, and
G.
Krinner
,
2013
:
An analysis of present and future seasonal Northern Hemisphere land snow cover simulated by CMIP5 coupled climate models
.
Cryosphere
,
7
,
67
80
, doi:.
Burnash
,
R. J. C.
,
R. L.
Ferral
, and
R. A.
McQuire
,
1973
: A generalized streamflow simulation system—Conceptual modeling for digital computers. Joint Federal and State River Forecast Center, U.S. National Weather Service and California Department of Water Resources, Sacramento, California, 204 pp.
Burnett
,
A. W.
,
M. E.
Kirby
,
H. T.
Mullins
, and
W. P.
Patterson
,
2003
:
Increasing Great Lake–effect snowfall during the twentieth century: A regional response to global warming?
J. Climate
,
16
,
3535
3542
, doi:.
Carroll
,
C.
,
2007
:
Interacting effects of climate change, landscape conversion, and harvest on carnivore populations at the range margin: Marten and lynx in the northern Appalachians
.
Conserv. Biol.
,
21
,
1092
1104
, doi:.
Changnon
,
S. A.
,
1969
: Climatology of severe winter storms in Illinois. Illinois State Water Survey, Bull. 53, State of Illinois, Department of Registration and Education, Urbana, Illinois, 45 pp.
Changnon
,
S. A.
,
2007
:
Catastrophic winter storms: An escalating problem
.
Climatic Change
,
84
,
131
139
, doi:.
Changnon
,
S. A.
,
D.
Changnon
, and
T. R.
Karl
,
2006
:
Temporal and spatial characteristics of snowstorms in the contiguous United States
.
J. Appl. Meteor. Climatol.
,
45
,
1141
1156
, doi:.
Christensen
,
J. H.
, and Coauthors
,
2007
:
Regional climate projections
. Climate Change 2007: The Physical Science Basis, S. Solomon et al., Eds., Cambridge University Press, 847–940.
CMAP
, cited
2013
: Climate adaptation guidebook for municipalities in the Chicago region. Chicago Metropolitan Agency for Planning. [Available online at https://www.cmap.illinois.gov/livability/sustainability-climate-change/climate-adaptation-toolkit.]
Cox
,
H. J.
, and
J. H.
Armington
,
1914
: The Weather and Climate of Chicago. University of Chicago Press, 375 pp.
Daly
,
C.
,
2006
:
Guidelines for assessing the suitability of spatial climate data sets
.
Int. J. Climatol.
,
26
,
707
721
, doi:.
Daly
,
C.
,
W. P.
Gibson
,
G. H.
Taylor
,
M. K.
Doggett
, and
J. I.
Smith
,
2007
:
Observer bias in daily precipitation measurements at United States cooperative network stations
.
Bull. Amer. Meteor. Soc.
,
88
,
899
912
, doi:.
Daly
,
C.
, and Coauthors
,
2008
:
Physiographically sensitive mapping of climatological temperature and precipitation across the conterminous United States
.
Int. J. Climatol.
,
28
,
2031
2064
, doi:.
Delworth
,
T. L.
, and Coauthors
,
2006
:
GFDL’s CM2 global coupled climate models. Part I: Formulation and simulation characteristics
.
J. Climate
,
19
,
643
674
, doi:.
Déqué
,
M.
,
C.
Dreveton
,
A.
Braun
, and
D.
Cariolle
,
1994
:
The ARPEGE–IFS atmosphere model: A contribution to the French community climate modelling
.
Climate Dyn.
,
10
,
249
266
, doi:.
Deser
,
C.
,
R.
Knutti
,
S.
Solomon
, and
A. S.
Phillips
,
2012a
:
Communication of the role of natural variability in future North American climate
.
Nat. Climate Change
,
2
,
775
779
, doi:.
Deser
,
C.
,
A. S.
Phillips
,
V.
Bourdette
, and
H.
Teng
,
2012b
:
Uncertainty in climate change projections: The role of internal variability
.
Climate Dyn.
,
38
,
527
546
, doi:.
Deser
,
C.
,
A. S.
Phillips
,
M. A.
Alexander
, and
B. V.
Smoliak
,
2014
:
Projecting North American climate over the next 50 years: Uncertainty due to internal variability
.
J. Climate
,
27
,
2271
2296
, doi:.
Draper
,
S. E.
, and
J. E.
Kundell
,
2007
:
Impact of climate change on transboundary water sharing
.
J. Water Resour. Plann. Manage.
,
133
,
405
415
, doi:.
Edwards
,
R. Y.
,
1956
:
Snow depths and ungulate abundance in the mountains of western Canada
.
J. Wildl. Manage.
,
20
,
159
168
, doi:.
Eisenberg
,
D.
, and
K. E.
Warner
,
2005
:
Effects of snowfalls on motor vehicle collisions, injuries, and fatalities
.
Amer. J. Public Health
,
95
,
120
124
, doi:.
Flato
,
G. M.
,
G. J.
Boer
,
W. G.
Lee
,
N. A.
McFarlane
,
D.
Ramsden
,
M. C.
Reader
, and
A. J.
Weaver
,
2000
:
The Canadian Centre for Climate Modelling and Analysis global coupled model and its climate
.
Climate Dyn.
,
16
,
451
467
, doi:.
Franz
,
K. J.
,
T. S.
Hogue
, and
S.
Sorooshian
,
2008a
:
Snow model verification using ensemble prediction and operational benchmarks
.
J. Hydrometeor.
,
9
,
1402
1415
, doi:.
Franz
,
K. J.
,
T. S.
Hogue
, and
S.
Sorooshian
,
2008b
:
Operational snow modeling: Addressing the challenges of an energy balance model for National Weather Service forecasts
.
J. Hydrol.
,
360
,
48
66
, doi:.
Gordon
,
H. B.
, and
S. P.
O’Farrell
,
1997
:
Transient climate change in the CSIRO coupled model with dynamic sea ice
.
Mon. Wea. Rev.
,
125
,
875
907
, doi:.
Gordon
,
H. B.
, and Coauthors
,
2002
: The CSIRO Mk3 Climate System Model. CSIRO Atmospheric Research Tech. Rep. 60, 130 pp.
Groisman
,
P. Ya.
,
R. W.
Knight
,
T. R.
Karl
, and
R. R.
Heim
Jr.
,
1993
: Inferences of the North American snowfall and snow cover with recent global temperature changes. Proc. Snow Watch ’92 Int. Workshop, Niagara-on-the-Lake, ON, Canada, Canadian Climate Centre, 44–51.
Groisman
,
P. Ya.
,
R. W.
Knight
, and
T. R.
Karl
,
2001
:
Heavy precipitation and high streamflow in the contiguous United States: Trends in the twentieth century
.
Bull. Amer. Meteor. Soc.
,
82
,
219
246
, doi:.
Gutowski
,
W. J.
, and Coauthors
,
2008
: Causes of observed changes in extremes and projections of future changes. Weather and Climate Extremes in a Changing Climate—Regions of Focus: North America, Hawaii, Caribbean, and U.S. Pacific Islands, T. R. Karl et al., Eds., U.S. Climate Change Science Program and the Subcommittee on Global Change Research, 81–116.
Hanson
,
R. T.
,
L. E.
Flint
,
A. L.
Flint
,
M. D.
Dettinger
,
C. C.
Faunt
,
D.
Cayan
, and
W.
Schmid
,
2012
:
A method for physically based model analysis of conjunctive use in response to potential climate changes
.
Water Resour. Res.
,
48
,
W00L08
, doi:.
Held
,
I. M.
, and
B. J.
Soden
,
2000
:
Water vapor feedback and global warming
.
Annu. Rev. Energy Environ.
,
25
,
441
475
, doi:.
Hoving
,
C. L.
,
D. J.
Harrison
,
W. B.
Krohn
,
R. A.
Joseph
, and
M.
O’Brien
,
2005
:
Broad-scale predictors of Canada lynx occurrence in eastern North America
.
J. Wildl. Manage.
,
69
,
739
751
, doi:.
Kalnay
,
E.
, and Coauthors
,
1996
:
The NCEP/NCAR 40-Year Reanalysis Project
.
Bull. Amer. Meteor. Soc.
,
77
,
437
471
, doi:.
Kapnick
,
S.
, and
T.
Delworth
,
2013
:
Controls of global snow under a changed climate
.
J. Climate
,
26
,
5537
5562
, doi:.
Kitoh
,
A.
,
A.
Noda
,
Y.
Nikaidou
,
T.
Ose
, and
T.
Tokioka
,
1995
:
AMIP simulations of the MRI GCM
.
Pap. Meteor. Geophys.
,
45
,
121
148
, doi:.
Knapp
,
K. K.
,
D.
Kroeger
, and
K.
Giese
,
2000
: Mobility and safety impacts of winter storm events in a freeway environment: Final report. Center for Transportation Research and Education, Iowa State University, 76 pp.
Krasting
,
J.
,
A.
Broccoli
,
K.
Dixon
, and
J.
Lanzante
,
2013
:
Future changes in Northern Hemisphere snowfall
.
J. Climate
,
26
,
7813
7828
, doi:.
Krohn
,
W. B.
,
C. L.
Hoving
,
D. J.
Harrison
,
D. M.
Phillips
, and
H. C.
Frost
,
2006
: Martes foot-loading and snowfall patterns in eastern North America: Implications to broad-scale distributions and interactions of mesocarnivores. Martens and Fishers (Martes) in Human-Altered Environments: An International Perspective, D. J. Harrison, A. K. Fuller, and G. Proulx, Eds., Springer, 115–132.
Kunkel
,
K. E.
,
N. E.
Wescott
, and
D. A. R.
Kristovich
,
2000
: Climate change and lake-effect snow. Preparing for a Changing Climate: The Potential Consequences of Climate Variability and Change, P. J. Sousounis, and J. M. Bisanz, Eds., U.S. EPA, Office of Research and Development Global Change Research Program, 25–28.
Kutikoff
,
S. L.
,
2013
: A climatology of lake-effect snowfall and evaluation of the Cobb method for the Great Lakes region. M.S. thesis, Earth and Atmospheric Sciences, University of Nebraska, 62 pp.
Laird
,
N.
,
R.
Sobash
, and
N.
Hodas
,
2009
:
The frequency and characteristics of lake-effect precipitation events associated with the New York State Finger Lakes
.
J. Appl. Meteor. Climatol.
,
48
,
873
886
, doi:.
Leathers
,
D. J.
, and
A. W.
Ellis
,
1996
:
Synoptic mechanisms associated with snowfall increases to the lee of Lakes Erie and Ontario
.
Int. J. Climatol.
,
16
,
1117
1135
, doi:.
Legutke
,
S.
, and
R.
Voss
,
1999
: The Hamburg atmosphere–ocean coupled circulation model ECHO-G. Tech. Rep. 18, German Climate Computer Centre (DKRZ), Hamburg, 62 pp.
Lemke
,
P.
, and Coauthors
,
2007
: Observations: Changes in snow, ice and frozen ground. Climate Change 2007: The Physical Science Basis, S. Solomon et al., Eds., Cambridge University Press, 337–383.
Lettenmaier
,
D.
, and
T.
Gan
,
1990
:
Hydrologic sensitivities of the Sacramento–San Joaquin River basin, California, to global warming
.
Water Resour. Res.
,
26
,
69
86
, doi:.
López-Moreno
,
J. I.
,
S.
Goyette
,
S. M.
Vicente-Serrano
, and
M.
Beniston
,
2011
:
Effects of climate change on the intensity and frequency of heavy snowfall events in the Pyrenees
.
Climatic Change
,
105
,
489
508
, doi:.
Mahanama
,
S.
,
B.
Livneh
,
R.
Koster
,
D.
Lettenmaier
, and
R.
Reichle
,
2012
:
Soil moisture, snow, and seasonal streamflow forecasts in the United States
.
J. Hydrometeor.
,
13
,
189
203
, doi:.
Manabe
,
S.
,
R. T.
Wetherald
, and
R. J.
Stouffer
,
1981
:
Summer dryness due to an increase of atmospheric CO2 concentration
.
Climatic Change
,
3
,
347
386
, doi:.
Mastin
,
M. C.
,
K. J.
Chase
, and
R. W.
Dudley
,
2011
:
Changes in spring snowpack for selected basins in the United States for different climate-change scenarios
.
Earth Interact.
,
15
, doi:.
Maurer
,
E. P.
,
L. D.
Brekke
, and
T.
Pruitt
,
2010
:
Contrasting lumped and distributed hydrology models for estimating climate change impacts on California watersheds
.
J. Amer. Water Resour. Assoc.
,
46
,
1024
1035
, doi:.
McKenney
,
D. W.
,
M. F.
Hutchinson
,
J. L.
Kesteven
, and
L. A.
Venier
,
2001
:
Canada’s plant hardiness zones revisited using modern climate interpolation techniques
.
Can. J. Plant Sci.
,
81
,
129
143
, doi:.
Mesinger
,
F.
, and Coauthors
,
2006
:
North American Regional Reanalysis
.
Bull. Amer. Meteor. Soc.
,
87
,
343
360
, doi:.
Miller
,
N. L.
,
K. E.
Bashford
, and
E.
Strem
,
2003
:
Potential impacts of climate change on California hydrology
.
J. Amer. Water Resour. Assoc.
, 39,
771
784
, doi:.
Mishra
,
V.
,
K. A.
Cherkauer
, and
S.
Shukla
,
2010
:
Assessment of drought due to historic climate variability and projected future climate change in the midwestern United States
.
J. Hydrometeor.
,
11
,
46
68
, doi:.
National Geophysical Data Center
, cited
2013
: 2-minute gridded global relief data (ETOPO2v2). NOAA/NESDIS/NGCD. [Available online at http://www.ngdc.noaa.gov/mgg/fliers/06mgg01.html.]
Norton
,
D. C.
, and
S. J.
Bolsenga
,
1993
:
Spatiotemporal trends in lake effect and continental snowfall in the Laurentian Great Lakes, 1951–1980
.
J. Climate
,
6
,
1943
1955
, doi:.
Notaro
,
M.
,
D.
Lorenz
,
D.
Vimont
,
S.
Vavrus
,
C.
Kucharik
, and
K.
Franz
,
2011
:
21st century Wisconsin snow projections based on an operational snow model driven by statistically downscaled climate data
.
J. Climatol.
,
31
,
1615
1633
, doi:.
Notaro
,
M.
,
A.
Zarrin
,
S.
Vavrus
, and
V.
Bennington
,
2013
:
Simulation of heavy lake-effect snowstorms across the Great Lakes basin by RegCM4: Synoptic climatology and variability
.
Mon. Wea. Rev.
,
141
,
1990
2014
, doi:.
O’Gorman
,
P. A.
, and
T.
Schneider
,
2009
:
The physical basis for increases in precipitation extremes in simulations of 21st-century climate change
.
Proc. Natl. Acad. Sci. USA
,
106
,
14 773
14 777
, doi:.
Ohmura
,
A.
,
2001
:
Physical basis for the temperature-based melt-index method
.
J. Appl. Meteor.
,
40
,
753
761
, doi:.
Peacock
,
S.
,
2012
:
Projected twenty-first-century changes in temperature, precipitation, and snow cover over North America in CCSM4
.
J. Climate
,
25
,
4405
4429
, doi:.
Räisänen
,
J.
,
2008
:
Warmer climate: Less or more snow?
Climate Dyn.
,
30
,
307
319
, doi:.
Raleigh
,
M. S.
, and
J. D.
Lundquist
,
2012
:
Comparing and combining SWE estimates from the SNOW-17 model using PRISM and SWE reconstruction
.
Water Resour. Res.
,
48
,
W01506
, doi:.
Rasmussen
,
R. M.
,
A.
Crook
, and
C.
Kessinger
,
1993
:
Snow-band formation and evolution during the 15 November 1987 aircraft accident at Denver airport
.
Wea. Forecasting
,
8
,
453
480
, doi:.
Robinson
,
P. J.
,
1989
:
The influence of weather on flight operations at the Atlanta Hartsfield International Airport
.
Wea. Forecasting
,
4
,
461
468
, doi:.
Roeckner
,
E.
, and Coauthors
,
2003
: The atmospheric general circulation model ECHAM5. Part I: Model description. Max Planck Institute for Meteorology Rep. 349, 127 pp.
Roesch
,
A.
,
2006
:
Evaluation of surface albedo and snow cover in AR4 coupled climate models
.
J. Geophys. Res.
,
111
,
D15111
, doi:.
Rogot
,
E.
, and
S. J.
Padgett
,
1976
:
Associations of coronary and stroke mortality with temperature and snowfall in selected areas of the United States, 1962–1966
.
Amer. J. Epidemiol.
,
103
,
565
575
.
Rooney
,
J. F.
,
1967
:
Urban snow hazard in the United States—Appraisal of disruption
.
Geogr. Rev.
,
57
,
538
559
, doi:.
Schmidt
,
G. A.
, and Coauthors
,
2006
:
Present-day simulations using GISS ModelE: Comparison to in-situ, satellite, and reanalysis data
.
J. Climate
,
19
,
153
192
, doi:.
Schummer
,
M. L.
,
R. M.
Kaminski
,
A. H.
Raedeke
, and
D. A.
Graber
,
2010
:
Weather-related indices of autumn–winter dabbling duck abundance in middle North America
.
J. Wildl. Manage.
,
74
,
94
101
, doi:.
Scott
,
D.
,
J.
Dawson
, and
B.
Jones
,
2008
:
Climate change vulnerability of the US Northeast winter recreation-tourism sector
.
Mitig. Adapt. Strategies Global Change
,
13
,
577
596
, doi:.
Severinghaus
,
C. W.
,
1947
: Relationship of weather to winter mortality and population levels among deer in the Adirondack region of New York. Trans. N. Amer. Wildl. Conf.,12, 212–223.
Southwick Associates
,
2012
: Hunting in America: An economic force for conservation. Southwick Associates, produced for National Shooting Sports Foundation in partnership with the Association of Fish and Wildlife Agencies, 16 pp. [Available online at http://www.nssf.org/PDF/research/HuntingInAmerica_EconomicForceForConservation.pdf.]
Stacy
,
E. W.
,
1962
:
A generalization of the gamma distribution
.
Ann. Math. Stat.
,
33
,
1187
1192
, doi:.
U.S. Fish and Wildlife Service
, cited
2011
: National Survey of Fishing, Hunting, and Wildlife-Associated Recreation. U.S. Department of the Interior, U.S. Fish and Wildlife Service, and U.S. Department of Commerce, U.S. Census Bureau. [Available online at http://www.census.gov/prod/www/fishing.html.]
Verme
,
L. J.
,
1968
:
An index of winter weather severity for northern deer
.
J. Wildl. Manage.
,
32
,
566
574
, doi:.
Vicuna
,
S.
,
R.
Leonardson
,
M.
Hanemann
,
L.
Dale
, and
J.
Dracup
,
2008
:
Climate change impacts on high elevation hydropower generation in California’s Sierra Nevada: A case study on the Upper American River
.
Climatic Change
,
87
,
123
137
, doi:.
Westerling
,
A.
,
H.
Hidalgo
,
D.
Cayan
, and
T.
Swetnam
,
2006
:
Warming and earlier spring increases western U.S. forest wildfire activity
.
Science
,
313
,
940
943
, doi:.
WICCI
, cited
2011
: Wisconsin’s changing climate: Impacts and adaptation. Wisconsin Initiative on Climate Change Impacts. [Available online at http://www.wicci.wisc.edu/.]
Wood
,
A. W.
,
L. R.
Leung
,
V.
Srifhar
, and
D. P.
Lettenmaier
,
2004
:
Hydrologic implications of dynamical and statistical approaches to downscaling climate model outputs
.
Climatic Change
,
62
,
189
216
, doi:.
Zappa
,
M.
,
F.
Pos
,
U.
Strasser
,
P.
Warmerdam
, and
J.
Gurtz
,
2003
:
Seasonal water balance of an Alpine catchment as evaluated by different methods for spatially distributed snowmelt modeling
.
Nord. Hydrol.
,
34
,
179
202
.
Zhang
,
Y.
,
S.
Wang
,
A. G.
Barr
, and
T. A.
Black
,
2008
:
Impact of snow cover on soil temperature and its simulation in a boreal aspen forest
.
Cold Reg. Sci. Technol.
,
52
,
355
370
, doi:.

Footnotes

*

Nelson Institute Center for Climatic Research Contribution 1183.

+

Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JCLI-D-13-00520.s1.

#

Current affiliation: Department of Biological Sciences, State University of New York at Oswego, Oswego, New York.

Supplemental Material