The sensitivity of snowpack sublimation estimates to instrument and measurement uncertainty perturbed in a Monte Carlo framework

Abstract

The bulk aerodynamic flux equation is often used to estimate snowpack sublimation since it requires meteorological measurements at only one height above the snow surface. However, to date the uncertainty of these estimates and the individual input variables and input parameters uncertainty have not been quantified. We modeled sublimation for three (average snowpack in 2005, deep snowpack in 2011, and shallow snowpack in 2012) different water years (October 1 to September 30) at West Glacier Lake watershed within the Glacier Lakes Ecosystem Experiments Site in Wyoming. We performed a Monte Carlo analysis to evaluate the sensitivity of modeled sublimation to uncertainties of the input variables and parameters from the bulk aerodynamic flux equation. Input variable time series were uniformly adjusted by a normally distributed random variable with a standard deviation given as follows: 1) the manufacturer’s stated instrument accuracy of 0.3°C for temperature (T), 0.3 m/s for wind speed (U_z), 2% for relative humidity (RH), and 1 mb for pressure (P); 2) 0.0093 m for the aerodynamic roughness length (z₀) based on z₀ profiles calculations from multiple heights; and 3) 0.08 m for measurement height (z). Often z is held constant; here we used a constant z compared to the ground surface, and subsequently altered z to account for the change in snow depth (d_s). The most important source of uncertainty was z₀, then RH. Accounting for measurement height as it changed due to snowpack accumulation/ablation was also relevant for deeper snow. Snow surface sublimation uncertainties, from this study, are in the range of 1% to 29% for individual input parameter perturbations. The mean cumulative uncertainty was 41% for the three water years with 55%, 37%, and 32% occurring for the wet, average, and low water years. The top three variables (z varying with d_s, z₀, and RH) accounted for 74% to 84% of the cumulative sublimation uncertainty.

Keywords

snow / sublimation / uncertainty / aerodynamic methods

Introduction

In mountainous regions, such as the western United States, most of the annual precipitation falls as snow and is stored in high-elevation mountain snowpacks. Snowpack sublimation is an important hydrologic process which can account for significant water losses to the atmosphere (^{Hood et al., 1999}; Liston and Elder, 2006a, b; ^{Molotch et al., 2007}; ^{Sexstone et al., 2016}). However, the amount of water that is exchanged between seasonal snowpacks and the atmosphere through sublimation is still poorly understood (^{Lang, 1981}; ^{Hood et al., 1999}; ^{Molotch et al., 2007}). The large degree of uncertainty can have significant consequences on hydrologic studies, water supply forecasting, and water supply modelling.

Sublimation losses from the snowpack have been estimated for various environments and can constitute a significant component of the water balance, with net sublimation losses estimated between 10%–35% of the seasonal snow accumulation, specifically: 12%–33% in the Canadian prairies (^{Woo et al., 2000}), 19% in the Wyoming Rocky Mountains (^{Hultstrand, 2006}), 15% in the Colorado Rocky Mountains (^{Hood et al., 1999}), and 18% in the Sierra Nevada Mountains (^{Kattelmann and Elder, 1991}). Over a 40-day period in the Colorado Rocky Mountains, total snowpack sublimation was greater than measured precipitation (^{Molotch et al., 2007}).

Sublimation occurs more readily under certain weather conditions, such as low relative humidity (RH) and increased wind speed (U_z). Vapour pressure (e) gradients (Δe) between the snowpack and the atmosphere, snow surface roughness length (z₀), wind speed, and atmospheric stability all have a significant contribution to sublimation magnitude and direction (upward as a loss or downward as vapour deposition). Methods for estimating sublimation from a snowpack include the bulk aerodynamic flux (BF) calculation, Bowen ratio-energy balance, snow evaporation pans, and aerodynamic profile (AP) methods. Newer techniques that directly measure atmospheric flux (eddy covariance, EC) have also been tested for use in snowpack sublimation monitoring in alpine and sub-alpine environments (^{Hood et al., 1999}; ^{Molotch et al., 2007}; ^{Marks et al., 2008}; ^{Reba et al., 2012}; ^{Sexstone et al., 2016}).

