Estimating degree-day factors of snow based on energy flux components

Meltwater from mountainous catchments dominated by snow and ice is a valuable source of fresh water in many regions. At mid-latitudes, seasonal snow cover and glaciers act like a natural reservoir by storing precipitation during winter and releasing it in spring and summer. Snowmelt is usually modelled either by energy balance or by temperature-index approaches. The energy balance approach is process-based and more sophisticated but requires extensive input data, while the temperature-index approach uses the degree-day factor (DDF) as a key parameter to estimate melt of snow and ice merely from air temperature. Despite its simplicity, the temperature-index approach has proved to be a powerful tool for simulating the melt process especially in large and data-scarce catchments.

The present study attempts to quantify the effects of spatial, temporal, and climatic conditions on the DDF of snow in order to gain a better understanding of which influencing factors are decisive under which conditions. The analysis is based on the individual energy flux components; however, formulas for estimating the DDF are presented to account for situations where observed data are limited. A detailed comparison between field-derived and estimated DDF values yields a fair agreement with bias = 0.14 mm ^∘C⁻¹ d⁻¹ and root mean square error (RMSE) = 1.12 mm ^∘C⁻¹ d⁻¹.

The analysis of the energy balance processes controlling snowmelt indicates that cloud cover and snow albedo under clear sky are the most decisive factors for estimating the DDF of snow. The results of this study further underline that the DDF changes as the melt season progresses and thus also with altitude, since melting conditions arrive later at higher elevations. A brief analysis of the DDF under the influence of climate change shows that the DDFs are expected to decrease when comparing periods of similar degree days, as melt will occur earlier in the year when solar radiation is lower, and albedo is then likely to be higher. Therefore, the DDF cannot be treated as a constant parameter especially when using temperature-index models for forecasting present or predicting future water availability.

1 Introduction

Meltwater from snow- and ice-dominated mountainous basins is a main source of fresh water in many regions (Immerzeel et al., 2020). Seasonal snow cover and glaciers act as natural reservoirs which significantly affect catchment hydrology by temporarily storing and releasing water on various timescales (Jansson et al., 2003). In such river basins, snow and glacier melt runoff modelling is a valuable tool when predicting downstream river flow regimes, as well as when assessing the changes in the cryosphere associated with climate change (Hock, 2003; Huss and Hock, 2018). Therefore, a more accurate quantification of the snowmelt and ice melt processes and related parameters is the key to a successful runoff modelling of present and future water availability.

Two different approaches are common in snowmelt modelling. The energy balance approach is process-based but data-intensive, since melt is deduced from the balance of incoming and outgoing energy components (Braithwaite, 1995; Arendt and Sharp, 1999). On the contrary, temperature-index models, also called degree-day models, merely use the air temperature as an index to assess melt rates (Martinec, 1975; Bergström, 1976; Quick and Pipes, 1977; DeWalle and Rango, 2008). The degree-day approach is very common and popular since air temperature is an excellent surrogate variable for the energy available in near-surface atmosphere that governs the snowmelt process (Lang and Braun, 1990). The relationship between temperature and melt is defined by the degree-day factor (DDF) (Zingg, 1951; Braithwaite, 2008), which is the amount of melt that occurs per unit positive degree day (Braithwaite, 1995; Kayastha et al., 2003; Martinec et al., 2008). There are different methods by which the DDF can be determined, e.g. by measurements using ablation stakes (Zhang et al., 2006; Muhammad et al., 2020), by using snow lysimetric outflows (Kustas et al., 1994), by estimating daily changes in the snow water equivalent (Martinec, 1960; Rango and Martinec, 1979, 1995; Kane et al., 1997), or by using satellite-based snow cover data (Asaoka and Kominami, 2013; He et al., 2014).

The DDF is usually treated as a decisive parameter subject to model calibrations because sufficient direct observations are typically lacking in large catchments. Most commonly, for calibrating the DDF, runoff is used (Hinzman and Kane, 1991; Klok et al., 2001; Luo et al., 2013; Bogacki and Ismail, 2016). However, it is also important to note that the calibration of the DDF using runoff can be significantly affected by other model parameters due to their interdependency (Gafurov, 2010; He et al., 2014). Researchers also select DDF directly from other studies; hence, the spatial transferability is not always good (e.g. Carenzo et al., 2009; Wheler, 2009). Despite its simplicity, this approach has proved to be a powerful tool for simulating the complex melt processes especially in large and data-scarce catchments (Zhang et al., 2006; Immerzeel et al., 2009; Tahir et al., 2011; Lutz et al., 2016).

Figure 1 Location of the Brunnenkopfhütte automatic snow and weather station in the Dreisäulerbach catchment – German Alps.

Extensive research has been devoted to the enhancement of the original degree-day approach. Braun (1984) introduced the temperature- and wind-index method by the inclusion of a wind-dependent scaling factor. A hybrid approach, which combines both temperature-index and energy balance methods, was introduced by Anderson (1973). Oerlemans (2001) introduced a simplified energy balance approach. Hock (1999) attempted to improve the simple temperature-index model by adding a term to consider potential incoming direct solar radiation for clear-sky conditions. The potential clear-sky solar radiation is calculated as a function of the position of the Sun, geographic location, and a constant atmospheric transmissivity (Hock and Noetzli, 1997; Hock, 1999). This model is comparable with the data requirements of a simple degree-day model. Pellicciotti et al. (2005) considered the net shortwave radiation instead of just incoming shortwave radiation by including snow albedo in their proposed degree-day model. Although all these enhancements focus on adding more physical foundation to the original degree-day method, the classical approach is still more popular because of its simplicity and merely dependence on air temperatures.

A weakness of the degree-day approach is the fact that it works well over longer time periods (e.g. 10 d, monthly, seasonal) but with increasing temporal resolution, in particular for sub-daily time steps, the accuracy decreases (Lang, 1986; Hock, 1999). In addition, the spatial variability of melt rates is not modelled accurately as the DDFs are usually considered invariant in space. However, melt rates can be subject to substantial small-scale variations, particularly in high-mountain regions due to topography (Hock, 1999). For example, topographic features (e.g. topographic shading, aspect, and slope angles) including altitude of a basin can influence the spatial energy conditions for snowmelt and lead to significant variations of the DDF (Hock, 2003; Marsh et al., 2012; Bormann et al., 2014). Under otherwise similar conditions, DDFs are expected to increase with (i) increasing elevation, (ii) increasing direct solar radiation input, and (iii) decreasing albedo (Hock, 2003).

Obviously, the DDF cannot be treated as a constant parameter as it varies due to the changes in the physical properties of the snowpack over the snowmelt season (Rango and Martinec, 1995; Prasad and Roy, 2005; Shea et al., 2009; Martinec et al., 2008; Ismail et al., 2015; Kayastha and Kayastha, 2020). The spatio-temporal variation in the DDF (Zhang et al., 2006; Asaoka and Kominami, 2013) not only affects the accuracy of snowmelt and ice melt modelling (Quick and Pipes, 1977; Braun et al., 1993; Schreider et al., 1997) but is also key to estimating heterogeneity of the snowmelt regime (Hock, 1999, 2003; DeWalle and Rango, 2008; Braithwaite, 2008; Schmid et al., 2012). Since melt depends on energy balance processes and topographic settings, changes in DDFs are a result of energy components that vary with different climatic conditions (Ambach, 1985; Braithwaite, 1995). Another topic that needs attention is the stationarity of the DDF under climate change (Matthews and Hodgkins, 2016). Moreover, it is important to consider the over-sensitivity of temperature-index models, often used by large-scale studies, to future warming (Bolibar et al., 2022; Vincent and Thibert, 2022). Future water availability under climate change scenarios is typically modelled with DDFs calibrated for the present climate, which increases the parametric uncertainty introduced by the hydrological models (Lutz et al., 2016; Ismail and Bogacki, 2018; Hasson et al., 2019; Ismail et al., 2020).

In order to allow for a more-process-based estimate of the DDF, the present study attempts to quantify the contribution of each energy balance component to snowmelt and subsequently to the overall DDF. Considering that degree-day models are typically utilised in large catchments with data-scarce conditions, we estimate energy balance components by formulas with minimum data requirement following the approach by Walter et al. (2005). Based on these formulas, the DDF contribution corresponding to the respective energy components is quantified in tables and graphs for common snowmelt conditions, which can be used for a rapid appraisal. The presented approach is open in the sense that if for any of the energy balance components observed data are available, or more sophisticated models are desired, these can easily replace each of the presented approximations.

The objective of this study is not to incorporate an energy-balance-based DDF approach into temperature-index models but rather to gain a quantitative insight on how different factors affect the DDF. Through this approach, we aim to obtain a good estimate and realistic limits for calibration of the DDF parameter to predict changes during the melt season for forecasting purposes or to study the effects of climate change.

Figure 2 Automatic snow and weather station at Brunnenkopfhütte, 1602 m a.s.l. (image credit: Wolfgang Bogacki).

2 Test site and datasets

2.1 Test site

The test site is located in the Dreisäulerbach catchment, which is a part of the Isar river system and lies in the sub-alpine region of Bavaria in the Ammergau Alps (Ammergauer Alpen), Germany. The Dreisäulerbach catchment approximately lies between latitudes 47^∘34^′55^′′–47^∘35^′05^′′ N and longitudes 10^∘56^′40^′′–10^∘57^′07^′′ E. It covers an area of about 2.3 km² and has a mean hypsometric elevation of just over 1200 m a.s.l. (Fig. 1). The elevation ranges from about 950 m a.s.l. at the Linderhof gauging station up to 1768 m a.s.l. at the Hennenkopf station.

