Snow

The creation of saturation conditions necessary for the formation of water droplets or ice particles occur mainly through convection or updraft, cyclonic cooling induced by circulation, frontal or non-frontal lifting of warm air, or orographic cooling by mountain barriers. Snow forms primarily through heterogeneous nucleation. This process involves air that is saturated with a temperature below 0 °C. Water vapor condenses and solidifies, or vapor is deposited on nuclei, which grow into ice and snow crystals. These freezing nuclei may be clay mineral dust (kaolinite, e.g., becomes active at −9 °C), aerosols, pollutants, ice crystal splinters from clouds above, or artificial seeding agents (solid CO₂ or “dry ice,” silver iodide, or urea). The crystals may continue growing through interactions between crystals (crystal aggregation) or with supercooled water droplets, a process called riming (the capture of supercooled cloud droplets by snow crystals) to form snow pellets and/or snowflakes (Reference MosimannMosimann et al., 1993). The minimum size of ice crystals involved in riming is ~60 μm diameter for hexagonal plates and 30 μm width and 60 μm length for columnar ice crystals (Reference ÁvilaÁvila et al., 2009). Under extremely low temperatures (below −40 °C), ice particles can also be formed by the spontaneous freezing of water molecules, which is called “homogeneous nucleation.” Homogeneous nucleation of water droplets occurs at −40 °C; at −10 °C approximately 1/10⁶ drops freeze and at −30 °C about 1/10³ drops freeze.

Ice crystal shapes are hexagonal in form from 0 °C to −80 °C and cubic form from −80 °C to −130 °C. The reason is that an oxygen molecule is tetrahedral; two together form a hexagon or tetrahedra offset by 60° form a cubic crystal. A cubic crystal will transform to a hexagon if warmed but not vice versa. Crystal types have a dependence on temperature and saturation vapor pressure over ice. Under various combinations of temperature and supersaturation conditions with respect to ice, a wide range of snowflakes/pellets results (Figure 2.1). In general, as the temperature decreases, plates → needles → prisms. They can be broadly classified as dendritic and sector plates that involve crystal growth on the a-axis (horizontal), or columns (prisms and needles) that involve growth on the c-axis (vertical) (Figure 2.1). Reference MasonMason (1994) suggests that transitions between crystal types in clouds can lead to more effective release of precipitation through the formation of precipitation elements that have a better chance of surviving below-cloud-base evaporation.

Figure 2.1 Types of snow crystals resulting from various combination of temperature and supersaturation.

(D. Kline, after Kobayashi, 1961)

Snowfall

Whenever snow crystals grow to a size when gravitational pull exceeds the buoyancy effect of air, snowfall occurs. Snowfall typically reaches the ground when the freezing level is not higher than about 250 m above the surface, and the surface air temperature averages ≤1.2 °C. Snow may fall as snowflakes, snow grains (the solid equivalent of drizzle; white, opaque ice particles ≤1 mm in diameter), or graupel (snow pellets of opaque conical or rounded ice particles 2–5 mm in diameter formed by aggregation).

Snowflakes

Snowflakes can be classified into many types (Reference Grey, Prowse and MaidmentGrey and Prowse, 1993; Reference Sturm, Holmgren and ListonSturm et al., 1995). Snowflakes form through the growth of ice crystals by the accretion of water vapor and by their aggregation in branched clusters. The saturation vapor pressure is lower over an ice surface than a water surface, reaching a maximum difference of 0.12 mb at −12 °C. As a result, in a mixed phase cloud, supercooled water droplets tend to evaporate and vapor is deposited onto ice crystals. This is known as the Bergeron–Findeisen process after its discoverers. Snowflakes grow in small cap clouds over elevated terrain when ice crystals falling from an upper cloud layer seed them. This is known as the seeder–feeder mechanism (Reference BarryBarry, 2008, p. 273). Ice crystals may float in the atmosphere as “diamond dust” when the air temperature is ≤−40 °C. The design and variations of snowflakes are way beyond human imaginations, as some examples in Figure 2.2 that show needle, sheath, and varieties of stellar crystals with plates, dendritic, and sector-like branches. Bentley, who was born in 1865, even believed that no two snowflakes are exactly alike (Reference TeelTeel, 1994).

Figure 2.2 Examples of snowflakes classified according to Reference Magono and LeeMagono and Lee (1966): 1. Needle, 2. Sheath, 3. Stellar crystal, 4. Stellar crystal with sector-like ends, 5. Stellar crystal with plates at ends, 6. Crystal with broad branches, 7. Plate, 8. Plate with simple extension, 9. Plate with sector-like ends, 10. Rimed plate with sector-like ends, 11. Hexagonal plate with dendritic extensions, 12. Plate with dendritic extensions, 13. Dendritic crystal, 14. Dendritic crystal with sector-like ends, 15. Rimed stellar crystal with plates at ends, and 16. Stellar crystal with dendrites

Depth hoar

Other than in permafrost areas (high latitudes or high elevations in middle latitudes), the ground is mostly warm or near freezing when the ground is snow covered. This is true even when the air is very cold, because snow is a good insulator. Therefore, there will usually be liquid water in the snowpack, and it is common for the snow near the ground to remain damp for most of the winter. Depth hoar forms at the base of a snowpack, as a result of large temperature gradients between the warm ground and the cold snow surface, when rising water vapor freezes onto existing snow crystals. It usually requires a thin snowpack combined with a clear sky or low air temperature, and it grows best at snow temperatures from −2 °C to −15 °C. Therefore, the occurrence of depth hoar is common in high Arctic regions such as Alaska, the Northwest Territory, Nunavut, and northern Siberia (Reference DerksenDerksen et al., 2009). Depth hoar constitutes about 20% of all snow layers in the high Arctic (46% in sub-Arctic), an average grain size of 6.5 mm (long-axis) and 2 mm (short-axis), and about 0.23 gm cm⁻³ in density (Reference DerksenDerksen et al., 2014).

Depth hoar consists of sparkly, large-grained, faceted, cup-shaped ice crystals up to 10 mm in diameter. Beginning and intermediate facets are 1–3 mm square; advanced facets can be cup-shaped 4–10 mm in size. Larger-grained depth hoar is more persistent and can last for weeks. Depth hoar is strong in compression but not so in shear, and hence often behaves like a stack of champagne glasses; it can fail in the form of collapsing layers, or in shear, with fractures often propagating long distances and around corners. Even though the stability of the snowpack depends on the cohesion between layers of snow and meteorological factors, almost all catastrophic avalanches, which involve the entire season’s snow cover, fail on depth hoar layers (Reference TremperTremper, 2008). Seismic sensors have been used for the remote detection of snow avalanches for they can also estimate avalanche velocity, size, and type (Reference LacroixLacroix et al., 2012).

Introduction

Snow is an integral component of the global climate system because of its linkages and its feedback between surface energy, moisture fluxes, clouds, precipitation, hydrology, and atmospheric circulation (Reference King, Armstrong and BrunKing et al., 2008). It is the most spatially extensive and seasonally variable component of the global cryosphere (see Table 1.1). On an average, snow covers almost 50% of the Northern Hemisphere’s land surface in late January, with an August minimum of about 1%. In addition, there is perennial snow cover over the Antarctic ice sheet (12 million km²) and at the higher elevations of the Greenland Ice Sheet (about 0.6 million km²) (Figure 2.3).

Figure 2.3 Averaged monthly snow cover area of Northern Hemisphere in (×10⁶) km² calculated from weekly snow cover extent maps produced primarily from daily visible satellite imagery of NOAA-AVHRR by the Rutgers Global Snow Lab

Since snow produces substantial changes in the surface characteristics, and the atmosphere is sensitive to physical changes of the earth surface, its presence over large areas of the Earth for at least part of the year exerts an important influence on the climate, both locally and globally. The best-known effect involves the albedo–temperature positive feedback, whereby an expanded (reduced) snow cover increases (decreases) the reflection of incoming solar radiation, reducing (increasing) the temperature, and thereby encouraging an expansion (reduction) of the snow cover. Fresh snow has a spectrally integrated albedo of 0.8–0.9, making it the most reflective natural surface. This value decreases with age to 0.4–0.7 as the snow density increases through settling and snow metamorphism and is reduced still further by impurities in or on the snow (e.g., mineral dust, soot, aerosols, biogenic matter) (see Figure 2.4). The cooling effect of snow cover is illustrated by the example that, in the Upper Midwest of the United States, winter months with snow cover are about 5–7 °C colder than the same months without snow cover. Snow, a poor conductor of heat, also insulates the soil surface and sea ice. Therefore, a better knowledge of the snow cover and its properties over large regions will lead to a better understanding of our climate.

Figure 2.4 Field measurements of broadband albedo at Mammoth Mountain in the Sierra Nevada for (a) 2003 and (b) 2005, showing albedo in the visible, near-infrared and all solar radiation

(adapted from Reference DozierDozier et al., 2009)

Snow is held in cold storage until there is sufficient energy to melt it to water or to sublimate it to water vapor. The storage of water in the seasonal snow cover introduces into the hydrological cycle an important delay of weeks to months, causing a peak in the annual runoff in spring and early summer when the river water is agriculturally more valuable. It is highly beneficial to be able to estimate the amount and timing of release of this stored precipitation to spring runoff, which allows a better management of water resources for irrigation and hydroelectric production planning. The dynamics of water storage in seasonal snowpack is also critical to the effective management of water resources globally. Snow water accumulated in winter in the Arctic river basins is critical for the springtime snowmelt, and the freshwater from its river systems accounts for about 50% of the net flux of freshwater into the Arctic Ocean (Reference Barry, Serreze and LewisBarry and Serreze, 2000), which is a large percentage when compared to the freshwater inputs to the tropical oceans, where freshwater input is dominated by direct precipitation. Frozen soil affects the snowmelt runoff and soil hydrology by reducing the soil permeability. Runoff affects ocean salinity and sea ice conditions (Reference PetersonPeterson et al., 2002), and the degree of surface freshening can affect the global thermohaline circulation (Reference BroeckerBroecker, 1997). In future, we expect even larger freshwater input from the Arctic as the global hydrologic cycle accelerates, higher precipitation is expected in high latitudes and higher runoff from Arctic river basins (Reference Nummelin, Ilicak, Li and SmedsrudNummerlin et al., 2016).

Snowpacks affect energy and water exchanges, and so both snow cover and SWE are important climatic and hydrologic variables. In particular, snow controls the climate and hydrology of the cryosphere and higher latitude regions significantly, and the amount and distribution of snow is affected by the climate and vegetation types. In the Canadian Prairies, mixed precipitation can occur within a certain range of temperature (Reference KienzleKienzle, 2008), but on a whole approximately one-third of its annual precipitation occurs as snowfall, and the shallow snow cover generates as much as 80% of the annual surface runoff. In the Colorado Rockies, the Sierra Nevada of California, and the Cascade Mountains of Washington, snowmelt can account for up to 65–80% of the annual water supply (Reference SerrezeSerreze et al., 1999), but shifts in snowmelt runoff due to climate warming is expected to affect water supply in snow-dominated river basins (Reference FosterFritz et al., 2011).

The snow covers of North America (NA) and Eurasia change seasonally, in accordance with the position of the Sun that shines directly at the Tropic of Cancer in the Northern Hemisphere (NH) on June 21 (summer solstice) and then moves southward, reaching the Tropic of Capricorn in the Southern Hemisphere on December 22 (winter solstice), before moving northward for the next 6 months; the cycle repeats itself on an annual time scale. The extent of snow cover in the NH lands reaches an average maximum of about 46.8 × 10⁶ km² in January and February, and an average minimum of about 3.4 × 10⁶ km² in August (see Figure 2.5) (Reference RopelewskiRopelweski, 1989; Reference Brown, Armstrong, Armstrong and BrunBrown and Armstrong, 2008; Reference RobinsonRobinson, 2008), which constitute 8% and about 0.5% of the Earth’s surface, respectively. From 1966 to 2008, the maximum January snow cover of NH ranged from as low as 42 × 10⁶ km² (in 1982) to as high as 50.1 × 10⁶ km² (in 2008) (GSL, 2008). For 1966–2008, the mean annual NH snow extent was 25.5 × 10⁶ km² (Reference RobinsonRobinson, 2008). The 1967–2016 mean annual NH snow extent was lower after mid-1980s relative to what it was before the mid-1980s. Even though the linear trend estimated was − 25,000 km² per year, this is due to a step-like drop in the mid-1980s, for the decline of NH snow cover is nonlinear, and above-average snow cover was even observed in 2015–2016 (Reference ConnollyConnolly et al., 2019).