Sublimation losses from the snowpack are typically calculated from a mass transfer equation, as per ^{Dingman (2002)} and ^{Fassnacht (2004)}. The latent heat flux (Q_E in kJ·s^–1·m^–2) is equal to the product of the latent heat of sublimation (L_S as 2838 kJ/kg at a temperature (T) of 0°C; ^{Datt, 2011}) and the rate of latent mass transfer (E in kg·s^–1·m^–2). The most common method for estimating snowpack sublimation is measuring snowpack Q_E using the BF equation (^{Moore, 1983}). This method has the advantage of requiring meteorological measurements at only one height above the snow surface. However, a primary assumption applied to the BF method is that the snow surface temperature follows the air temperature for the estimation of saturation vapour pressure (e_sat). This assumption is often inaccurate at temperatures colder than 0°C (^{Raleigh et al., 2013}) and can lead to an over-estimation of sublimation (^{Bernier and Edwards, 1989}; ^{Marks et al., 2008}). In addition, this method also assumes the snow surface is saturated with respect to ice or water (i.e., 100% relative humidity), which may not always be the case (^{Box and Steffen, 2001}). With the surface temperature and humidity assumptions, the estimated sublimation loss (upward flux) is a function of the difference in vapour pressure between the measurement height and the surface (^{Fassnacht, 2004}) and can never be downward, such as in the form of frost deposition (^{Sexstone et al., 2016}). In reality, the snow surface temperature is limited to 0°C but the air temperature can be warmer and vapour deposition can occur depending on the relative humidity of the air, i.e., the vapour gradient.

A more accurate method for calculating the snowpack Q_E and snowpack turbulent fluxes as compared to the BF method is the AP method which requires the measurement of U_z, T, and RH at multiple heights above the snowpack (^{Cline, 1997a}; ^{Hood et al., 1999}; ^{Sexstone et al., 2016}). Both the BF and AP methods require an estimate of z₀ to define the U_z profile. Using z₀ as a variable rather than a constant parameter will alter the computed snowpack sublimation loss estimates (^{Fassnacht, 2010}); z₀ has been seen to vary by almost three orders of magnitude (2×10^–5 to 9.76×10^–3 m) (^{Brock et al., 2006}). Previous research has used z₀ as a parameter with values of 1×10^–3 m (^{Hultstrand, 2006}), 5×10^–2 m (^{Fassnacht, 2004}) and 5×10^–4 m (^{Box and Steffen, 2001}).

The EC method is considered the most accurate/direct method for calculating snowpack Q_E (^{Molotch et al., 2007}; ^{Sexstone et al., 2016}). A two-tower approach (two eddy covariance sensors) 35 miles north of Bangor, ME was used to test the EC method uncertainties in Q_E measurements, results from the study state that a Q_E measurement uncertainty of 0.005 kJ·s^–1·m^–2 was reported for an entire calendar year (^{Hollinger and Richardson, 2005}). This method is considered a direct atmospheric flux measurement technique to determine the vertical turbulent fluxes within the atmospheric boundary layer. While EC systems are fairly robust, the EC procedure requires adequate site conditions, such as long fetches, and a high frequency sonic anemometer that can be cost and energy prohibitive. Sublimation measurements that require extensive meteorological measurements and equipment for the EC and AP methods are relatively limited in complex mountain regions (^{Sexstone et al., 2016}), typically limited to research facilities, which makes the BF calculation the most common method used to estimate snowpack sublimation at local and regional scales.

The objectives of this research for one seasonally snow-covered alpine environment are (i) to quantify the sensitivity of sublimation estimates using the BF method from the uncertainty of the input measurements; (ii) to quantify snow sublimation uncertainty as a function of peak snow water equivalent and total precipitation; and (iii) to provide guidance to account for instrumentation errors and what input variables need the greatest attention while quantifying snowpack sublimation. This study provides an evaluation of the sensitivity of snow sublimation calculations and possible uncertainties which can improve our understanding of water resources, water supply forecasting, and water supply modelling.

Study site

For this study, sublimation was estimated for the West Glacier Lake watershed (WGLW) within the Snowy Range Mountains, Wyoming (41°22′30″N latitude and 106°15′30″W longitude) (Fig. 1). WGLW is part of the US Forest Service’s Glacier Lakes Ecosystem Experiments Site (GLEES) developed to conduct research on the effects of atmospheric deposition on alpine and subalpine ecosystems (^{Musselman, 1994}). Approximately 575 ha in size, GLEES consists of three small watersheds, beneath a northeast-southwest ridge, and WGLW ranges in elevation from 3277 m at the lake outlet to 3493 m at the top. Mean annual temperature is –1°C at the outlet and –2.5°C at the top of the basin (^{Korfmacher and Hultstrand, 2006}). Mean annual precipitation is 1200 mm, with approximately 75% to 85% falling as snow, which remains from late November to early June (^{Wooldridge et al., 1996}; ^{Korfmacher and Hultstrand, 2006}). This region is dominated by high westerly winds up to 26 m/s with an average of 8 m/s (^{Korfmacher and Hultstrand, 2006}). These climatic conditions combine to create an environment where snow accumulation, snow redistribution, and snowpack sublimation can have significant impacts on the watershed hydrology.

Methodology

Bulk aerodynamic flux

The meteorological variables needed for the BF equation are T (in °C), RH (in %), U_z (in m/s), and station pressure (P in mb). For the BF method, the latent mass flux is calculated as:

(1)E={\displaystyle \frac{\mathrm{0.622}\cdot {\rho }_a}{P\cdot {\varphi }_m\cdot {\varphi }_v}}\cdot U_Z\cdot {\displaystyle \frac{k^{\mathrm{2}}(q_a-q_o)}{\ln {\left({\displaystyle \frac{Z_a+Z_d}{Z_o}}\right)}^{\mathrm{2}}}},