The area is mostly characterised by south-facing slopes but also contains northern slopes in southern parts of the catchment. The catchment is densely forested, which during the winter season is fully snow covered. The mean annual temperature in the observation period (i.e. November 2016–May 2021) is about 5.8 ^∘C, and the long-term (i.e. 1961–2018) mean annual precipitation at the Ettal-Linderhof station of the Water Science Service Bavaria is reported to be 1676 mm (Kopp et al., 2019).

In order to observe the seasonal snow dynamics, snow measurement instruments in addition to a standard meteorological station have been installed at the Brunnenkopfhütte test site at an elevation of 1602 m a.s.l. (see Fig. 2). The installed station has various sensors that measure temperature, pressure, wind, solar radiation (incoming, outgoing), snow depth, snow scale, a snowpack analyser, and a pluviometer. Table 1 summarises the observed monthly meteorological data at Brunnenkopfhütte station. Figure 3 presents the observed snow water equivalent (SWE) at the test site.

Table 1 Observed monthly average meteorological data at Brunnenkopfhütte (November 2016–May 2021).

T_a: air temperature; P: precipitation; u: wind speed; RH: relative humidity; A: albedo (only considered when ground is snow covered); K_T: clearness index; SR_in: incoming shortwave radiation.

Figure 3 Observed SWE (mm) at the Brunnenkopfhütte snow station (period: winter 2016/2017–2020/2021).

2.2 Datasets

The present study utilises three different datasets. Data sources and the aim of using these datasets are described below.

1. We use observed hydro-meteorological datasets from a test site (i.e. Brunnenkopfhütte) with the aim to show how the DDF of snow can be estimated for a specific site under naturally varying hydro-meteorological conditions.

2. In order to demonstrate the variation of the DDF over time, location, and altitude, as well as its significance for temperature-index modelling, we use elevation-zone temperature data of the upper Jhelum basin from a previous study (Bogacki and Ismail, 2016).

3. In the discussion section (Sect. 5), we perform a brief analysis in order to show the influence of climate change on the DDF in poorly monitored regions. In this specific analysis, projected changes in temperature are based on a previous study (Ismail et al., 2020). These projected changes in temperature are the median of four global climate models (GFDL-ESM2M, HadGEM2-ES, IPSL-CM5A-LR, and MIROC5) that are driven by two representative concentration pathways (RCP2.6 and RCP8.5). These data are provided by the Inter-Sectoral Impact Model Intercomparison Project (ISIMIP) (Hempel et al., 2013; Frieler et al., 2017).

3 Materials and methods

The primary objective of this paper is to analyse the contribution of individual energy balance components to snowmelt, in order to better understand and predict how the lumped degree-day factor will vary with the season, latitude, altitude, and the actual meteorological conditions. In addition, we want to demonstrate following the approach of Walter et al. (2005) how these energy balance components can be estimated with minimal data requirements, as limited data availability is the major reason to apply temperature-index models.

3.1 Degree-day factor

The basic formulation of the degree-day method calculates daily snowmelt depth M (mm) by multiplying the number of degree days T_DD (^∘C d) with the degree-day factor DDF (mm ^∘C⁻¹ d⁻¹) (Zingg, 1951; Braithwaite, 1995; Rango and Martinec, 1995).

$$\begin{array}{}\text{(1)}& M=\mathrm{DDF}\times T_{\mathrm{DD}}\end{array}$$

Degree days T_DD are only defined if a characteristic air temperature lies above a reference temperature T₀; otherwise, T_DD is set to 0 ^∘C d. Typically, the freezing point T₀ = 0 ^∘C is chosen as reference temperature (DeWalle and Rango, 2008). Depending on the availability of temperature data, the characteristic air temperature is usually calculated as the mean of maximum and minimum daily air temperatures (Braithwaite, 1995) or the mean of hourly observations (Rango and Martinec, 1995; DeWalle and Rango, 2008). Other approaches like daily maximum temperature (Bagchi, 1983), integrating the positive part of a diurnal cycle (Ismail et al., 2015), or averaging the positive degree-day sum of m daily observations (Braithwaite and Hughes, 2022) are also common.

By a simple re-arrangement of Eq. (1) to

$$\begin{array}{}\text{(2)}& \mathrm{DDF}={\displaystyle \frac{M}{T_{\mathrm{DD}}}},\end{array}$$

the DDF can be calculated for given degree days T_DD, if the daily melt depth M is known either by observation or by calculation. Likewise, the portion of the degree-day factor DDF_i associated with the melt depth M_i, which relates to any of the individual energy balance components (see Eq. 4), can be determined.

The energy needed to melt ice at 0 ^∘C into liquid water at 0 ^∘C is defined by the latent heat of fusion of ice (333.55 kJ kg⁻¹). Thus the melt depth M_i caused by an energy flux Q_i (W m⁻²) over a certain time period Δt (s) can be calculated from the relation (USACE, 1998; Hock, 2005)

$$\begin{array}{}\text{(3)}& M_i={\displaystyle \frac{Q_i}{\mathit{\lambda }{\mathit{\rho }}_{\mathrm{w}}}}\mathrm{\Delta }t\cong \mathrm{3.00}\times {\mathrm{10}}^{-\mathrm{6}}{Q}_i\mathrm{\Delta }t,\end{array}$$

where ρ_w is the density of water at 0 ^∘C (999.84 kg m⁻³). In the context of degree-day factor models, the time period Δt is usually taken as 1 d = 86 400 s, though some authors (Hock, 1999; McGinn, 2012) have calculated degree-day factors also for sub-daily, e.g. hourly, periods. According to the relation given in Eq. (3), an energy flux of 1 W m⁻² for a single day will result in a melt depth of 0.26 mm.

3.2 Energy balance

The energy flux available for snowmelt Q_M can be calculated from the balance of energy fluxes entering or leaving the snowpack and the change in the internal energy stored in the snowpack ΔQ (e.g. USACE, 1998):

$$\begin{array}{}\text{(4)}& Q_{\mathrm{M}}=Q_{\mathrm{S}}+Q_{\mathrm{L}}+Q_{\mathrm{H}}+Q_{\mathrm{E}}+Q_{\mathrm{G}}+Q_{\mathrm{P}}-\mathrm{\Delta }Q,\end{array}$$

where Q_S and Q_L are the net shortwave and longwave radiation, Q_H is the sensible heat, Q_E is the latent energy of condensation or vaporisation, Q_G is the heat conduction from the ground, and Q_P is the energy contained in precipitation (all terms in W m⁻²).

In the following sections, the individual components of the energy balance are discussed in more detail. As heat conduction from the ground into the snowpack is small and can be in general neglected except when first snow falls on warm ground (Anderson, 2006), the component Q_G is not considered in the further analysis.

3.2.1 Shortwave radiation

Shortwave radiation emitted from the Sun is usually the largest source of energy input to the snowpack. The net energy flux Q_S (W m⁻²) entering the snowpack by absorption of shortwave radiation is

$$\begin{array}{}\text{(5)}& Q_{\mathrm{S}}=\left(\mathrm{1}-A\right)S_i,\end{array}$$

where A is the snow albedo (–) and S_i the incident solar radiation (W m⁻²) on the snow surface. A widely used approach (Masters, 2004; Yang and Koike, 2005; Badescu, 2008) to determine the incident solar radiation on Earth's surface is the introduction of a clearness index K_T (–)

$$\begin{array}{}\text{(6)}& S_i=K_T{S}_{\mathrm{0}},\end{array}$$

where S₀ is the mean daily potential extraterrestrial solar radiation (W m⁻²) that would insolate a horizontal surface on the Earth's ground if no atmosphere would be present.

Potential insolation at the top of atmosphere

The potential insolation, which is only dependent on the changing position of the Sun during the year in relation to the geographic location of the incident point on the Earth's surface, can be calculated from the equation (Masters, 2004)

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{S}_{\mathrm{0}}& =G_{\mathrm{S}}{\displaystyle \frac{\mathrm{1}}{d_{\mathrm{r}}^{\mathrm{2}}}}{\displaystyle \frac{\mathrm{1}}{\mathit{\pi }}}\left(\cos \left(\mathrm{\varnothing }\right)\cos \left(\mathit{\delta }\right)\cos \left({\mathit{\omega }}_{\mathrm{s}}\right)\right.\\ & \left.+{\mathit{\omega }}_{\mathrm{s}}\sin \left(\mathrm{\varnothing }\right)\sin \left(\mathit{\delta }\right)\right),\end{array}\end{array}$$

where G_S is the solar constant (W m⁻²), d_r the relative distance of the Earth to the Sun (–), ϕ the geographic latitude (rad) of the incident point, δ the solar declination (rad), and ω_s the sunrise hour angle (rad). The solar constant G_S is slightly varying with the occurrence of so-called sunspots. Measurements by Kopp and Lean (2011) indicate a present value of about 1361 W m⁻².

Both Sun position variables, the relative distance of the Earth to the Sun and the solar declination, can be calculated quite exactly by rigorous astronomical algorithms (Meeus, 1991; Reda and Andreas, 2004), but for non-astronomical purposes, more simple formulas are sufficiently accurate. The relative distance of the Earth to the Sun, which varies throughout the year due to the elliptical orbit of the Earth, can be approximated by (Masters, 2004)

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{1}}{d_{\mathrm{r}}^{\mathrm{2}}}}\approx \mathrm{1}+\mathrm{0.034}\cos \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }J}{\mathrm{365}}}\right),\end{array}$$

