The generation and release of meltwater during the spring snowmelt season can be delayed because of spring storm episodes with snow accumulation and/or sustained subfreezing temperatures. The delayed release of snowmelt often extends beyond the particular storm event because of changes to the internal state of energy in the snowpack that prevents transmission of meltwater. Following a storm, two energy deficits internal to the snowpack must be overcome before surface melt can drain and exit the snowpack: 1) cold content created by heat lost during the episode must be removed and 2) dry pore space must be filled with liquid water to residual saturation. This study investigates the role of these two processes in spring snowmelt following past storm episodes in western Montana. The analysis addresses ~10 yr of historical snowpack and air temperature data from 33 stations in the Columbia River headwaters. Results indicate that the addition of pore space has a greater impact on delaying snowmelt than does the addition of cold content, with snow accumulation events responsible for 86% of the collective energy deficit imposed on the snowpack during storm episodes. Nearly all refreezing events occur within one month of peak snowpack, but accumulation events are common up to 50 days later. Under standardized conditions representing clear weather during the spring season, these energy deficits could all be overcome in a matter of hours.
As atmospheric conditions warm during the spring snowmelt season, the energy gains by the snowpack surface outweigh the energy losses. Heat is transferred through the snowpack by diffusive processes and by advection of latent heat, as meltwater generated at the surface percolates downward either along a uniform wetting front or as discrete “pipes,” and refreezes (e.g., Marsh and Woo 1985). The added heat initially removes cold content and drives snow grain metamorphism. Once the snow becomes isothermal at the melting temperature, additional heat serves to melt ice and fill grain pores with liquid water. Pores must fill with liquid water to the residual saturation (irreducible water content) before percolating meltwater can overcome surface tension forces and be transmitted through the snowpack (Colbeck 1976). The residual saturation varies by texture, grain size, and grain shape (Colbeck et al. 1990), with normal ranges from 2.5% to 8.5% by volume (Colbeck 1974, 1976; Coleou and Lesaffre 1998; Lemmela 1973; Singh et al. 1997). When sufficient energy has been added to the snowpack so that it is both isothermal at 0°C and its pore space has been filled to residual saturation, the snow is considered “ripe” and capable of transmitting meltwater and producing runoff (e.g., Lundquist et al. 2004).
During the spring season, the melting snowpack is often subjected to a continuum of snow accumulation and/or sustained subfreezing temperature events that inhibit or stop the generation and release of snowmelt (e.g., Fig. 1). These late season events can have the immediate impact of reducing energy input to the snowpack, thereby shutting down snowmelt. However, these events can also have lasting impacts caused by the alteration of two aspects of the snowpack’s energy state: 1) dry pore space can be added in a layer of new snow pore space that must be filled to residual saturation (Anderson 1976; Colbeck 1976) before surface melt can drain through the snowpack and 2) the snow can lose energy and generate cold content that must be overcome before meltwater can be generated and transmitted through the snowpack. Another important change caused by spring snow accumulation events is a substantial increase to snow albedo caused by new snow on the surface. While this effect impacts the ability of the snowpack to absorb solar radiation and therefore plays a key role in snowmelt generation, this alone is not considered a change to the internal energy state of the snowpack.
Spring weather events can create a spring energy hurdle (SEH) resulting from the combination of the energy barriers created by the residual saturation constraint Qrs and added cold content Qcc. These energy deficits can be described in megajoules per square meter (MJ m−2). Before snowmelt can again be generated and transmitted through the snowpack, sufficient energy must be added to the snowpack to overcome the sum of these individual hurdles. The rate at which the SEH is overcome is dependent on the ensuing energy balance (dictated by snow properties, site properties, meteorology, etc.) and the magnitude of the SEH.
The timing of snowmelt runoff is influenced by SEH events of varying frequency and magnitude. While the physical processes governing these energy changes are well understood, their frequency and magnitude in the natural world has not been explicitly documented, particularly not in western Montana. The potential exists for SEH events to create substantial delays in the timing of spring runoff, but their role in snowmelt timing is unknown. This limits our ability to interpret historical snowmelt timing and changes therein and to fully understand and project snowmelt timing for purposes of flood forecasting and water resources. Here, we investigate historical SEH events occurring in the western Montana snowpack with the purpose of elucidating current processes and the sensitivity of such events to climate parameters. The aim of this research is to 1) quantify the frequency and magnitude distribution of historical SEH events in western Montana and 2) partition the relative contribution of added pore space and cold content to the distribution of observed SEH events. These results will be used to interpret the importance of SEH events in delaying the timing of snowmelt runoff.
To investigate the role of SEH events in delaying snowmelt in western Montana, we make simplified calculations of the magnitude of historical events using time series of air temperature, snowpack density, and snowfall. In both refreezing and accumulation scenarios, Qrs and Qcc are calculated separately and then summed to create the total energy hurdle Qh so that
Accumulation events primarily add Qrs, but they can also add a Qcc component if the snow falls at subfreezing temperature. Similarly, freezing events add Qcc, but they also add a Qrs component due to freezing of pore water. Before meltwater transmission can resume following an event, energy input to the snowpack must exceed Qh.
Because of highly limited historical data for Montana’s snowpack, we use a method to approximate these quantities that is based on simplified physics but is feasible given the data constraints explained below. In cases where the snowpack is well instrumented, estimates using methods that capture the physics in more detail are certainly possible. However, as no prior work has addressed this topic, our analysis provides valuable new insight into changes in the internal energy state of the snowpack caused by spring storm events despite the limitations.
Data for our analysis are supplied by 33 Snowpack Telemetry (SNOTEL) sites that are operated and maintained by the National Resource Conservation Service (NRCS). All sites used in this study are located west of the Continental Divide in Montana (Fig. 2) and ranged in elevation from 1295 to 2515 m. Western Montana is characterized by a large relief where elevations range from 560 m in the lowest valleys to roughly 3000 m at the highest peaks. As meteorological stations in the mountains are extremely rare, the SNOTEL network is the primary source of both snow and atmospheric data in this region. The distribution of SNOTEL stations in western Montana is generally representative of moderate elevations where there is, on average, one station for every 2000 km2 of mountainous terrain. Further, because of their scarcity at high elevations, SNOTEL sites poorly sample Montana’s high-elevation snowpack (Gillan et al. 2010). Hence, while our sampling includes every available site suitable for analysis, we note that data are sparse and potentially biased with respect to their spatial distribution and sampling of the mountain snowpack.
SNOTEL sites measure three important snowpack variables, which include snow water equivalent (SWE) observed by a pressure-sensing snow pillow, snow depth observed with a sonic distance ranger, and atmospheric air temperature observed with a thermistor. Many SNOTEL stations have SWE records beginning in the 1970s and temperature records starting in the 1980s. Sonic rangers were not installed until more recently, appearing at most sites in the early 2000s (Table 1). Our methods described below for estimating SEH values are highly dependent on calculations of bulk snow density, derived by dividing SWE by depth. Hence, our dataset is limited to roughly the last 10 yr (dependent upon individual sites), when depth data exist. We obtained temperature, SWE, and depth data at 3-h intervals directly from the NRCS (Bozeman, Montana, office). With 33 SNOTEL stations in western Montana, 313 station years of data were analyzed.
a. Cold content
The net energy change in energy of the snowpack, , is a function of multiple energy fluxes, such that
where Qsn is the net shortwave radiation flux, Qln is the net longwave radiation flux, Qh is the sensible heat flux, Qe is the latent heat flux, Qg is the ground heat flux, and Qp is heat advected from liquid water (Gray and Male 1981). Each of the heat fluxes varies over time and may dominate the heat balance under particular conditions. While the terms in Eq. (2) are well understood, confidently implementing observationally based energy balance calculations can be problematic because of numerous difficulties with accurately measuring the energy fluxes into and out of the snowpack over time (e.g., Helgason and Pomeroy 2012). In our case, the necessary meteorological data for determining historical SEH magnitudes from full energy balances do not exist for the 33 sites across western Montana. We are therefore forced to approximate the cold content generated during storm events using a scheme based purely on temperature observations.
Our approach here is to model the change in cold content of a snowpack during a sustained cold period by solving the heat equation (Paterson 1994)
where T is the snow temperature at depth z, k is the thermal diffusivity of snow, ρs is the density of snow, ρ is the density of water, F is the latent heat of fusion, ci is the specific heat capacity of ice, and S′ is the volumetric snow liquid water content at depth z. On the right-hand side of the expression, the first term accounts for heat transfer through conduction and the second term describes latent heat transfer in the snowpack. Latent energy is released when liquid water in the snowpack freezes. When there is no liquid water present, the second term on the right-hand side is zero and conduction is the primary method of heat transfer. The sum of these energy components describes the change in the temperature of the snowpack over time.
An analytic solution to Eq. (3) requires detailed knowledge of both the initial conditions of the snowpack and boundary conditions; notably, the net energy balance at the snow surface. Considering the difficulties of computing an accurate energy balance described above, we employ a numerical approximation [Eq. (4)] first adopted by Bengtsson (1982), where
The approximation yields the refreezing depth zf of the liquid content in ripe snow and added cold content as a function of snowpack conditions, time, and the overlying air temperature time series, where Ta is air temperature (°C) over the period 0 to time τ and Tm is the mean of the negative air temperature (subfreezing) over the same period. The parameter a represents the snowpack conditions:
Bengtsson (1982) showed that Eq. (4) provides reasonable results at diurnal time scales when refreezing times are relatively short (i.e., 12 h) and refreezing depths are relatively small (i.e., 200 mm). Extending the analysis, we find that this simplification is also appropriate for larger refreezing depths and longer time frames (Fig. 3). Assuming the refreezing is driven by the conduction of heat, the approximation only fails in the case of unrealistically low temperatures that are required to create a very large refreezing depth over relatively small time frames. For example, to reach a freezing depth of 34 cm in 24 h, an unrealistic mean hourly temperature of −217°C is required.
The cold content represents the negative energy over the depth of refreezing and is also approximated from Eq. (4), which assumes that average snowpack temperature in the refrozen depth is proportional to the average subfreezing temperature of the overlying air. Similarly, just the linear portion of this approximation is considered:
Our simplified approach to calculating heat exchange does not explicitly represent each term in the energy balance. Of note is that longwave radiation is not explicitly accounted for and is often the dominant mechanism for heat flux over a snowpack (e.g., Cline 1997; Marks and Dozier 1992). However, longwave radiation flux and turbulent heat flux do have a temperature dependency, which is the physical basis for the high success of degree-day snowmelt modeling schemes (Ohmura 2001). Indeed, the widely applied SNOW-17 model also employs an even simpler approximation of heat storage and exchange based entirely on temperature scaling (Anderson 1976). The impact of our approximation methods on our results may be inconsequential in some cases, but may be substantial in some circumstances. We address this in more detail later in section 4.
b. Residual pore space
The second component of the energy hurdle,
results when empty pore space is created in the upper layers of the snowpack when either liquid water refreezes or new dry pore space is added during snow accumulation events. In the case of accumulation events, the accumulation depth za is used in place of the refreezing depth zf in Eqs. (6) and (7). When liquid water refreezes in the snowpack, however, porosity is reduced, which is not explicitly accounted for in our simplified methods. The bulk snow density calculated from SWE and depth data is used to derive both the residual saturation rs (Coleou and Lesaffre 1998) and thermal diffusivity (Sturm et al. 1997). In addition, for accumulation events a new snow density ρn is calculated and used over the depth of new snow instead of the bulk snow density ρs.
c. Late season storm events
Our goal is to capture SEH events resulting from multiday storm episodes consistent with frontal passage, such as cold fronts that bring wet and/or cold conditions that delay spring snowmelt. Hence, we do not explicitly consider the diurnal events associated with the daily sun cycle. Our events are defined by new snowfall and/or cold temperatures that did not exceed 0°C for 12 or more hours. Although arbitrary, this 12-h threshold limited the investigation to punctuated cold and accumulation events taking place beyond the scale of diurnal freeze–thaw cycles (which are also important to controlling snowmelt timing, but beyond the scope of this paper). Thus, spring cold and accumulation events and their corresponding energy hurdle constituents (Qh, Qrs, and Qcc) were determined for every year of every station by considering the following criteria.
Strictly cold events are defined by having negligible snowfall (<1 cm), more than 1 day of consecutive subfreezing temperatures, and daytime temperatures that do not exceed 0°C for 12 or more hours. The termination of cold events is marked when temperatures exceed 0°C for 12 or more hours.
The beginning of a late season snowfall event is marked as the start of accumulating snowfall.
The end of a late season snowfall event is marked differently in two possible cases. In the first, if temperatures exceed 0°C for 12 or more hours, the end of the event is marked as the maximum of the rising limb of the snow depth curve. In the second, if temperatures do not exceed 0°C for 12 or more hours, cold content continues to be added to the newly accumulated snow, and the end of the event is marked when temperatures subsequently exceed 0°C for 12 or more hours.
In the case of a snowfall event where subfreezing temperatures do not occur, Qcc = 0 and Qh = Qrs.
d. Delay times
The hydrologic significance of SEH events is the extent to which they hinder the snowmelt generation and transmission process and ultimately delay runoff. Ideally, the melt generation delay resulting from SEH events could be obtained directly from changes in SWE recorded at SNOTEL sites. A time series of d(SWE)/dt does demonstrate change in water equivalent mass from a site; however, it does not necessarily reveal information about internal snowpack energy states or the generation of melt. The liquid–solid fractions in the snowpack are often reproportioned with no change in mass on the pillow, or the pillow mass often shows reductions during periods of net energy loss of the snowpack due to draining of pore water at the base. Further, meltwater can move laterally along ice layers and grain-size boundaries (Marsh and Woo 1985; Sommerfeld et al. 1994; Williams et al. 1999), leading to pillow mass changes that are not directly proportional to surface melt generation. Hence, delay in melt runoff conditions caused by SEH events cannot always be confidently back calculated from SWE records.
Our goal here is to investigate the extent to which multiday storm events alter the energy state of an initially melting snowpack and the potential of such changes to influence the timing of snowmelt runoff. The time to overcome an SEH event not only depends upon the magnitude of the energy hurdle, but also on the rate of energy input related to weather conditions following the formation of the hurdle. In other words, the time to balance the heat deficit depends not only on the size of the deficit, but also on the rate of added heat. We therefore focus our attention on normalized delays such that all SEH events are overcome using a common benchmark rate of energy input. This enables isolation and comparison of processes under circumstances where all energy balance terms have been calculated.
To compare different SEH events under a common rate of energy input with known energy balance terms, we calculate hypothetical delays using a common benchmark rate representing clear-sky and relatively warm spring melt season conditions. We generated a time series of net energy balance at the snow surface from a model constrained by field measurements collected at the Twin Lakes SNOTEL site (Fig. 4) during the 2011 melt season. We installed a meteorological station to measure wind speed, wind direction, relative humidity, air temperature, precipitation, snow depth, net radiation, and incoming and reflected shortwave radiation at 5-min intervals throughout the winter and spring of 2011. To compute the energy balance at the snow surface, we used SNOWPACK, a one-dimensional mass and energy transfer model designed for avalanche and hydrologic forecasting (Bartelt and Lehning 2002). The benchmark period was selected as a 10-day period with high atmospheric pressure and clear-sky conditions. This period had strong diurnal swings with high rates of energy input during peak daylight hours associated with warm and clear weather and little change in conditions from day to day. The peak rates of energy input were similar to those measured at a site in Colorado during similar conditions (Cline 1997). This 10-day time series served as a forcing function on SEH events to compute the hypothetical benchmark delays, comparable between all SEH events.
A total of 678 SEH events occurred across western Montana during the ~10-yr study interval. Of the 313 station years analyzed, 254 station years (81%) experienced SEH events. Accumulation events far outnumbered refreezing events, with 528 accumulation events and 150 refreezing events. The maximum accumulation event deposited 1.2 m of snow and produced a SEH of 17.5 MJ m−2. Four days of subfreezing temperatures (minimum temperature of −8°C) resulted in the largest refreezing hurdle of 3.9 MJ m−2. Accumulation events were not only more frequent but were responsible for 86% of the total energy deficit imposed on the snowpack by late season weather events (Table 2). The quantity Qrs comprised the majority of the hurdle of both refreezing and accumulation SEH events. Cold content accrued during accumulation events accounted for just 3% of the total energy deficits of accumulation events.
Precipitation is the dominant factor in determining SEH magnitude throughout the spring melt season (Fig. 5). Further, as the melt season progresses and mean atmospheric conditions become warmer, cold content becomes less important to generating SEH events. The bulk of refreezing events occur within 1 month of peak SWE, while accumulation events persist well over 2 months after peak SWE. Accumulation and refreezing events occurred as late as 23 June and 7 June (serial day 174 and 158).
The total distribution of accumulation SEH events across western Montana can differ significantly from year to year, reflecting the different character of spring melt seasons (Fig. 6). As an example, the magnitude of SEH events in western Montana is compared between 2009 and 2010 (Fig. 7). The 2009 distribution is right skewed and characterized by many small events and increasingly fewer larger events. However, the 2010 distribution demonstrates a slight shift toward a higher frequency of relatively larger events. In 2010, there are obvious large outliers at roughly 10 MJ m−2, but there are also a greater number of hurdles in the 2–5 MJ m−2 range and fewer hurdles less than 2 MJ m−2.
The magnitude and frequency of both accumulation and refreezing SEH events also demonstrates variability between sites. The Flattop Mountain SNOTEL site, for example, experienced 20 accumulation events and 1 refreezing event from 2002 to 2009 (Fig. 8). Many of these accumulation events are among the largest observed in the distribution. In contrast, the Copper Camp SNOTEL site resides at a similar elevation and observed 21 accumulation events and 9 refreezing events from 2003 to 2011. While both of these sites received accumulation events at a similar frequency, the Copper Camp SNOTEL site received far more refreezing events and some of the largest ones in the distribution. No statistically significant patterns based on site location or site elevation could be identified.
The delay between SEH formation and reestablishment of melt runoff conditions ranges from hours to at most 2 days when energy hurdles are subjected to the benchmark time series of energy input. Thus, all SEH events are overcome relatively quickly under the benchmark’s representation of warm and clear-sky spring conditions. The energy uptake under these conditions is in fact so great relative to the size of hurdles that a real-world delay period would be to a large part dictated by the point within the diurnal cycle in which the recovery began. The maximum melt delay from a refreezing event was 0.17 days, and 100% of refreezing SEH events were overcome in less than 1 day (Fig. 9). Delays from accumulation events were longer, with the maximum melt delay from an accumulation event at 2.1 days. However, 98% of accumulation SEH events were overcome in less than 1 day with the benchmark conditions.
The temperature-based approach to computing SEH magnitude used here is highly simplified and does not explicitly represent all physical processes controlling heat exchanges in the snowpack. This likely has a small effect on SEH values derived from snow accumulation events, where the temperature of the snowfall is close to the air temperature and the pore space is the primary cause of energy deficit. Ohmura (2001) shows that surface temperature does proxy outgoing longwave radiation since the snow emissivity is strongly scaled to lower-atmosphere temperature. Regardless, our methods likely underestimate SEH magnitude during periods when heat loss through longwave radiation is the dominate energy flux and the snow temperature is much colder than the near-surface air temperature. Nevertheless, if we were to double our estimates of the cold content magnitudes, accumulation events would remain the prevailing mechanism in the record for generating SEH events.
Our analysis indicates that historical SEH events in western Montana are dominated (86%) by accumulation events adding pore space to the snowpack. New snow accumulation creates larger change in the snowpack energy state simply because melting snow to fill pore spaces requires far more energy than eliminating relatively small amounts of cold content. A phase change of 1 kg of solid ice to liquid water requires nearly 164 times more energy than raising the same kilogram of ice by 1 K. Hence, the climate of western Montana is far more likely to produce bigger energy hurdles from snowfall than from extended periods with deep freezing. This is likely the case during the snowmelt season at most mountain areas of the western United States, where snowfall rates are high and temperatures are relatively mild. While spring cold events do occasionally occur, the dry/cold conditions typical of an arctic climate are likely required for the cold content component to have a greater impact on snowpack energy hurdles than snowfall.
The importance of precipitation events during spring is exacerbated by the increased albedo associated with new snowfall. While an albedo increase does not change the internal energy state of the snowpack, it does affect the ability of the snowpack to absorb solar heat (e.g., Waliser et al. 2011). Radiative heat transfer can sometimes be a particularly important mechanism for gaining heat during spring snowmelt (e.g., Marks and Dozier 1992; Cline 1997). Albedo diminishes with time, but on a time scale of days under warm spring temperatures (e.g., Robinson and Kukla 1984). Hence, this is a compounding effect to the energy hurdle created by the new snowfall. Nevertheless, if the rate of heat transfer in our benchmark example were reduced by 50%, the time to overcome the energy hurdles would typically remain on the order of hours and not days or weeks.
To explore the sensitivity of the Qrs component of accumulation SEH events to snowfall depth and density, we calculate a family of SEH events with a fixed Qcc = 0 (Fig. 10). From density we infer porosity, which is directly related to residual saturation such that more meltwater is needed to meet residual saturation requirements of lower density snow (Coleou and Lesaffre 1998). We then approximate corresponding melt generation delays using variations of the benchmark energy input rate. We do not account for densification of the snowpack over time, which would tend to reduce the delays in our lower density scenarios.
The time between a snowfall event and melt runoff is dependent upon both the magnitude of the snowfall and the rate of energy input following the event. Using our benchmark input rate, which represents warm and clear-sky conditions, a new snowfall depth of at least 0.4 m is required to delay melt beyond one day. While such snowfall events have occurred in western Montana during the spring snowmelt season, they are relatively rare at and below the elevations of SNOTEL sites. However, actual energy-input rates can be expected to be lower than our clear-sky benchmark in many cases, because mountain weather conditions are not necessarily clear sky. A delay of 24 h can result from 0.2 m of new snow depth if the energy input is reduced by 50% from our benchmark. Snow density also influences delays, with lower densities leading to longer melt delays.
The sensitivity of SEH magnitude to snowfall amount and density implies that the delays in the onset of melt following a storm episode are highly variable across the mountain landscape. Snow accumulation in mountainous terrain displays large spatial heterogeneity (e.g., Clark et al. 2011; Elder et al. 1991; Neumann et al. 2006). Further, snowfall generally increases with elevation, with substantial basin SWE accumulated at the highest elevations above instrumented sites in some cases (Gillan et al. 2010; Rice et al. 2011). Energy hurdles and the subsequent delay in establishment of runoff conditions should therefore increase with elevation from both greater snowfall and decreased energy input following snowfall.
The release of snowmelt from the base of the snowpack is commonly interrupted during the spring melt season in western Montana because of multiday episodes of cold weather and/or snow accumulation. After an event concludes, melt release is further delayed until cold content and residual saturation constraints within the snowpack are satisfied. These “energy hurdle” events occurred in over 80% of 313 SNOTEL station years analyzed. The largest accumulation event was 17.5 MJ m−2, while the largest refreezing SEH event was 3.9 MJ m−2. Accumulation SEH events not only occur more frequently (528 accumulation events versus 150 refreezing events) but constitute 86% of the collective energy deficit imposed on the snowpack by late season storm episodes. Accumulation hurdles depend on both snowfall amount and on snowfall density, as the pore space component is directly proportional to the density-derived residual saturation.
The time scale to reestablish runoff conditions after formation of the observed energy hurdles is typically hours, and at most 2 days, under relatively warm clear-sky conditions. Snowfall events are more important to delaying the establishment of runoff conditions than refreezing events because the energy deficits created in the snowpack are larger and because the increased snow albedo reduces the rate of solar heat gain. However, the heat input to the snowpack during spring season has the potential to be large relative to the size of energy hurdles, so even the largest hurdles could be overcome relatively quickly. Thus, we find that in western Montana during the time period studied, the antecedent snowpack conditions after a spring storm event can be easily overcome when subsequent weather conditions permit.
This research was funded by grants from Inland Northwest Research Alliance and the U.S. National Science Foundation (EAR-0609570). We thank three anonymous reviewers for detailed and constructive comments, one of whom was particularly generous with their time.