where {\rho }_a is the density of air (kg/m³), {\varphi }_m is the stability function for momentum calculated as a function of the Richardson number (Ri, unitless), {\varphi }_v is the stability function for water vapour calculated as a function of the Ri, k is von Karman’s constant (0.4), q_o is the specific humidity (kg/kg) at the surface of the snow, q_a is the specific humidity at measurement height Z_a (in m), and Z_d is the zero-plane displacement (in m). A value of zero is used for Z_d.

The q_o at each level in the profile is determined by (^{Saucier, 1983}):

(2)q_o={\displaystyle \frac{\mathrm{0.622}e}{P-\mathrm{0.378}e}},

where e is the vapour pressure (in mb), calculated from the equation:

(3)e={\displaystyle \frac{e_{sat}\left(RH\right)}{\mathrm{100}}},

where e_sat is the saturation vapour pressure over water (in mb), estimated from the equation:

(4)e_{sat}=\mathrm{6.11}*\exp \left({\displaystyle \frac{\mathrm{17.3}*T}{T+\mathrm{237.3}}}\right),

with T in °C. Stability functions are calculated as a function of the Ri as described by ^{Ohmura (1981)}. The Richardson number is determined by:

(5)Ri={\displaystyle \frac{g}{\overline{\theta }}}\left({\displaystyle \frac{\delta \theta /\delta z}{{(\delta U/\delta z)}^{\mathrm{2}}}}\right),

where g is the acceleration due to gravity (9.81 m/s²) and \overline{\theta } is the mean potential temperature of the levels (in °C) (^{Andreas, 2002}). The stability factors ({\varphi }_m and {\varphi }_v) are estimated as a function of the calculated Ri, based on the value of Ri, as per ^{Cline (1997a)}.

Air density is calculated from the equation:

(6){\rho }_a={\displaystyle \frac{P}{R\cdot T}},

where P is air pressure (mb *100), R is the specific gas constant (287.05 J/kg/K), and T is the air temperature in degrees Kelvin.

Data

The GLEES maintains an 18-meter tower equipped with standard meteorological sensors located at 3286 m elevation between east and west Glacier Lakes (Fig. 1). T, RH, U_z, wind Direction (W_d), solar radiation (Q_h), and soil temperature (T_soil) are measured every 15 minutes. For this study, we utilized the quality controlled hourly meteorological data (T, RH, U_z, P) from the GLEES tower for the water years (October 1 to September 30) 2005, 2011, and 2012. These years were selected because they represent an average snow season (2005), a wet snow season with cooler T, higher RH and above average precipitation (2011), and a drier snow season that melted out early with warmer T, lower RH, and below average precipitation (2012) (Fig. 2 and Table 1). Snow depth (d_s) and snow water equivalent (SWE) were obtained from the Natural Resources Conservation Service Snow Telemetry (SNOTEL) site Brooklyn Lake, Wyoming located approximately 1 km to the southeast at an elevation of 3115 m (Fig. 2). Snow depth data were only available since water year 2004. The SNOTEL data were used to determine the snow-cover period for the three water years and the snowpack and average meteorological conditions were computed (Table 1).

The z₀ parameter was empirically derived from a nearby research site, Niwot Ridge Subnivean Lab Colorado that measures meteorological variables at multiple heights allowing for AP calculations. Niwot Ridge is similar in elevation and meteorology as West Glacier Lake region and is assumed that the Niwot Ridge z₀ parameter was transposable to West Glacier Lake. The z₀ parameter estimates were computed during the snowcover season (^{Cline, 1997b}) using U_z measurements at two different heights, z_a and z_b:

(7)z_o=\exp {\displaystyle \frac{(U_a*LN(z_b)-U_b*LN(z_a\left)\right)}{(z_a-z_b)}}.

An average z₀ value of 0.0043 m and a standard deviation to 0.0013 m were computed from sub-hourly wind measurements that were at 0.5 and 2.0 m above the snowpack for water years 1994 and 1995. These are similar to values (0.0022 m to 0.0050 m) presented in ^{Brock et al. (2006)}.

Monte Carlo simulations