where J is the day number, with J = 1 on 1 January. The solar declination can be obtained from the sinusoidal relationship

$$\begin{array}{}\text{(9)}& \mathit{\delta }\approx \mathrm{0.409}\sin \left({\displaystyle \frac{\mathrm{2}\mathit{\pi }}{\mathrm{365}}}\left(J-\mathrm{81}\right)\right),\end{array}$$

which puts the spring equinox on day J = 81. Knowing the solar declination δ, the sunrise hour angle ω_s can be calculated from

$$\begin{array}{}\text{(10)}& \cos {\mathit{\omega }}_{\mathrm{s}}=-\tan \left(\mathrm{\varnothing }\right)\tan \left(\mathit{\delta }\right).\end{array}$$

On the Northern Hemisphere, the maximum extraterrestrial radiation occurs at the summer solstice with a fairly identical mean daily energy flux of about 480 W m⁻² over latitudes 30–60^∘ N, as the Sun's lower-altitude angle at higher latitudes is compensated for by longer daylight hours. On the contrary, minimum extraterrestrial radiation at the winter solstice varies strongly with latitude, e.g. 227 W m⁻² at 30^∘ and only 24 W m⁻² at 60^∘ N.

Clearness index

When the solar radiation passes through the atmosphere, it is partly scattered and absorbed. While even on a clear day only about 75 % of the incoming radiation reaches the ground, by far the largest reflection is caused by clouds. A vast number of solar radiation models exist that parameterise this effect, which is denoted as clearness index K_T or atmospheric transmissivity τ, as a function of meteorological variables. For a review see e.g. Evrendilek and Ertekin (2008), Ahmad and Tiwari (2011), or Ekici (2019).

A fundamental and widely used solar radiation model that is proposed in the context of evapotranspiration calculations (Allen et al., 1998) is the Ångström–Prescott model, which relates the clearness index to the relative sunshine duration

$$\begin{array}{}\text{(11)}& K_T={\displaystyle \frac{S_i}{S_{\mathrm{0}}}}=a+b{\displaystyle \frac{n}{N}},\end{array}$$

where n is the actual and N the maximal possible duration of sunshine (h), where the latter can be calculated from the sunrise hour angle ω_s,

$$\begin{array}{}\text{(12)}& N={\displaystyle \frac{\mathrm{24}}{\mathit{\pi }}}{\mathit{\omega }}_{\mathrm{s}}.\end{array}$$

The parameters a and b in Eq. (11) are regression parameters that usually have to be fitted to observed global radiation. If no actual solar radiation data are available, the values a = 0.25 and b = 0.50 are recommended (Allen et al., 1998). Though the Ångström–Prescott model has the disadvantage that the parameters have to be fitted and the actual duration of sunshine has to be observed, it has the benefit that both parameters allow for a direct physical interpretation. The parameter a represents the clearness index K_T on overcast days (n = 0), while their sum a+b gives the clearness index on clear days (n=N).

Commonly, in remote mountainous regions, only temperature data are available, due to which another group of solar radiation models is typically utilised. In these approaches the difference between daily maximum and minimum air temperature ΔT (^∘C) is used as a proxy for cloud cover, as clear-sky conditions result in a higher temperature amplitude between day and night compared to overcast conditions. Typical models are the exponential approach proposed by Bristow and Campbell (1984) and its later modifications or the simple empirical equation by Hargreaves and Samani (1982).

$$\begin{array}{}\text{(13)}& K_T=k_{\mathrm{H}}\sqrt{\mathrm{\Delta }T},\end{array}$$

with the empirical coefficient k_H = 0.16 for inland and k_H = 0.19 for coastal locations. Since the influence of cloud cover on the clearness index and thus on the DDF can be illustrated much more directly by Ångström–Prescott type models, this model type is used in the paper.

It is obvious that the attenuation of extraterrestrial solar radiation is a function of the distance the rays have to travel through the atmosphere, as absorption and scattering occurs all along the way. Several solar radiation models consider altitude as a variable, for which the models below were calibrated including high-altitude stations and are of Ångström–Prescott type; thus the altitude effects can be compared directly.

Jin et al. (2005):

$$\begin{array}{}\text{(14)}& {\displaystyle }\text{(a)}{K}_T=\left(\mathrm{0.0855}+\mathrm{0.0020}\mathrm{\varnothing }+\mathrm{0.030}z\right)+\mathrm{0.5654}{\displaystyle \frac{n}{N}},\text{(15)}& {\displaystyle }\begin{array}{rl}\text{(b)}{K}_T& =\left(\mathrm{0.1094}+\mathrm{0.0014}\mathrm{\varnothing }+\mathrm{0.0212}z\right)\\ & +\left(\mathrm{0.5176}+\mathrm{0.0012}\mathrm{\varnothing }+\mathrm{0.0150}z\right){\displaystyle \frac{n}{N}}.\end{array}\end{array}$$

Rensheng et al. (2006):

$$\begin{array}{}\text{(16)}& K_T=\left(\mathrm{0.122}+\mathrm{0.001}\mathrm{\varnothing }+\mathrm{0.0257}z\right)+\mathrm{0.543}{\displaystyle \frac{n}{N}}.\end{array}$$

Liu et al. (2019):

$$\begin{array}{}\text{(17)}& K_T=\left(\mathrm{0.1755}+\mathrm{0.0136}z\right)+\left(\mathrm{0.5414}+\mathrm{0.0117}z\right){\displaystyle \frac{n}{N}}.\end{array}$$

For all models, z is the altitude (km) and ϕ the latitude (degree). In Eqs. (14) and (16) only the parameter a of Eq. (11) is a function of altitude z, while in Eqs. (15) and (17) also the parameter b is dependent on z.

In order to evaluate the altitude effect separately from other parameters, the clearness index K_T is split into two components:

$$\begin{array}{}\text{(18)}& K_T=K_{T_{\mathrm{0}}}\cdot K_z,\end{array}$$

where $K_{T_{\mathrm{0}}}$ is the clearness index at z = 0 m a.s.l. and K_z is a clearness altitude factor (–), which represents the increase in K_T with altitude relative to $K_{T_{\mathrm{0}}}$. At sea level, K_z = 1 for all models and all values of relative sunshine duration $n/N$. Though the clearness altitude factors K_z obtained from Eqs. (14)–(17) are different for each model, they all show a linear increase with altitude, the slope of which depends on the cloudiness (see Fig. 4). It should be noted that although K_z is higher for overcast than for clear-sky conditions, the absolute increase in the clearness index K_T with altitude is higher under clear-sky conditions because of the higher base value $K_{T_{\mathrm{0}}}$. For example, at z = 2000 m a.s.l., the Jin (b) (Eq. 15) model has a clearness altitude factor K_z of 1.27 for overcast and 1.10 for clear-sky conditions. However, when multiplying by the respective clearness factors at sea level $K_{T_{\mathrm{0}}}$ of 0.15 and 0.72, the resulting clearness indices K_T at z= 2000 m a.s.l. increase by 0.04 under overcast and 0.07 under clear-sky conditions to an absolute value of 0.19 and 0.79, respectively.

Figure 4 Clearness altitude factor K_z for latitude 45^∘ and different altitudes, based on different models presented in Eqs. (14)–(17) (i.e. Jin (a), Jin (b), Rensheng, and Liu).

Albedo

While the albedo of fresh snow is well above 0.9 (Hock, 2005), indicating that most of the shortwave radiation is reflected, it may drop significantly within a few days due to snow metamorphism. Well-aged snow generally has an albedo in the range of 0.4–0.5 (Anderson, 2006). Snow albedo is primarily dependent on the grain size of the snow crystals near the surface but also on aerosols in the snow and dust deposits. Respective snow albedo models are proposed by Wiscombe and Warren (1980), Warren and Wiscombe (1980), or Amaral et al. (2017). However, because of their data requirements, surrogate exponential decay models as formulated by USACE (1956) are commonly in use, which assumes the decrease in albedo as a function of time after the last significant snowfall. For example, Walter et al. (2005) use the empirical relationship

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{A}_n& =\mathrm{0.35}-\left(\mathrm{0.35}-A_{max}\right)\\ & \exp {\left[-\left(\mathrm{0.177}+\ln {\left({\displaystyle \frac{A_{max}-\mathrm{0.35}}{A_{n-\mathrm{1}}-\mathrm{0.35}}}\right)}^{\mathrm{2.16}}\right)\right]}^{\mathrm{0.46}},\end{array}\end{array}$$

where A_n−1 is the albedo of the previous day and A_max is the maximum albedo (∼ 0.95) of fresh snow. Following Eq. (19), the snow albedo will decrease from 0.95 to 0.52 after 10 d and to 0.43 after 30 d if no new snowfall occurs.

3.2.2 Longwave radiation

The net longwave radiation flux over the snow surface Q_L (W m⁻²) is the balance between incoming longwave radiation that is emitted by the atmosphere Q_L,in (W m⁻²) and outgoing radiation from the snowpack Q_L,out (W m⁻²).

$$\begin{array}{}\text{(20)}& Q_{\mathrm{L}}=Q_{\mathrm{L},\mathrm{in}}-Q_{\mathrm{L},\mathrm{out}}\end{array}$$

Longwave radiation is a function of the temperature of the emitting body and can be calculated with the Stefan–Boltzmann law:

