Arctic snow presence, absence, properties, and water amount are key components of Earth’s changing climate system that incur far-reaching physical and biological ramifications. Recent dataset and modeling developments permit relatively high-resolution (10-km horizontal grid; 3-h time step) pan-Arctic snow estimates for 1979–2009. Using MicroMet and SnowModel in conjunction with land cover, topography, and 30 years of the NASA Modern-Era Retrospective Analysis for Research and Applications (MERRA) atmospheric reanalysis data, a distributed snow-related dataset was created including air temperature, snow precipitation, snow-season timing and length, maximum snow water equivalent (SWE) depth, average snow density, snow sublimation, and rain-on-snow events. Regional variability is a dominant feature of the modeled snow-property trends. Both positive and negative regional trends are distributed throughout the pan-Arctic domain, featuring, for example, spatially distinct areas of increasing and decreasing SWE or snow season length. In spite of strong regional variability, the data clearly show a general snow decrease throughout the Arctic: maximum winter SWE has decreased, snow-cover onset is later, the snow-free date in spring is earlier, and snow-cover duration has decreased. The domain-averaged air temperature trend when snow was on the ground was 0.17°C decade−1 with minimum and maximum regional trends of −0.55° and 0.78°C decade−1, respectively. The trends for total number of snow days in a year averaged −2.49 days decade−1 with minimum and maximum regional trends of −17.21 and 7.19 days decade−1, respectively. The average trend for peak SWE in a snow season was −0.17 cm decade−1 with minimum and maximum regional trends of −2.50 and 5.70 cm decade−1, respectively.
Ample evidence indicates the Arctic is changing. Available long-term temperature observations show warming trends of variable strength throughout the Arctic (Serreze et al. 2000; Overland et al. 2004; Chapin et al. 2005; Hinzman et al. 2005; White et al. 2007) and model simulations of future scenarios point to a warmer Arctic overall, especially with continued summer sea ice loss (Solomon et al. 2007; Holland et al. 2010). Arctic permafrost temperatures, monitored from boreholes 10–20 m deep, increased 2°C in the last 20–30 yr (Romanovsky et al. 2010). The Arctic Ocean’s summer sea ice extent continues to shrink and the coverage and thickness of multiyear ice is in marked decline (Serreze and Francis 2006; Stroeve et al. 2007; Gerland et al. 2008; Holland et al. 2010).
Arctic terrestrial precipitation trends are inherently more difficult to detect given snowfall measurement challenges with sparsely distributed observations, rare long-term records, chronic station record discontinuities, variable gauge designs, low precipitation amounts, and high winds (Adam and Lettenmaier 2003; Yang et al. 2005; Bonsal and Kochtubajda 2009; Turner and Overland 2009). These factors make point observations of snowfall especially problematic in terms of broader representation and for identification of long-term trends, yet valiant attempts in identifying trends have been made nonetheless. New et al. (2001) reported a 0.32 cm decade−1 increase in precipitation (from 1901 to 1998) using station data from 60° to 80°N latitude. Hinzman et al. (2005) highlighted a handful of mostly nonsignificant and slight increases or decreases in long-term Arctic precipitation records; two stations with significant trends registered precipitation increases. Unfortunately, given the uncertainties and spatial variability of snowfall, the distribution of snow gauge observations does not lend itself to coarser-scale extrapolations of precipitation trends beyond isolated landscapes of intensive study.
More frequently, snow-cover data (mostly satellite visible data) are used to identify changes in the arrival and longevity of terrestrial snow (Frei and Robinson 1999; Serreze et al. 2000; Dye 2002; Stone et al. 2002; Brown et al. 2007; Brown and Mote 2009; Zhao and Fernandes 2009; Brown et al. 2010; Choi et al. 2010). With terrestrial snow, the presence or absence frequently serves as a surrogate measure of snow and cryosphere change, since visible wavelengths cannot be used to estimate the amount of water present. These approaches are particularly adept at capturing changes in snow arrival, departure, and duration. Unfortunately, changes in snowmelt processes and their important energy feedback consequences (Chapin et al. 2005; Euskirchen et al. 2009), and changes from years to decades in snow water equivalent (SWE) depth distributions, are still lacking. While substantial strides are being made in SWE algorithms using passive microwave data (Armstrong and Brodzik 2001; Wulder et al. 2007; Derksen et al. 2010), uncertainties with algorithm applications, snow properties, land cover, and coarse-scale measurement issues persist, limiting the confidence and applicability of these data for trend analyses.
An alternative method of estimating Arctic snow properties (e.g., cover, SWE, duration) and their changes over time is through models. General circulation models (GCMs) have been used to address Arctic climate change and precipitation questions for past and future conditions (Räisänen 2008; Walsh et al. 2008; Finnis et al. 2009a,b); these studies generally find higher temperatures lead to increases in Arctic precipitation. Efforts have also been made to link land surface hydrology models to reanalyses or GCM data, to specifically address Arctic (Pohl et al. 2007; Slater et al. 2007) or relatively high-resolution (0.5°) global (Adam et al. 2009) snow processes.
Inasmuch as these coarse-scale modeling approaches can answer critical climate-related questions, one substantial deficiency is their resolution. In most cases, the scales at which GCMs and other global models operate (0.5°–2.5°; Shukla et al. 2010) are too coarse to capture key snow processes, heterogeneities, and land-cover and snow interactions (Liston 2004; Liston and Hiemstra 2011). Modeling capabilities have grown, along with improvements in computing power and the emergence of relatively high-resolution topographic, land cover, and meteorological data products. Together, these tools and datasets can be combined to provide a reasonable facsimile of cryospheric processes, and allow improved understanding of the specific implications of climate changes related to snow. In addition, this can now be done at much higher resolution than previously possible.
The purpose of this paper is to perform and analyze the spatial and temporal evolution of snow, snow processes, and snow characteristics in high northern latitudes (north of approximately 55°N), at the highest achievable spatial and temporal resolution; all in an effort to understand regional variations in key climate-relevant snow-related features. This is accomplished by driving a local- to regional-scale meteorological and snow-evolution modeling system with 3-hourly, ⅔° longitude by ½° latitude gridded atmospheric reanalysis datasets. The resulting snow-related modeling and analysis datasets span 30 yr (1 August 1979–31 July 2009), covering a pan-Arctic domain with a 10-km grid increment and 3-h time step.
2. Model description
To quantify spatial and temporal variations in Arctic-system snow properties and characteristics, we performed model simulations using SnowModel (Liston and Elder 2006a), a spatially distributed snow-evolution modeling system designed for application in all landscapes, climates, and conditions where snow occurs. It is an aggregation of four submodels: EnBal (Liston 1995; Liston et al. 1999) calculates surface energy exchanges and snowmelt; SnowPack (Liston and Hall 1995; G. E. Liston and S. H. Mernild 2011, unpublished manuscript) simulates snow depth and water-equivalent evolution; SnowTran-3D (Liston and Sturm 1998; Liston et al. 2007) accounts for snow redistribution by wind; and SnowAssim (Liston and Hiemstra 2008) is available to assimilate field and remote sensing datasets.
SnowModel is designed to run on grid increments of 1–200 m and temporal increments of 10 min to 1 day. It can be applied using much larger grid increments (up to tens of kilometers) if the inherent loss in high-resolution (subgrid) information (Liston 2004) is acceptable. Processes simulated by SnowModel include snow precipitation; blowing-snow redistribution and sublimation; interception, unloading, and sublimation within forest canopies; snow-density evolution; and snowpack ripening and melt. SnowModel incorporates first-order physics required to simulate snow evolution within each of the global snow classes [i.e., ice, tundra, taiga, warm forest (or alpine), prairie, maritime, and ephemeral] defined by Sturm et al. (1995) and G. E. Liston and M. Sturm (2011, unpublished manuscript). Required SnowModel inputs include temporally variant precipitation, wind speed and direction, air temperature, and relative humidity obtained from meteorological stations and/or an atmospheric model located within or near the simulation domain. Spatially distributed, time-invariant topography and land cover are also necessary.
Meteorological forcings required by SnowModel are provided by MicroMet (Liston and Elder 2006b), a quasi–physically based, high-resolution (e.g., 1-m to 10-km horizontal grid increment), meteorological distribution model. MicroMet is a data assimilation and interpolation model that utilizes meteorological station datasets and/or gridded atmospheric model or (re)analyses datasets. MicroMet minimally requires near-surface air temperature, relative humidity, wind speed and direction, and precipitation data. The model uses known relationships among meteorological variables and the surrounding landscape (primarily topography) to distribute those variables over any given landscape in physically plausible and computationally efficient ways (Liston and Elder 2006b). MicroMet performs two kinds of adjustments to the meteorological data: 1) all available data fields, at a given time, are spatially interpolated over the domain; and 2) physically based submodels are applied to each MicroMet variable to quantify topographic, elevation, and vegetation effects at any given point in space and time. At each time step, MicroMet simulates and distributes air temperature, relative humidity, wind speed, wind direction, incoming solar radiation, incoming longwave radiation, surface pressure, and precipitation, and makes them accessible to SnowModel.
MicroMet and SnowModel constitute a physically based modeling system that creates value-added snow information (e.g., snow depth, snow density, snowmelt rate, snow thermal properties, snow-cover duration, and sublimation) from basic meteorological variables (e.g., air temperature, humidity, precipitation, and wind). The products yielded are based on our physical understanding of snow-evolution processes and features, and their interactions with the atmosphere and surrounding land surface. MicroMet and SnowModel have been used to distribute observed and modeled meteorological variables and evolve snow distributions over complex terrain in Colorado, Wyoming, Idaho, Oregon, Alaska, Arctic Canada, Siberia, Japan, Tibet, Chile, Germany, Austria, Norway, Greenland, and Antarctica as part of a wide variety of terrestrial modeling studies (e.g., Liston and Sturm 1998, 2002; Greene et al. 1999; Liston et al. 2000, 2002, 2007, 2008; Prasad et al. 2001; Hiemstra et al. 2002, 2006; Hasholt et al. 2003; Bruland et al. 2004; Liston and Winther 2005; Mernild et al. 2006, 2008, 2009, 2010; Liston and Hiemstra 2008, 2011; Mernild and Liston 2010).
3. Model simulation
a. Model configuration and simulation domain
Snow evolution and surface energy fluxes were simulated using MicroMet and SnowModel for the 30-yr period, 1 August 1979–31 July 2009. The simulation covered a 7250 km by 7250 km domain centered on the North Pole (Fig. 1). This domain encompasses the majority of the Arctic system, defined to be the northern region of Earth where energy and moisture interact with midlatitudes (Roberts et al. 2010). The simulation domain incorporates many of the common definitions of the terrestrial Arctic system: land north of the Arctic Circle; the majority of land north of the July 10°C surface air temperature isotherm and the annual-average 0°C surface air temperature isotherm that circle the North Pole; the southern boundary of land draining into the northern high-latitude oceans (the simulation domain does not quite reach this in a couple locations); and land north of the northern treeline and the tundra-taiga ecotone. The model simulation was performed using a 10-km horizontal grid increment (525 625 grid cells) and 3-h time step. This relatively fine grid increment allows improved representation of snow-evolution processes associated with topographic and land-cover variations. Most (86%) of the domain is above 55°N latitude, but the corners extend as far south as 43°N latitude. Because blowing snow does not typically move across 10-km grid cells into adjacent cells (in the natural system, wind-transported snow particles are typically either captured in a topographic drift trap, captured in vegetation protruding above the snow surface, or they sublimate away completely before they travel 10 km), SnowTran-3D and the associated snow-transport processes, including blowing snow sublimation, were not included in the simulation. The lack of quality pan-Arctic snow data also precluded application of SnowAssim in these simulations.
Topographic data used in the model simulation were obtained from the National Oceanic and Atmospheric Administration (NOAA) Global Land One-km Base Elevation Project (GLOBE; Hastings et al. 1999), which provided 1-km digital elevation model (DEM) data that were resampled to 10 km. The land cover distribution used in the simulation was a hybrid dataset created primarily from 300-m Global Land Cover (GlobCover; more information available online at http://ionia1.esrin.esa.int/) data augmented with the Circumpolar Arctic Vegetation Map (CAVM; CAVM Team 2003). The CAVM was utilized to correct GlobCover’s misclassified snow/ice (areas that should have been barren or rock) at high northern latitudes (>82°N) in northern Canada and Greenland. The resulting hybrid dataset was resampled to 10 km and reclassified into SnowModel land cover classes (Liston and Elder 2006a).
Topography of the study area ranges from sea level to over 5000 m, and land cover includes bare, wetlands, tundra, shrubs, deciduous and coniferous trees, glaciers, and ice sheets (Fig. 1). The 30-arc-s (~1 km) resolution version (G. E. Liston and M. Sturm 2011, unpublished manuscript) of Sturm et al.’s (1995) global snow classification, regridded to the pan-Arctic SnowModel simulation grid (Fig. 2), shows the land area of this domain contains 10% ice (i.e., glaciers and ice sheets), 39% tundra, 34% taiga, 7% warm forests (alpine in Sturm et al. 1995), 9% prairie, and 1% maritime snow classes. These snow classes take into account the wind, precipitation, and temperature regimes these snow covers evolve within, combining to yield unique combinations of typical depths, densities, layering, crystal morphology, and grain characteristics.
Climatologically, air temperatures over the entire domain are typically well below 0°C for much of the fall, winter, and spring months, with the northernmost regions spending considerable time in winter darkness. Near-surface temperature inversions are common throughout the snow-covered Arctic (e.g., Mernild and Liston 2010). The associated thermal stability inhibits vertical mixing and produces variable local and regional climates. This, in combination with local and regional terrain influences, can produce local and regional meteorological and snow conditions that are much more complex than coarse-scale patterns (e.g., Lynch et al. 2001; Liston and Sturm 2002; Taras et al. 2002). The fall, winter, and spring snow seasons in this domain can range from less than 30 days to over 300 days each year. Snow can begin accumulating as early as 1 September, and it can be 1 July before it melts completely.
b. Meteorological forcing
Atmospheric forcing data were provided by the National Aeronautics and Space Administration (NASA) Modern Era Retrospective-Analysis for Research and Applications (MERRA) products (Bosilovich 2008; Bosilovich et al. 2008, 2011; Cullather and Bosilovich 2011; Rienecker et al. 2011; Robertson et al. 2011). This reanalysis program has the specific goal of improving the representation of water cycle processes and features within the analyses while taking advantage of modern satellite era datasets. MERRA covers the period 1979–present, on a ⅔° longitude by ½° latitude global grid. Surface atmospheric forcing variables are available hourly. The analysis assimilates a wide range of satellite observations in addition to more conventional radiosonde, dropsonde, aircraft, and surface observations. Bosilovich et al. (2008) analyzed precipitation outputs from an early version of the MERRA reanalysis system, and concluded the MERRA precipitation fields were an improvement over the previous generations of reanalyses.
In preparation for the model simulation, hourly MERRA 10-m air temperature, specific humidity, and u and υ wind components, and surface pressure and precipitation variables were aggregated to 3-hourly values. This was done to improve computational efficiency while still resolving the diurnal cycle and the associated energy-related processes. MicroMet then used these to create the 3-hourly, 10-km atmospheric forcing distributions required by SnowModel (air temperature, relative humidity, wind speed and direction, precipitation, and incoming solar and longwave radiation); see Liston and Elder (2006b) for additional details. Water-equivalent precipitation was provided from MERRA, and MicroMet’s temperature threshold parameterization was used to define whether rain or snow fell on each model grid cell. The 10-km atmospheric fields were ingested by SnowModel to simulate the time evolution and spatial distribution of water and energy fluxes and states. Simulated variables included the following: surface (skin) temperature, albedo, outgoing longwave radiation, latent heat flux, sensible heat flux, liquid precipitation, solid precipitation, snowmelt, sublimation, snowmelt runoff, snow depth, snow density, and snow water equivalent. In addition, we generated secondary products such as the timing and distribution of rain-on-snow events, changes in snow and growing season lengths, hydrologic budgets, winter soil microbial activity, changes in snow thermal characteristics, and changes in surface energy exchanges.
c. Verification datasets
Three independent datasets that cover similar spatial and/or temporal domains to the simulation presented herein are available for model verification. Unfortunately, these datasets do not possess an identical spatial or temporal resolution to the SnowModel products, but they remain useful for comparisons and identifying large differences. First, NOAA data describing snow-free or snow-covered land conditions are available weekly over the Northern Hemisphere on a polar grid having increments ranging in size from 125 to 200 km (Robinson et al. 1993; Choi et al. 2010). These datasets are produced by trained observers interpreting visible-band satellite imagery. They suffer from difficulties distinguishing between clouds and snow, snow obscured beneath forest canopies, the inability to identify subpixel snow features, and they require daylight. The somewhat subjective nature of this dataset makes trend analyses more uncertain, but statistically significant trends have been detected (Choi et al. 2010). These data were obtained from the National Snow and Ice Data Center (NSIDC) for the period 1979–2007, on the 25-km by 25-km Equal-Area Scalable Earth Grid (EASE) Grid (Armstrong and Brodzik 2007).
Second, the U.S. Air Force Environmental Technology Application Center (USAF/ETAC) developed a monthly snow-depth climatology (one depth distribution for each of the 12 months) on a global, 1° longitude by 1° latitude grid (Foster and Davy 1988). These distributions were created by gathering snow depth data from many worldwide sources, plotting them on monthly hemispheric maps, and manually analyzing them based on a confidence level assigned to each input dataset. See Brown and Frei (2007) for a discussion of the numerous dataset limitations. This dataset was obtained from the NASA Goddard Distributed Active Archive Center (DAAC).
Finally, global, monthly satellite-derived SWE data for 1979–2007 were obtained from NSIDC. These data are also available on the 25-km EASE-Grid, and are derived from Scanning Multichannel Microwave Radiometer (SMMR) and selected Special Sensor Microwave Imagers (SSM/I; Armstrong et al. 2007). See König et al. (2001) for a summary of the various limitations of these sensors for measuring snow, and the shortcomings of the associated SWE products.
Other potential verification datasets are also less than ideal. Numerous GCM simulations have output snow data (e.g., Frei et al. 2003; Meehl et al. 2007). For example, the Meehl et al. (2007) multimodel dataset includes 15 GCMs that output SWE and/or snow cover (Brown and Mote 2009). Unfortunately, the spatial resolution of these models is typically greater than 1.5° longitude by 1.5° latitude, and the temporal output resolution is often monthly. Available station data have also been used to create continental- and global-scale gridded snow datasets (e.g., Brasnett 1999). Brown et al. (2003) used station data to generate monthly mean SWE and snow depth distributions over North America on a 0.3° longitude–latitude grid. Unfortunately, for the model simulations described herein, there are few stations north of 55°N latitude. At finer spatial scales, numerous snow-related field experiments exist (e.g., Sturm and Liston 2003; and the numerous Arctic snow-related field studies listed in Liston 2004; Kohler et al. 2006; Sturm et al. 2010). These datasets typically do not have the pan-Arctic distribution of interest in this study (they commonly cover less than 1% of the simulation domain, and typically span years instead of decades). Many of these have been used as part of MicroMet and SnowModel development exercises (see section 2).
4. Model results
The 3-hourly MicroMet- and SnowModel-simulated atmospheric and snow data were aggregated (averaged or summed, depending on the variable) to daily values over the 30-yr simulation period for spatial and trend analyses. In addition, for each year, averages (of variables like air temperature), sums (of variables like solid precipitation), snow-onset and snow-free dates, maximum SWE depth during the snow year, and other secondary variables were calculated. In these calculations a “snow” year was assumed to be 1 August–31 July, to allow deep, maritime snow to melt, and to simulate early-season Arctic snow precipitation and accumulation events. Because of the focus here on land-based seasonal snow evolution, any grid cells containing glaciers, ice sheets, and oceans were omitted from calculations and analyses.
With our focus on trends, we performed only limited comparisons among the SnowModel simulation and verification datasets described above. Complete SnowModel authentication studies have been made for numerous areas around the globe (see section 2), and we are confident that, assuming the MERRA forcing data are appropriate (e.g., Bosilovich et al. 2008, 2011; Cullather and Bosilovich 2011; Rienecker et al. 2011; Robertson et al. 2011), simulated snow fields will be a reasonable representation of nature. As part of our analyses, we also provide comparisons with previously published Arctic-relevant results for variables of interest. In addition, the MERRA data archive includes MERRA SWE distributions. In general, because of the strong control the MERRA precipitation forcing fields have over the SWE distributions and trends, we find the SnowModel results have added higher-resolution spatial information, and the coarser-resolution MERRA SWE fields and trend patterns are similar to those produced by SnowModel.
The 30-yr-average, 10-m air temperature for days with snow on the ground, is shown in Fig. 3a. The linear trends in annual averages for each grid cell are provided in Fig. 3b. The spatial variability in linear trends across this Arctic-system domain is considerable, with variations in magnitude and changes in sign occurring over distances of a few hundred kilometers. As general examples of positive (increasing with time) and negative (decreasing with time) temperature trends, annual values during the snow-covered period, averaged over 250 km by 250 km regions in Alaska and Canada are presented in Figs. 3c,d, respectively, along with the associated linear trend lines. These two regions were chosen because they typically represent opposite trends for the analyzed variables.
The domain-average linear trend for air temperature with snow on the ground was 0.17°C decade−1 (p < 0.10, where p is the level of significance; Fig. 3; Table 1). In addition, a region was defined to be a 250 km by 250 km area that was free of ice or ocean grid cells (such as the boxes in Figs. 3a,b). To find the region within the simulation domain that had the minimum average trend, and the region that had the maximum average trend, a 25 by 25 gridcell “box” was swept over the simulation domain, incrementing the box position by one grid cell in x and/or y, while calculating the average trend over each box position. Box positions that included ice and/or ocean grid cells were rejected, and the minimum and maximum regional (i.e., box) trends were identified and saved. For air temperature with snow on the ground, minimum (−0.55°C decade−1; p < 0.05) and maximum (0.78°C decade−1; p < 0.01) regional trends over the simulation domain were calculated (Table 1). These regions corresponded to the minimum and maximum colors/values and patterns shown in Fig. 3b (exact positions not shown).
Positive trends of this variable covered 73% of the simulation domain (Fig. 3b; Table 1). At the coarsest scale, the overall domain-averaged temperature trend was small since it comprised larger-magnitude temperature shifts in contrasting directions occurring across the Arctic system. Also note that because the snow-covered season changed throughout the simulation, the time period over which the temperature averaging occurred also varied from one year to the next. Because of this, snow-on-the-ground air temperature trends do not just reflect changes in air temperature; they also include information regarding changes in snow-season timing and length. The domain-average, annual air temperature trend (i.e., with and without snow on the ground) was 0.38°C decade−1 (p < 0.01), and positive trends covered 99% of the domain (Fig. 4; Table 1).
In addition to air temperature, precipitation is a key climate system variable. Figure 5a displays the 30-yr-average annual solid precipitation for each grid cell. The areas with dry, continental climates typically have annual water-equivalent precipitation totals of 10–30 cm, while maritime, coastal climates can have annual solid precipitation amounts in excess of 75 cm. The linear trends in solid precipitation are plotted in Fig. 5b. Again, considerable spatial variability, and both positive and negative trends are found (Figs. 5b–d). When averaged over the simulation domain, the solid precipitation linear trend almost goes to 0 (−0.02 cm decade−1; not statistically significant), while the region minimum (−3.03 cm decade−1; p < 0.01) and maximum (8.00 cm decade−1; p < 0.01) were of considerably greater magnitudes (Table 1). Negative solid precipitation trends covered 64% of the domain (Fig. 5b; Table 1). The fraction of annual precipitation falling as snow, and the associated trends, are shown in Fig. 6 and Table 1. Averaged over the temporal and spatial domains (Fig. 6a), 47% of the precipitation fell as snow. The domain averaged linear trend was −0.78% decade−1 (p < 0.01), with a regional minimum of −5.41% decade−1 (p < 0.01) and a regional maximum of 2.53% decade−1 (p < 0.10; Table 1). The higher latitudes of the simulation domain showed a clear and widespread decrease in the amount of snowfall relative to total yearly precipitation (Fig. 6b).
Given the important role snow cover plays in Arctic surface energy and moisture budgets, and other aspects of the Arctic climate system (e.g., Serreze and Barry 2005; McGuire et al. 2006; Euskirchen et al. 2007), quantifying changes and variations in snow-cover duration, timing, and spatial patterns are essential for a comprehensive understanding of high-latitude climate changes. The core snow season was defined to be the longest continuous period (in days) with snow cover for each year, for each grid cell (Fig. 7a, inset). Using the model datasets and this definition, the core season snow-onset date (typically in the fall, but could be in winter), core snow-free date (typically in the spring, but could be winter or summer), length of core snow season, date of first snow accumulation (snow that lasted at least 24 h), date of last snow on the ground, and total number of days with snow on the ground were calculated. Over the 30-yr period, this domain averaged 219 ±50 days with snow on the ground, 214 ±52 days of which were during the core snow season (Fig. 7a). The 28-yr (1979–2007), average number of days with snow on the ground from the NOAA snow-cover dataset is given in Fig. 8. Even with the much coarser resolution, the general patterns and magnitudes are similar to the SnowModel simulation, with a domain average of 215 ±57 days with snow on the ground. The distribution of 30-yr linear trends are shown in Figs. 7b–d. Throughout much of the Arctic, there was a decrease in the snow-cover duration, with a regional peak of −17.0 days decade−1 (p < 0.01) (negative trends covered 75% of the domain) and a domain average of −2.6 days decade−1 (p < 0.01), but there were also regions of increased snow duration, with a maximum of 8.1 days decade−1 (p < 0.10; Table 1). The trends in total number of days with snow on the ground (both the red and blue periods shown in Fig. 7a, inset) were similar to those of the core snow period (Table 1). The Alaska region, shown in Fig. 7c, had a 15-day snow-cover decrease over the 30-yr period. Trends in snow-cover duration are strongly controlled by the combination of snow precipitation inputs (Fig. 5), air temperature (Fig. 3), and solar-radiation-related ablation processes.
The snow-onset date in the fall (Fig. 9a) typically occurred later in the year over the 30-yr period (Fig. 9b), with a domain average of 1.3 days decade−1 (p < 0.01; Table 1). Regional extremes occurred that ranged from a decrease of 10.8 days decade−1 (p < 0.01) to an increase of 14.1 days decade−1 (p < 0.01; Table 1). Positive trends (snow onset later in the year) of this variable covered 65% of the domain (Table 1). The snow-free date in the spring (Fig. 10a) shifted to be earlier in the year (Fig. 10b), with a domain average of −1.3 days decade−1 (p < 0.01; Table 1). Regional extremes ranged from −9.9 days decade−1 (p < 0.10) to 3.7 days decade−1 (p < 0.05), with negative trends covering 80% of the domain (Figs. 10b–d; Table 1).
From a regional hydrologic perspective, moisture contained within the winter snowpack represents water storage that is made available to the climate system as a liquid when it melts in the spring. Throughout the Arctic, this spring melt is typically the largest single hydrologic event of the year, leading to a snowmelt discharge hydrograph that contains as much as 80% of the total annual runoff from many Arctic drainage basins (e.g., McNamara et al. 1998; Yang et al. 2003, 2004). This moisture storage is captured by the peak SWE during the snow season (Fig. 11a). For comparison, Fig. 12a presents the USAF/ETAC peak SWE distribution. To generate this display, the original USAF/ETAC snow-depth climatology data were converted to SWE using the Sturm et al. (1995) and G. E. Liston and M. Sturm (2011, unpublished manuscript) snow classification (Fig. 2), while using the snow densities in each class and each month defined by Sturm et al. (2010). Then, for each grid cell, the maximum SWE value was selected out of the 12 months of data. Figure 12b presents the 28-yr-average (1979–2007) SWE from the SMMR and SSM/I datasets. A comparison of these two datasets, with the SnowModel simulation, reveals obvious remote sensing data limitations with relatively warm coastal accumulations (e.g., the mountains of western British Columbia and Norway), and lack of spatial structure associated with land cover, continentality, and minor mountain ranges. Trends in SWE (Figs. 11b–d) are associated with the corresponding trends in snow precipitation (Figs. 5b–d). Negative trends in peak SWE covered 61% of the simulation domain (Table 1).
Snow density evolves as a function of snowpack temperature (and the associated air temperature and surface energy balance), wind speed (wind breaks up the snow particles and allows them to pack more tightly together), compaction (due to snow overburden pressure), and temperature and vapor pressure gradients within the snowpack (e.g., Liston et al. 2007). Changes in snow density represent an integrated measure of the “snow climate” the snowpack evolves within. Because of the relatively large gridcell increment used in this model simulation, blowing snow processes were not included in the simulations. As a consequence, we estimate snow density values simulated by the model (Fig. 13a) in the nonforested areas of the simulation domain (Fig. 1b) are approximately 50 kg m−3 lower than those found in nature, based on observations presented in Sturm et al. (2010). The snow density trends (Figs. 13b–d) averaged over the domain are slightly positive (Table 1), while regional trends had minimums and maximums as large as −21.0 kg m−3 decade−1 (p < 0.01) and 15.4 kg m−3 decade−1 (p < 0.01), respectively (Table 1).
Total snow-season, static-surface sublimation is presented in Fig. 14. In the natural system, sublimation can occur from the static snow surface, from forest–canopy snow interception, and from blowing snow particles. Static-surface sublimation of snow on the ground (in nature and SnowModel) depends on air temperature, the air’s moisture deficit, wind speed, and the other components of the surface energy balance. Because this simulation did not include blowing-snow processes, sublimation of wind-transported snow particles are not included in the Fig. 14 totals; we therefore expect nonforested area sublimation totals presented here are underestimated by an unknown amount. In previous studies, snow sublimation has been found to be an important component of the Arctic moisture balance, representing between 10% and 50% of the total winter precipitation, and blowing-snow sublimation can be a key component of that balance (e.g., Liston and Sturm 1998, 2002, 2004). Static-surface sublimation trends (Figs. 14b–d; Table 1) showed considerable spatial variability, with 77% of the domain having negative trends.
Ice-crust formation, resulting from rain-on-snow (ROS) events, can have considerable consequences for animals living in snow-covered areas, because these crusts limit access to winter forage (Aanes et al. 2000; Putkonen et al. 2009). Following Rennert et al. (2009), ROS events were defined as a minimum of 3-mm MicroMet-simulated of liquid precipitation falling on a minimum of 5 mm of SnowModel-simulated SWE depth (Fig. 15a). The spatial variation of ROS trends across the Arctic system domain is considerable (Figs. 15b–d).
As a whole, the domain-averaged temperature trend during snow-covered periods was 0.17°C decade−1, and 73% and 27% of the domain showed a positive or negative trend, respectively (Table 1). Since our model simulations were performed at a relatively high resolution, marked differences in domain elevation and land-cover patterns increased the spatial heterogeneity compared with coarser-resolution simulations. As described elsewhere (Turner and Overland 2009), the temperature change trend pattern is expected to be heterogeneous with some regions cooling and others warming over the same time period.
The largest negative air temperature trends were in the Yamalo-Nenets and Khanty-Mansi regions (i.e., the Ural Mountains and the West Siberian Plain) in the Russian Federation; the Russian Amur region and northeast China; northern Kamchatka and southern Chukotka Russian Federation; Southwest Alaska; and the northern Canadian Archipelago (Fig. 3). The largest warming trends were in Scandinavia; vast areas of Northwest Territories, Nunavut, and Quebec, Canada; and Greenland. The pan-Arctic climate record air temperature anomaly patterns identified by Overland et al. (2004) for the coincident 1979–2002 record also appear in Fig. 3, showing fidelity with meteorological station forcing within the MERRA dataset. Likewise, 1991–2005 temperature difference patterns reported for Canada’s Mackenzie River Delta (Bonsal and Kochtubajda 2009) are in general agreement with our SnowModel results. Furthermore, the 1979–2009 0.3°C decline in temperature associated with the Alaska Box (Fig. 3c) is corroborated with an identical observed winter temperature decline during 1977–2005 (Shulski and Wendler 2007).
Snow precipitation is distributed largely as expected (Figs. 5 and 11); higher elevations and coastal ranges adjacent to warmer ocean waters have increased precipitation, while interior continental regions and much of the Arctic Ocean margins are drier. Annual snow precipitation patterns appear superficially similar (considering vast scale differences) to 1979–93 40-yr European Centre for Medium-Range Weather Forecasts (ECMWF) Re-Analysis (ERA-40) data (Serreze et al. 2005). Pan-Arctic domain-averaged snow precipitation decreased slightly at −0.02 cm decade−1 (Table 1), which contrasts with the 1979–95 ~0.3 cm decade−1 increase in fall, winter, and spring precipitation (55°–85°N) reported by Serreze et al. (2000), and the ~0.3 cm decade−1 annual precipitation increase for (60°–80°N) during 1979–98 reported by New et al. (2001). With respect to longer trends, Ye (2001) reported a winter precipitation increase of 0.4 cm decade−1 over northern central Eurasia for 1926–93. Rawlins et al. (2006) found a decreasing snowfall trend of −0.3 cm decade−1 for the former Soviet Union during 1956–99.
In terms of specific sites with long-term records, Hinzman et al. (2005) reported annual precipitation trends for Barrow, Alaska (−1.29 cm decade−1); Fairbanks, Alaska (0.14 cm decade−1); Fort McMurray, Canada (2.60 cm decade−1); Alert, Canada (0.12 cm decade−1); Yakutsk, Russia (0.50 cm decade−1); and Tiksi, Russia (0.06 cm decade−1). The only significant trends were associated with Yakutsk and Fort McMurray records. With the exceptions of Fairbanks and Fort McMurray, these general trends were visible in the SnowModel simulation (Fig. 5). A later examination of Fairbanks’ annual precipitation record during 1916–2006 showed a nonsignificant trend of −0.54 cm decade−1, with the strongest declines occurring with winter and spring precipitation (Wendler and Shulski 2009). The negative Fairbanks trend more closely resembles SnowModel’s −1.33 cm decade−1 winter precipitation estimate for the Fairbanks-area (central Alaska) regional box (Fig. 5c).
The driving atmospheric forcing dataset and associated reanalysis system (MERRA in our case) strongly influences SnowModel’s simulated values and biases (Adam and Lettenmaier 2003; Yang et al. 2005; Drobot et al. 2006; Bosilovich et al. 2008, 2011; Walsh et al. 2008). Reanalysis data possess uncertainties associated with assimilated observational datasets, temporal data discontinuities, model physics, and assimilation methods and programs (Serreze et al. 2003; Bosilovich et al. 2008, 2011). As part of previous studies, it is clear that SnowModel outputs are governed in large part by the air temperature and precipitation forcing used; no amount of skill in MicroMet and SnowModel can compensate for temperature and/or precipitation input deficiencies or biases. As a consequence, simulated snow characteristics and trends were dictated by MERRA forcings and are susceptible to biases masquerading as trends that may be associated with changes in data streams and observational inputs [see Robertson et al. (2011) for a discussion of MERRA sensitivity to changes in the observing system]. Substantial effort has been dedicated to producing and evaluating reanalyses and their biases (see Bosilovich et al. 2008, 2011; Cullather and Bosilovich 2011), with precipitation being a key variable of interest.
MERRA data were selected for this project because they represented improvements over previous generations of reanalyses (Bosilovich et al. 2008, 2011; Cullather and Bosilovich 2011; Rienecker et al. 2011; more information available online at http://gmao.gsfc.nasa.gov/merra/). MERRA uses the Goddard Earth Observing System Data Assimilation System, version 5 (GEOS-5; Rienecker et al. 2008). A thorough evaluation (Bosilovich et al. 2008) of NASA’s Goddard Earth Observing System, version 4 (GEOS-4), 60°–90°N precipitation for all Januarys (1979–2005) showed the closest agreement between GEOS-4 and Global Precipitation Climatology Project (GPCP) observations (Adler et al. 2003) compared with other reanalyses [e.g., the Japanese 25-year reanalysis (JRA-25; Onogi et al. 2007), ERA-40 (Uppala et al. 2005), the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR; Kalnay et al. 1996), and the NCEP Department of Energy (DOE; Kanamitsu et al. 2002)]. Furthermore, the MERRA atlas (available online at http://gmao.gsfc.nasa.gov/ref/merra/atlas/) has updated comparisons with observational data (GPCP) and coincident reanalyses [i.e., JRA-25, NCEP-DOE, ERA-40, and the Climate Forecast System Reanalysis (CFSR; Saha et al. 2010)]. With annual precipitation during 1979–2009, MERRA shows remarkably similar spatial trends with other reanalyses in comparisons with GPCP data. MERRA tends to closely match observed data (GPCP) with exceptions of the Russian and Scandinavian Arctic, southern Greenland, and Beringia, where MERRA slightly and chronically underestimates precipitation (0–2.5 mm day−1). Localized overestimates (0–2.5 mm day−1) of precipitation are found in the higher elevations of North America (e.g., Pacific coastal mountains and the northern Rockies) and central Asia. Identical seasonal precipitation trends are apparent in September–November and December–February comparisons, while agreement increases in March–May comparisons.
Many have demonstrated that Arctic and Northern Hemisphere snow duration is shorter than in the past (Robinson and Dewey 1990; Groisman et al. 1994; Frei and Robinson 1999; Brown 2000; Serreze et al. 2000; Dye 2002; Groisman et al. 2006; Brown et al. 2007; Turner et al. 2007; Brown and Mote 2009; Brown et al. 2010; McCabe and Wolock 2010; Choi et al. 2010), and SnowModel results offer additional confirmation of this trend at a higher resolution (Fig. 7). Averaged over the simulation domain, snow arrives 1.29 days decade−1 later in fall (Fig. 9; Table 1) and departs 1.28 days decade−1 earlier in spring (Fig. 10; Table 1), with the exception of coastal high-elevation mountain systems where already deep snows are getting deeper (Fig. 11). Areas with a normally shallow snow cover have abbreviated snow-covered periods (e.g., Alaska’s North Slope), while deeper snow areas are augmented and have a longer snow-covered season (Fig. 7; Räisänen 2008). Our distributed snow-onset and snow-free estimates (Figs. 9–10) are largely comparable with Dye’s (2002) 1972–2000 calculated trends of 0.4–3.6 day decade−1 later fall snow arrival and −3.2 to −5 day decade−1 earlier snow-free spring conditions. Spatially, SnowModel snow-cover duration (Fig. 7) duplicates remote sensing analyses for Canada (Brown et al. 2007) and much coarser-resolution GCM simulation and NOAA data (Brown and Mote 2009). Furthermore, SnowModel data show a remarkable spatial agreement with coincident (1972–2008) NOAA weekly snow-cover data trends (Brown et al. 2010) with one notable exception: the NOAA dataset shows western North America coastal snow having a declining trend while SnowModel’s trend is lengthened (Figs. 7, 9, and 10), but other trends are largely duplicated. Snow-cover duration is important for a number of reasons. Altered snow regimes can produce substantial differences in surface energy balance reflected in air temperature records, especially in spring when solar radiation is intense (Chapin et al. 2005). In addition, snow cover plays an important role in governing ground climate below the snow (Taras et al. 2002; Sturm et al. 2005; Lawrence and Slater 2010).
Tracking snow density changes through time (Fig. 13) and estimates of sublimation (Fig. 14) are novel contributions of this modeling effort. Density values are not normally reported in most modeling studies, although efforts are being made to describe functional snow density classes using observations (Sturm et al. 2010). Sublimation changes through time could have important ramifications for water budgets, since sublimation values are likely changing (Fig. 14) alongside precipitation amounts (Fig. 5).
Numerous studies have suggested the number of ROS events (Fig. 15) will increase under a warming climate (e.g., Putkonen and Roe 2003; Rennert et al. 2009). Extensive winter mortality in reindeer, caribou, and muskoxen populations can result from ice-crust formation that prevents access to winter food. For example, in October 2003 an extreme ROS event killed approximately 20 000 muskoxen on Banks Island in western Arctic Canada (Putkonen et al. 2009). As another dramatic example, on Svalbard, Norway, during winter 1993/94, rain, followed by below-freezing temperatures, produced a 10-cm-thick ice crust that led to a 78% decrease in the reindeer population (360 individuals in 1993 to 78 individuals in 1994; Aanes et al. 2000).
To help understand the physics behind our model results, and as part of our model output data analyses, we performed numerous regressions among the model output variables. The only clear relationship was between solid (snow) precipitation and peak SWE. All other attempts at finding simple relationships among basic atmospheric variables and the more complex snow-property variables (e.g., between air temperature and snow-free date) failed. Unfortunately, even the strong relationship between snowfall and peak SWE is partially an artifact of the modeling tools and spatial resolution used, and is not completely representative of the natural system; in tundra environments (Fig. 2) local snow accumulations are largely controlled by wind redistribution of snow, and large-scale precipitation patterns may have little in common with local patterns (e.g., Liston and Sturm 2002; Liston 2004).
While there is something attractive about defining simple relationships among first-order atmospheric forcings, like air temperature and precipitation, and snow features and characteristics, such relationships rarely exist. For example, at first glance, using air temperature to define melt rates and snow-free dates may have some merit. But on further consideration, in reality it is necessary to know: 1) SWE (this defines how much snow must be melted), 2) incoming solar radiation (this defines the melt rate), and 3) whether the snow surface (skin) temperature is at the melting point (this defines whether melting will occur). While air temperature may be a reasonable proxy for skin temperature (as opposed to doing a full energy balance to obtain the skin temperature), incoming solar radiation is the energy that melts high-latitude snow covers (an order of magnitude more than that provided by the sensible heat flux associated with air temperature). Another possibility is that air temperature could be used to define snow-onset dates. Without more information this again fails, because information is needed regarding the timing of precipitation relative to the (near or below freezing) air temperatures, and whether snow accumulation melts or remains on the ground in subsequent days. As a third example, attempts to define discrete snow-related events (e.g., ROS) using general indices (e.g., monthly mean air temperature) are destined to fail because these indices are unrelated to the physical processes associated with the event of interest. And finally, clearly defined and simple relationships between air temperature and solid precipitation would make a welcome contribution to our understanding of snow-evolution changes. Unfortunately, the natural system is not so simple, and precipitation mechanisms depend on many factors in addition to air temperature, such as surface fluxes, moisture advection, cloud microphysical processes, and orographic lifting of air masses. Ultimately, we concluded that simple relationships among basic meteorological variables and snow-evolution features do not generally exist, and reproducing the complexity of the natural system requires physically based modeling systems capable of accounting for threshold and other nonlinear aspects of snow interactions with the environment.
This model simulation and the associated analyses have not addressed how these snow-related features may change in the future. Numerous other studies have examined this (e.g., Hosaka et al. 2005; Pohl et al. 2007; Räisänen 2008; Adam et al. 2009; Brown and Mote 2009; Lawrence and Slater 2010), and the broad consensus is that the general decrease in Arctic snow and snow-related features and characteristics identified in this study will continue well into the future.
This model simulation is not without its limitations. The SnowModel snow-evolution simulation assumed one-way atmospheric forcing, where the atmospheric conditions were prescribed at each time step without regard for the snow conditions (or lack of snow conditions) at the ground surface. In the natural system, atmospheric variables like air temperature and humidity would be modified depending on the surface state, such as snow-covered or snow-free conditions, a dry or melting snowpack, and/or protruding or buried vegetation. These feedbacks were not included in the simulation presented herein. These kinds of interactions are included in regional and global climate models, and 30-yr simulations with these models at 10-km resolution are expected to be possible in the near future (Bromwich et al. 2010; Hines et al. 2010; Shukla et al. 2010). In addition, MicroMet was used to downscale the MERRA atmospheric forcing data from its ⅔° longitude by ½° latitude grid to the 10-km SnowModel simulation grid. An atmospheric model with more physics and dynamics (a regional or global climate model) would be expected to create improved downscaled temperature and precipitation fields. Again, simulations with such a modeling system are expected to soon be possible (e.g., Bromwich et al. 2010; Hines et al. 2010). Also, as noted previously, because of the relatively large grid increment, this simulation did not include blowing snow processes. This influences the simulated snow depth, snow density, SWE, sublimation distributions, and other associated aspects of the snow-cover evolution physics in all nonforested areas of the domain (Fig. 1). A solution to this limitation would be to implement a subgrid blowing snow parameterization to account for the relevant processes and interactions (e.g., Bowling et al. 2004). Alternatively, SnowModel, with SnowTran-3D (Liston and Sturm 1998; Liston et al. 2007), could be run at much higher resolution (e.g., ≤250-m grid increment) to define the blowing snow sublimation fluxes. An additional limitation is the assumption of static vegetation distributions in the simulation. Numerous studies have observed changes, for example, in Arctic tree line location and shrub distributions (see Sturm et al. 2005).
High temporal and spatial resolution snow data represent a critical deficiency in current Arctic monitoring and modeling efforts. The SnowModel pan-Arctic simulation dataset offers a rich look at snow climatology and properties during 1979–2009 at unprecedented spatial (i.e., 10 km) and temporal (i.e., 3 hourly) resolution. The data were realistically distributed in time and space, and air temperature and snow-onset and departure trends largely concur with previous lower-resolution climate studies. The merging of state-of-the-art atmospheric forcing datasets (i.e., MERRA) with relatively high-resolution modeling tools (i.e., MicroMet and SnowModel) allowed a detailed mapping of spatial variations in pan-Arctic snow trends. These trends exhibit strong regional variations, which are attributed to a combination of spatial variations in atmospheric forcing and spatial variations in snow-evolution physics. The nonlinear interactions among snow accumulation and ablation processes created heterogeneity far beyond that seen in the basic atmospheric forcing variables (cf. Figs. 3a and 5a with 11a). MicroMet and SnowModel can be thought of as detailed process models that take our understanding of snow physics and dynamics, and convert basic meteorology such as air temperature, humidity, precipitation, wind, and radiation, into the evolution of complex snow variables such as depth, density, and sublimation.
The relatively high spatial resolution of this Arctic dataset allows important insights into the regional distributions of snow-related features. A pervasive characteristic of the simulated snow fields is strong regional variability. Throughout the Arctic, regions of positive and negative trends are the rule rather than the exception. Positive snow-season air temperature trends covered 73% of the simulation domain, negative solid precipitation trends covered 64% of the domain, the number of days with snow on the ground decreased for over 76% of the domain, the snow-onset date was later in the year for 65% of the area, the snow-free date was earlier for 80% of the domain, and the maximum SWE decreased for over 61% of the simulation domain.
Almost without exception, the domain-averaged 30-yr trends indicate decreasing snow throughout the Arctic: the number of days in the core snow season and the total number of days with snow cover has decreased over the last 30 yr. The onset of snow in the fall occurs later in the year, and the snow-free date occurs earlier. The maximum SWE found during the snow season is decreasing, average snow density is increasing, and the number of rain-on-snow events is increasing. All of these are associated with increasing Arctic snow-season temperatures.
Overall the Arctic has warmed over the last 30 yr, with some notable but small regional exceptions: this trend is predicted to continue and perhaps accelerate (e.g., Chapman and Walsh 2007; Turner and Overland 2009). While climate models suggest winter precipitation will rise (e.g., Meehl et al. 2007; Räisänen 2008; Deser et al. 2010), in fact over the past 30 yr, coastal winter precipitation has risen, but in continental areas the signal is spatially variable. This spatial pattern has been loosely translated into a similar pattern for snow-season length, where the overwhelming signal (75% of the pan-Arctic) is seeing shorter snow seasons, driven by later trends in snow onset and earlier snow-free dates. Peak winter SWE is paramount for ground/permafrost insulation and runoff; SWE trends were related closely to winter precipitation trends.
This work and its resultant datasets have implications for future avenues of investigation of snow–climate interactions. Snow has a clear influence on land surface hydrology, ground temperatures, and permafrost (Bartlett et al. 2005; Zhang 2005; Lawrence and Slater 2010), that can be better quantified with the SnowModel dataset. Given the importance of snow and cryospheric processes on ecosystem structure and function (Chapin et al. 2005; Clein et al. 2007; Euskirchen et al. 2007; Post et al. 2009), this snow-properties dataset represents another leap forward toward a more explicit understanding of links among snow, landscapes, and a changing climate. For example, snow distribution patterns and snow-free duration may well be associated with observed changes in tundra normalized difference vegetation index (NDVI) or phenology (Bunn and Goetz 2006; Verbyla 2008; Bhatt et al. 2010; de Beurs and Henebry 2010). In addition, it is likely that soil temperatures and microbial processes are also changing in light of the altered Arctic snow regime (e.g., Sturm et al. 2001, 2005; Liston et al. 2002).
With a duration averaging 219 ±50 days yr−1, the pan-Arctic snow cover influences numerous climate-related interactions among the atmosphere, hydrosphere, biosphere, and other aspects of the cryosphere. Our current understanding and remote sensing and model-based representations of these interactions still suffer weaknesses associated with spatial resolution, lack of interactive physics and dynamics, inadequate observational datasets, and incomplete insight into the critical linkages and feedbacks among key processes. Future research efforts will fill these gaps and improve our understanding of the role snow plays throughout the Arctic. This model simulation and its associated analyses contribute to the next generation of Arctic-system understanding, where high-resolution datasets and information play a key role.
This work was supported by NSF Grants 0629279 and 0632133 and NASA Grant NNX08AI03G. The authors thank Matthew Sturm for his insightful contributions to the direction of this research, and Skip Walker and Jed Kaplan for their land-cover mapping discussions and advice. Andy Bunn, Matt Callihan, Kelly Elder, Sebastian Mernild, David Robinson, Drew Slater, John Strack, Sveta Stuefer, Matthew Sturm, Amy Tidwell, and two anonymous reviewers provided helpful comments on the manuscript. We would also like to thank the NASA Goddard Earth Sciences (GES) Data and Information Services Center (DISC) and Global Modeling and Assimilation Office (GMAO) for providing the MERRA datasets, NSIDC for providing the EASE-Grid weekly snow-cover extent and monthly SWE data, and NASA for providing the USAF/ETAC snow depth climatology data.