To evaluate the sensitivity of simulated sublimation using the BF method to uncertainties of the input variables and parameters data, we performed a Monte Carlo analysis. Monte Carlo methods utilize computational algorithms to model the probability of different outcomes in a process that cannot easily be predicted due to the intervention of random variables and/or uncertainty (^{Hastings, 1970}). Only one variable or parameter was perturbed at a time to determine which one was most sensitive in the sublimation calculation; we did not examine joint-uncertainties (^{Graham et al., 2010}; ^{Sexstone et al., 2016}). Cumulative sublimation uncertainty was quantified by the addition of individual variable/parameter uncertainties (^{Bliss et al., 2011}). Seven numerical experiments were performed (Table 2) by individually changing four meteorological variables (T, RH, U_z, P) that had an hourly time step, two parameters (z₀ and z) that are usually assumed to be constant, and one parameter (z) that was used as a variable. For each hourly time step, the variable or parameter was adjusted using a random number that was selected from a normal distribution with a mean of zero and a standard deviation based on the instrumentation measurement error range which was set to the manufacturer’s stated instrument accuracy (Table 2). The perturbed RH values were constrained to a maximum of 100%. For z₀, a standard deviation was calculated from field data, as presented above (Eq. (7)). For the Monte Carlo analysis, sublimation was computed for each year 1000 times, each of the 1000 times using the randomly selected perturbation to the individual input variable or parameter to assess the sensitivity of the sources of uncertainty. The measurement height changes as snow accumulates or ablates, but z is often held constant (^{Fassnacht, 2004}). Therefore, z was used as constant of 3.0 m with a standard deviation of 0.08 m based on stated instrument accuracy (Judd snow depth sensor) (^{Ryan et al., 2008}). The value of z was also adjusted to account for the change in d_s (Fig. 2); the Monte Carlo approach was not used in this last test (Table 2).

Surface temperature estimates

We evaluated the use of dewpoint temperature to represent the snow surface temperature, as shown in ^{Raleigh et al. (2013)}. Sublimation was computed using air temperature and relative humidity to estimate the dewpoint temperature (Eq. (4)), without perturbations. These sublimation estimates were compared to the base case scenarios for the three years.

Results

For the three different snow years, the estimated sublimation for the non-perturbed base case (Table 1 and lines in Fig. 2(c)) was inversely related with the amount of snow (Figs. 2(a) and 2(b)), which is in part a function of the other meteorological conditions. The simulated sublimation was somewhat consistent over the winter of 2011, as seen by approximately constant slope of cumulative sublimation (Fig. 2(c)). In 2012, the change in cumulative sublimation rate decreased after February while in 2005 it started to increase then (Fig. 2(c)).

Using the Monte Carlo approach with input variable/parameter perturbations yielded a range of results (Table 3), which were consistent for some variables, such as RH, U_z, and z₀. The range of variability from the perturbations became larger in mid-February 2011, mid-March 2012, and early April 2011 (Fig. 2(c)). The shape of the z₀ simulation was similar for each year (Fig. 3 and the standard deviation in Table 2) and coefficient of variation was the same (Table 3). Sublimation sensitivity based on the perturbed input variables/parameters with the Monte Carlo simulations show consistent inter-annual variability in the simulations for RH, U_z, and z₀, but not for T or the two z tests (Fig. 4 and Table 3). Since the instrument accuracy of P was so high (~0.1%), the P perturbations has a negligible effect on sublimation estimates.

Calculated sublimation was most sensitive to the perturbations of z₀, RH, z (in 2011 due to deep snow; Fig. 4(b)), U_z and finally T (Table 3 and Fig. 4). Temperature perturbation yielded the most noticeable sublimation variability in 2005 (Fig. 4), about four to five as much as in the other years in the other years likely since 2005 had the most freezing days. However, this variation in temperature was small compared to the other years. When the snowpack is shallower, such as in 2005 and 2012 (Fig. 2(a)), z has less of an influence on the sensitivity of the calculated sublimation (Figs. 3(a) and 3(c)).

The range of variability in calculated sublimation using the Monte Carlo approach was most for 2011 and least for 2005 (Fig. 2(c)). However, the variation in simulated sublimation (Table 4) was much greater for the lower snow years when taken as a percentage of peak SWE (Fig. 2(b)) and annual precipitation (Table 1). The range of simulated sublimation to peak SWE and to annual total precipitation was greater by a factor of two for both 2005 and 2012 (Table 4). More importantly, sublimation was estimated to be at least half of the peak SWE in 2012, and a third of peak SWE in 2005 (Table 4). The maximum simulated sublimation was more than two-thirds for 2005 and almost 100% for 2012 (Table 4). These values can be considered scaled when subsequently compared to the annual precipitation total, since peak SWE was 67% of annual precipitation in both 2005 and 2012; it was 83% in 2011 (Table 1). For the base case, sublimation was computed to be 37%, 21%, and 52% of the annual precipitation in 2005, 2011, and 2012, respectively (Table 4).

The average daytime dewpoint temperature has a negative bias (^{Raleigh et al., 2013}). This use of dewpoint temperature (T_d) for the surface temperature (T_s) yielded a negative cumulative sublimation estimate throughout the winter for each year (Fig. 5), which is not correct. Interestingly, this assumption yielded the most deposition (negative cumulative sublimation in Fig. 5) and were somewhat of a translation of the base case (Fig. 2(c)). The computed mean hourly sublimation rate is also negative every hour and smaller in magnitude than when the surface air temperature is used, as stated above (Fig. 6). Therefore, using dewpoint temperature for surface temperature was not evaluated further.