$$\begin{array}{}\text{(21)}& L=\mathit{\epsilon }\mathit{\sigma }{T}^{\mathrm{4}},\end{array}$$

where L is the radiative flux (W m⁻²), ε and T are the emissivity (–) and the absolute temperature (K) of the emitting body, respectively, and σ is the Stefan-Boltzmann constant ($\mathrm{5.67}\times {\mathrm{10}}^{-\mathrm{8}}$ W m⁻² K⁻⁴).

In particular, fresh snow is nearly a perfect blackbody with respect to longwave radiation with a high emissivity of 0.99 (Warren, 1982; USACE, 1998; Anderson, 2006). For old snow, Brutsaert (1982) gives an emissivity value of 0.97. Given a melting snowpack having a surface temperature of 0 ^∘C, the outgoing energy flux can be taken as constant with $Q_{\mathrm{L},\mathrm{out}}\sim \mathrm{310}$ W m⁻².

For the atmospheric longwave radiation, usually the air temperature T_a (K) is used in Eq. (21). However, while the snowpack longwave emissivity is virtually constant, the emissivity of the atmosphere is highly variable. Typical values under clear-sky conditions range from 0.6–0.8, primarily depending on air temperature and humidity (Anderson, 2006), whereas for overcast conditions it can be close to 1.0.

A number of empirical and more physically based approaches exist to estimate atmospheric longwave emissivity from standard meteorological data (see Hock, 2005, for a discussion). For clear-sky conditions, Brutsaert (1975) developed a theoretically based formula depending on air temperature and vapour pressure measured at screen level:

$$\begin{array}{}\text{(22)}& {\mathit{\epsilon }}_{\mathrm{ac}}=\mathrm{1.24}{\left({\displaystyle \frac{p_{\mathrm{v}}}{T_{\mathrm{a}}}}\right)}^{\frac{\mathrm{1}}{\mathrm{7}}},\end{array}$$

where ε_ac is the clear-sky longwave emissivity (–), p_v the actual vapour pressure (hPa), and T_a the air temperature (K). Later, Brutsaert reconciled Eq. (22) with an empirical approach proposed by Swinbank (1963),

$$\begin{array}{}\text{(23)}& {\mathit{\epsilon }}_{\mathrm{ac}}=\mathrm{9.2}\times {\mathrm{10}}^{-\mathrm{6}}{T}_{\mathrm{a}}^{\mathrm{2}},\end{array}$$

that considers the strong correlation between vapour pressure and air temperature; thus, only air temperature is needed as input variable. Using the above relation, at an air temperature of 10 ^∘C the atmospheric longwave radiation flux into the snowpack amounts to Q_L,in = 281 W m⁻² under clear-sky conditions, which is less than the outgoing flux of 310 W m⁻²; i.e. the snowpack will lose energy in this situation.

The variability of atmospheric emissivity due to cloud cover, which increases the longwave emissivity, is significantly higher than variations under clear-sky conditions. Monteith and Unsworth (2013) give the simple linear relationship

$$\begin{array}{}\text{(24)}& {\mathit{\epsilon }}_{\mathrm{a}}=\left(\mathrm{1}-\mathrm{0.84}c\right){\mathit{\epsilon }}_{\mathrm{ac}}+\mathrm{0.84}c,\end{array}$$

where ε_a is the atmospheric longwave emissivity, c is the fraction of cloud cover (–), and ε_ac is calculated by Eqs. (22) or (23). For overcast conditions and an air temperature of 10 ^∘C, Eq. (24) yields an atmospheric emissivity of 0.96, which results in an atmospheric longwave radiation flux of Q_L,in = 351 W m⁻² and thus a positive flux of Q_L = 41 W m⁻² into the snowpack.

Although cloud cover is difficult to parameterise, as clouds can be highly variable in space and time and their effects on radiation depend on the different cloud genera, a strong correlation between cloud cover and sunshine duration is obvious. Doorenbos and Pruitt (1977) give a tabulated relation between cloudiness c and relative sunshine hours $n/N$ (see Eq. 11), which can be fitted by the quadratic regression

$$\begin{array}{}\text{(25)}& c=\mathrm{1}-\mathrm{0.5544}{\displaystyle \frac{n}{N}}-\mathrm{0.5483}{\left({\displaystyle \frac{n}{N}}\right)}^{\mathrm{2}}.\end{array}$$

Nevertheless, in simple sky models usually a linear relation between cloudiness and relative sunshine hours is applied as a first approximation (e.g. Brutsaert, 1982; Annandale et al., 2002; Pelkowski, 2009), which, as Badescu and Paulescu (2011) showed by using probability distributions to develop relations between cloudiness and relative sunshine hours, is a first good estimate.

3.2.3 Sensible heat exchange

Sensible heat exchange describes the energy flux due to temperature differences between the air and the snow surface while air is permanently exchanged by wind turbulence. A frequent approach to parameterise turbulent heat transfer is the aerodynamic method, which explicitly includes wind speed as a variable (Braithwaite et al., 1998; Lehning et al., 2002; Hock, 2005):

$$\begin{array}{}\text{(26)}& Q_{\mathrm{H}}={\mathit{\rho }}_{\mathrm{a}}{c}_p{C}_{\mathrm{H}}u\left(T_{\mathrm{a}}-T_{\mathrm{s}}\right),\end{array}$$

where ρ_a is the air density (kg m⁻³), c_p the specific (isobaric) heat capacity of air (1006 J kg^−1 ∘C⁻¹), C_H the exchange coefficient for sensible heat (–), u the mean wind speed (m s⁻¹), T_a the air temperature (^∘C), and T_s the temperature at the snow surface (^∘C).

The density of air ρ_a is a function of atmospheric pressure, air temperature, and humidity:

$$\begin{array}{}\text{(27)}& {\mathit{\rho }}_{\mathrm{a}}={\displaystyle \frac{M_{\mathrm{d}}\left[p-\left(\mathrm{1}-e\right)p_{\mathrm{v}}\right]}{R{T}_{\mathrm{a}}}},\end{array}$$

where p is the atmospheric pressure (Pa), p_v the vapour pressure (Pa) (see Eq. 32), T_a the air temperature (K), M_d the molar mass of dry air (0.02897 kg mol⁻¹), R the universal gas constant (8.31446 J mol⁻¹ K⁻¹), and e the ratio of molar weights of water and dry air equal to 0.622. At usual air temperatures humidity has only a minor effect on the air density.

The decrease in atmospheric pressure with altitude z (m a.s.l.) can be estimated by the isothermal barometric formula

$$\begin{array}{}\text{(28)}& p\left(z\right)=p_{\mathrm{0}}\exp \left(-{\displaystyle \frac{g{M}_{\mathrm{d}}}{R{T}_{\mathrm{a}}}}z\right),\end{array}$$

where p₀ is the atmospheric pressure at sea level (Pa) and g the gravitational acceleration (m s⁻²). At an air temperature of 0 ^∘C and a standard atmospheric pressure at sea level of 101.325 kPa, the air density is 1.29 kg m⁻³, while for example at an altitude of 2000 m a.s.l. the atmospheric pressure reduces to 78.9 kPa and the air density becomes 1.01 kg m⁻³.

The exchange coefficient C_H can be approximated with (Campbell and Norman, 1998)

$$\begin{array}{}\text{(29)}& C_{\mathrm{H}}={\displaystyle \frac{k^{\mathrm{2}}}{\ln \left(\frac{z_u}{z_{\mathrm{m}}}\right)\ln \left(\frac{z_T}{z_{\mathrm{h}}}\right)}},\end{array}$$

where k is the von Kármán constant 0.41 (–), z_u and z_T the height of wind and temperature observation above the snow surface (m), z_m the momentum roughness parameter, and z_h the heat roughness parameter. For a snow surface, the roughness parameters are given by Walter et al. (2005) as z_m∼0.001 m and z_h∼0.0002 m.

As it can be seen from Eq. (26), the sensible heat component depends mainly on wind speed and temperature. During stable clear weather periods with typically light winds, the turbulent exchange is smaller on average than the radiation components. For example, a wind speed of 1 m s⁻¹ and an air temperature of 5 ^∘C will result in a sensible heat flux of about 15.5 W m⁻². However, at warm rain events or at föhn conditions with strong warm winds, turbulent exchange can significantly contribute to the melting process. For example, a föhn event of 14 h duration on 8 December 2006 at Altdorf (Switzerland, 440 m a.s.l.) with an average air temperature of about 16 ^∘C, average relative humidity of 37 %, and average wind speed of 14.6 m s⁻¹ resulted in a mean sensible heat flux of about 700 W m⁻² during that duration.

3.2.4 Latent energy of condensation or vaporisation

The latent energy exchange reflects the phase change of water vapour at the snow surface, either by condensation of vapour contained in the air or by vaporisation of snow. Thus, it can either warm or cool the snowpack (Harpold and Brooks, 2018). The energy flux is dependent on the vapour gradient between the air and the snow surface and is, like the sensible heat exchange, a turbulent process that increases with the wind speed. Thus, the aerodynamic formulation is analogous to Eq. (26):

$$\begin{array}{}\text{(30)}& Q_{\mathrm{E}}={\mathit{\rho }}_{\mathrm{a}}{\mathit{\lambda }}_{\mathrm{v}}{C}_{\mathrm{E}}u\left(q_{\mathrm{a}}-q_{\mathrm{s}}\right),\end{array}$$