Figure 2.5 Seasonal variation in the mean monthly snow and sea ice cover extent for January, April, and July over Northern Hemisphere using data of NSIDC over 1967–2005 for snow and 1979–2005 for ice; for January and July over Antarctic/Southern Hemisphere over 1987–2002 for snow and 1979–2003 for ice (Reference MaurerMaurer, 2007) by Lambert Azimuthal Equal-Area (http://nsidc.org/data/atlas) projection; January 31, 2008 snow and ice chart of NH

adapted from NOAA-AVHRR image of NOAA (http://wattsupwiththat.com/2008/02/09/jan08-northern-hemisphere-snow-cover-largest-since-1966/)

In the NH, most mid-summer snow cover is found over the Greenland and some parts of the Canadian High Arctic (Figure 2.5a), while about 60% of winter snow cover is found over Eurasia and 40% over Canada and the upper portion of the United States, sometimes down to latitude 30° N (Figure 2.5c). Figure 2.6 is a composite monthly NOAA-AVHRR image of NA that shows large seasonal variations in snow cover between the four seasons. In contrast, in South America, there is only a small area covered with snow in July.

Figure 2.6 Seasonal variation in the mean monthly snow cover extent for (a) July, (b) October, (c) January, and (d) April over North America computed from snow charts derived from weekly visible satellite images of NOAA-AVHRR over 1972–1993 (www.tor.ec.gc.ca/CRYSYS/cry-edu.htm); (e) Northern Hemisphere snow and forest covers for January, 2005 computed from the NSIDC Equal-Area Scalable Earth Grid (EASE-Grid) snow cover product (Armstrong and Brodzik, 2005) and the University of Maryland global land cover classification

(Reference Hanson and HookeHansen et al., 2000) (taken from Reference Rutter and EsseryRutter et al., 2009)

Snow cover is observed in situ at hydrometeorological stations, from daily depth measurements, (monthly) snow courses, and in special automated networks such as about 730 SNOwpack TELemetry (SNOTEL) automated systems of snow pressure pillows, sonic snow depth sensors, precipitation gauges, and temperature sensors distributed across the United States. The extent of snow cover is also observed and mapped daily (since June 1999) over the NH from the operational satellites of the NOAA, United States.

Canada has extensive in situ snow depth and snow course networks that are a valuable database for monitoring cryospheric changes and for validating satellite data such as those shown in Figures 2.5 and 2.6. However, most of the field observations are concentrated in southern latitudes and lower elevations, where the majority of the population lives. At many northern sites, manned stations have been replaced by automatic weather station (AWS) that use acoustic sounders to measure the height of the snow surface.

Besides seasonal variability, snow cover is subject to inter-annual fluctuations but only about 40% of these have been found to be associated with continental to hemispheric scale forcing (Reference Robinson, Frei and SerrezeRobinson et al., 1995), and the rest could be partly attributed to regional forcings or “coherent” regions. Regarding the boreal winter variability, Reference Saito and CohenSaito and Cohen (2003) noted that snow cover at continental scale varies similarly as the inter-annual to inter-decadal oscillations of an internal mode of the atmosphere, but it leads the atmosphere by several months through their mutual oscillations. By principal component analysis (PCA) and composite analysis, Reference Frei and RobinsonFrei and Robinson (1999) found that over western NA, snow cover extent (SCE) is associated with the longitudinal North American ridge, the Pacific North America (PNA) index, while over eastern NA, it is associated with the meridional oscillation of the 500-mb geopotential height, the North Atlantic Oscillation (NAO), and the teleconnection patterns are coupled to tropospheric variability during autumn and winter. Reference Gobena and GanGobena and Gan (2006) found during El Niño winters, the southeasterly flow of warm dry Pacific air and the northwesterly flow of cool dry Arctic air will be the dominant flow over western Canada and Pacific Northwest (PNW) of the United States, giving rise to drier climate (less snowfall) over these regions. On the other hand, La Niña winters are associated with an erosion of the western Canadian ridge and strengthening of the Pacific Westerly, giving rise to greater moisture supply and so more winter snowpack in western Canada and the PNW of the United States. For Europe, Reference Henderson and LeathersHenderson and Leathers (2010) found that large (small) snow extent is associated with negative (positive) 850 hPa zonal wind anomalies, negative (positive) phase of the North Atlantic Oscillation, negative (positive) 1,000–500 hPa thickness anomalies, and generally positive (negative) Northern European precipitation anomalies. Besides solar radiation, snowpacks are related to surface air temperature, precipitation, storm tracks, and mid-tropospheric geopotential heights at 500 mb. Reference Frei and RobinsonFrei and Robinson (1999) postulate that snow extent, by exerting an influence on lower tropospheric dynamics (e.g., air temperature), could even modulate atmospheric circulations.

Reference BrownBrown (2000) observed some decline in NH snow cover in recent decades, but the declines are not statistically significant. From 1972 to 2000, using weekly NH snow cover data of high latitude and high elevation areas derived from visible bands of NOAA satellite observations, Reference DyeDye (2002) found that the week of the last-observed snow cover in spring shifted earlier by 3–5 days per decade estimated from a linear regression analysis, and the duration of the snow-free period increased by 5–6 days per decade, primarily as a result of earlier snow cover disappearance in spring. Similarly, based on the 1966–2007 snow cover data of NOAA satellites and simulations from the Coupled Model Intercomparison Project Phase 3 Model (CMIP3), on the response of NH land area with seasonal snow cover to warming and increasing precipitation, Reference Brown and MoteBrown and Mote (2009) found the largest decrease in snow cover duration (SCD) was concentrated in zones where seasonal mean air temperatures were in the range of −5 to +5 °C, which extended around the mid-latitudinal coastal margins of the continents. Regional studies in the western United States (e.g., Reference Adam, Hamlet and LettenmaierAdam et al., 2009) show that losses of snowpack associated with warming trends have been ongoing since the mid-twentieth century, especially near boundaries of areas that currently experience substantial snowfall. These findings very likely reflect clear signals of human-induced impact on the climate shown by the changing snowpacks of NH and by the river flows of western United States (Reference Barnett, Pierce, Hidalgo, Bonfils, Santer, Das, Bala, Wood, Nozawa, Mirin, Cayan and DettingerBarnett et al., 2008), and Reference Brown and RobinsonBrown and Robinson (2011) detected significant contractions of snow cover in NH over 1922–2010.

According to the CMIP5 Report, observed SCE has decreased in the NH, especially in spring (Reference Vaughan and StockerVaughan et al., 2013). Satellite records indicate that over the period 1967–2012, the annual mean SCE negative trend is statistically significant, with the largest change of −53% (−40% to −66%) occurred in June. Over 1922–2012, the available snow cover data for March and April show a 7% decline and a strong negative (−0.76) correlation with the March–April 40° N to 60° N land temperature. The spring SCE in Arctic land areas north of 60° N has significantly declined since 1960s, estimated respectively at −3.5% in May and −13.4% in June per decade between 1981 and 2018, relative to the 1981–2010 mean, from multiple data sets (SROCC, 2019). From surface observations, satellite data, and model-based analyses, the snow cover duration has also become shorter, between −0.7 and −3.9 days per decade depending on region and time period, but all spring snow cover duration trends from all data sets are negative (Reference Bulygina, Groisman, Razuvaev and KorshunovaBulygina et al., 2011; Reference Liston and HiemstraListon and Hiemstra, 2011; Reference Estilow, Young and RobinsonEstilow et al., 2015). These same multisource data sets also identify reductions in autumn snow extent and duration at −0.6 to −1.4 days per decade (Reference Brown, Vikhamar-Schuler, Bulygina, Derksen, Loujus and MudrykBrown et al., 2017). The CMIP6 report will provide further updates on observed snow cover changes (Reference Eyring, Bony, Meehl, Senior, Stevens, Stouffer and TaylorEyring et al., 2016).

The mountain snow cover is characterized by a very strong interannual and decadal variability (Reference Mankin and DiffenbaughMankin and Diffenbaugh, 2015). Long-term in situ snow cover data are limited in some regions, particularly in High Mountains of Asia and Northern Asia. For key mountainous regions, Reference StewartStewart (2009) found that higher temperatures have decreased snowpack and resulted in earlier melt in spite of precipitation increases at mid-elevation regions but not at high-elevation regions, which remain well below freezing during winter. Mountain snow cover at lower elevations has generally declined in duration by about several days per decade mainly due to more precipitation falling as rain and to higher melt rate at most elevations, mostly due to increased air temperature (Reference Marty, Tilg and JonasMarty et al., 2017). The mean snow depth and SWE have declined since the mid-twentieth century, with regional variations but at higher elevation, snow cover trends are generally insignificant (SROCC of Reference Pörtner, Roberts, Masson-Delmotte, Zhai, Tignor, Poloczanska, Mintenbeck, Nicolai, Okem, Petzold, Rama and WeyerIPCC, 2019).

At lower elevations in the European Alps, Western North America, Himalaya, and subtropical Andes, the snow depth or mass is projected to decline by 25% between the recent past period (1986–2005) and the near future (2031–2050). By 2081–2100, reductions of up to 80% are expected under RCP8.5, 50% under RCP4.5, and 30% under RCP2.6 climate scenarios of CMIP5 (SROCC of Reference Pörtner, Roberts, Masson-Delmotte, Zhai, Tignor, Poloczanska, Mintenbeck, Nicolai, Okem, Petzold, Rama and WeyerIPCC, 2019, ch. 2). At higher elevations, projected reductions are smaller as temperature increases at higher elevations mainly affect the ablation component of snow mass evolution. The projected increase in wintertime snow accumulation may even result in a net increase in winter snow mass. All elevation levels and mountain regions are projected to exhibit sustained interannual variability of snow conditions throughout the twenty-first century.

Snow cover, depth distribution, and blowing snow

At continental scale or larger, snow cover distribution primarily depends on latitude and seasons (Figures 2.5 and 2.6). At the macro or regional scale, for areas up to 10⁶ km², and distances from 10 to 1,000 km, snow cover distribution depends on latitude, elevation, orography, and meteorological factors. For example, snowfall caused by orographic cooling tends to increase with a rise in elevation, and frontal activities involving cold fronts generally produce more intense snowfall over relatively smaller areas as against warm fronts that produce moderate or light snowfall over larger areas, because the former has relatively steep leading edge while the latter has mild leading edge. On the mesoscale, with distances of 100 m to 10 km, snow distribution depends on the blowing effect of wind, relief, and vegetation patterns, while on the microscale, 10–100 m, the influencing factors are more local. Over highly exposed terrain, the effects of meso- and microscale differences in vegetation and terrain features may produce wide variations in accumulation patterns and snow depths.

Blowing snow occurs when the force of wind exceeds the shear strength of the snowpack surface that resists snow particles to move. Blowing snow increases with wind speeds and the amount of snowfall but decreases with increasing surface roughness. The effects of wind on the accumulation and distribution of a snowpack are most pronounced in open environments, for example, the Canadian Prairies or Siberian steppes, with three modes of snow particle movement: snow particles begin in motion by creeping or rolling on snowpack surface, then by saltation or bouncing when wind speed increases, and finally in turbulent diffusion or snow particles suspended in the air under high wind speed. These three modes of transport typically occur less than 1 cm above ground under a low wind speed U < 5 m s⁻¹, between 1 and 10 cm for U = 5–10 m s⁻¹, and between 1 and 100 m for U > 10 m s⁻¹, respectively.

Based on wind tunnel studies with surface wind speeds of up to 40 m s⁻¹, Reference DyuninDyunin et al. (1977) argued that saltation accounts for most drifting snow at all conceivable wind speeds. However, Reference Budd, Dingle and RadokBudd et al. (1964) found that turbulent suspension was the primary mechanism from snowdrift studies at Byrd Station, Antarctica. Suspension increases at about U⁴, whereas saltation increases linearly with U at high wind speeds (at which most transport occurs), so suspension dominates the overall effect of wind (Reference Pomeroy and GrayPomeroy and Gray, 1990; Pomeroy, p.c. December 2009). At low wind speeds, saltation is the dominant process.

Blowing snow is important in open environments, especially for high elevation alpine areas above treeline, in the Prairies, and in the tundra of the North American Arctic and Siberia. In these regions snow depth variation depends mainly on terrain features because without the hindering effect of vegetation cover, wind causes snow drift and redistribution to smooth topography, so that mountain tops and plateau tend to have thin snowpack as snow tends to be blown to valleys and low-lying areas which as a result tend to have relatively thick snowpack. In the coastal tundra and open subarctic forest near Churchill, Manitoba of Canada, Reference Kershaw and McCullochKershaw and McCulloch (2007) found that snowpack characteristics measured from 2002 to 2004 also depend on vegetation characteristics, ecosystems, and associated micro-climates. Ecosystems that dominate the circumpolar north are such as wetland, black and white spruce forest, burned forest, forest-tundra transition, and tundra. In lower latitudes as the forest canopy density generally increases, higher snow accumulation has been found in forests of medium density (25–40%) than large open areas because of reduced wind effects by densely forested areas (e.g., Reference VeatchVeacth et al., 2009). Although forest structure and canopy interception, and local terrain characteristics will influence snow retention at local scales, Reference LundquistLundquist et al. (2013) show that where the mean DJF temperatures exceed −1 °C, forests with lower total canopy cover are likely to enhance snow retention by minimizing mid-winter and early spring melt.

Reference Gordon, Savelyev and TaylorGordon et al. (2009) developed a camera system to measure the relative blowing snow density profile near the snow surface in Churchill, Manitoba, and Franklin Bay, Northwest Territory. Within the saltation layer, they found that the observed vertical profile of mass density is proportional to exp(−0.61z/H), where H, the average height of the saltating particles, varies from 1.0 to 10.4 mm, while z, the extent of the saltation layer, varies from 17 to over 85 mm. At greater heights, z > 0.2 m, the blowing snow density varies according to a power law (ρ_s∝z^−ɤ), with a negative exponent 0.5 < γ < 3. Between these saltation and suspension regions, results suggest that the blowing snow density decreases following a power law with an exponent possibly as high as γ≈8.

Modeling blowing snow

Reference Pomeroy, Gray and LandinePomeroy et al. (1993) developed the first comprehensive blowing snow model for the prairies environment. It estimates saltation, suspension and sublimation using readily available meteorological data. They show that within the first 300 m of fetch, transport removes 38–85% of the annual snowfall. However, beyond 1 km of fetch, sublimation losses from blowing snow dominate over transport losses. In Saskatchewan, sublimation losses are 44–74% of annual snowfall over a 4-km fetch. Subsequently, Reference PomeroyPomeroy (2000) showed that the ratio of snow removed and sublimated by blowing snow to that transported at prairies (arctic) sites was 2:1 (1:1), respectively.

Reference Essery, Long and PomeroyEssery et al. (1999) developed a distributed model of blowing snow transport and sublimation to consider physically based treatments of blowing snow and wind over complex terrain for an Arctic tundra basin. By considering sublimation, which typically removes 15–45% of the seasonal snow cover, the model is able to reproduce the distributions of snow mass, classified by vegetation type and landform, which they approximated with lognormal distributions. The representation used for the downwind development of blowing snow with changes in wind speed and surface characteristics is shown to have a moderating influence on snow redistribution. Spatial fields of snow depth have power spectra in one and two dimensions that occur in two frequency intervals separated by a scale break between 7 and 45 m (Reference Trujillo, Ramirez and ElderTrujillo et al., 2007). The break in scaling is controlled by the spatial distribution of vegetation height when wind redistribution is minimal and by the interaction of the wind with surface concavities and vegetation when wind redistribution is dominant.

Reference Liston and SturmListon and Sturm (1998) developed a SnowTran-3D that simulates wind-driven snow-depth evolution over topographically variable terrain and tested it in an arctic–tundra landscape, for it is generally applicable to treeless areas characterized by strong winds, below-freezing temperatures, and solid precipitation. Reference LawlerListon et al. (2007) extended the SnowTran-3D to version 2.0 that simulates wind-related snow distributions over the range of topographic and climatic environments globally. This version includes three primary enhancements to the original model: (1) an improved wind sub-model, (2) a two-layer sub-model describing the spatial and temporal evolution of friction velocity that must be exceeded to transport snow, and (3) a 3-D, equilibrium-drift profile sub-model that forces snow accumulations to duplicate observed drift profiles.

In mountainous regions, wind plays a prominent role in determining snow accumulation patterns and turbulent heat exchanges, strongly affecting the timing and magnitude of snowmelt runoff. Reference Winstral and MarksWinstral and Marks (2002) use digital terrain analysis to quantify aspects of the upwind topography related to wind shelter and exposure. They develop a distributed time-series of snow accumulation rates and wind speeds to force a distributed snow model. Terrain parameters were used to distribute rates of snow accumulation and wind speeds at an hourly time step for input to ISNOBAL, an energy and mass balance snow model that accurately modeled the observed snow distribution (including the formation of drifts and scoured wind-exposed ridges) and snowmelt runoff. In contrast, ISNOBAL forced with spatially constant accumulation rates and wind speeds taken from the sheltered meteorological site at Reynolds Mountain in southwest Idaho, a typical snow-monitoring site, overestimated peak snowmelt runoff, and underestimated snowmelt inputs prior to the peak runoff.

Remote sensing of SWE

Because snowpack can attenuate gamma radiation, over a thousand flight lines have been conducted mainly in USA and Canada to collect airborne gamma data for SWE survey during winters and distributed by the NOHRSC of NOAA-USA (www.nohrsc.noaa.gov/snowsurvey). However, as airborne data, gamma radiation data are not as readily available as space borne data and also because gamma radiation is attenuated by water in all phases, the effect of soil moisture has to be accounted for to avoid inaccurate estimation of SWE.

Since 1978, spatially distributed SWE data have been retrieved from the brightness temperature (TB) in Kelvin from passive microwave remote sensing sensors such as the Television Infrared Observation Satellite (TIROS-N) launched in 1978 (Reference Dozier, Schneider and McGinnisDozier et al., 1981), and Scanning Multichannel Microwave Radiometer (SMMR) flew on NASA’s Nimbus 7 from October 25, 1978 to August 20, 1987, and the Special Sensor Microwave Imager (SSM/I) mounted on the Defense Meteorological Satellite Program satellites since September 7, 1987. Since May 2002, TBs retrieved from the Advanced Microwave Scanning Radiometer-EOS (AMSR-E) sensor aboard the Aqua satellite have been used to estimate SWE. The AMSR-E sensor uses one of the largest ever microwave antennas to detect faint microwave emissions from the Earth’s surface. AMSR-E produces global and continuous daily SWE and snow depth data sets at 25-km spatial resolution. The frequencies (resolution) of SSM/I are 6.6 (150 km), 19, 22, 37 (25 km) to 85 (12.5 km) GHz, while that of AMSR-E are 6.9 (50 km) to 89 GHz (5 km) resolution. These TBs are provided in the form of horizontal (H) and vertical (V) polarizations.

The basis of this approach is that microwaves are not sensitive to water vapor or liquid water in the troposphere, especially at long wavelength, and microwave radiation emitted by the land surface is attenuated by snow cover. However, the effects of any liquid water due to snowmelt, the obscuring of the snow surface by the vegetation canopy, changes in the grain size of the snow, and terrain irregularities greatly complicate such determinations. Moreover, the typical satellite footprint is of the order of 12–25 km, meaning that the signal is a complex spatial average and hard to relate to point measurements. Nevertheless, such methods are being used for operational SWE mapping over the high plains and prairies of North America (Reference Goodison, Walker, Thirkettle, Kite and WankiewiczGoodison et al., 1990; Reference Singh and GanSingh and Gan, 2000, Reference Singh, Gan and Gobena2009), and other parts of NH (Reference TakalaTakala et al., 2011), unless the snow packs are wet, reasonably accurate SWE can mostly be estimated from the passive microwave data of SMMR, SSM/I, and AMSR-E. Passive microwave data represent an important resource for monitoring trends and spatial variability of terrestrial snowpacks (Reference Wulder, Nelson, Derksen and SeemannWulder et al., 2007).

Because of volumetric scattering (Figure 2.13), the dominant loss mechanism for microwave radiation greater than 15 GHz incident on a snowpack, it is possible to relate empirically the brightness temperature (TB) of certain frequencies and polarizations (H or V) to the SWE of the snow packs (e.g., Reference Singh and GanSingh and Gan, 2000; Reference Armstrong and BrodzikArmstrong et al., 2001; Reference DerksenDerksen et al., 2009). The volumetric scattering of snow on the incident microwave radiation depends on the snow grain size, snow density, depth, SWE, temperature, degree of metamorphism, nocturnal crust development, ice lenses, and other factors (Reference MätzlerMätzler, 1994). In theory, the greater the depth of snow, the lower should be the TB if climatic and snowpack conditions remain the same, but some studies have indicated that TB can increase with depth beyond a certain snow depth (Reference Mätzler, Schanda and WoodMätzler et al., 1982; Reference SchandaSchanda, 1983). The snow scattering property is affected by factors such as snow metamorphism, which dictates the internal structure of snow as it ages, the multiple melt-freeze cycles, which together contribute to the complicated physical processes in the formation of snow structure; and snow redistribution that depend mainly on wind, terrain features, and land use (Reference Rott, Nagler, Choudhury, Kerr, Njoki and PampaloniRott and Nagler, 1995; Reference Rosenfeld and GrodyRosenfeld and Grody, 2000).

Figure 2.13 In the volumetric scattering of snow, some fractions of the radiation can enter beyond the boundary into the lower, layered snowpacks. The reflected radiation will consist of that reflected at the top planar boundary, and that which is first partially transmitted into the second medium and reflected at the second planar boundary and partially transmitted back to the first medium, and so on

Reference ShiShi et al. (2009) used SNTHERM to simulate the snow profiles of the 1200-km transect of snow stratigraphy measured in the Alaska region, called the Snow Science Traverse-Alaska Region (SnowSTAR2002), and the Microwave emission Model for Layered Snowpacks (MEMLS) to simulate the 19- and 37-GHz TB for both SNTHERM-MEMLS and SnowSTAR2002-MEMLS cases, which matched well with the passive microwave data. They concluded that the simulation of snow microphysical profiles is a viable strategy for retrieving SWE from passive microwave data.

Over the past several decades, various snow retrieval algorithms for passive microwave data have been developed and applied (e.g., Reference ChangChang et al., 1982; Reference Hallikainen and JolmaHallikainen and Jolma, 1992; Reference Walker and GoodisonWalker and Goodison, 1993; Reference FerraroFerraro et al., 1994; Reference GanGan, 1996; Reference WilsonWilson et al., 1999; Reference Singh and GanSingh and Gan, 2000; Reference Gan, Kalinga and SinghGan et al., 2009). In general, microwave-derived maps tend to underestimate snow extent during fall and early winter, due to a weak signal from shallow and intermittent snow cover, and the underestimation can be as much as 20%, decreasing to a few percent during mid-winter and spring. Conversely, a thinner atmosphere between the snow-covered surface and satellite could lead to over-estimating SCE, as is the case for the Tibetan Plateau 3,000–5,000 m above sea level (Reference SavoieSavoie et al., 2009). SWE products for NH from the National Snow and Ice Data Center (NSIDC) and GlobSnow from the European Space Agency (ESA) have been developed using different algorithms. By comparing the data with historical snow depth measurements obtained from ~7,400 meteorological stations across the NH, Reference Liu, Li, Huang and TianLiu et al. (2014) show that for SWEs between 30 and 200 mm, GlobSnow products agree better with ground measurements than NSIDC products which tend to suffer from “microwave saturation,” underestimating SWEs over 120 mm. However, GlobSnow products tend to overestimate SWE less than 30 mm marginally more often than that of NSIDC products. Using the new GlobSnow 3.0 dataset, Pulliainen et al. (2020) show that the 1980–2018 annual maximum snow mass in the NH was about 3,062 ± 35 gigatons with different continental trends over this study period, e.g., snow mass decreased by 46 gigatons per decade across NA but had a negligible trend across Eurasia.

Linear and nonlinear regression algorithms to estimate SWE

The 1979–1987 SMMR SWE data provided by the NSIDC (National Snow and Ice Data Center) were retrieved from Equation 2.8,

$SWE\left(mm\right)=4.77\left(T_{B H19}-T_{B H37}\right)$(2.8)

where TB is in K. The SSM/I SWE data of NSIDC were retrieved using Equation 2.9, and false SWE signals from deserts are filtered (Reference Armstrong and BrodzikArmstrong and Brodzik, 2002). To ensure that snow packs are detectable by passive microwave data, SWE less than 7.5 mm is considered unreliable and set to zero.

$SWE\left(mm\right)=4.77\left(T_{B H19}-T_{B H37}-5\right)$(2.9)

The daily SWE is adjusted for surface forest cover (A_F in percent) using the BU-MODIS (NSIDC, 2005) land cover data so that,

$SWE\left(mm\right)=\frac{SWE}{\left(1-A_F\right)}$(2.10)

Any pixels with forest cover higher than 50% are set to the 50% threshold, so that the forest correction by Equation (2.10) is limited to a maximum factor of 2, given that microwave data can only detect snowpack properties at canopy densities less than 60–70% (Reference Pulliainen, Koskinen and HallikainenPulliainen et al., 2001). Using measurements from snow course transects in the former Soviet Union, Reference Armstrong and BrodzikArmstrong and Brodzik (2002) reported a general tendency for nearly all of the algorithms tested to underestimate SWE, especially when the forest-cover density exceeded 30–40%. Ideally, specific SWE retrieval algorithms should be derived for different vegetation types such as boreal forest or tundra (Reference Derksen, Walker, Goodison and StrappDerksen et al., 2005) but that is only feasible with adequate field campaigns which are expensive and time consuming.

Reference HallikainenHallikainen (1989) developed an algorithm similar to Equation 2.9 to account for the surface effect of land cover when there is no snow. Reference Goodison, Walker, Thirkettle, Kite and WankiewiczGoodison et al. (1990) used an algorithm similar to Equation (2.9) but vertically polarized data TB _V19 and TB_V37, to model the snow of Canadian prairies. Reference DerksenDerksen et al. (2009) found TB _V37 of AMSR-E to be the appropriate polarization for remote sensing of the SWE of tundra snowpack in the Arctic.

The retrieval of SWE from passive microwave data based on Equations 2.8–2.10 applied in a continental framework is expected to have varying accuracy, depending on the density and types of vegetation, and snowpack characteristics, frozen and unfrozen water, north- or south-facing slopes, topographic variability, mountains versus flat plains, and so on (Reference ChangChang et al., 1997; Reference TaitTait, 1998; Reference Goita, Walker and GoodisonGoita et al., 2003; Reference TakalaTakala et al., 2011), and the possible effect of wet snow albeit even though NSIDC selectively used “cold” pass daily TB to prepare these monthly SWE data. Lake ice could also cause low SWE values and the presence of snowmelt, ice lenses, surface hoar, depth hoar and very deep snow packs could confound the interpretation of TB for SWE. In view of the aforementioned limitations, SWE data derived from passive microwave data should be more dependable in the Canadian Prairies dominated by grassland than in the Canadian Arctic with countless frozen water bodies or in forested, mountainous areas such as the Rocky Mountains.

Besides the effect of forest cover, various screening criteria have been proposed to eliminate data affected by depth hoar, wet snow, and water bodies. Reference ChangChang et al. (1982) indicated the possibility of discriminating the effect of depth hoar and the underlying ground condition (frozen or unfrozen) using the polarization factor, p-factor = (V37 − H37)/(V37 + H37). Reference Singh and GanSingh and Gan (2000) used a p-factor >0.026 to eliminate the SSM/I data that were affected by depth hoar and the presence of a water body of significant size in the vicinity of the footprint that has the effect of causing an underestimation of predicted SWE because of the high dielectric constant of the water body (or the presence of water in the snowpack due to above-freezing temperatures). Other recommended screening criteria are such as TB _V37 < 250 K, (TB _V19 − TB _V37) ≥ 9 K, TB _V37 − TB _H37 ≥ 10 K, TB _V37 > 225 K (e.g., Reference Goodison, Walker, Choudhury, Kerr, Njoku and PampaloniGoodison and Walker, 1994).

The traditional TB difference between 19 and 37 GHz has been shown to be inappropriate for the lake-rich environment in the Arctic, and retrieving tundra SWE can be challenging because the SWE versus 37-GHz TB relationship could have a reversed slope that occurs beyond a theoretical limit of approximately 120-mm SWE (Reference de la MareDerksen et al., 2009). In the northern environment, where a footprint could include both frozen water bodies (A_W) and tundra (A_TUNDRA), Reference GanGan (1996) assumed microwave emission of the former to be related to air temperature (T_a),

$SWE=K_1\left(A_{TUNDRA}\right)\left(T_{B H19}-T_{B H37}\right)+K_2\left(A_W\right)\left(T_a\right)+K_3$(2.11)

Unlike Equations (2.8) and (2.9), coefficients K₁, K₂, and K₃ for Equation (2.11) will be region or basin dependent.

The above algorithms assume relatively simple TB-SWE relationships, even though TB is influenced by many snow parameters, and so more complicated relationships have been developed, such as Equation (2.12) that Reference Singh and GanSingh and Gan (2000) applied to the Red River Basin of North Dakota and Minnesota,

$SWE\left(mm\right)=K_4\left(T_{B V19}-T_{B H37}\right)+K_5\left(AMSL\right)+K_6\left(1-A_F\right)+K_7\left(1-A_W\right)+K_{8}TPW$(2.12)

where AMSL is the average elevation above mean sea level and TPW is the total precipitable water. Again, the coefficients are site specific. For the Red River Basin, by first screening SSM/I data to remove data representing wet snow, and calibrate the screened SSM/I data (dry snow) with corresponding airborne SWE data, they found K₄ = 0.2357, K₅ = 0.0064, K₆ = 4.0399, K₇ = 20.0287, and K₈ = 1.0825. Reference Singh and GanSingh and Gan (2000) also developed another algorithm similar to Equation (2.11) that was based on the surface/ground TB converted from the at-satellite TB by applying the atmospheric attenuation model of Reference ChoudhuryChoudary (1993) which accounted for the attenuation of atmosphere water vapor (based on TPW in cm) on the satellite TB data.

Above are parametric regressions relating TB to SWE. The non-parametric approach has also being used, such as the Projection Pursuit Regression (PPR) of Reference Friedman and StuetzleFriedman and Stuetzle (1981) by Reference Gan, Kalinga and SinghGan et al. (2009). PPR models the response variable as a sum of functions of linear combinations of predictor variables. Suppose y and x’s denote response and predictor vectors respectively, PPR finds the number of terms M_o, direction vectors (ɑ₁, ɑ₂, . ., ɑ_Mo), and nonlinear transformations (Φ₁, Φ₂, . ., Φ_Mo) as shown in Equation (2.13),

(2.13)$\hat{y} \approx \overline{y}+\sum _{m=1}^{Mo}{\beta }_m{\varphi }_m\left({\alpha }_m^{T}x\right).$

Through minimizing the expected distance or mean square error between y, which is the observed SWE, and $\hat{y}$ , the estimated SWE, using Equation (2.14), the model parameters β_m (the response linear combinations), ɑ_m (the direction vectors), Φ_m (the predictor functions), for m = 1, 2, …, M_o are obtained.

$L_2\left(\beta ,\alpha ,\varphi ,x,y\right)=E{\left[y – ŷ\right]}^2$(2.14)

A successful application of PPR lies in selecting an optimum number of terms, M_o, determined by trial and error, often by starting the algorithm with a large M_o and then decreasing M_o such that the increase in accuracy due to an additional term is not worth the increased complexity (Reference FriedmanFriedman, 1985). The optimum M_o is determined in terms of the fraction of variance it cannot explain (Reference Friedman and StuetzleFriedman and Stuetzle, 1981). From Equation (2.14), this unexplained variance, U is given as

$U=\frac{L_2(\beta ,\alpha ,\varphi ,x,y)}{\text{Var}\left(y\right)}.$(2.15)

Artificial neural network (ANN) algorithms to estimate SWE

ANN has been widely used in many fields of research because of its ability to model nonlinear and poorly understood systems with their inherent non-linearity and complex internal structure. Other than some drawbacks such as being classified as black box models, the problem of over-fitting and tedious training, ANN can approximate almost any function. Reference Gan, Kalinga and SinghGan et al. (2009) used the Modified Counter Propagation Network (MCPN) to model the SWE for the Red River basin of North Dakota and Minnesota from SSM/I data. The MCPN that makes use of the self-organizing feature map (SOFM) learning algorithm (Reference Kohonen, Oja, Simula and KangasKohonen et al., 1996), consists of an interconnected network of three layers, namely, the input, hidden, and output layers.

The unsupervised clustering procedure of SOFM performs the input-hidden layer transformation (SD_i → IP_j), the non-linear part of the Input/Output mapping. The training of SOFM is carried out by computing the distance d_j between the normalized input vector (the input snow data SD_i, given in Table 2.2) and the weighting vector w_ji as

$d_j=\sqrt{\left[\sum _{i=1}^{n_O}{\left({\text{SD}}_i-w_{ji}\right)}^2\right]}j\text{=1,}\text{.}\text{.}\text{.,}{n}_1$(2.16)

where n_o is the number of input variables, and n₁ is the number of hidden nodes, among which the winning node (c) has the smallest d_j (d_c = min (d_j), j = I_c). To complete the SOFM training, the updating of weights w_ji is performed only for the hidden nodes in the neighborhood Λ_c surrounding the winner node as

$\begin{array}{c}{w}_{ji}\left(t\right)=w_{ji}\left(t-1\right)+\eta \left(t\right)\left({\text{SD}}_i-w_{ji}\left(t-1\right)\right),\text{for}j\in {\Lambda }_c\left(t\right)\\ i=1,2,\dots ,n_0,0<\eta \left(t\right)<1\\ w_{ji}\left(t\right)=w_{ji}(t-1)\text{otherwise,}\end{array}$(2.17)

with t being the iteration counter for the training process, and η(t) is the learning rate which together with Λ_c(t) are decreased iteratively from their initial settings of η₀= 0.2 ~0.5 and Λ₀ = n₁/2. Before the hidden-output layer transformation, the intermediate output parameters, IP_j, corresponding to the input vector, SD_i, are computed as

$\begin{array}{c}{\text{IP}}_j \text{=1–}{d}_j\text{for}j\in \text{Ω}\\ {\text{IP}}_j \text{=0 otherwise,}\end{array}$(2.18)

where Ω the size of hidden nodes centered on the neighborhood of I_c, should be equal to or greater than that of SD_i. The training of the weights (M_kj) required for the hidden-output transformation (IP_j → SWE_k) is performed by a simple recursive gradient search. In the neighborhood of Ω surrounding the active node, these weights are recursively updated as

$v_{kj}\left(t\right)=v_{kj}(t-1)+\beta \left(T{\text{SWE}}_k-{\text{SWE}}_k\right)y_j,\text{for}j\in \Omega ,k=c,$(2.19)

where β is the learning step size (0 ≤ β ≤ 1), and Tswe_k is the observed (target) SWE and SWE_k is its estimated value after each iteration. Adjustments are made to M_kj to obtain their final values, which together with IP_j are used to compute the output, SWE_k

$\begin{array}{c}{\text{SWE}}_k=\sum _j{M}_{kj}{\text{IP}}_i\text{for}j\in \Omega ,k \text{=} c\\ \\ {\text{SWE}}_k=\varnothing \text{for}k\ne c.\end{array}$(2.20)

At 25 × 25 km resolution, we expect the micro-scale spatial variability of snowpack to be mostly averaged out, and so we can mostly expect to detect meso- to macro-scale variability of snowpack from the above data. Figure 2.14 shows the monthly distributions of SWE in NH for November 2003–February 2004 derived from SSM/I (left) and AMSR-E (right), which are very similar to each other with limited differences. As expected, such SWE data are subjected to errors, particularly SWE values from mountainous areas with large topographic variability due to a possible mixture of deep snow on north-facing slopes but almost snow free on south-facing slopes, or forested areas because of mixed microwave emissions from trees, snow canopy and ground surface. Figure 2.15 shows some obvious differences between the Northern Hemisphere SWE map with enlarged area in Northern Europe and Western Russia for October 24–31, 2003, derived from (a) the NSIDC Global Monthly EASE-Grid SWE Climatology, (b) with additional snow-covered area (red) determined from at least 25% of component MODIS CMG pixels indicating snow cover, and (c) by combining 89 GHz data with MODIS snow-covered area (selected day with maximum difference in [37 − 85] GHz) (Reference Armstrong, Brodzik, Savoie and KnowlesArmstrong et al., 2006). Figure 2.15c shows an improvement in shallow snow SWE estimates for in passive microwave imagery, wet areas containing melting snow or wet snow packs could return low or zero SWE values. Reference Singh and GanSingh and Gan (2000) employed screening procedures to eliminate potentially erroneous SWE data in their analysis.

Figure 2.14 Northern Hemisphere monthly average snow water equivalent (SWE) derived from SSM/I (left) and AMSR-E (right), November 2003–February 2004

(Courtesy of Brodzik, M. J., NSIDC/CIRES)

Figure 2.15 (a) shows Northern hemisphere snow water equivalent (SWE) map derived from NSIDC Global Monthly EASE-Grid SWE Climatology, with enlarged area in Northern Europe and Western Russia for October 24–31, 2003; (b) shows additional snow-covered area (red), determined from at least 25% of component MODIS CMG pixels indicating snow cover; (c) shows improvement in shallow snow SWE estimates by combining 89 GHz data with MODIS snow-covered area (select day w/ max diff 37 − 85 GHz)

(Courtesy of Brodzik, M. J., NSIDC/CIRES)

Reference Gan, Kalinga and SinghGan et al. (2009) retrieved the SWE for the Red River basin of North Dakota and Minnesota using SSM/I, physiographic and atmospheric data by MCPN, a Projection Pursuit Regression (PPR) and a nonlinear regression. They used the airborne gamma-ray measurements of SWE as the observed SWE. They screened the SSM/I data for the presence of wet snow, large water bodies like lakes and rivers, and depth hoar. They found MCPN to produce encouraging results in both calibration and validation stages (R² was about 0.9 for both calibration (C) and validation (V)), better than PPR (R² was 0.86 for C and 0.62 for V), which in turn was better than the multivariate nonlinear regression at the calibration stage (R² was 0.78 for C and 0.71 for V). MCPN is probably better than the linear and nonlinear regression counterparts because of its parallel computing structure and its ability to learn and generalize information from the SWE-SSM/I relationships.

Reference TakalaTakala et al. (2011) developed an algorithm assimilating synoptic snow depth data of weather stations with passive microwave radiometer data to produce a 30-year time-series of seasonal SWE for the NH, which was validated against independent SWE data from Russia, the former Soviet Union, Finland and Canada. The results show an overall strong performance in the SWE retrieved with root mean square errors below 40 mm when SWE <150 mm, but the retrieval uncertainty increases when SWE is above 150 mm. The SWE estimates obtained by this assimilation approach are better than SWE obtained from typical stand-alone retrieval algorithms shown above.

Active microwave data

In contrast to passive microwave data, active microwave data acquired by synthetic-aperture radar (SAR) sensors can attain much higher spatial resolutions. SWE can also be estimated using active microwave data such as the SIR-C/X-SAR data by Reference Shi and DozierShi and Dozier (2000), or possible future data from the dual-frequency (X-/Ku-band) synthetic aperture radar (SAR) satellite proposed by the Cold Land Processes Working Group (CLPWG) of NASA. Active and passive microwave sensors have different sensitivities to the same snowpack properties, and so they can provide complementary information. With improved resolution (about 10 km), passive microwave measurements can help supporting the high-resolution radar measurements, make linkages across process scales, and help relating future, dual-frequency SAR missions such as the Snow and Cold Land Process (SCLP) mission to the long-term microwave record of snow.

Empirical snowmelt-runoff models

Empirical or statistical models (Equation 2.8) such as linear (all γ_i = 1, i = 1, 2, …, n) or nonlinear (at least one or more of γ_i ≠ 1, i = 1, 2, …, n) regressions based on SWE measured using snow pillows or snow course samplers at selected sites and/or baseflow measured for winter months such as November to March, as independent variables, for example, B₁, B₂, …, B_n, have been popularly used to forecast spring snowmelt of river basins. These models assume that under average snowpack conditions, empirical relationships derived between spring snowmelt runoff and measured snowpack data in the past are applicable in future years.

$Q=a_o+a_1{B}_1{}^{\gamma }1+a_2{B}_2{}^{\gamma }2+\dots \dots .+a_n{B}_n{}^{\gamma }n.$(2.21)

Another statistical method for predicting spring snowmelt runoff is to relate SCA with spring snowmelt.

Degree-day or temperature index (TI) method

Many operational snowmelt runoff models use the degree day (temperature index) approach that estimates the snowmelt rate M (mm d⁻¹) as the difference between the mean daily air temperature (T_a) and a melt-threshold or base temperature (T_thm) adjusted by some optimized melt factors m_f in mm d⁻¹ °C⁻¹, and a depletion curve that relates how of the original snow cover remains versus mean areal extent of snow cover (e.g., WMO, 1986). m_f depends on vegetation types, the slope and aspect of the land surface, percent of snow cover, time of the year and the climatic regime (Reference Frank and LeeFrank and Lee, 1966). So in applying Equation (2.22) to each basin zone, a different m_f is usually used to reflect the vegetation characteristics and climate of each zone. The TI method has been widely used in operational snowmelt forecasting because by adjusting the degree-day, (T_a − T_thm), with an appropriate m_f, it can approximately represent the daily energy supply for melting snowpack on a regional basis.

$M=m_f\left(T_a-T_{thm}\right).$(2.22)

In addition, m_f has been found to vary seasonally, and tends to increase as the season progresses because of a decrease in snow albedo, and an increase in incoming solar radiation; in some cases the melt factor is allowed to vary through the melt season, such as a form of sinusoidal function in Equation 2.23,

$m_f=\frac{m_{f\max }+m_{fmin}}{2}+\sin \left(\frac{2\pi n}{366}\right)\left[\frac{m_{f\max }-m_{fmin}}{2}\right],$(2.23)

where m_fmax and m_fmin are the maximum and minimum melt rate, and n is the Julian date.

Even though the degree-day approach may work well in mid-latitude or temperate environments when there is a strong correlation between T_a and the dominant energy responsible for snowmelt, it does not adequately account for many climatic factors related to snowmelt. For example, Reference Male and GrangerMale and Granger (1981) showed that in open, non-forested areas the short wave radiation exchange is a dominant melt-producing energy flux, but it generally does not correlate well with air temperature. Moreover, the consideration of a uniform snow accumulation in the whole elevation range followed by a snowmelt process based on an assumed areal depletion curve is hardly imaginable in a mountain basin (Reference MartinecMartinec, 1980). Some newer models use satellite images to update its areal snow cover distribution and do not solely rely on an areal depletion curve. Examples of more widely used degree-day models are such as the NWSRFS (National Weather Service River Forecast System) (Reference AndersonAnderson, 1976), SRM (Snowmelt Runoff Model) (Reference Martinec, Rango, Roberts, Baumgartner and ApflMartinec et al., 1998), HBV (Swedish Meteorological and Hydrological Model) (Reference Bergström and SinghBergström, 1995; Reference Seibert and VisSeibert and Vis, 2012), and others.

Modified degree-day or modified temperature index (MTI) method

Reference Singh, Gan and GobenaSingh et al. (2005) introduced a modified degree-day method that includes the near-surface soil temperature (T_g) as an additional predictor,

$m=m_f\left(m_{rf}\right)\left(T_r-T_{thm}\right)$(2.24)

where m_rf is an adjustment factor for m_f so as to better capture the onset of initial snowmelt and is estimated from

$m_{rf}={\left[{\beta }_1+{\beta }_2\left({\tan }^{-1}{T}_g+{\beta }_3\right)\right]}^{\psi },$(2.25)

where β₁, β₂, β_3, and ψ are model parameters derived through calibration (see Reference Singh, Gan and GobenaSingh et al., 2005). As a tangent function, m_rf is relatively small when T_g < 0 °C, and reaches an upper limit of one when T_g ≥ 0 °C. The effect of m_rf is “felt” mostly during the onset of snowmelt because its value approaches 1.0 when T_g approaches 0 °C. The desired rate of change of m_rf can be achieved by adjusting the parameter ψ. T_r is a reference temperature computed as a weighted average of T_g and T_a,

$T_r=\chi T_a+(1-\chi )T_g,$(2.26)

where χ is also a model parameter. The rationale for the proposed modification comes from past studies that pointed out the importance of T_g as an indicator of spring snowmelt (Reference Woo and ValverdeWoo and Valverde, 1982), as well as the analysis of data observed at the Paddle River Basin (PRB) in central Alberta by Reference Singh, Gan and GobenaSingh et al. (2005). Analysis of hourly data observed for 6 years during the spring season (March 1–April 30) at PRB shows that there is a significant correlation between net radiation (R_n) and T_g at daily time step. Moreover, the data revealed that Pearson’s correlation coefficient between T_g and R_n (ranging from 0.62 to 0.89) was mostly higher than that between T_a and R_n (ranging from 0.47 to 0.87). Since R_n generally dominates the energy balance during spring snowmelt in the Canadian Prairies (Reference ShookShook, 1995), adding T_g as another predictor should improve the performance of Equation 2.22.

Modeling snowmelt by energy balance method (EBM)

A more physically based snowmelt modeling approach is an one-dimensional (1-D) mass and energy balance model, such as the US Army Cold Regions Research and Engineering Laboratory (CRREL) model, SNTHERM developed by Reference JordanJordan (1991) for predicting snowpack properties and temperature profiles (Figure 2.16). SNTHERM calculates energy exchange at the surface and bottom of the snowpack, grain growth, densification and settlement, melting and liquid water flow, heat conduction and vapor diffusion, and solar absorption. The model’s surface boundary conditions require incoming solar and longwave radiation; wind speed, air temperature and humidity at some reference height; and precipitation. Lower boundary conditions include soil textural properties (currently clay or sand used as defaults), wetness and temperature profile.

Figure 2.16 Energy fluxes considered in SNTHERM, a physically based 1-D energy balance model for snow and soil by Reference JordanJordon (1991)

(figure taken from SNTHERM fact sheet)

SNTHERM accounts for changes in albedo due to grain growth, sun angle and cloud cover but it does not account for the decrease in effective albedo when the snow depth is shallow and when radiation penetrates through the snowpack to the underlying soil. It is expected that the snow albedo decreases exponentially to the soil albedo when the radiation penetrates through the snowpack to the underlying soil, and in a forest environment, the accumulation of forest litter could reduce snow albedo. In the energy balance approach, which is mostly considered in the vertical direction only, both the energy content of the snowpack plus a soil layer underneath that interacts thermally with the snowpack should be considered. This procedure provides a simple approximation of the effects of frozen ground, or snow falling on warm ground. The model output provides snow depth, profiles of snow temperature, water content, density, grain size, and surface fluxes of sensible heat and evaporation. Reference ShiShi et al. (2009) used SNTHERM to simulate the snow profiles of a 1,200-km transect of snow stratigraphy measured in Alaska. SNTHERM has been extensively tested for the prediction of temperature, snow depth, SWE, and snow-covered sea ice in sites such as the California Sierra Nevada (Reference MarksMark, 1988), Greenland (Reference Rowe, Kuiven and JordanRowe et al., 1995), the Canadian Boreal forest (Reference HardyHardy et al., 1998), the Arctic (Reference Jordan, Andreas and MakshtasJordan et al., 1999), Antarctic (Reference Andreas, Jordan and MakshtasAndreas et al., 2005), southern Finland (Reference Koivusalo and BurgesKoivusalo and Burges, 1996), and other Nordic environments (Reference KotlyakovLanglois et al., 2009).

A comparative study of three snow models with different complexities was carried out by Reference JinJin et al. (1999) to assess how a physically detailed snow model can improve snow modeling within general circulation models. The three models were (a) SNTHERM; (b) a simplified three-layer model, Snow–Atmosphere–Soil Transfer (SAST), which includes only the ice and liquid-water phases; and (c) the snow submodel of the Biosphere-Atmosphere-Transfer-Scheme (BATS), which calculates snowmelt from the energy budget and snow temperature by the force–restore method. SNTHERM gave the best match to observations with the SAST simulation being close. BATS captured the major processes in the upper layers of a snowpack where solar radiation is the main energy source and gave satisfactory seasonal results.

CROCUS is a model of the Centre d’Etudes de la Neige, Grenoble (Reference BrunBrun et al., 1992). It is a 1-D physical model that determines mass and energy balance for a snow cover and is used for operational avalanche forecasting. The snow cover is represented as a pile of layers parallel to the ground. Energy exchanges are projected orthogonally to the slope. The model describes the evolution of the internal state of the snow cover as a function of meteorological conditions. The variables describing the snow cover are temperature, density, liquid water content, and snow type of each layer. To match the natural layers, the thickness and number of layers are adjusted by the model. The model simulates the heat conduction, melting/refreezing of snow layers, settlement, metamorphism, and percolation. It simulates dry and wet snow metamorphism with experimental laws derived from laboratory data. Snow grains are characterized by their size and type. This allows an accurate albedo of the snow cover to be calculated.

Reference Bartelt and LehningBartelt and Lehning (2002) presented a 1-D physical model of the snowpack (SNOWPACK) with equations for heat transfer, water transport, vapor diffusion, and mechanical deformation. In their model, new snow, snow drift, and ablation are considered, and the snow layers are treated in terms of height, density, and microstructure (grain size, shale, and bonding). The model is used for avalanche warnings in Switzerland. Reference LangloisLanglois et al. (2009) simulated the SWE in southern Quebec using SNOWPACK, CROCUS, and SNTHERM models driven by local and reanalysis meteorological data, which are in agreement with ground measurements of SWE in winter seasons of 2004/2005, 2005/2006, and 2007/2008, with correlation coefficients ranging between 0.72 and 0.99.

One-dimensional vertical EBM

A 1-D energy balance snowmelt model developed by Reference Singh, Gan and GobenaSingh et al. (2009), which they refer to as SDSM-EBM, is briefly described. The model incorporates processes for snow interception by forest canopy, separate snowpack energy and mass balance for open and forested areas, separate water balance for liquid and ice phases, snow sublimation, compaction, refreezing, and so on. The snow interception capacity at different levels of canopy is estimated as a function of LAI, forest types (Reference Hardy and Hansen-BristowHardy and Hansen-Bistow, 1990), tree species, and prevailing forest structure (Reference Golding and SwansonGolding and Swanson, 1986). Reference PomeroyHedstrom and Pomeroy (1998) tested the snow interception model for the Canadian southern boreal forest, which Reference Singh, Gan and GobenaSingh et al. (2009) implemented in SDSM-EBM.

The transfer of energy at the snow surface and the snow/soil interface governs the snowmelt. The amount of energy available for melting snow is determined from the 1-D, energy equation (Equation 2.27) applied to a control volume of snow having upper and lower interfaces with air and ground, respectively. During melting, the snowpack is isothermal at 0 °C ( $T_{sp}^{t+1}$ = 0 °C) and the heat for snowmelt, Q_o, can be calculated as

$Q_o=Q_n\pm Q_h\pm Q_e\pm Q_g+Q_v+Q_f-Q_{cc},$(2.27)

where

• Q_n = Net radiation (short-wave and long-wave) absorbed by the snow,

• Q_h = Convective or turbulent sensible heat flux exchanged at the surface due to a difference in temperature at the snow-air interface,

• Q_e = Convective or turbulent latent heat exchanged at the surface due to vapor movement as a result of a difference in vapor pressure (evaporation, sublimation, condensation) at the snow-air interface,

• Q_g = Ground heat flux at snow-interface by conduction,

• Q_v = Advective heat of precipitation,

• Q_f = Energy released by freezing of liquid water content, and

• Q_cc = Cold content of snowpack in the previous time step.

Part of Q_o may be used to overcome the cold content of the snowpack, Q_cc (see Equation 2.37). All fluxes are computed in terms of W m⁻² or older units such as langley per min. The net shortwave radiation (Q_sn = Q_si − Q_so) can either be measured by a net radiometer or by two pyrometers measuring the incoming short wave radiation (Q_si) and the reflected short wave radiation (Q_so) which depends on the albedo of the snow-covered area, and the net longwave radiation (Q_ln = Q_li − Q_lo) can be similarly measured using two pyrgeometers with one inverted. By knowing the air temperature and the surface temperature of snowpack, and the emissivity of the atmospheric medium, it is also possible to estimate Q_li and Q_lo using the Stefan-Boltzmann’s constant and the blackbody theory.

$Q_n=Q_{sn}+Q_{ln}$(2.28)

Q_s can also be converted to Q_sn [= Q_s(1 − α)] as a function of the areal albedo of a partially ablated snow cover (α), which can be taken as the larger of that retrieved from satellite images such as that of NOAA-AVHRR or MODIS, or computed from

$\alpha ={\alpha }_{sn}{A}_{sn}+{\alpha }_g(1-A_{sn})$(2.29)

where α_sn is the snow albedo that can be estimated from an albedo decay function (Reference Riley, Israelsen and EgglestonRiley et al., 1972), α_g is the albedo of the ground surface and A_sn is the fraction of snow-covered area, which can be tracked using either a linear, $A_{sn}=\frac{SD}{h}$ , or a non-linear, parameterizing the SWE probability distribution based on the basin’s topography and vegetation in a dimensionless form, $A_{sn}\left(D\right)=\oint \left(\frac{D}{D_{\max }}\right)$ , depletion curve. SD (D) is the snow depth (SWE), and h (D_max) is the depth (SWE) at which snow cover is complete (or the depth (SWE) below which bare patches start to emerge). In computing A_sn, $\frac{D}{D_{max}}$ is replaced by $\frac{SD}{h}$ . From snow data collected at the Paddle River Basin (PRB), Reference Singh, Gan and GobenaSingh et al. (2009) found that h ranges between 0.07 and 0.1 m as the cutoff for a partial snow cover. In forest-covered areas, α_sn can be further modified to account for the effect of litter fall fraction (Reference HardyHardy et al., 1998).

Q_s can also be estimated empirically (e.g., Reference BrasBras, 1990) but it is subjected to errors because of possible cloud cover effects and air pollution problems. Instead of measuring snow albedo, it can also be estimated from the age of the last snowfall since α_sn declines with age and α_sn ranges from about 0.9 or higher for freshly fallen snow to less than 0.4 for shallow, dirty, and wet snow. Snowpack albedo has also being related to the cumulative maximum air temperature (US Army Corps of Engineers, 1956), snow grain size, or snow surface density which tends to increase with time.

The other important components of energy balance are Q_h, Q_e and Q_v. Q_g, usually being the smallest, is often ignored. The sensible heat, Q_h, is due to the turbulent flux of energy exchanged at the snow surface, as the result of a difference in temperature between air (T_a) and snow surface (T_s), in contrast to Q_e that is caused by the vapor pressure difference between air and the snow surface. As a turbulent heat transport through convection (like Q_e), Q_h is affected by wind, such that

$Q_h=l_v{k}_2\left[\frac{P_a}{P_o}\right]{\overline{u}}_b(T_a-T_s){\left(z_a{z}_b\right)}^{-}{}^{1/6},$(2.30)

where k₂ = 0.00357 cm m^1/3 hr day^{−1 o}C⁻¹ km⁻¹, P_o (P_a) is the standard (actual) atmospheric pressure, ${\overline{u}}_b$ is the wind speed (m s⁻¹), z_a and z_b = heights where air temperature and wind speed are measured respectively, and l_v ≈ 250 kJ kg⁻¹. A simpler way to estimate Q_h transported through convection is

$Q_h=C_h{\rho }_a{C}_p{\overline{u}}_b(T_a-T_s),$(2.31)

where ρ_a is the air density, C_p is the specific heat of air (1,004 J kg^{−1 o}C⁻¹), and C_h is a bulk transfer coefficient that depends on displacement, roughness height, and atmospheric stability. An even simpler expression to estimate Q_h is

$Q_h=B_h{\overline{u}}_b(T_a–T_s),$(2.32)

where B_h is the bulk transfer coefficient that replaces C_hρ_aC_p given in Equation (2.31). In SDSM-EBM, T_s is computed using one of three simplified heat flow models, namely the force-restore method, the surface conductance method, and the Kondo and Yamazaki method (Reference Singh and GanSingh and Gan, 2005).

The latent heat flux (Q_e) is the sum of surface sublimation/condensation (Q_e,surf) and blowing snow sublimation (Q_e,bss). When water vapor is transported to the snow surface, it changes phase to either liquid or solid, releasing Q_e,surf at the snow surface,

$Q_{e,surf}=l_{e}E=C_e{\rho }_a\left(\frac{0.622l_e}{P_o}\right){\overline{u}}_b\left(e_a-e_s\right),$(2.33)

where E is the rate of vapor transfer, C_e the bulk transfer coefficient for water vapor, l_e is the latent heat of sublimation (2,836 kJ kg⁻¹) or the latent heat of vaporization (2,501 kJ kg⁻¹), e_a and e_s are the actual and snow surface vapor pressure (assumed saturated at T_s) in Pascal, respectively, and P_o is the standard atmospheric pressure at sea level (≈101.33 kPa). e_s is estimated as

$e_s=611\exp \left(\frac{17.27T_s}{237.3+T_s}\right).$(2.34)

The bulk transfer coefficient under neutral condition is computed from Reference BrutsaertBrutsaert (1982),

$C_e=\frac{{\kappa }^2}{{\left[\ln \left(\right(z_r-d_0)/z_0)\right]}^2},$(2.35)

where κ is the von Karman constant, z_r is the reference height, d₀ is the zero displacement height and z₀ is the roughness height. The model can be run assuming neutral conditions or one of three options, namely that of Reference Price and DunnePrice and Dunne (1976), Reference LouisLouis (1979), and Reference MorrisMorris (1989) can be used to adjust C_e.

If the precipitation is rainfall, there is positive advective heat Q_v that contributes toward melting the snowpack, Q_o, given as

$Q_v=c_s{\rho }_{w}RT,$(2.36)

where R = intensity of rainfall in cm hr⁻¹, T = temperature of rainfall in °C, and c_s = specific heat (2.093 kJ kg^{−1 o}C⁻¹ for snow, 4.186 kJ kg^{−1 o}C⁻¹ for water).

The cold content of a snowpack is the heat required to raise the temperature of the snowpack to 0 °C, if the temperature of the snowpack, T_p, which can be different from the snowpack surface temperature T_s, is less than 0 °C.

$Q_{cc} \approx -{\rho }_p{\text{c}}_{\text{s}}D{T}_p,$(2.37)

where ρ_p is the snowpack density, c_s is the snowpack specific heat, D is the snowpack depth in cm, and T_p is the snowpack temperature in °C. If D_cc is the cold content expressed in terms of the equivalent depth of ice at 0 °C that can be melted to water at 0 °C,

(2.38)$D_{cc}=\frac{Q_{cc}}{-{\rho }_w{l}_{fs}}=\frac{{\rho }_p{C}_{s}D{T}_p}{{\rho }_w{l}_{fs}}.$

In computing Q_cc, the heat capacity of the entrapped air is neglected, and the snowpack temperature T_p may be assumed as T_s if no information is available. In other words, T_p is assumed to be independent of depth below snow surface (z).

After overcoming the cold content of the snowpack, Q_cc, the depth of snowmelt, M, due to Q_o acting for a time interval Δt, is then

$M=\frac{(Q_o\Delta t-Q_{cc})}{{\rho }_w{l}_{fs}\theta },$(2.39)

where ρ_w is the density of water and θ is the thermal quality of snow, which is the fraction of ice in a unit mass of wet snow, or the ratio of the heat necessary to produce a given amount of water from snow to the amount of heat needed to produce the same quantity of melt from pure ice at 0 °C. θ usually ranges between 0.95 and 0.97.

The total energy needed to melt a snowpack is Q_oΔt = ρ_pDl_ms + Q_cc, where l_ms is the actual latent heat to melt the snowpack of depth D. If Q is the amount of energy needed to produce the same amount of melt from pure ice at 0 °C, such that Q = ρ_pDl_fs then by definition, the thermal quality of snow, θ is given by

$\theta =\frac{l_{ms}}{l_{fs}}-\frac{c_s{T}_p}{l_{fs}}$(2.40)

In theory θ can also be estimated if we know the liquid water content (W_liq) of a snowpack since by knowing that we can calculate l_ms (Equation 2.7), which will be easier than trying to measure l_ms. Among various methods available to determine W_liq of a snowpack, a popular approach is to measure W_liq by a time-domain reflectometry because of the large difference in dielectric constant between water and ice (Reference SchneebeliSchneebeli et al., 1998).

Given the insulating effect of vegetation cover, Q_g, the exchange of energy between the snowpack and the underlying ground by conduction is often ignored except where the ground is frozen, or in tundra. If Q_g needs to be accounted for, it can be computed as

$Q_g={\lambda }_g\left(\frac{\partial T_g}{\partial z}\right),$(2.41)

where λ_g, the thermal conductivity of soil, is about 0.4–2.1 W m⁻¹ °C⁻¹ for unfrozen silt and clay, and T_g is the ground temperature that changes with elevation z.

Snowpack water balance

The water balance equations can be expressed in terms of water and ice at both canopy and ground levels as

(2.42)${\rho }_w{c}_s{W}^{t+\Delta t}{T}_{sp}^{t+\Delta t}={\rho }_w{c}_s{W}^t{T}_{sp}^t+(Q_n+Q_h+Q_e+Q_g+Q_p)+Q_o.$

Here, W^t+Δt accounts for both the addition of precipitation (P_r or P_s) during the time step and the change in water and ice mass due to Q_e (sublimation or freezing) depending on whether $T_{sp}^t$ is isothermal at zero or less than zero,

$W_{liq}^{t+\Delta t}=\{\begin{array}{c}{W}_{liq}^t+P_r+\frac{Q_e}{{\rho }_w{l}_v}if{T}_{sp}^t=0\\ W_{liq}^t+P_{s}if{T}_{sp}^t<0\end{array}$(2.43)

$W_{ice}^{t+\Delta t}=\{\begin{array}{c}{W}_{ice}^t+P_{s}if{T}_{sp}^t=0\\ W_{ice}^t+P_s+\frac{Q_e}{{\rho }_w{l}_{fs}}if{T}_{sp}^t<0\end{array}.$(2.44)

The net energy exchange in the snowpack (Q*) is then equal to

$Q^{*}=(Q_n+Q_h+Q_e+Q_g+Q_p).$(2.45)

If Q* < 0, the snowpack is losing energy to the atmosphere (cooling), and some liquid water (if available) may be re-frozen. The amount of energy released to the snowpack (positive value) by re-freezing liquid water is given by

$Q_o=min(-Q^{*},{\rho }_w{l}_{fs}{W}_{liq}^{t+\Delta t}).$(2.46)

The resulting changes in the liquid and ice phases are given by

$W_{liq}^{t+\Delta t}=W_{liq}^{t+\Delta t}-\frac{Q_o}{{-}_w{l}_{fs}}$(2.47)

$W_{ice}^{t+\Delta t}=W_{ice}^{t+\Delta t}+\frac{Q_o}{{-}_w{l}_{fs}}$(2.48)

$W^{t+\Delta t}=W_{liq}^{t+\Delta t}+W_{ice}^{t+\Delta t}.$(2.49)

The negative snowpack temperature, $T_{sp}^{t+\Delta t}$ (associated with its cold content), is then updated from Equation (2.42). If Q* > 0, the snowpack is gaining energy from the atmosphere (heating), and in the process the negative $T_{sp}^{t+\Delta t}$ will increase until it just reaches the isothermal condition ( $T_{sp}^{t+\Delta t}$ → 0). When $T_{sp}^{t+\Delta t}$ becomes positive, it is set equal to zero and Q_o is computed by Equation (2.46) and applied to Equation (2.47) and Equation (2.48) to compute the new liquid and ice components of SWE.

At each time step, the compaction of snowpack, S_comp is based on the present snowpack density ρ_sp (=W^t+Δt / SD^t+Δt), maximum ρ_s,max, and a settlement constant, c_s (Reference Riley, Israelsen and EgglestonRiley et al., 1972), as

(2.50)$S_{comp}=S{D}^{t+\Delta t}{c}_s\left(1-\frac{{\rho }_{sp}}{{\rho }_{s,\max }}\right).$

The depth of snowpack after compaction is the difference between SD and S_comp. ρ_s,max and c_s are manually calibrated such that the simulated SD matches the corresponding snow course data for a given land-use. The effect of rain compaction on snow is also based on Equation (2.50), where SD^t replaces SD^t+Δt when precipitation is in the form of rain.

During melt (positive Q_o) and $T_{sp}^{t+\Delta t}$ is isothermal at 0 °C, water is removed as meltwater (M_ij) when the liquid phase increases beyond the current liquid water holding capacity (LWHC) of the snowpack at the expense of the ice phase, or M_ij is held within the pack when snowmelt first appears at the bottom of the snowpack.

$M_{ij}=W_{liq}^{t+\Delta t}-\left(\text{LWHC}\right)W^{t+\Delta t}$(2.51)

where i is the sub-basin number and j is the land use type, and recommended values for LWHC, a function of snowpack properties and the presence of depth hoar, are 0.02W to 0.05W (US Army Corps of Engineers, 1956), and 0.05W for ρ_sp < 400 kg m⁻³ (Reference Riley, Israelsen and EgglestonRiley et al., 1972). As meltwater contributes runoff at the bottom of snowpack, the new $W_{liq}^{t+\Delta t}$ ,

$W_{liq}^{t+\Delta t}=W_{liq}^{t+\Delta t}-M_{ij}.$(2.52)

The final SWE is computed from Equation (2.49). Routing the meltwater through the snowpack is usually neglected because the routing time for moderately deep snow covers is usually less than the hourly time step of the 1-D energy balance snowmelt model.

Reference Singh, Gan and GobenaSingh et al. (2009) tested the standard and modified temperature index methods, and the 1-D EBM on a small watershed called the Paddle River Basin (PRB) (53º 52′ N, 115º 32′ W), a sub-basin of the Mackenzie River basin, in a semi-distributed approach whereby PRB was divided into units that have similar drainage patterns identified from digital elevation model data. This semi-distributed approach provides a trade-off between modeling resolution, complexity and data availability (Reference Biftu and GanBiftu and Gan, 2001) so that the snowmelt (M_i) for sub-basin i and at each time step is the sum of melt from each land cover, weighted by its corresponding areal fraction ϕ_j as

$M_i=\sum _{j=1}^n{\varphi }_j{m}_{ij},$(2.53)

where n is the total number of land cover classes considered. They found that both EBM and MTI models show good agreements between their simulated and observed values of PRB of 2006/2007 based on ρ_s,max of 250 kg m⁻³ for the deciduous forest (DF) areas (Figure 2.17a) and the coniferous forest (CF) (Figure 2.17b), respectively. Apparently better simulations of the snow depth could be achieved by varying ρ_s,max with time, which is in line with the anticipated increase in snow density with time. ρ_s,max continues to change due to interaction of the snowpack with freshly fallen snow and settlement, and usually reaches the highest value at the end of snow accumulation period. When compared to the snow pillow data collected, EBM reproduced the snow accumulation between January 1, 2007 and March 13, 2007 almost perfectly. After March 13, the model overestimated snow accumulation and similarly lagged the ablation by about 2 days. On the other hand, MTI slightly overestimated snow accumulation up to the beginning of March and then reproduced the late season accumulation and ablation almost perfectly. Similarly, the time series of SWE estimated from the daily snow depth averaged from hourly automated snow depth measurements and the mean measured ρ_s from bi-weekly snow surveys conducted for the BERMS (Boreal Ecosystem Research and Monitoring Sites) sites located within the southern boreal forest of Saskatchewan for 2006/2007 shown in Figure 2.17c, and for 2007/2008 in Figure 2.17d, respectively, are compared with SWE retrieved by the GlobSnow project (Reference LuojusLuojus et al. 2010). Despite of snow interactions with mixed vegetation and relatively dense forest, the BERMS data generally compare well with the GlobSnow data.

Figure 2.17 Comparison of observed SWE and snow depth (SD) to those obtained by EBM (dashed line) and MTI (solid line) of a sub-basin of PRB for 2006/2007 for (a) the deciduous forest (DF) and (b) the coniferous forest (CF) with daily snow pillow data, respectively; and comparison of observed SWE for the BERMS sites located within the southern boreal forest of Saskatchewan for (c) 2006/2007 and (d) for 2007/2008, respectively with the GlobSnow data

(taken from Reference TakalaTakala et al., 2011)

As an energy balance model, EBM accounts for the Q_cc explicitly but as an empirical model, MTI does not account for the Q_cc of the snowpack and so occasional warm air masses occurring in the winter may cause MTI to over-simulate snowmelt runoff. This is the limitation of the temperature index (TI) approach, where without considering the effect of T_g, the model performance could be rather poor in both the calibration and the validation stages, as shown in Figure 2.17c, when the effect of T_g in the MTI model (Equation 2.23) was completely ignored by setting χ to 1 and ψ to 0, making it a standard TI model (Equation 2.21). In order to appreciate better the improvement achieved by introducing T_g and also to perform a fair comparison between TI and MTI, Reference Singh, Gan and GobenaSingh et al. (2009) re-calibrated the m_f of TI (Equation 2.21) based on T_a only. It was necessary to use an artificially low m_f for TI to perform well at the calibration stage (figure not shown), for example, m_f of 0.04 mm hr⁻¹ °C⁻¹ for DF, which is a very low melt factor. However, the model performance decreased considerably when these m_f values were used for all validation stages (Figure 2.17c). Apparently in a Prairie environment where the seasonal snow cover is shallow to moderately deep, where T_g is found to have a fairly strong correlation with net radiation (Q_n), and the onset of major snowmelt usually happens when T_g approaches 0 °C, using both T_g and T_a in a temperature index approach (e.g., MTI) should generally lead to more accurate results than using T_a alone, which is much less data demanding compared to the EBM model.

Reference Singh, Gan and GobenaSingh et al. (2009) further assessed the three snowmelt models in terms of hourly simulated runoff at the outlet of the PRB within the framework of the semi-distributed hydrologic model DPHM-RS (Reference Biftu and GanBiftu and Gan, 2001). At the calibration stage (figure not shown), EBM (R² = 0.85, and RMSE = 1.01) and MTI (R² = 0.79, RMSE = 1.24) performed reasonably well even though in terms of the timing of peak flows, both EBM and MTI show a tendency to somewhat lag the peak flows observed. At the validation stage, they found the results to be less satisfactory (R² = 0.5 for both EBM and MTI) partly because of errors in the lapse rate and gradient used to distribute point temperature and precipitation measurements to sub-basins, errors in the rating curve-discharge relationships due to icing, and the “regulatory” effects of beaver dams in PRB’s streamflow.

Even though the degree-day or temperature index (TI) approach has been a popular tool in modeling the spring snowmelt in temperate climates, it seems that either the data intensive, energy balance or the modified TI approach of Reference Singh, Gan and GobenaSingh et al. (2009) should be applied to watersheds of northern climates. More extensive research on modeling the snowmelt processes in the cryosphere using the MTI approach should be done since the energy balance approach is generally not practical because its data requirements are usually not met except in watersheds chosen for intensive research studies.

Most energy based snowmelt models consider energy balance in the vertical direction only which may not be sufficient in a mountainous watershed where the effects of terrain features can be important. A snowpack may receive radiation reflected from a surrounding slope or shielded from the incoming radiation by adjacent terrains. For example, if the diffused radiation from the atmosphere to a horizontal plane is Q_D, the diffused radiation received by a slope of inclination angle Ź from the atmosphere will be 0.5Q_Dcos²Ź where 0.5cos²Ź is the sky-view factor (Reference Dewalle and RangoDeWalle and Rango, 2008). If there is an adjacent terrain of albedo α, the incoming shortwave radiation reflected by the adjacent terrain and received by this slope will be Q_siα(1 − 0.5cos²Ź).

The effects of forest cover on snowmelt processes can be approximately accounted for by estimating the amount of shortwave radiation penetrated through the canopy, which, similar to the interception of snow by canopy, is a function of the LAI. If the incoming short wave radiation is Q_si, then the amount of radiation reaching the snowpack becomes Q_siexp(−ξLAI) where ξ is the extinction coefficient for shortwave radiation in a forest. Values of ξ depend on the forest types, for example, Reference Baldocchi, Matt, Hutchison and McMillenBaldocchi et al. (1984); Reference Chen, Rich, Gower, Norman and PlummerChen et al., (1997). The incoming longwave radiation, Q_li, for a snowpack under forest cover will be the weighted sum of incoming longwave radiation that penetrated through the overhead canopy and the longwave radiation emitted by the forest.

Two-dimensional energy balance approach

The energy balance model presented above and most if not all energy balance models do not directly account for the effects of the distribution of two-dimensional (2-D) patchy snowcover on the local advective energy of melt processes. This limits their application in a shallow snow cover (<60 cm) environment particularly in late melting periods when beside radiative fluxes, turbulent fluxes should also be accounted for. It has been found that the maximum snowmelt rate occurs when the land is only partially snow-covered, and often when it is slightly less than 60% (Reference Schmidt, Grise and PaceShook and Grey, 1997). Because of lower albedo, the bare ground absorbs a larger amount of solar radiation than the adjacent snow patches. The energy imbalance induces an advective, turbulent transfer of latent and sensible heat from the bare ground to snow patches, enhancing the melt rate. Since advective melting is the greatest along the leading edge of a snowfield, under constant climatic conditions, the melt rate of a patchy snow cover should be related to the perimeter of the patches (Reference ShookShook, 1993), or linearly relate the effective patch size for a regular pattern to the average patch size of a complex snow-cover pattern (Reference EsseryEssery et al., 2006).

Besides radiative, sensible and latent heat fluxes, there can be local advection (turbulent energy) due to large area of bare patches within a snow field, which may significantly alter the energy balance, and becomes increasingly important to melt as the snow cover dwindles. The proportion between radiation and turbulent energy sources to melt depend on the size of the snowfield. The smaller snowpatches will be dominated by turbulent melt throughout the season, or until they disappear. For large snowfields, melting is dominated by radiative melt early in the season and turbulent melt late in the season as they decrease in area (Reference Shook and GrayShook and Gray, 1997). Near the leading edge of an alpine snowfield, Reference Olyphant and IsardOlyphant and Isard (1988) found that advection may contribute more than 30 MJ m⁻² d⁻¹ of melt energy on a very windy day and more than 12 MJ m⁻² d⁻¹ on a relatively windless day. They also found that the corresponding advective energy at 1 km from the leading edge decreases to 5 and 2 MJ m⁻² d⁻¹ on windy and windless days, respectively. Further, the effects of wind on a snowpack are more pronounced on the windward than on the leeward slopes.

The sensible heat and latent heat fluxes are primarily related to wind speed, atmospheric stability, temperature and vapor pressure of air and snow surfaces, when they are considered as 1-D flux elements. However, near the edge of patchy snow cover, such simplification may not give good estimation and thus requires the use of a 2-D model that considers the development of the boundary layer beginning at the edge of the snowpack. A 2-D turbulent diffusion model, proposed by Reference WeismanWeisman (1977), considers this aspect of energy exchange when air flows over a snow cover. The model assumes a steady turbulent flow of warm moist air mass moving from a homogeneous surface onto a ripe snowpack (isothermal at 0 °C), no change in albedo over the snowpack and so no change of the radiation terms from point to point. Also the model assumes that albedo does not vary over the snowpack, and so Q_n remains constant. This assumption implies that only Q_e and Q_h vary over the snowpack. The model, therefore, considers the 2-D aspect of Q_e and Q_h, and quantify the snowmelt M′, due to condensation or sublimation and sensible heat, as,

$M\text{'}=Q_e+Q_h.$(2.54)

The steady-state equations that describe the mean airflow over the snow are continuity, conservation of momentum, sensible heat and water vapor.

$\text{Continuity}:\frac{\partial u}{\partial x}+\frac{\partial w}{\partial z}=0$(2.55)

$\text{X momentum}:u\frac{\partial u}{\partial x}+w\frac{\partial w}{\partial z}=\frac{1}{\rho }\frac{\partial \tau }{\partial z}$(2.56)

$\text{Sensible heat}:u\frac{\partial T}{\partial x}+w\frac{\partial T}{\partial z}=\frac{1}{\rho c_p}\frac{\partial H}{\partial z}$(2.57)

$\text{Vapor}:u\frac{\partial q}{\partial x}+w\frac{\partial q}{\partial z}=\frac{1}{\rho }\frac{\partial V}{\partial z},$(2.58)

where u, w are mean wind components in the x and z directions, T is temperature, q is specific humidity, ρ is air density, c_p is specific heat of air at constant pressure, and $\tau$ , H and V are the turbulent fluxes of momentum, sensible heat, and water vapor, respectively. The molecular diffusion, lateral and forward turbulent diffusion, and the pressure term in x-momentum equation have all been neglected. The problem of airflow over a sudden change in surface temperature and humidity has been solved using the mixing length theory. A vapor diffusion equation is included in the set of conservation equations and a vapor buoyancy term is included in the stability length.

Reference WeismanWeisman (1977) found that a stable atmospheric condition dampens the turbulent diffusion of fluxes, and the melt rate is at a maximum near the leading edge and decreases by one-third approximately 15–20 m from the leading edge. He provides an approximation for the advection of energy that relates the dimensionless melt energy at a point ( ${\hat{M}}^{\prime }$ ) and total average melt over the snowpack ${\hat{M}}^{\prime }$ as a function of dimensionless horizontal downwind distance ( ${\hat{x}}^{}$ ) and dimensionless snowpack fetch ${\hat{x}}^{}$ , respectively.

${\hat{M}}^{\prime }=a{\hat{x}}^{-b}$(2.59)

${\hat{M}}^{\prime }=c{\hat{x}}_o^{-d}.$(2.60)

The constants a, b, c and d depend on the stability parameters, A_*, which is associated with temperature change, and B_*, which is associated with specific humidity change. The snow surface temperature T_s can be retrieved from NOAA-AVHRR, MODIS or Landsat-TM data. Once the dimensionless melt is known, it can be converted to dimensional melt energy M′ by multiplying ${\hat{M}}^{\prime }$ by the energy gradient at the leading edge of the snowpack. Tables 2.2 and 2.3 list some typical values. The above equation underestimates the melt flux for values of $\hat{x}$ less than 10⁴, which comes out to be about 25 m from the leading edge for an average roughness of snow (0.002 m).

The stability parameters, A^* is associated with temperature change and B^* is associated with specific humidity change, given as

${\text{A}}^{*}=\frac{kg{z}_o(T_o-T_{osoil}^{})}{u_{*a}^2{T}_{osoil}^{}}$(2.61)

(2.62)${\text{B}}^{*}=-0.61\frac{kg{z}_o(q_o-q_{osoil}^{})}{u_{*a}^2},$

where k = von Karman constant, 0.4, z_o roughness height (m), $u_{*a}^{}$ = frictional velocity upwind of the leading edge of sow (m s⁻¹), T_o = temperature of snow surface (0 °C), $T_{osoil}^{}$ = temperature of soil surface upwind of leading edge (0 °C), q_o= surface specific humidity at snow surface temperature (assume saturation) and $q_{osoil}^{}$ = surface specific humidity at the soil surface temperature, upwind of leading edge. The knowledge of soil surface temperature upwind of the leading edge of the snow patch can also be obtained from the empirical equations for vapor pressure deficit and psychrometric equations.

Once the dimensionless melt is known, it can be converted to dimensional melt energy M′ by multiplying ${\hat{M}}^{\prime }$ by the energy gradient at the leading edge of the snowpack.

$M\text{'}=\hat{M}\text{'}\left[\rho u_{*a}^{}{c}_p\right(T_o-T_{osoil}^{})+\rho u_{*a}^{}{l}_v(q_o-q_{osoil}^{}\left)\right],$(2.63)

where ρ is air density (kg m⁻³), and c_p is specific heat of air (J kg⁻¹ °C⁻¹).

Reference Olyphant and IsardOlyphant and Isard (1988) found that Equation (2.27) does not provide a good estimate near the leading edges of discontinuous melting snow surfaces where the advective heat contributes greatly to the energy balance of late-lying snow, and they introduced a modified version of Reference WeismanWeisman (1977)’s 2-D boundary layer model. The model pertains to a unit width of snow surface aligned with the prevailing wind. The horizontal pressure gradient is assumed to be negligible, as is the forward component of turbulent diffusion. Time derivatives have been added to the momentum and heat flux equations in an effort to simulate time-dependent responses to the diurnal variation of ambient forcing. The relevant boundary layer equations include one continuity and three conservation equations (momentum, sensible heat and latent heat). The estimation of fetch length, which fluctuates with the wind direction, and the irregular shape of snow patches can be attempted empirically. So far, there have been limited attempts at 2-D snowmelt modeling, which is predominantly 1-D.

Intercomparison of snowmelt models

Since the intercomparison project by the World Meteorological Organization (WMO) on snowmelt runoff models conducted in 1986, some other model intercomparison studies have been conducted, ranging from small scale, individual studies, to international effort, of which the latest and most comprehensive is the SnowMIP2 (Reference Rutter and EsseryRutter et al., 2009). Reference Singh, Gan and GobenaSingh et al. (2009), for example, compared three semi-distributed snowmelt models, namely, a temperature-Index or degree-day model based on air temperature (T_a), a modified temperature index model based on T_a and near ground-surface temperature (T_g), and an energy-balance model that considers liquid and ice phases separately. For the Canadian Prairies where snowpack is shallow to moderately deep, and winter is relatively severe, apparently the modified temperature index model based on T_a and T_g can perform comparably well with the more complex energy balance model, because the advantage of using both T_a and T_g is partly attributed to T_g showing a stronger correlation with solar radiation than T_a during the spring snowmelt season, and partly to the onset of major snowmelt which usually happens when T_g approaches 0 °C.

In the SnowMIP2 project, thirty-three snowpack models of varying complexity were compared across a wide range of hydrometeorological and forest canopy conditions, for up to two winter snow seasons applied to five northern sites in the NH: Alptal of Switzerland; Boreal Ecosystem Research and Monitoring Sites (BERMS) of Canada; Fraser of USA; Hitsujigaoka of Japan and Hyytiȁlȁ of Finland (see Figure 2.11). As expected, Reference Rutter and EsseryRutter et al. (2009) concluded that it is easier to model SWE at open than at forested sites, where precipitation phase and duration of air temperatures >0 °C are key factors to the divergence and convergence of modeled estimates of the sub-canopy snowpack, and more consistent results are found between open than forested sites (Figure 2.18). Furthermore, a model that estimates snowpack accurately at a forest site may not do well at an open site and vice versa. Calibrated models at forest sites perform better than uncalibrated models, though the benefits of calibration may not translate to subsequent years, nor to open conditions, which is expected since forests have large influences on snow dynamics and many recent models have included vegetation canopy (e.g., Reference Niu and YangNiu and Yang, 2004; Reference Bartlett, MacKay and VerseghyBartlett et al., 2006).

Reference ChenChen et al. (2014) evaluated the snowpack simulated by six LSMs against one year of SWE data at 112 SNOw TELemetry (SNOTEL) sites in the Colorado River Headwaters region and 4-year flux tower data at two AmeriFlux sites. While all models captured the key characteristics of the seasonal SWE evolution well at all SNOTEL sites, no single model was the best capturing the combined features of the peak SWE, the timing of peak SWE, and the length of snow season. The models responded differently to different forest coverage, and deficient in the treatment of snow albedo and its cascading effects on surface energy deficit, surface temperature, stability correction, and turbulent fluxes, for all models substantially overestimated (underestimated) radiative flux (heat flux). There are also significant intermodel differences in modeling snowmelt and sublimation efficiency, for models with high snow accumulation and melt rates were able to reproduce the observed seasonal evolution of SWE. It seems the parameterization of cascading effects of snow albedo and below-canopy turbulence and radiation transfer is critical to simulate the SWE and the winter land-atmosphere interactions correctly.

Figure 2.18 Observed SWE (black dots) and the average of thirty-three modeled estimates for (a) (gray line) of forest site at Hyytiȁlȁ (61° 51′ N, 24° 17′ E) of Finland (taken from Reference Rutter and EsseryRutter et al., 2009), and (b) (green line) of open and forest site at Alptal (47° 03′ N, 8° 48′ E) of Switzerland

(taken from Reference EsseryEssery et al., 2009)

An integrated approach to modeling snow accumulation and ablation processes

With the availability of geographical information systems (GIS), digital terrain elevation (DTED) data, snow products such as that of NOHRSC, re-analysis data, spatially distributed remotely sensed (RS) data to augment our limited ground-based, point observations, and an exponential growth in computing power, it will be desirable to integrate distributed, physically based snow accumulation and ablation processes with RS data, and ground measurements of snowpack to better model snowmelt processes. The idea is to progress from empiricism (e.g., degree-day method) to a discipline of applied science, and to model hydrological processes from measurable causative factors. The building and applications of distributed snow models depend largely on our ability to retrieve useful snow accumulation and snowmelt information of reasonable resolution from RS data in order to augment limited ground measurements.

Ideally, such distributed models should be developed in the direction of sub-grid parameterizations rather than the traditional quest for refining the resolution of small-scale parameterizations. This means finding a trade-off between the resolution of processes modeled, the types of data available and the information contained in the data, and the accuracy required. Otherwise, such models may be difficult to apply because of excessive data demand and the difficulty in obtaining the parameters required at all grid elements. On the other hand, some distributed snow models are based on rectangular or square grid elements of constant size without considering sub-grid parameterization, irrespective of the terrain features. This means that process descriptions may become “artificial,” since nature does not behave as a system of symmetrical grids placed side by side.

There have also been developments in land surface schemes using existing Surface Vegetation Atmosphere Transfer Schemes (SVATS), such as the Land Data Assimilation

Systems (LDAS) (http://ldas.gsfc.nasa.gov/) which are forced with gauged precipitation, RS data, radar precipitation, and output from numerical weather prediction models. Then in situ or remotely sensed measurements of LDAS storages (such as snow), water and energy fluxes will be used to validate and constrain the LDAS predictions using certain data assimilation techniques. However, one of the possible drawbacks of LDAS is the discretization of grids symmetrically without including sub-grid parameterization. Even with subgrid parameterizations in distributed models, ideally it may be more desirable to develop semi-distributed models that discretize river basins according to terrain characteristics, and that are designed to model snow, water, and energy dynamics with practical details under the forcing of fluxes and the influence of terrain and vegetation characteristics, in a framework that mimics nature as much as possible.

Snowmelt in Greenland and Antarctic

About 20 years of SSM/I passive microwave data has been used to monitor the regional snowmelt of Greenland since 1988 (Reference TedescoTedesco, 2007; Reference TedescoTedesco et al., 2008). Defining a melting index (MI) as melting area multiplied by melting days, 2005 was the year with the highest MI, followed by 2002, 1998, 2004, and 2007. In 2007, areas higher than 2,000 m in southern Greenland experienced about 30 more days of melting than the study period of 1988–2007, and 53% higher MI than the average, even though overall Greenland had a MI about 20% higher than the average. In contrast, in 2008 northern Greenland experienced record melting, and they attributed record melting of snow in different parts of Greenland to warmer surface temperature.

Snow cover extent

Based on the 12-month running means, the monthly SCE over NH lands (including Greenland), Eurasia and NA for about 40 years (1966–2006) mainly show negative snow cover anomalies since the late 1980s, with the exception of 1996 and 2002 for NH and Eurasia, and 1997 for Eurasia and NA (Reference RobinsonRobinson, 2008). The SCE over NA has been shown to increase during autumn and early winter (November–January), but decrease over early spring (Reference Frei and RobinsonFrei and Robinson, 1999; Reference Dyer and MoteDyer and Mote, 2006), which indicates an earlier onset of the spring snowmelt. Reference Burn and KokeljCallaghan et al. (2011b) found Arctic snow cover showing different regional responses to climate warming and increasing winter precipitation that has characterized the Arctic climate in recent decades, with the largest and most rapid decreases in SWE and snow cover duration (SCD) observed over maritime regions of the Arctic with the highest precipitation amounts. Eurasian Arctic region has experienced larger declines in SCD (12.6 days), compared to the North American Arctic region (6.2 days) between 1982 and 2011 (Reference Barichivich, Briffa and MyneniBarichivich et al., 2013). Further, the North American sector of the Arctic exhibiting decreases in SCE and snow depth when in situ observations became available since 1950s, while widespread decreases in SCE are only apparent over Eurasia in the Arctic after 1980s, but snow depths are increasing in many regions of Eurasia. Warming and more frequent winter thaws are contributing to changes in snowpack structure.

Even though the SCE in January 2008 of 50.1 × 10⁶ km² was the largest recorded SCE for January in the last 40 years, the 2008 annual SCE over NH lands averaged 24.4 × 10⁶ km², which is 1.1 million km² less than the 39-year average, putting 2008 as the fourth least extensive snow cover on record, and the 12-month running means ran below the long-term average throughout 2008, following the generally lower than average extents in 2007. The SCE in January has been declining consistently since the notable minimum in the late 1980s and early 1990s over Eurasia, even though in North America, the SCE anomalies have somewhat rebounded over the course of 2008 from a 2005–2007 minimum. From a snow model that generated monthly SCE of NH from 1905 to 2002, Reference McCabe and WolockMcCabe and Wolock (2010) found a substantial decrease in the March SCE of NH since 1970s attributed to an increase in the mean winter temperature and a contraction of the circumpolar vortex and a poleward movement of storm tracks, which resulted in decreased precipitation (and snow) in the low- to mid-latitudes and an increase in precipitation (and snow) in high latitudes.

Are these generally negative snow cover anomalies since the late 1980s an indication of the warming trends detected for North America, Europe, and Asia in NH? From the Mann–Kendall’s trend analysis of SCE of NH over 1972–2006, Reference Déry and BrownDery and Brown (2007) found significant declines in SCE during spring over North America and Eurasia, with lesser declines during winter and evidence of a poleward amplification of decreasing SCE trends during spring over Eurasia and North America, which is consistent with an enhanced snow-albedo feedback over northern latitudes that act to reinforce an initial anomaly in the cryospheric system.

In recent years, the fraction of precipitation in mid- and higher-latitude regions falling as snow has decreased, the fall accumulation has been occurring later while the spring ablation earlier, and the northern hemisphere’s (NH) SCE and snow depth have decreased, which would have significant socioeconomic consequences as more than one-sixth of the world’s population depend on meltwater for their water supply. Climate warming has been the dominant factor in the observed decrease in winter SCE and earlier spring ablation over North America (NA) (Reference Dyer and MoteDyer and Mote, 2006), given air temperature anomalies over NH’s mid-latitude land areas has been shown to explain about 50% of the observed variability in SCE (Reference Brown and RobinsonBrown and Robinson, 2011).

Numerous observations have shown that the SCE of NH has been decreasing rapidly in recent decades, especially in spring, and in warmer locations where SCE are closely linked to temperature variability, for example, March–April SCE of NH is decreasing at 3.4%±1.1% per decade (1979–2005) (Reference Brown and RobinsonBrown and Robinson, 2011; Reference Hernández-Henríquez, Déry and DerksenHernández-Henríquez et al., 2015). Over 1967–2012, the annual mean SCE has significantly decreased, with the largest change, −53% [−40 to −66%], occurred in June. Over 1922–2012, a much smaller 7% [4.5–9.5%] decline for March–April was found. Majority of climate stations show decreasing trends, and stations at lower elevation or higher average temperature were the most liable to show a decrease.

The satellite-based visible SCE of National Oceanic and Atmospheric Administration (NOAA) climate data record (CDR) developed from weekly Northern Hemisphere SCE data span from October 4, 1966 to the present. Satellite data incorporated into the NH SCE CDR at various starting dates include platforms of the Environmental Science Services Administration (ESSA) series, NOAA Polar-orbiting Operational Environmental Satellites (POES), NOAA Geostationary Operational Environmental Satellites (GOES), the Defense Meteorological Satellite Program (DMSP) series, the Meteosat series, Geostationary Meteorological Satellites (GMS), Aqua, Terra, and the Multi-functional Transport Satellites (MTSAT) (Reference Estilow, Young and RobinsonEstilow et al., 2015). Starting in June 1999 the CDR is completely derived using the Interactive Multisensor Snow and Ice Mapping System (IMS) involving a diverse set of products: satellite imagery, snow and ice analysis maps, National Centers for Environmental Prediction (NCEP) model data, and surface observations, and improved discrimination between snow and cloud covers. With an 89 × 89 Cartesian grid laid over an NH polar stereographic projection and a cell size of 190.6 km × 190.6 km, cell areas range from ~10,700 km² near the Equator to ~41,800 km² near the pole.

Despite of many up and down cycles in snowfall that had occurred over the NH in last five decades, Figure 2.19 of the Rutgers University Snow Lab based on this CDR data set, the NH SCE exhibits an overall increasing trend since the late 1980s, in contrast to the overall decreasing trend observed between late 1960s and 1980s. The increasing trend in the SCE of the NH observed in the last three decades are attributed to increasing SCE in the Fall and winter seasons, especially the Fall season, which together had offset the overall decreasing trend in the spring SCE of the NH shown in Figure 2.19d. Trend estimates for annual precipitation of the 1951–2008 time series of CRU TS 3.10.01 (updated from Reference Mitchell and JonesMitchell and Jones, 2005) are 5.82 ± 2.72 and 1.13 ± 2.01, for GHCN V2 (updated through 2011; Reference VoseVose et al., 1992) are 4.52 ± 2.64 and 1.39 ± 1.98, and for GPCC V6 (Reference BeckerBecker et al., 2013) are 2.69 ± 2.54 and 1.50 ± 1.93 mm yr⁻¹, for both latitudinal bands 30° N–60° N and 60° N–90° N, respectively (Reference HartmannHartmann et al., 2013). Even though confidence in precipitation change is medium for the years after 1950 because of insufficient data, the non-significant positive annual precipitation trends over both latitudinal bands of the NH should have contributed to increasing trends in the Fall and Winter SCE of NH of the last five decades. However, likely because of more rainfall than snowfall in spring and earlier onset of spring snowmelt, the Spring SCE of NH generally show decreasing trends despite of positive annual precipitation trends at latitudinal bands north of 30° N.

Figure 2.19 (a) Twelve-month running anomalies of monthly snow extent plotted on the seventh month using values from November 1966 to December 2018 for NH (https://climate.rutgers.edu/snowcover/chart_anom.php?ui_set=0&ui_region=nhland&ui), and seasonal snow cover extent of NH for (b) Fall, (c) Winter, and (d) Spring. The mean 1981–2010 SCE for the NH in Fall, Winter, and Spring are 18,912, 45,480, and 29,788 km², respectively

Snow water equivalent

SWE of snowpacks, the depth of liquid water produced by a snowpack after complete melting, can be estimated from in situ depth (d) and density (ρ_s) in kg m⁻³ measurements, SWE = d(ρ_s/ρ_w), where the density of water ρ_w is about 1,000 kg m⁻³. So far passive microwave radiometry is the only spaceborne technique proven useful in extracting SWE data (Reference Walker and SilisWalker and Silis, 2002). SWE has been estimated from the brightness temperature (TB) of microwave frequency and polarization (H or V) of the TIROS-N sensor, the Scanning Multichannel Microwave Radiometer (SMMR) from the Nimbus 7 satellite with data between October 25, 1978 to August 20, 1987; the Special Sensor Microwave/Imager (SSM/I) data from the Defense Meteorological Satellite Program (DMSP) satellite with data available since September 7, 1987 (Reference Armstrong, Brodzik, Savoie and KnowlesArmstrong et al., 2006), and succeeded by the Special Sensor Microwave Imager/Sounder (SSMIS) in 2005, and the Advanced Microwave Scanning Radiometer‐Earth Observing System (AMSR-E) data from the Aqua satellite with data available between 2002 and 2011 or the Advanced Microwave Scanning Radiometer 2 (AMSR2) on the JAXA GCOM W1 spacecraft, launched May 18, 2012 and is still operating. Limitations to such SWE data are a coarse resolution (typically 25 km), underestimated in forests which mask the snow cover and vegetation suppresses the scattering signal, as does liquid water (wet snow). From retrievals of SMMR-SWE over Canada, Reference Dong, Walker and HouserDong et al. (2005) concluded that for practical applications, the uncertainty of such SWE data is acceptable for regions more than 200 km away from large open water bodies, daytime air temperatures lower than 2 °C, and SWE values are above 100 mm.

Reference MeisterMote et al. (2005) detected widespread declines in modeled April 1 SWE over the Western United States and Canada between 1950 and 1997. Reference Gan, Barry, Gobena and RajagopalanGan et al. (2013) detected significant decreasing trends in SWE in NA over 1979–2007, which are more extensive in Canada than in the United States where such decreasing trends are mainly found along the Rocky Mountains. However, because of uncertainties associated with SWE data retrieved from passive microwave data, results obtained for mountainous or forested areas of NA, the tundra in Arctic Canada with many frozen lakes or snow packs that consist of depth hoar and wind slab, should be treated with caution. To assess the possible impact of climatic change to the snowpack of North America, Reference Gan, Barry, Gobena and RajagopalanGan et al. (2013) analyzed the trends of temperature and precipitation data and then the SWE-air temperature and SWE-precipitation relationships. Using the gridded, 2-m air temperature data of the North American Regional Reanalysis (NARR) (Reference MesingerMesinger et al., 2006) and that of the University of Delaware (Reference Willmott and RobesonWillmott and Robeson, 1995), they detected statistically significant temperature trends for 1979–2007 (January–April) mainly in the southern states of USA. They also detected significant decreasing precipitation trends from the University of Delaware and the NARR data but there are limited agreements between the locations of significant decreasing precipitation trends and significant decreasing SWE trends. However, extensive areas of negative correlations between SWE and temperature exist both across the United States and Canada, and the distribution of these areas of negative correlation closely follows the areas of the decreasing trends detected from the SWE data, which for Canada is mainly east of the Canadian Rocky Mountains, while for the United States is mainly on the American Rocky. Apparently, extensive decreasing trends in SWE data of passive microwave detected in Canada and parts of USA are caused more by increasing temperatures than by decreasing precipitation. Higher air temperature means more rainfall and less snowfall especially in areas and seasons where the average air temperature is close to 0 °C, earlier onset of spring snowmelt, and generally less snowpack.

From SWE estimated from satellite measurements for October to April over 1979–2007 in NA, Reference Gan, Barry, Gobena and RajagopalanGan et al. (2013) found that about 25–30% of the pixels covered with snow showed statistically significant negative trends, compared with 5–10% that showed increasing trends. The overall mean trend magnitudes estimated are about −0.4 to −0.5 mm a⁻¹ which means the overall SWE of NA had decreased by about 10–13 mm or possibly more, or in terms of snow depth could range from about 4 to 13 cm (depending on the snow density) over 1979–2007. The detected changes should have significant impacts on the spring snowmelt of, say, the Canadian Prairies and the Washington Cascades. Using a combination of satellite data, in situ data and land surface model simulations, Reference Park, Yabuki and OhataPark et al. (2012) found negative trends in snow depth more pronounced in North America than in Eurasia between 1948 and 2006. Reference Dyer and MoteDyer and Mote (2006) found significant (>1 cm yr⁻¹) negative regional snow depth trends in central Canada along a southeast line from the Yukon Territory to the Great Lakes using in situ snow depth data for 1960–2000 for North America.

Global Cryosphere Watch (2015) has compiled a comprehensive snow data set from passive-microwave, in situ, and analysis/reanalysis products for non-alpine regions of the NH. The data set at 25 km resolution is available at daily, weekly and monthly time steps since 1979. By combining different data sources, GlobSnow is probably better than SWE data derived from stand-alone remote sensing data sets (Reference Frei, Tedesco, Lee, Foster, Hall, Kelly and RobinsonFrei et al., 2012).

Based on the non-parametric Mann–Kendall test at 0.05 significant levels, more snow-covered pixels of the monthly GlobSnow SWE data set of NH for November–April (1988–2017) show decreasing than increasing trends (Figure 2.20). From the total number of snow-covered pixels analyzed, up to 15.5% (7.7%) of the pixels show statistically significant decreasing (increasing) trends. December has the largest SCE and the greatest percentage of statistically significant decreasing trends, of which the majority are located north of 55° latitude. April exhibits the greatest percentage of statistically significant positive trends and most of these pixels are located in Asia.

Figure 2.20 Northern Hemisphere monthly snow cover for December–March of GlobSnow SWE data set of 1988–2017 that exhibit increasing (red) and decreasing (black) trends statistically significant at the 0.05 level using Mann–Kendall test

Figure 2.21 shows spatial distributions of grids of the Northern Hemisphere with significant SWE (5-day averages) trends at p < 0.05 significant level over 1988–2017 using the Mann–Kendall test. There are more statistically significant negative than positive SWE trends especially across Canada, the high Arctic, and Europe, with scatterered positive trends in Russia.

Figure 2.21 Spatial distributions of grids of the Northern Hemisphere GlobSnow data set with significant SWE (5-day averages) trends at p < 0.05 significant level over December–March of 1988–2017 using the Mann–Kendall test, with more statistically significant negative than positive SWE trends especially across Canada, the high Arctic, and Europe, while scatted positive trends in Russia

As expected, based on the Mann–Kendall’s test applied to the ERA-Interim reanalysis data, Figure 2.22 shows predominantly warming trends over the Northern Hemisphere in the last 30 years (1988–2017). On the other hand, some scattered cooling trends are still detected in the northern Pacific Ocean in January and March, over some eastern parts of Siberia in February, and a small part of the north Atlantic Ocean in March.

Figure 2.22 The spatial distribution of grids with statistically significant (p < 0.05) temperature Pentad (5-day averages) trends by the Mann–Kendall test in the Northern Hemisphere based on the 1988–2017 ERA-Interim reanalysis temperature data

The Spearman’s rank correlation between SWE and air temperature across the Northern Hemisphere is predominantly negative, which is expected because under the effect of climate warming, there should be an overall less snowpack partly attributed to more rainfall than snowfall, if the total precipitation will remain unchanged or even higher (Figure 2.23).

Figure 2.23 Spatial distributions of grids of the Northern Hemisphere with significant Spearman’s rank correlation between SWE and air temperature over 1988–2017 at p < 0.05 significant level

The Self Organizing Map (SOM) (Kohonen, 2001) is an artificial neural network developed to uncover an underlying structure in a data set through an unsupervised learning to produce a 2-D array of nodes, called a map, that are organized in such a way that similar items are placed close to each other. By applying the SOM and K-means clustering analysis, Northern Hemisphere grids were grouped into 20 different clusters based on the SWE niveograph of November–March for 1988–2017, for example, one grid can be assigned to a different cluster in different years. Figure 2.24 shows that grids in snow-dominated areas of Canada, Siberia, and northern Europe are generally assigned to more clusters (up to 20) than other regions of NH that also experience snowfall.

Figure 2.24 (a) Number of different clusters derived from the SOM and K-means clustering analysis that each grid in Northern Hemisphere is assigned over 1988–2017, and an example of (b) the SWE in boxplots for all grids of NH assigned to Cluster #5

Figure 2.25 shows some limited statistically significant Spearman’s rank correlation between SWE of the Northern Hemisphere and large-scale climate anomalies such as AO, PNA, NAO, and PDO, respectively. There are more negative (black) than positive (red) correlation especially between PNA and SWE, which means less snowpack during the warm than the cold phase of PNA. For example, in western Canada, Reference Gan, Gobena and WangGan et al. (2007) show that El Niño (La Niña) leads to a 14% decrease (20% increase) in the mean winter precipitation, strong positive (negative) PNA leads to a 12% decrease (9% increase) in mean winter precipitation, and PDO is associated with an 8% decrease (9% increase) in mean winter precipitation. The detected teleconnections could occur at interannual or interdecadal levels depending on the climate anomaly, and their strength changes in time and space.

Figure 2.25 The spatial distribution of grids with significant (p < 0.05) positive (red) and negative (black) Spearman’s rank correlation between cluster areas, AO, PNA, NAO, and PDO, respectively

Climate indices significantly correlated with any of the first four PCs of the SWE of NH at 0.05 significance level are shown in Table 2.4. Apparently AO, NAO and PNA exert a stronger influence on the snowpack of the NH between January and April. PC1 is more strongly correlated with AO and NAO, while PDO and PNA are more correlated with PC2 to PC4 in some months. Therefore, some of the observed changes in SWE are partly attributed to the influence of large-scale climate anomalies but climate warming has undoubtedly played a bigger role, especially in the detected negative trends of SWE of NH.

Reference Wu, Che, Li, Wang and YangWu et al. (2018) found slower snowmelt rates over the entire NH over 1980–2017, but with higher ablation rates in locations with large SWE, because of the reduction of SWE in deep snowpack regions. However, moderate and high snow ablation rates show a decreasing trend. Reference Kapnick and HallKapnick and Hall (2012) detected significant losses of spring mountain snowpack at western United States over past several decades. For Canada, there has been extensive decreasing snow depths and snow cover duration and extent since the mid-1970s, with the largest declines in western Canada and proportionally greater changes later in winter and spring (Reference DeBeer, Wheater, Carey and ChunDeBeer et al., 2016). Reference Berghuijs, Woods and HrachowitzBerghuijs et al. (2014) show that in conterminous USA, catchments with a higher fraction of precipitation falling as snow tends to have higher mean streamflow, which is likely to decrease in catchments that experience significant reductions in the fraction of precipitation falling as snow because of a warmer climate.

Most studies also show negative trends in snow depth and snow duration over past decades in the mountain cryosphere of Europe (Reference Beniston, Farinotti, Stoffel, Andreassen, Coppola and EckertBeniston et al., 2018), with less pronounced changes at high elevations (Reference Terzago, Fratianni and CremoniniTerzago et al., 2013). Spring SWE shows a decreasing trend in Alps (Reference Marty, Tilg and JonasMarty et al., 2017). However, there were positive trends of maximum snow depth and SWE in higher and colder parts of the Fennoscandian Mountains although it turns out to be negative trends in recent years (Reference Kivinen and RasmusKivinen and Rasmus, 2015). Reference Matti, Dahlke, Dieppois, Lawler and LyonMatti et al. (2017) show that the flood seasonality for snowmelt‐dominated Scandinavian catchments have changed over the twentieth century with statistically significant decreasing (increasing) trends in summer (winter and spring) maximum and mean daily flows in some catchments. Changes in annual flood occurrences generally point toward a shift from snowmelt‐dominated to rainfall‐dominated flow regimes in some regions, with flood peaks due to spring snowmelt consistently occurring earlier, often at the expense of the much needed summer streamflow for growing crops (Reference Kerkhoven and GanKerkhoven and Gan, 2011).

There are other cryospheric observations also support the effect of a warmer climate in the recent past. For example, the number of snow days in Switzerland decreased abruptly at the end of the 1980s. Reference MartyMarty (2008) found that records at 34 long-term stations between 200 m and 1,800 m asl for 1948–2007 show an unprecedented series of low winter snow in the last 20 years. The abrupt change in 1988 resulted in a loss of 20–60% of the total number of snow days with no clear trend since then. The decrease is shown to be correlated with an increase in winter temperatures.

SCE over North America has been shown to increase during autumn and early winter (November-January), although it also decreased in early spring (Reference Frei and RobinsonFrei and Robinson, 1999; Reference Dyer and MoteDyer and Mote, 2006), which indicates an earlier onset of the spring snowmelt. Some parts of the Canadian Arctic stretching from western Hudson’s Bay to the North Slope of Alaska shows a persistent pattern of high SWE (Reference Armstrong and BrodzikArmstrong and Brodzik, 2002; Reference Andreadis and LettenmaierAndreadi and Lettenmaier, 2006). However, as noted earlier, SWE data derived from passive microwave for the tundra in Arctic Canada should be treated with caution (Reference DerksenDerksen et al., 2009). Less snowfall has been observed in the lower Missouri River Basin (Reference BensonBerger et al., 2002) and in New England (Reference HuntingtonHuntington et al., 2004).

Reference MoteMote et al. (2005) found evidence of declining snowpack in the western United States, except in the southern Sierra Nevada. Reference Bedford and DouglassBedford and Douglass (2008) analyzed daily SWE from 28 SNOTEL stations in the Great Salt Lake Basin for 1982–2007. They found an advance of about 15 days in the date of peak SWE as well as (less robust) evidence of a decrease in peak and April 1 SWE amounts.

Through three spring indicators – lilacs, honeysuckles, and streamflow – Reference CayanCayan et al. (2001) found earlier onset of the spring season by up to three weeks in the western North America since the 1970s. By simulating the snow energy balance to climatic changes projected by nine regional climate models to the end of the twenty-first century in the Pyrenees, Moreno et al. (2008) concluded that the most significant changes to future snowpack processes are related to temperature. Reference TedescoTedesco (2007) concluded that record melting of snow in southern Greenland in 2007 was caused by higher surface temperature. The above documented observations about the cryosphere agree with the enhanced increasing trends of surface air temperature anomaly changes over other continents, especially enhanced temperature trends occurred in the last several decades since 1960s reported by IPCC (2013) (see Figure 2.26 adapted from Reference Jones, Stott and ChristidisJones et al., 2013). In Figure 2.26, multi-model means are shown as thick lines, and 5–95% ranges shown as thin light lines, and as black lines for HadCRUT4. Mean temperatures are shown for Antarctica and six continental regions. Temperatures are shown with respect to 1880–1919 for all regions apart from Antarctica where temperatures are shown with respect to 1950–2010.

Figure 2.26 Global, land, ocean, and continental annual mean temperatures anomaly with respect to 1880–1919 for CMIP3 and CMIP5 historical (red) and historicalNat (blue) simulations, and Hadley Centre/Climatic Research Unit gridded surface temperature data set 4 (HadCRUT4, black). Weighted model means shown as thick dark lines and 5–95% ranges shown as shaded areas

(taken from Reference Jones, Stott and ChristidisJones et al., 2013)