Discussion

The three snow seasons represent a range of snow conditions within the available period of record. Given the large variability in the amount of snowfall (Fig. 2 and Table 1), the importance of sublimation to the seasonal water balance varied widely (^{Hultstrand, 2006}; Table 4). There were large variations in calculated sublimation were due to z₀ (Fig. 3) and there is much uncertainty in estimating the value of z₀ (^{Andreas, 2002}; Fassnacht et al., 2015), as it varies spatially and temporally (^{Brock et al., 2006}). Only for deep snow (2011; Fig. 2(a)) was perturbation in z relevant (Fig. 4(c)). When using the BF method, the actual measurement height should be estimated from snow depth (e.g., Fig. 2(a); ^{Fassnacht, 2010}).

There was some sensitivity to U_z, but much less to T (mostly in 2005), and essentially none to P (Fig. 4). The calculated sublimation sensitivity due to RH and U_z is based on instrumentation accuracy. More advanced methods, such as the EC method (^{Box and Steffen, 2001}; ^{Reba et al., 2012}; ^{Sexstone et al., 2016}), could improve accuracy. However, such more accurate sensors may not be as robust in the field (^{Sexstone et al., 2016}).

The mean BF cumulative sublimation uncertainty was 41% for the three water years with 55%, 37%, and 32% occurring for the wet, average, and low water years (Table 3). The mean uncertainty for each input parameter in order from largest to smallest uncertainty is z varying with d_s (15.8%), z₀ (9.6%), RH (7.1%), U_z (4.9%), z at constant height (2.0%), and T (1.9%) (Table 3). The uncertainty for top three variables account for 74% to 84% of the cumulative sublimation uncertainty (Table 3).

Latent heat flux estimates over snow covered surfaces are highly correlated for the EC and BF methods, but with some discrepancies (e.g., ^{Marks et al., 2008}), such as a positive bias in estimating sublimation using the BF method (^{Box and Steffen, 2001}; ^{Marks et al., 2008}). The BF sublimation estimates calculated here (Table 1 and Fig. 2(c)) are larger than those presented in the literature (e.g., ^{Hood et al., 1999}). The EC and AP methods were shown to compare well (^{Reba et al., 2012}), and both are considered better than BF (^{Sexstone et al., 2016}). However, the BF method provided reasonable snowpack sublimation estimates when EC instrumentation were not available (^{Sexstone et al., 2016}), which is often the case (^{Fassnacht, 2004}).

For the average and deep snowpack, there was limited inter-annual year variability whereas the drier and shallower snowpack had larger inter-annual variability of total modeled snow sublimation from the snow surface using the base case (Fig. 2(c)). This was also observed by ^{Sexstone et al. (2016)} and ^{Reba et al. (2012)} suggesting snow sublimation is largely dependent on the amount of snowfall. ^{Sexstone et al. (2016)} stated that snow sublimation rates do not scale with snowpack depth or SWE, i.e., during a low snow season a larger percentage of the snowpack is lost to sublimation and hence less total snow is available for melt. The results of this study confirm the statements made by ^{Reba et al. (2012)} and ^{Sexstone et al. (2016)}. It should be noted that blowing snow (^{Fassnacht, 2004}) and blowing snow sublimation were not considered (^{Sexstone et al., 2018}). The variability in the BF-based sublimation estimates were a function of the amount of snow (Fig. 2(c)), specifically how the snow depth is represented (Fig. 4).

The contribution from turbulent transfer to the ablation of the snow cover has been shown to vary between 5% and 60% (^{Stewart, 1982}; Kattelman and Elder, 1991; ^{Cline, 1997a}, ^b; ^{Hunsaker et al., 2012}), with sublimation losses ranging from 20% of the snowpack in the alpine of the California Sierra Nevada (^{Marks and Dozier, 1992}; ^{Marks et al.,1992}, ¹⁹⁹⁸) to 45% of the snowpack in a subalpine forest in the Rocky Mountains (^{Hood et al., 1999}; ^{Pomeroy and Essery, 1999}; ^{Marks et al., 2001}; ^{Pomeroy et al., 2006}; ^{Molotch et al., 2007}; ^{Reba et al., 2012}; ^{Sexstone et al., 2016}). This analysis shows sublimation losses generally greater than this range (21% to 52% of the snowpack for the base case in Table 4), but it was largely dependent on the amount of snow that accumulated (Fig. 2(a)). It should be noted that the reported percentages could be lower if direct measurements of SWE and d_s were made in the GLEES basin where a deeper snowpack maybe present, but this would also Change z.