where λ_v is the latent heat of vaporisation of water at 0 ^∘C (2.501×10⁶ J kg⁻¹); C_E the exchange coefficient for latent heat (–), which is assumed to be equal to the exchange coefficient for sensible heat C_H; q_a the specific humidity of the air (–); and q_s the specific humidity at the snow surface (–).

The specific humidity q_a can be derived from measurements of relative humidity or dew point temperature. In cases where such data are not available, Walter et al. (2005) approximate the dew point temperature by the minimum daily temperature. For any air temperature T (^∘C), the saturation vapour pressure p_s (Pa) can be calculated by an empirical expression known as the Magnus–Tetens equation in the general form (Lawrence, 2005)

$$\begin{array}{}\text{(31)}& p_{\mathrm{s}}=C{e}^{\frac{AT}{B+T}},\end{array}$$

where A, B, and C are coefficients after Allen et al. (1998): A = 17.2694, B = 237.3 ^∘C, and C = 610.78 Pa. At the snow surface, according to Lehning et al. (2002) the air temperature can be assumed equal to the snow surface temperature, and Eq. (31) is applied with coefficients for saturation vapour pressure over ice: A = 21.8746, B = 265.5 ^∘C, C = 610.78 Pa (Murray, 1967). At a temperature of 0 ^∘C, both coefficient sets yield the same saturation vapour pressure of p_s = 611 Pa.

Knowing the relative humidity ψ (–) and the saturation vapour pressure p_s at a given air temperature, the actual vapour pressure p_v (Pa) can be calculated through the relation

$$\begin{array}{}\text{(32)}& p_{\mathrm{v}}=\mathit{\psi }{p}_{\mathrm{s}}\end{array}$$

and subsequently the respective specific humidity by

$$\begin{array}{}\text{(33)}& q={\displaystyle \frac{e{p}_{\mathrm{v}}}{p-\left(\mathrm{1}-e\right)p_{\mathrm{v}}}}\approx {\displaystyle \frac{e}{p}}{p}_{\mathrm{v}},\end{array}$$

with p the atmospheric pressure (Pa) and e the ratio of molar weights of water and dry air equal to 0.622 as in Eq. (27). Assuming melting conditions with a snow temperature T_s = 0 ^∘C and saturated vapour conditions, the vapour pressure at the snow surface is $p_{\mathrm{v},\mathrm{snow}}=p_{\mathrm{s}}$ (0 ^∘C) = 611 Pa. While at positive air temperatures the sensible heat flux is always warming the snowpack, the latent heat flux can cool the snow by vaporisation if the relative humidity of the air is low. Even when assuming a relative humidity of 100 %, the latent heat flux into the snowpack will be comparatively small if wind speed is low, e.g. about 13 W m⁻² at an air temperature of 5 ^∘C and a wind speed of 1 m s⁻¹.

3.2.5 Precipitation heat

The heat transfer into the snowpack by lowering the rain's temperature, which is usually assumed to be equal to the air temperature T_a (^∘C), down to the freezing point at 0 ^∘C can be estimated as

$$\begin{array}{}\text{(34)}& Q_{\mathrm{P}}=c_{\mathrm{w}}P{T}_{\mathrm{a}},\end{array}$$

where c_w is the specific heat capacity of water (4.2 kJ kg^−1 ∘C⁻¹) and P is the daily rainfall depth (kg m⁻² d⁻¹). The energy input from precipitation is usually quite small, and even during extreme weather conditions, like heavy warm rain storms with temperatures of 15 ^∘C and a precipitation depth of 50 mm, that may occur, e.g. during early winter in the Alps, the mean daily energy flux from rain would be a moderate 36.5 W m⁻².

3.2.6 Change in internal energy

The rate of change in the energy stored in the snowpack ΔQ (W m⁻²) represents the internal energy gains and losses due to changes in the snowpack's temperature profile and due to phase changes, i.e. melting of the ice portion or refreezing of liquid water in the snowpack. Until the snowpack temperature is isothermal at 0 ^∘C, any melt produced in the surface layer that exceeds the liquid water holding capacity of the porous snow matrix will percolate downward and will be captured and refrozen in colder lower layers. This internal mass and energy transport process absorbs at least parts of the incoming energy, which reduces the energy available for melt and will thus reduce the actual DDF.

Under data-scarce conditions and particularly when only daily data are available, it is difficult to properly quantify the change in the internal energy of the snowpack (see discussion in Sect. 5.2.1). Therefore, in our study we focus on melt periods when the snowpack is “ripe”, i.e. when the snowpack temperature is isothermal at 0 ^∘C and the residual volumetric water content of about 8 % (Lehning et al., 2002) is filled with liquid water. This assumption is not a limitation when analysing the contribution of each individual energy flux component towards a resulting DDF as presented in following sections, but the additional energy needed for warming the snowpack has to be taken into account when estimating the total DDF if a snowpack is not ripe (see Fig. 11).

4 Results

In this section, the contribution of each energy flux component Q_i to the lumped daily DDF is presented. For this purpose, the respective melt depth M_i is calculated according to Eq. (3) and further converted into the corresponding degree-day factor component DDF_i using Eq. (2). For the following exemplary calculations, the air temperature is assumed to stay always above 0 ^∘C; thus, degree days T_DD (^∘C d) in Eq. (2) have the same numerical value as the daily average air temperature T_a (^∘C) used in the calculation of several energy flux components.

Besides demonstrating the dependency of the DDF components on decisive parameters of the energy flux components, the presented tables in the Supplement (“S” in numbering of figures and tables means that they are in the Supplement) and graphs in this section, which are based on the relationships given in Sect. 3, can be used to estimate the DDF component values if observed data are either not available or not sufficient for more sophisticated approaches. It should be noted that parameters are normalised where applicable, i.e. set to hypothetical values like clearness index K_T = 1 or wind speed u = 1.0 m s⁻¹; thus, final DDF values can be obtained by multiplying the given figures by the actual values of those parameters. Furthermore, all results are based on the assumption that the snowpack is isothermal at 0 ^∘C and in a fully ripe state.

4.1 Shortwave radiation component – DDF_S

Shortwave-radiation-induced melt is usually considered the largest DDF component, especially at higher elevations as well as under dry climates. The net energy flux Q_S is calculated using Eq. (5), which consists of three factors: (a) latitude, (b) albedo, and (c) clearness index K_T. The dependency of DDF_S on these factors is demonstrated in Fig. 5 for the period between winter solstice (21 December) and summer solstice (21 June). As shortwave radiation is independent of air temperature and hence of degree days, the corresponding melt is divided by a hypothetical degree-day value of 1 ^∘C d to arrive at the presented DDF_S values. In the event of actually higher degree days, the given DDF_S values have to be divided accordingly.

Figure 5 Variation of solar-radiation-based DDF_S for a degree-day value of 1 ^∘C d for (a) different latitudes under constant snow albedo and clearness index, (b) snow albedos under constant latitude and clearness index, and (c) different clearness indices under constant latitude and snow albedo. The latitude = 48^∘ corresponds to the location of the Brunnenkopfhütte test site.

Figure 5a shows the variation of DDF_S depending on latitude for the range 30–60^∘ N, while albedo (A = 0) and clearness index (K_T = 1) are set constant. Obviously, there is a significant difference in DDF_S for different latitudes around the winter solstice due to solar inclination, making latitude the predominant factor for DDF_S at this time of the year. However, around the summer solstice, DDF_S has nearly the same value at different latitudes because the lower solar angle at higher latitudes is counterweighted by a larger hour angle, i.e. longer sunlight hours. Thus, with the progress of the melting season, the factors albedo and clearness index become more important than latitude.

Figure 5b shows the influence of albedo on the DDF_S at a given latitude (Brunnenkopfhütte test site – latitude 48^∘) and normalised constant clearness index (K_T = 1). Snow albedo varies between 0.9–0.4, covering the range between fresh and well-aged snow. As expected, the influence of albedo increases with increasing incoming solar radiation towards the summer solstice. A good estimate of albedo is therefore much more important when the snowmelt season progresses than in early spring. If for example the same degree-day value of 10 ^∘C d is assumed on 21 March and on 21 May, the difference in DDF_S between fresh (A = 0.9) and aged (A = 0.4) snow would be 0.8 and 4.6 mm ^∘C⁻¹ d⁻¹ in March compared to 1.2 and 7.1 mm ^∘C⁻¹ d⁻¹ in May, respectively.

The dependency of DDF_S on the clearness index K_T is shown in Fig. 5c. In line with Eq. (6), DDF_S values under clear sky (K_T = 0.75) are always higher than under overcast conditions (K_T = 0.25). Similar to albedo, the influence of the clearness index becomes more pronounced, with increasing solar angle when the snowmelt season progresses.

The influence of altitude on DDF_S in terms of increasing K_T values can be assessed by multiplying a clearness index $K_{T_{\mathrm{0}}}$ at sea level, which may be obtained by any of the numerous solar radiation models, with a clearness altitude factor K_z (see Eq. 18). Figure 4 shows the range of clearness altitude factors for latitude 45^∘ derived from Eqs. (14)–(17). All K_z values show a linear increase with altitude, with the slope depending on cloudiness. It should be noted that although the increase in K_z relative to $K_{T_{\mathrm{0}}}$ is higher under overcast than under clear-sky conditions, the absolute increase in the clearness index K_T with altitude is larger for clear-sky conditions (see Sect. 3.2.1). When using the intersection of all models and sky conditions, which is indicated by the dark grey area in Fig. 4, in order to get one overall rough estimate of K_z for all conditions, the clearness altitude factor and thus the resulting DDF_S are found to increase by about 6.4 % per each 1000 m of altitude.

4.2 Longwave radiation component – DDF_L

The net longwave energy flux Q_L is calculated using Eq. (21), in which the outgoing radiation from the snowpack can be assumed to be constant. Thus, the contribution of the longwave radiation component DDF_L is mainly dependent on the air temperature and the emissivity of the atmosphere (e.g. in particularly cloudiness conditions). Figure 6 and Table S1 present the DDF_L as a function of degree days T_DD and cloudiness. For a wide range of degree days, especially in conjunction with low cloudiness, the outgoing longwave energy flux is higher than the incoming, resulting in a theoretically negative degree-day factor that will reduce the total DDF. This means that the DDF_L component under clear-sky conditions is usually contributing to a cooling of the snowpack rather than to melting. Under overcast conditions, the DDF_L is relatively constant around 1 mm ^∘C⁻¹ d⁻¹ with a maximum value of 1.3 mm ^∘C⁻¹ d⁻¹ at T_DD = 20 ^∘C d. Although this contribution to the total DDF is small compared to the shortwave radiation component DDF_S, it can be of importance at the onset of snowmelt in early spring, when the solar radiation is still low and the albedo of fresh snow is high.

Figure 6 Longwave radiation component (DDF_L) for selected cloudiness (%) and degree days (^∘C d).

4.3 Sensible heat component – DDF_H

The sensible heat flux Q_H as given by Eq. (26) is mainly proportional to wind speed and the temperature difference between the air and the snow surface. Furthermore, air density, besides its dependency on temperature, is a function of relative humidity and atmospheric pressure, and thus of altitude (Eqs. 27 and 28). Since the influence of the relative humidity on air density is negligible, a relative humidity of RH = 0 % is assumed in the below analysis on the response of the DDF_H to changes in temperature and degree days, wind speed, and altitude. It should be noted that this analysis assumes typical melt conditions with a snowpack temperature of T_s = 0 ^∘C and positive air temperature, whereas negative air temperature would lead to a negative sensible heat flux resulting in a cooling of the snowpack and a decrease in total DDF.

Figure 7 and Table S2 show the variation in DDF_H depending on altitude and degree days, while the wind speed is assumed to be constant at u=1 m s⁻¹. The latter allows the DDF_H to be easily calculated for any other wind speed by multiplying the given value by the actual wind speed. The DDF_H principally decreases with altitude, with less pronounced differences due to temperature at higher altitudes. However, the most important factor is the wind speed. For example, with a daily average wind speed of u = 1.0 m s⁻¹, e.g. at sea level z = 0 m a.s.l., the DDF_H only decreases from 0.806 to 0.781 mm ^∘C⁻¹ d⁻¹ when degree days increase from 1 to 10 ^∘C d (see Fig. 7 and Table S2). In contrast, for a degree day of 1 ^∘C d, the DDF_H increases proportionally from 0.806 to 8.061 mm ^∘C⁻¹ d⁻¹ when wind speed increases from 1 to 10 m s⁻¹. Thus, wind speed is a decisive variable when estimating the DDF_H.

If wind speed observations are not available, they may be roughly estimated based on the topographic and climate characteristics of the study area. However, average values may not represent the actual wind conditions and thus DDF_H on a certain day. While for example the geometric mean of observed daily wind speed at the Brunnenkopfhütte station is about 0.8 m s⁻¹, resulting in a DDF_H of approx. 0.7 mm ^∘C⁻¹ d⁻¹, the maximum daily average wind speed is about 4.5 m s⁻¹, which increases DDF_H to approx. 3.9 mm ^∘C⁻¹ d⁻¹.

Figure 7 Variation of sensible heat component (DDF_H) at different altitude based on different degree days, RH = 0 %, and u = 1 m s⁻¹.

4.4 Latent heat component – DDF_E

The latent heat flux Q_E approximated by an aerodynamic model as in Eq. (30) indicates that the latent heat component DDF_E is mainly dependent on the humidity gradient near the snow surface and on the wind speed. Additionally, altitude has an influence, as the air density decreases with altitude. Figure 8 and Table S3 give the resulting DDF_E as a function of degree days for different values of relative humidity and at daily average wind of u = 1.0 m s⁻¹, whereas air density values are assumed at an elevation of 0 m a.s.l. In line with the sensible heat component DDF_H, DDF_E for any other wind speed can be obtained by multiplication by the actual value. For relative humidity < 30 % the DDF_E is negative over the whole range of degree days; hence, the latent heat component will reduce the total DDF under these conditions. Even if the air is humid and warm, contribution of latent heat is moderate, e.g. DDF_E=1.0 mm ^∘C⁻¹ d⁻¹ at a relative humidity of 100 % and T_DD = 20 ^∘C d.

Figure 8 Latent heat component (DDF_E) for selected relative humidity (%), degree days (^∘C d), and u = 1 (m s⁻¹), These values are for u=1 m s⁻¹; for a different wind speed these values can be multiplied for the desired wind speed. Air density values are assumed at an elevation of 0 m a.s.l.

Figure 9 Variation of latent heat component (DDF_E) depending on altitude for different relative humidity values, degree days, and u = 1 m s⁻¹.

Figure 9 shows the combined effect of altitude, relative humidity, and temperature on DDF_E. At a high relative humidity (e.g. RH = 100 %), similar to the DDF_H the DDF_E values principally decrease with altitude, with less pronounced differences due to temperature at higher altitudes. At lower relative humidity (e.g. RH = 50 %), the altitude effect is less noticeable, and at low temperatures even a reversal of the effect can be observed. Thus, altitude reduces the positive DDF_E associated with high humidity, while it also reduces the cooling effect of a negative latent heat flux, which is associated with low humidity and lower air temperature.

As the above analysis shows, humidity is the main variable influencing the DDF_E. In general, humid air will promote condensation at a cooler snow surface, which releases latent energy and contributes to a positive DDF, while dry air will promote evaporation and sublimation from the snow surface, which abstracts energy from the snowpack. Thus, mainly depending on the humidity of the air, the latent heat energy flux is usually a heat sink, while only in the event of high humidity in conjunction with higher temperature it becomes a heat source to the snowpack. Especially in spring, when relative humidity is comparatively low in middle and northern latitudes, large parts of the incoming solar radiation can be consumed by evaporation from the snow surface, significantly reducing the energy available for melt and thus reducing the corresponding DDFs (Lang and Braun, 1990; Zhang et al., 2006).

4.5 Precipitation heat component – DDF_P

Rainfall can affect the snowpack energy budget by adding sensible heat due to warm rain and by releasing latent heat if the rain refreezes in the snowpack (DeWalle and Rango, 2008). The latter effect is not considered in this study, as the snowpack is assumed at a 0 ^∘C melting condition. Given that the precipitation heat Q_P is linearly dependent on air temperature in Eq. (34), a division by respective degree days makes DDF_P independent of temperature and proportional to rainfall, resulting in a DDF_P = 0.0125 mm ^∘C⁻¹ d⁻¹ for a precipitation depth of 1 mm d⁻¹. DDF_P for any other precipitation can be obtained by respective multiplication. The exemplary values in Table 2 show, however, that the contribution of the precipitation heat component DDF_P is modest compared to other DDF components. Even high rainfall of 50 mm d⁻¹ would release only a small amount of sensible heat, resulting in a DDF_P of 0.6 mm ^∘C⁻¹ d⁻¹.

Table 2 Precipitation heat component (DDF_P) (mm ^∘C⁻¹ d⁻¹) for selected precipitation (mm d⁻¹).

5 Discussion

While the previous section focuses on the characteristics of each individual energy-flux-based DDF component, this section mainly discusses the influence of spatial, seasonal, or meteorological conditions on the overall DDF. The discussion section bifurcates into two sub-sections: (i) influence of individual factors on the DDF, such as latitude, altitude, albedo, season, and rain on snow events, and (ii) application of energy-flux-based DDF estimates, discussing how this value can be estimated for a temperature-index model by using different available datasets and applied under varying meteorological and climate change conditions.

5.1 Influence of individual factors on the DDF

In this section all conclusions are under the assumption that the snowpack is isothermal at T_s=0 ^∘C and in ripe condition; hence, all net incoming energy is available for melt and contributes to the total DDF. Apart from the discussed variables, we assumed the standard values of u = 1 m s⁻¹, RH = 70 %, A=0.5, and P = 0 mm and typical melt conditions of T_DD = 5 ^∘C d, unless otherwise stated.

5.1.1 Influence of latitude

While topographic factors like slope, aspect, or shading in mountainous regions result in a high local variability of melt conditions, larger-scale regional patterns of DDFs (e.g. a dependency on latitude) could not be detected in a data review by Hock (2003). This observation is supported by a brief analysis of the effect of latitude below, where the DDF is compared not on the same date but on the same degree days. As an illustrative example, typical melt conditions of T_DD = 5 ^∘C d at a latitude of about 35^∘ N in the upper Jhelum catchment (Bogacki and Ismail, 2016) are compared to similar conditions at a latitude of 48^∘ N (Brunnenkopfhütte, 1602 m a.s.l.). As zone-wise temperature data (see Sect. 2.2) indicate, in the upper Jhelum catchment at an elevation zone of 1500–2000 m a.s.l., melting conditions usually occur around mid-February, while at Brunnenkopfhütte comparable degree days are obtained about 1 month later in mid-March. Figure 10a compares the energy-flux-based DDF components at both latitudes. The decisive solar radiation component is very similar at the two locations, both under clear-sky and overcast conditions; thus the total DDF is virtually identical at both latitudes. Therefore, at least in moderate latitudes and when compared under similar melt conditions, no significant effect of latitude on DDF could be found.