The understanding of sublimation is not well established at regional macroscales, where land surface models (LSM) have been used to simulate cold season processes. LSM have a coarse resolution (5–30 km), usually larger than the scale of many cold season processes, which restricts how a LSM represents the variability of elevation, vegetation, and meteorology in complex terrain, due in part to a lack of reliable data. Model generated net sublimation has shown large variability based on the model and domain, small scale variability, and snowpack model algorithms (^{Pan et al., 2003}; ^{Sheffield et al., 2003}; ^{Bowling et al., 2004}; ^{Reba et al., 2012}; ^{Svoma, 2016}). For example, snowpack sublimation as a percent of precipitation was estimated from 0–15% using various LSMs for a grassland catchment Russia (^{Slater et al., 2001}), 10%–35% in Alaska using the Variable Infiltration Capacity (VIC) model (^{Bowling et al., 2004}), 8%–20% using various models over the Pacific Northwest of the United States (^{Sheffield et al., 2003}), and from 0–4% in low valleys to 20%–30% in the high mountains of the Upper Colorado River Basin (UCRB) using SnowModel (Liston and Elder, 2006b). Sublimation in isolated areas of the UCRB has been modeled to exceed 30% of the annual precipitation (^{Phillips, 2013}). This study focused on site specific calculations representing a small localized region, and showed some large interannual variability (Tables 3 and 4, Figs. 1, 2, and 3). Thus, without accurate measurements of sublimation at the localized scale, it is more difficult to evaluate large scale estimates of sublimation (^{Svoma, 2016}). Advancing the understanding of how localized variability affects large scale sublimation can still be achieved through careful model experiments, especially if the sensitivity of results is examined through variable/parameter adjustment and model simulations (^{Strasser et al., 2011}; ^{Phillips, 2013}).

One of the simplistic assumptions of applying the BF method is that the snow surface temperature tracks the air temperature measurement, and is thus likely a major uncertainty in the experimental setup and can lead to substantial overestimations of sublimation (^{Bernier and Edwards, 1989}; ^{Marks et al., 2008}). Recent experimental studies have used outgoing longwave radiation measurements to measure snow surface temperature and highlight that the snow surface tends to be consistently colder than the air temperature, especially at night (^{Reba et al., 2012}; ^{Sexstone et al., 2016}). Most models that utilize the BF equation do not assume that snow surface temperature tracks air temperature but rather solve the snow energy balance equation for snow surface temperature (e.g., ^{Liston and Elder, 2006a,b}). We tested using dewpoint temperature to estimate the snow surface temperature and this yielded a negative (downward) sublimation (Figs. 5 and 6). Downward sublimation or deposition as frost does occur some nights and into the early morning, so an improvement is necessary to the assumption for the snow surface temperature.

Conclusions

Surface sublimation sensitivity was evaluated and quantified in an alpine environment based on the BF method over three water years (average snow season, a wet snow season, and a low snow season). The magnitude and range of perturbed snow sublimation estimates show considerable uncertainty with perturbed input variables. Of the factors affecting the calculated snow sublimation estimates, z₀ and RH are the most significant. For deep snow conditions, where the distance between the instrument and the snow surface can be small, z is important. Wind speed uncertainty caused sensitivity in the sublimation estimates. Temperature perturbation only yielded noticeable sublimation variability in 2005.

Sublimation calculations are derived from a number of variables and parameters, many of which have rather high degrees of uncertainty. As a result, snowpack sublimation is, often reported as a single number, but can be better characterized as a range of values. Snow surface sublimation uncertainties, from this study, are in the range of 1% to 29% for individual input parameter perturbations, with the top three variables (z₀, RH, and z) accounting for 74% to 84% of the cumulative sublimation uncertainty. Surface sublimation uncertainties from this study provide a means to properly account for instrumentation errors and what variables need the greatest attention while performing snowpack sublimation estimates in high elevation regions.

References

[1] Andreas E L (2002). Parameterizing scalar transfer over snow and ice: a review. J Hydrometeorol, 3(4): 417–432

[2] Bernier P, Edwards G C (1989). Differences between air and snow surface temperatures during snow evaporation. In: Proceedings of Eastern Snow Conference, 46: 117–120

[3] Bliss A, Cuffey K, Kavanaugh J (2011). Sublimation and surface energy budget of Taylor Glacier, Antarctica. J Glaciol, 57(204): 684–696

[4] Bowling L C, Pomeroy J W, Lettenmaier D P (2004). Parameterization of blowing-snow sublimation in a macroscale hydrology model. J Hydrometeorol, 5(5): 745–762

[5] Box J E, Steffen K (2001). Sublimation on the Greenland ice sheet from automated weather station observations. J Geophys Res, 106(D24): 33965–33981

[6] Brock B W, Willis I C, Sharp M J (2006). Measurement and parameterization of aerodynamic roughness length variations at Haut Glacier d’Arolla, Switzerland. J Glaciol, 52(177): 281–297

[7] Cline D (1997a). Snow surface energy exchanges and snowmelt at a continental, midlatitude alpine site. Water Resour Res, 33(4): 689–701