Figure 10 Influence of (a) latitude, (b) altitude, (c) albedo, (d) season, and (e) rain on snow events on the DDF under clear-sky and overcast conditions. DDF components: shortwave radiation (DDF_S), longwave radiation (DDF_L), sensible heat (DDF_H), latent heat (DDF_E), precipitation heat (DDF_P) and increase in shortwave component due to the increase in the clearness index with altitude (DDF_A).

5.1.2 Influence of altitude

Contrary to the compensating effect in the case of latitude, the delayed onset of snowmelt due to altitude influences the DDF noticeably, which becomes important in temperature-index models where calculation is usually based on elevation bands. In order to demonstrate the influence of altitude on the DDF, two elevation zones with an altitude of 1500–2000 and 3500–4000 m a.s.l., respectively, are compared at 35^∘ latitude in the upper Jhelum catchment. As indicated above, typical melt conditions of T_DD = 5 ^∘C d occur at 1500–2000 m a.s.l. usually around mid-February, while at 3500–4000 m a.s.l. similar degree days are obtained for mid-May. The resulting DDFs (see Fig. 10b) show a significant difference, both under clear-sky conditions and under overcast conditions, because of the different input in solar radiation caused by the alteration in solar angle between February and May. Figure 10b shows an additional term DDF_A on top of the solar radiation component that represents the increase in incoming solar radiation due to the clearness altitude factor, which takes into account the increase in the clearness index with altitude (see Sect. 3.2.1). Averaging the factors proposed by different solar radiation models (see Fig. 4) results in an additional component DDF_A of 0.4 and 1.4 mm ^∘C⁻¹ d⁻¹ under clear-sky conditions and of 0.5 and 1.6 mm ^∘C⁻¹ d⁻¹ under overcast conditions at 1500–2000 and 3500–4000 m a.s.l., respectively.

While snow albedo is assumed constant at 0.5 in Fig. 10b, taking into consideration the decrease in albedo as the snow ages (see Table 1, e.g. A=0.74 in February and A=0.42 in May) results in a more pronounced difference with altitude, i.e. a DDF of 0.3 compared to 10.5 mm ^∘C⁻¹ d⁻¹ under clear-sky conditions and of 2.7 versus 7.3 mm ^∘C⁻¹ d⁻¹ under overcast conditions for the two altitudes, respectively.

The increase in DDF with increasing altitude is in line with previous studies (e.g. Hock, 2003; Kayastha and Kayastha, 2020). For example, in the Nepalese Himalayan region, seasonal-average DDF increases from 7.7 to 11.6 mm d^−1 ∘C⁻¹ with respect to altitude ranging 4900 to 5300 m a.s.l. (Kayastha et al., 2000), whereas Kayastha and Kayastha (2020) found that the model-calibrated range of the DDF in the central Himalayan basin varies between 7.0–9.0 mm d^−1 ∘C⁻¹ over an approximate altitude range of 4000–8000 m a.s.l. As Kayastha et al. (2000) pointed out, higher values of the DDF usually occur at very low temperatures since at higher altitudes the major driving factor for melt is the energy input by solar radiation.

5.1.3 Influence of albedo

As already discussed in the sections before, snow albedo is a critical parameter for the DDF since according to Eq. (5) the albedo directly controls the net solar radiation flux into the snowpack. While albedo of fresh snow is well above 0.9, hence reflecting most of the incoming shortwave radiation, it drops rapidly when larger grains form due to snow metamorphism. Figure 10c demonstrates the effect of ageing snow after a new snow event, when a simple exponential decay model as given in Eq. (19) is used and typical melting conditions T_DD = 5 ^∘C d are assumed. Since directly after a new snow event (Day = 0) the fresh snow albedo is high (A=0.95), the overall DDF is generally small. Under clear-sky conditions, if longwave radiation cooling is larger than net shortwave radiation flux, even a negative DDF value, i.e. no melt, may occur. If there is no new snow event in between, albedo will decrease following the exponential decay model to 0.52 after 10 d, resulting in a DDF of 5.8 mm ^∘C⁻¹ d⁻¹ under clear-sky conditions and 4.4 mm ^∘C⁻¹ d⁻¹ under overcast conditions. The increase in the DDF with exponential decay in albedo is in agreement with the findings of MacDougall et al. (2011), who found that the DDF is sensitive to albedo with values of > 4.0 mm d^−1 ∘C⁻¹ at an albedo of 0.6. As described qualitatively in the literature (e.g. Hock, 2003), under all sky conditions the DDF continuously increases with decreasing albedo, with the increase, however, being more pronounced under clear-sky than under overcast conditions.

5.1.4 Influence of season

Since the solar angle rises from its minimum at the winter solstice in December to its maximum on 21 June, the solar radiation component DDF_S increases during the snowmelt season, and thus the DDF is expected to increase. Figure 10d shows the influence of season on the DDF at the Brunnenkopfhütte test site during the melt period, assuming average degree days of 1, 4, and 7 ^∘C d in March, April, and May, respectively (see Table 1). Under clear-sky conditions, the total DDF increases from a negative value of −3.6 mm ^∘C⁻¹ d⁻¹ in March to 6.6 mm ^∘C⁻¹ d⁻¹ in May. Under overcast conditions, however, the DDF is virtually stable, ranging from 4.4 to 4.5 mm ^∘C⁻¹ d⁻¹ in the same period. The stability of the DDF under overcast conditions found in our study is in agreement with the study of Kayastha et al. (2000), where the DDFs observed during July–August are small compared to June because of prevailing cloud cover due to monsoon activity, which reduces the incoming shortwave radiation.

An evaluation of the individual DDF components shows that under clear-sky conditions the high impact of solar radiation in combination with low degree days at the onset of the snowmelt season is counterweighted by a strong negative longwave radiation component that decreases as the season progresses. Under overcast conditions, DDF_L is neutral or slightly positive, while the DDF_S component decreases because degree days rise faster than the input from solar radiation, which implies that sky conditions (i.e. overcast, and clear sky) are more decisive for an estimate of the DDF than the day of the year.

The effect of cloud cover is further amplified by the decrease in albedo while the melt season progresses, which becomes more significant under clear-sky conditions. In the present example (with the average monthly albedo as specified in Table 1) only 30 % of incoming solar radiation is contributing to melt in March, while it is about 60 % in May, enhancing the marked increase in the DDF under clear-sky conditions.

5.1.5 Influence of rain on snow events

In general, the precipitation heat component alone has only a minor effect on the DDF. However, in conjunction with certain weather conditions, like breaking in of warm and moist air, rain over snow events may lead to sudden melt and severe flooding.

Figure 10e shows the different DDF components resulting from a hypothetical rain over snow events assuming an air temperature of 15 ^∘C, a precipitation of 70 mm d⁻¹, a daily average wind speed of 10 m s⁻¹, a relative humidity of 100 %, and overcast conditions. Although the amount of precipitation is substantial and the rain temperature is comparatively high, the contribution of DDF_P is still modest. However, air temperature, relative humidity, and in particular wind speed associated with such events increase the sensible and latent heat components significantly. Thus, the resulting overall DDF is much higher than under usual melt conditions, which may lead to a considerable melt that adds to the runoff already caused by the heavy rain.

5.2 Application of energy-flux-based DDF estimates

5.2.1 DDF estimates under field conditions

In addition to the analysis of the influence of individual factors on the DDF, the dataset from the Brunnenkopfhütte test site is used to compare energy-flux-based estimates with field-derived DDFs in order to demonstrate how naturally varying meteorological conditions during the melt season, and in particular the cold content of the snowpack, affect the accuracy of DDF estimates.

For this purpose, daily melt was estimated from the daily difference of observed snow water equivalent during melt periods (see Fig. 3). Energy-flux-based melt was calculated by the formulas given in Sect. 3 using observed daily data from the Brunnenkopfhütte automatic snow and weather station (e.g. air temperature, wind speed; see Sect. 2.1) where applicable.

The daily degree-day sum is calculated from hourly air temperature data as proposed by Braithwaite and Hughes (2022). In general, operational degree-day models typically use constant degree-day factors for a certain time period (e.g. 10 d period). In this backdrop, both energy-flux-based and data-derived daily melt values were accumulated on a 10 d basis and divided by the degree days of the respective period. The 10 d averaging procedure also smooths daily noise in the observed data, in particular inaccuracies in the determination of daily melt and unrealistic DDF values because of daily temperature averages just above 0 ^∘C.

The comparison between field-derived and estimated (energy-flux-based) DDFs (see Fig. 11) yields a fair agreement with bias = 0.14 mm ^∘C⁻¹ d⁻¹ between estimated and field-derived values, and root mean square error (RMSE) = 1.12 mm ^∘C⁻¹ d⁻¹. Noteworthy in Fig. 11 are the new snow events, where the snowpack is no longer ripe and a certain amount of the incoming energy is needed to bring the snowpack back to ripe state, thus does not contributing to melt. For these events all estimated DDFs considerably overestimate the field-derived ones and were thus excluded from the calculation of the error metrics.

Figure 11 Comparison of field-derived vs. estimated (energy-flux-based) 10 d DDF for the Brunnenkopfhütte test site (period: November 2016–May 2021) – hollow points represent DDFs during periods with new snow events (new snow: precipitation ≥ 5 mm d⁻¹).

It is of interest to estimate DDFs also in cases where the snowpack is not ripe, e.g. because of new snow events or due to radiational cooling during clear cold nights. An approach to account for the snowpack's energy deficit, i.e. the energy needed to bring the snowpack temperature isothermal at 0 ^∘C, is the concept of cold content (Marks et al., 1999; Schaefli and Huss, 2011). The cold content is usually either estimated as a function of meteorological parameters or calculated by keeping track of the residuals of the snowpack energy balance (Jennings et al., 2018). For the latter, the SNOWPACK model (Lehning et al., 2002) is an excellent tool, which provides a highly detailed simulation of the vertical mass, energy, and besides other state variables the snow temperature distribution inside a snowpack. However, SNOWPACK requires a considerable number of meteorological input variables and preferably at least hourly observations, both of which are usually not available in the context where degree-day models are employed.

Especially suited for data-scarce conditions, Walter et al. (2005) apply a lumped approach that accounts for the cold content by changing the (isothermal) snowpack temperature depending on the daily net energy flux. When the incoming energy flux is sufficient to raise the snow temperature to 0 ^∘C or when it is already at 0 ^∘C the day before, all additional available energy produces melt. This approach, which does not need any additional data, however, seems to significantly overestimate the snowpack temperature particularly in situations with negative energy fluxes at night but a positive daily net balance, as a comparison with SNOWPACK simulations using data from Brunnenkopfhütte shows (see Sect. S2 in the Supplement). Therefore, an appropriate parameterisation of the cold content under limited data availability that would enable satisfactory estimates of DDFs in situations when the snowpack is not completely ripe remains subject to further research.

5.2.2 DDF estimates for temperature-index modelling

Snowmelt runoff models using the temperature-index approach have proven to be useful tools for simulation and forecasting in large snow- or glacier-dominated catchments, particularly in remote mountainous regions where data are usually scarce. A good estimate of the degree-day factor as the decisive model parameter is important either to stay in a realistic range when calibrating this parameter or in the case of forecasting when estimating its changes while the season progresses. In order to demonstrate the alteration of DDFs over time and altitude, energy-flux-based DDFs are estimated using 10 d average temperature (i.e. period 2000–2015) for the key elevation zones in the upper Jhelum catchment (Bogacki and Ismail, 2016). In the current example, the upper Jhelum catchment is discussed because of elevation-zone-wise data availability in comparison to test site where only point data are available (for more details, see Sect. 2). Because of the lack of data other than temperature and precipitation, prevailing conditions during the melt season are crudely approximated by the standard conditions used in this section, assuming persistent clear-sky conditions and albedo declining according to Eq. (19) after the last fresh snow event just before the beginning of the melting period.

Figure 12 shows the overall picture of how the DDF for snow will change over time and under climate change (i.e. present, RCP2.6, and RCP8.5; for DDF estimates under climate change, see Sect. 5.2.3). For example, Fig. 12a shows the development of DDFs in the elevation zones over time. As expected, melt starts earlier in lower-elevation zones and successively progresses to higher altitudes. Interestingly, the DDF in the first 10 d period of melting in each elevation zone increases with altitude. This is a combined effect of (i) higher solar radiation input and decreasing albedo while the season progresses and (ii) the onset of melt in higher-elevation zones starting at a lower degree-day threshold than in lower zones. In contrast to Fig. 10d, the DDF in Fig. 12 decreases continuously in all elevation zones in the subsequent melting periods since air temperature and thus degree days rise faster than melt. The range of DDFs for snow estimated by the energy flux components is in good agreement with earlier studies for the Himalayan region, e.g. 7.7–11.6 mm d^−1 ∘C⁻¹ (Kayastha et al., 2000), 5–9 mm d^−1 ∘C⁻¹ (Zhang et al., 2006), 5–7 mm d^−1 ∘C⁻¹ (Tahir et al., 2011), and 7.0–9.0 mm d^−1 ∘C⁻¹ (Kayastha and Kayastha, 2020).

Figure 12 (a) DDF estimates for a temperature-index modelling in the present climate; (b) influence of climate change – 2071–2100 under RCP2.6; (c) influence of climate change – 2071–2100 under RCP8.5.

5.2.3 DDF estimates under the influence of climate change

Climate change will ultimately influence snowmelt patterns depending on the projected changes in temperature and precipitation. In recent studies, usually model parameters including DDFs are considered as constant when assessing the climate change impact on future water availability from snow and glacier fed catchments (Lutz et al., 2016; Hasson et al., 2019; Ismail et al., 2020). However, due to the physical processes on which they depend, these parameters are subject to climate change. In this section, an attempt is made to estimate the influence of climate change on the DDFs in different elevation zones. For this analysis, results from ISIMIP data (see Sect. 2.2), which predict the temperature change for the period 2071–2100 to ΔT = 2.3 ^∘C under RCP2.6 and ΔT = 6.5 ^∘C under RCP8.5, are added to the temperatures in the present climate for each elevation zone.

The first effect to be observed in Fig. 12b and c is the common finding that snowmelt will start earlier under climate change as temperatures rise earlier above freezing. In addition, due to being earlier in the year, the DDFs in corresponding elevation zones are generally smaller compared to the current climate, though there are some outliers at the start of melting, due to division by low degree-day values. In the case of the pessimistic RCP8.5 scenario (Fig. 12c), a seasonal snow cover will no longer be established in the lowest elevation zone (i.e. 2500–3000 m a.s.l.) as air temperature at this elevation is projected to stay well above freezing throughout the winter. In general, the results of this brief analysis indicate that the DDFs are expected to decrease under the influence of climate change, as the snowmelt season will shift earlier in the year when solar radiation is small and snow albedo values are expected to be on the higher side. Musselman et al. (2017) highlight similar findings about slower snowmelt in a warmer world due to a shift of the snowmelt season to a time of lower available energy. These results may contribute the important aspect to the recent discussion on the linearity of temperature-index models used for glacier mass balance predictions (Bolibar et al., 2022; Vincent and Thibert, 2022) that the key impact of climate change on the DDF is the shift in the melt season. This effect should be larger on snowmelt because the entire melt season will be shifted to a time with lower solar radiation, while the period of glacial melt in the present climate will be preserved and will only start earlier and last longer under climate change.

6 Conclusions

Degree-day models are common and valuable tools for assessing present and future water availability in large snow- or glacier-melt-dominated basins, in particular when data are scarce like in the Hindu Kush–Karakoram–Himalayas mountain ranges. The present study attempts to quantify the effects of spatial, temporal, and climatic conditions on the degree-day factor (DDF), in order to gain a better understanding of which influencing factors are decisive under which conditions. While this analysis is physically based on the energy balance, formulas with minimum data requirement for estimating the DDFs are used to account for situations where observed data are limited. In addition, resulting tables (see Sect. S1 in the Supplement) and graphs for typical melt conditions are provided for a quick assessment. A comparison between field-derived and estimated DDFs at the Brunnenkopfhütte test site shows a fair agreement with bias = 0.14 mm ^∘C⁻¹ d⁻¹ and RMSE = 1.12 mm ^∘C⁻¹ d⁻¹ over periods without new snow events, since fresh snow increases the cold content of the snowpack and contradicts the condition of the snowpack being ripe and isothermal at 0 ^∘C. Further research is needed for cases where, under the constraint of limited data availability, also changes in the cold content of the snowpack are to be considered, with a specific focus on approaches that parameterise the diurnal dynamics of vertical temperature distribution in the snowpack.

Furthermore, the use of these DDF estimates directly as a model parameter or the incorporation of an energy-balance-based DDF approach into a degree-day model is not intended. One important aspect of temperature-index models is that the DDF is a lumped parameter, which is usually subject to calibration and accounts for uncertainties in different variables and parameters, i.e. temperature estimates, runoff coefficients, etc. Thus, the DDFs estimated by the energy balance approach are rather aimed to validate the results of parameter calibration and to highlight necessary adjustments due to climate change.

The analysis of the energy-balance processes controlling snowmelt indicates that cloud cover is the most decisive factor for the dynamics of the DDF. Under overcast conditions, the contribution of shortwave radiation is comparatively low, whereas the other components are in general small. Therefore, the total DDF value is not very high, and variations due to other factors are usually limited, apart from exceptional rainstorm events, for which energy balance models are the more suitable approach.

Under clear-sky conditions, on the other hand, shortwave radiation is the most prominent component contributing to melt. The increase in solar angle while the melt season progresses in combination with declining albedo and a decreasing cooling effect by the longwave radiation component along with increasing air temperature leads to a pronounced temporal dynamic in the DDF. Whereas incoming solar radiation and net longwave radiation can be determined fairly accurate by under clear-sky conditions, albedo becomes the crucial parameter for estimating the DDF, especially when new snow events occur during the melt period.

Clear-sky conditions promote the effect of increasing DDF with altitude if similar melting conditions are compared, since melting temperatures arrive later in the season at higher altitudes. The opposite effect can be observed with regard to climate change. Under higher temperatures at a given altitude, climate change will shift the snowmelt season earlier in the year. Consequently, when comparing periods of similar degree days, our study suggests DDFs are to decrease, since solar radiation is to generally decrease and albedo to typically increase.

Therefore, and as pointed out by many researchers, the DDF cannot be considered as a constant model parameter. Rather, its spatial and temporal variability must be taken into account, especially when using temperature-index models for forecasting present or predicting future snowpack and glacier changes, as well as the resulting water availability projections.