[8] Cline D (1997b). Effect of seasonality of snow accumulation and melt on snow surface energy exchanges at a continental Alpine site. J Appl Meteorol, 36(1): 32–51

[9] Datt P (2011). Latent Heat of Sublimation. In: Singh V P, Singh P, Haritashya U K, eds. Encyclopedia of Snow, Ice and Glaciers, Encyclopedia of Earth Sciences Series. Dordrecht: Springer

[10] Dingman S L (2002). Physical Hydrology. 2nd ed, Prentice Hall, Upper Saddle River, New Jersey, 646 pp

[11] Fassnacht S R (2004). Estimating alter-shielded gauge snowfall undercatch, snowpack sublimation, and blowing snow transport at six sites in the coterminous USA. Hydrol Processes, 18(18): 3481–3492

[12] Fassnacht S R (2010). Temporal changes in small scale snowpack surface roughness length for sublimation estimates in hydrological modeling. Cuadernos de Investigación Geográfica, 36(1): 43–57

[13] Fassnacht S R, Oprea I, Shipman P, Kirkpatrick J, Borleske J, Motta F, Kamin D J (2015). Geometric Methods in the Study of the Snow Surface Roughness. Proceedings of the 35th Annual American Geophysical Union Hydrology Days Conference. Fort Collins, Colorado, 41–50

[14] Graham C B, van Verseveld W, Barnard H R, McDonnell J J (2010). Estimating the deep seepage component of the hillslope and catchment water balance within a measurement uncertainty framework. Hydrol Processes, 24(25): 3631–3647

[15] Hastings W K (1970). Monte Carlo sampling methods using Markov chains and their applications. Biometrika, 57(1): 97–109

[16] Hollinger D Y, Richardson A D (2005). Uncertainty in eddy covariance measurements and its application to physiological models. Tree Physiol, 25(7): 873–885

[17] Hood E, Williams M, Cline D (1999). Sublimation from a seasonal snowpack at a continental mid-latitude alpine site. Hydrol Processes, 13(12–13): 1781–1797

[18] Hultstrand D M (2006). Geostatistical Methods for Estimating Snowmelt Contribution to the Seasonal Water Balance in an Alpine Watershed. Dissertation for Master degree. Watershed Science Program, Colorado State University, 130 pp

[19] Hunsaker C T, Whitaker T W, Bales R C (2012). Snowmelt runoff and water yield along elevation and temperature gradients in California’s Southern Sierra Nevada. J Am Water Resour Assoc, 48(4): 667–678

[20] Kattelmann R, Elder K (1991). Hydrologic characteristics and water balance of an alpine basin in the Sierra Nevada. Water Resour Res, 27(7): 1553–1562

[21] Korfmacher J L, Hultstrand D M (2006). Glacier Lakes Ecosystem Experiments Site (GLEES) hourly meteorological tower data: 1989–2005. U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station, Fort Collins, CO

[22] Lang H (1981). Is evaporation an important component in high alpine hydrology? Nord Hydrol, 12(4–5): 217–224

[23] Liston G E, Elder K (2006a). A meteorological distribution system for high-resolution terrestrial modeling (MicroMet). J Hydrometeorol, 7(2): 217–234

[24] Liston G E, Elder K (2006b). A distributed snow-evolution modeling system (snowmodel). J Hydrometeorol, 7(6): 1259–1276

[25] Marks D, Cooley K R, Robertson D C, Winstral A (2001). Long-term snow database, Reynolds Creek Experimental Watershed, Idaho, United States. Water Resour Res, 37(11): 2835–2838

[26] Marks D, Dozier J (1992). Climate and energy exchange at the snow surface in the alpine region of the Sierra Nevada. Water Resour Res, 28(11): 3043–3054

[27] Marks D, Dozier J, Davis R E (1992). Climate and energy exchange at the snow surface in the alpine region of the Sierra Nevada: meteorological measurements and monitoring. Water Resour Res, 28(11): 3029–3042

[28] Marks D, Kimball J, Tingey D, Link T (1998). The sensitivity of snowmelt processes to climate conditions and forest cover during rain-on-snow: a study of the 1996 Pacific Northwest flood. Hydrol Processes, 12(10–11): 1569–1587

[29] Marks D, Winstral A, Flerchinger G, Reba M, Pomeroy J, Link T, Elder K (2008). Comparing simulated and measured sensible and latent heat fluxes over snow under a pine canopy. J Hydrometeorol, 9(6): 1506–1522

[30] Molotch N P, Blanken P D, Williams M W, Turnipseed A A, Monson R K, Margulis S A (2007). Estimating sublimation of intercepted and sub-canopy snow using eddy covariance systems. Hydrol Processes, 21(12): 1567–1575

[31] Moore R D (1983). On the use of the bulk aerodynamic formulae over melting snow. Nord Hydrol, 14(4): 193–206

[32] Musselman R C (1994). The Glacier Lakes Ecosystem Experiments Site. Rocky Mountain Forest and Range Experiment Station, General Technical Report RM-249, 94 pp

[33] Ohmura A (1981). Climate and energy balance on arctic tundra, Axel Heiberg Island, Canadian Arctic Archipelago, Spring and Summer 1969, 1970 and 1972. Zürcher Geographische Schriften, 3, Dept. of Geography, Swiss Federal Institute of Technology (ETH), Zürich, 447 pp

[34] Pan M, Sheffield J, Wood E F, Mitchell K E, Houser P R, Schaake J C, Robock A, Lohmann D, Cosgrove B, Duan Q Y, Luo L F, Wayne Higgins R, Pinker R T, Dan Tarpley J (2003). Snow process modeling in the North American Land Data Assimilation System (NLDAS): 2. Evaluation of model-simulated snow water equivalent. J Geophys Res, 108(D22): 8850

[35] Phillips M (2013). Estimates of sublimation in the Upper Colorado River basin. Unpublished Master of Science Thesis, Department of Atmospheric Science, Colorado State University, 57pp

[36] Pomeroy J W, Bewley D, Essery R H, Hedstrom N, Link T, Granger R, Sicart J, Ellis C, Janowicz J R (2006). Shrub tundra snowmelt. Hydrol Processes, 20(4): 923–941

[37] Pomeroy J W, Essery R L H (1999). Turbulent fluxes during blowing snow: field tests of model sublimation predictions. Hydrol Processes, 13(18): 2963–2975

[38] Raleigh M S, Landry C C, Hayashi M, Quinton W L, Lundquist J D (2013). Approximating snow surface temperature from standard temperature and humidity data: new possibilities for snow model and remote sensing evaluation. Water Resour Res, 49(12): 8053–8069

[39] Reba M L, Pomeroy J, Marks D, Link T E (2012). Estimating surface sublimation losses from snowpacks in a mountain catchment using eddy covariance and turbulent transfer calculations. Hydrol Processes, 26(24): 3699–3711

[40] Ryan W A, Doesken N J, Fassnacht S R (2008). Evaluation of ultrasonic snow depth sensors for U.S. snow measurements. J Atmos Ocean Technol, 25(5): 667–684

[41] Saucier W J (1983). Principles of Meteorological Analysis. Mineola: Dover Earth Science, 464pp

[42] Sexstone G A, Clow D W, Fassnacht S R, Liston G E, Hiemstra C A, Knowles J F, Penn C A (2018). Snow sublimation in mountain environments and its sensitivity to forest disturbance and climate warming. Water Resour Res, 54(2): 1191–1211

[43] Sexstone G A, Clow D W, Stannard D I, Fassnacht S R (2016). Comparison of methods for quantifying surface sublimation over seasonally snow-covered terrain. Hydrol Processes, 30(19): 3373–3389

[44] Sheffield J, Pan M, Wood E F, Mitchell K E, Houser P R, Schaake J C, Robock A, Lohmann D, Cosgrove B, Duan Q Y, Luo L F, Wayne Higgins, Pinker R T, Dan Tarpley J, Ramsay B H (2003). Snow process modeling in the North American Land Data Assimilation System (NLDAS): 1. Evaluation of model-simulated snow cover extent. J Geophys Res, 108(D22): 8849

[45] Slater A G, Schlosser C A, Desborough C E, Pitman A J, Henderson-Sellers A, Robock A, Vinnikov K Y, Mitchell K, Boone A, Braden H, Chen F, Cox P M, de Rosnay P, Dickinson R E, Dickinson R E, Dai Y J, Duan Q, Entin J, Etchevers P, Gedney N, Gusev Y M, Habets F, Kim J, Koren V, Kowalczyk E A, Nasonova O N, Noilhan J, Schaake S, Shmakin A B, Smirnova T G, Verseghy D, Wetzel P, Xue Y, Yang Z L, Zeng Q (2001). The representation of snow in land surface schemes: results from PILPS 2(d). J Hydrometeorol, 2(1): 7–25

[46] Stewart B J (1982). Sensitivity and Significance of Turbulent Energy Exchange over an Alpine Snow Surface. Dissertation for Master degree. University of California, Santa Barbara, CA

[47] Strasser U, Warscher M, Liston G E (2011). Modeling snow canopy processes on an idealized mountain. J Hydrometeorol, 12(4): 663–677

[48] Svoma B M (2016). Difficulties in determining snowpack sublimation in complex terrain at the macroscale. Adv Meteorol, 2016: 9695757

[49] Woo M, Marsh P, Pomeroy J W (2000). Snow, frozen soils and permafrost hydrology in Canada, 1995–1998. Hydrol Processes, 14(9): 1591–1611

[50] Wooldridge G L, Musselman R C, Sommerfeld R A, Fox D G, Connell B H (1996). Mean wind patterns and snow depth in an alpine-subalpine ecosystem as measured by damage to coniferous trees. J Appl Ecol, 33(1): 100–108

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature