Journal cover Journal topic
Hydrology and Earth System Sciences An interactive open-access journal of the European Geosciences Union
Journal topic
Hydrol. Earth Syst. Sci., 23, 1–17, 2019
https://doi.org/10.5194/hess-23-1-2019
Hydrol. Earth Syst. Sci., 23, 1–17, 2019
https://doi.org/10.5194/hess-23-1-2019

Research article 02 Jan 2019

Research article | 02 Jan 2019

# A simple model for local-scale sensible and latent heat advection contributions to snowmelt

A simple model for local-scale sensible and latent heat advection contributions to snowmelt
Phillip Harder1, John W. Pomeroy1, and Warren D. Helgason1,2 Phillip Harder et al.

Abstract

1 Introduction

Sensible and latent turbulent heat fluxes contributing to snowmelt are complicated during snow-covered area (SCA) depletion by the lateral redistribution of energy from snow-free surfaces to snow. Unfortunately, many calculations of the snow surface energy balance have largely been limited to one-dimensional model frameworks (Brun et al., 1989; Gray and Landine, 1988; Jordan, 1991; Lehning et al., 1999; Marks et al., 1999) which simulate melt at points without considering variations in SCA. Despite the sophistication of these methods, they have neglected local-scale advection of energy.

To model advection to snow patches, Weisman (1977) applied mixing length theory, implicitly accounting for both latent heat advection (LEA) and sensible heat advection (HA). This work was proposed at a time when knowledge of the statistical properties of snow cover were inadequate for easy extension to natural snowpacks. Gray et al. (1986) noted that

The major obstacle to the development of an energy balance model for calculating melt quantities is the lack of reliable methods for evaluating the sensible heat flux. A priority research need is the development of “bulk methodologies” for calculating this term, especially for patchy, snow-cover conditions.

There remains a pressing need for an approach that can estimate areal average HA and LEA contributions during snowmelt that can easily integrate with existing one-dimensional snowmelt models. This work seeks to understand the implications of including local-scale HA and LEA with one-dimensional snowmelt models. To address this objective, this paper presents a simple and easily implementable HA and LEA model. The specific objectives are to validate the proposed model with observations of advection, to re-evaluate the scaling relationships of snow-cover geometry with current datasets of snow cover and to quantify the implications of including advection upon snowmelt.

2 Methodology

## 2.1 Model framework

Horizontal gradients of scalar properties are a first-order control on the advection flux. For snowmelt the gradients are conceptualized as snow-free surfaces upwind of a transition to a snow-covered surface. During melt periods upwind snow-free surfaces are typically comprised of dry soil and/or ponded water which correspond to warm dry and/or warm moist near-surface air properties, respectively. In contrast snow is commonly assumed to be ≤0C with saturated near-surface air (Fig. 1a). Conceptual air temperature and specific humidity profiles over snow, soil, and water surfaces are shown in Fig. 1b to articulate the atmospheric boundary layer dynamics observed by Granger et al. (2002, 2006) and Harder et al. (2017). Assuming the changes in profiles are solely due to exchange with the surface the magnitude and direction of the energy flux can be quantified by the integrated differences in profiles between the surface and the mixing height; the point above the surface where differences due to surface heterogeneity disappear with atmospheric mixing (Granger et al., 2002). When the upwind snow-free surface is warm the cooling of the air as it moves over the snow will lead to sensible heat advection to the snowpack and vice versa. Latent heat advection is dependent upon surface temperature as well as saturation. Thus, air from a dry soil may increase in humidity as it moves over snow, which induces greater sublimation and therefore a reduction in snowmelt energy (Harder et al., 2017). In contrast, a wet upwind condition will lead to a decrease in humidity as the air moves over the relatively drier snow due to condensation upon the snow surface, which imparts a release of latent heat or an increase in snowmelt energy (Harder et al., 2017).

Figure 1(a) Conceptual cross section of the advection process during snowmelt and (b) conceptual specific humidity and air temperature profiles between snow (0 C, 100 % RH), soil (6 C, 60 % RH) and water (1 C, 100 % RH) surfaces and the mixing height (3 C, RH of 60 %).

$\begin{array}{}\text{(1)}& {f}_{\mathrm{s}}=\begin{array}{ll}\mathrm{1}& \mathrm{SCA}>\mathrm{0.6},\\ \left(\frac{\mathrm{SCA}-\mathrm{0.4}}{\mathrm{0.2}}\right)& \mathrm{0.4}\le \mathrm{SCA}\le \mathrm{0.6},\\ \mathrm{0}& \mathrm{SCA}<\mathrm{0.4}.\end{array}\end{array}$

An fs of 1 implies the exchange of advection energy is limited by the snow-free patch areas and an fs of 0 implies the exchange of advection energy is limited by the snow patch areas. Conceptually early advection from snow-free patches will have a more effective energy exchange mechanism than later advection to isolated snow patches. The unstable temperature profile above a relatively rough warm snow-free surface patch will enhance exchange with the atmosphere, and therefore surrounding snow cover, per unit area of snow-free surface. In contrast, the stable temperature profiles above a cool and smooth isolated snow patch will limit energy exchange per unit area of snow surface. The stability influences upon surface exchange dynamics are implicitly accounted for in the parameterization of stability terms by Weisman (1977) and are expressed in Sect. 2.1.1.

Figure 2Conceptual model of advection dynamics for (a) the early melt period where energy is limited to what is transported out of soil (brown) patches to the surrounding snow (white), and for (b) the later melt period where snow patches remain and advection energy is limited to that exchanged over the discrete patches.

During snowmelt, meltwater may infiltrate into the frozen soil and any excess will pond prior to and during the runoff phase; these interactions will influence the near-surface humidity of the snow-free surface. Thus, LEA may enhance sublimation when the upwind surface is dry or condense and enhance melt when the upwind surface is wet (Harder et al., 2017). Any attempt to model advection must quantify the dynamic spatial properties of the snow and snow-free patch distributions, SCA, fractional water coverage of ponded water, and horizontal gradients of temperature and humidity between snow and snow-free surfaces. With quantification of these processes, existing simple advection parametrizations can be extended to calculate HA and LEA contributions to snowmelt in a manner that accounts for the dynamics of the driving variables and processes and still be easily implemented in snowmelt energy balance models. The SLHAM model quantifies the components of the conceptual model outlined in Figs. 1 and 2.

Table 1Parameterizations for extended GM2002.

esc: snow surface vapour pressure (kPa); Tsc: snow surface temperature (K); esf: snow-free surface vapour pressure (kPa); Tsf: snow-free surface temperature(K); g: acceleration due to gravity (9.81 m s−2); u: wind speed (m s−1); κ: von Kármán constant (0.4); u: friction velocity (m s−1); qsc: snow surface specific humidity (kg kg−1); z0s: snow surface roughness (0.005 m); qsf: snow-free surface specific humidity (kg kg−1); γ: psychrometric constant (kPa K−1).

### 2.1.1 Advection versus distance from surface transition

Granger et al. (2002) developed a simplified approach to estimate the advection over a surface transition from boundary layer integration. Advected energy, QA (W m−2), was presented as a power function of patch length, L (m), downwind of a surface transition as

$\begin{array}{}\text{(2)}& {Q}_{\mathrm{A}}\left(L\right)=a{L}^{b}.\end{array}$

The coefficient a (W m−2) scales with wind speed and the horizontal scalar gradient and the coefficient b (–) is a function of the Weisman (1977) stability parameters (W). Parametrizations for these coefficients vary for sensible (HA) and latent (LEA) heat advection and whether advection is from a snow-free patch or to a snow patch; parametrizations are summarized in Table 1. The GM2002 approach is restricted to considering HA contributions to snow. To extend this approach to LEA the a and b parameterizations of GM2002 were assumed to remain valid. The parameterization for coefficient a in the case of LEA was modified to use the surface vapour pressure gradient (kPa) with division by the psychrometric constant (γ, kPa K−1). This relates the horizontal water vapour gradient to be in terms of an equivalent temperature gradient; in the units of the original a parametrization. The coefficient b for LEA uses the humidity stability parameter of Weisman (1977) rather than the temperature stability parameter.

The humidity of the air at the surface interface is rarely observed but is needed to quantify the LEA term. The esc was estimated by assuming saturation at the Tsc. The esf is more challenging as it varies with the surface fraction of ponded water (Fwater, no unit) as

$\begin{array}{}\text{(3)}& {e}_{\mathrm{sf}}={F}_{\mathrm{water}}{e}_{\mathrm{wat}}+\left(\mathrm{1}-{F}_{\mathrm{water}}\right){e}_{\mathrm{soil}}.\end{array}$

The surface water vapour pressure for water surfaces (ewat, kPa) was estimated by assuming saturation at the surface temperature of the ponded water (Twat, K). Assuming negligible evaporation from dry soil surfaces during snowmelt, the surface water vapour of soil (esoil, kPa) can be taken to be the same as actual vapour pressure observed above the surface. The Tsf was also weighted by Fwater as

$\begin{array}{}\text{(4)}& {T}_{\mathrm{sf}}={F}_{\mathrm{water}}{T}_{\mathrm{wat}}+\left(\mathrm{1}-{F}_{\mathrm{water}}\right){T}_{\mathrm{soil}},\end{array}$

where Tsoil (K) is the dry soil surface temperature. The remaining uncertainties in applying this framework are the representation of the statistical distribution of L, and estimation of Fwater and SCA.

### 2.1.2 Fractional coverage of ponded water

To estimate Fwater, the meltwater in excess of frozen soil infiltration capacity was estimated using the parametric frozen soil infiltration equation of Gray et al. (2001). Gray et al. (2001) parameterized the maximum infiltration of the limited condition (INF, mm) as

$\begin{array}{}\text{(5)}& \mathrm{INF}=C{S}_{\mathrm{0}}^{\mathrm{2.92}}{\left(\mathrm{1}-{S}_{i}\right)}^{\mathrm{1.64}}{\left(\frac{\mathrm{273.15}-{T}_{\mathrm{si}}}{\mathrm{273.15}}\right)}^{-\mathrm{0.45}}{t}_{\mathrm{0}}^{\mathrm{0.44}},\end{array}$

where C (2.1, no unit) is a coefficient representing prairie soils, S0 (–) is a surface saturation (generally assumed to be 1), Si (–) is the antecedent soil saturation, Tsi (K) is the initial soil temperature, and t0 (h) is the infiltration opportunity time. The t0 term is estimated as the cumulative hours of active snowmelt over the course of the snowmelt period. Excess meltwater (Mexcess, mm) is calculated as

$\begin{array}{}\text{(6)}& {M}_{\mathrm{excess}}=\sum _{t=\mathrm{0}}^{i}{M}_{t}-{\mathrm{INF}}_{i},\end{array}$

where M (mm) is the snowmelt since the beginning of melt (t=0) to the present time step i.

To relate Mexcess to Fwater, an elevation profile of the microtopography must be known. For simplicity, the furrows that define the microtopography of an agricultural field were assumed to be represented by a half-period, trough to peak, of a sine curve (Fig. 3). Thus Fwater is given by the solution of

$\begin{array}{}\text{(7)}& {S}_{\mathrm{ret}}=\frac{\mathrm{1}}{\mathit{\pi }}\mathrm{sin}\left({F}_{\mathrm{water}}\mathit{\pi }\right)-{F}_{\mathrm{water}}\mathrm{cos}\left({F}_{\mathrm{water}}\mathit{\pi }\right),\end{array}$

where the ratio of filled detention storage (Sret, no unit) is determined from

$\begin{array}{}\text{(8)}& {S}_{\mathrm{ret}}=\frac{{M}_{\mathrm{excess}}}{{S}_{max}},\end{array}$

where a user-defined Smax (mm) is the maximum detention storage of the surface. Any Mexcess that is greater than Smax is removed as runoff and thereafter unavailable to future infiltration.

Figure 3Conceptual water area–volume relationship diagram where a cross section of land surface microtopography (brown is soil and blue is water) is assumed to follow a sinusoidal profile.

### 2.1.3 Snow-covered area

The SCA constrains the overall exchange of energy between the snow surface and the atmosphere. Essery and Pomeroy (2004) developed an SCA parameterization from the closed form fit to the parametric SCA curve produced by homogeneous melt of a log-normal SWE distribution,

$\begin{array}{}\text{(9)}& \mathrm{SCA}=\mathrm{tanh}\left(\mathrm{1.26}\frac{\mathrm{SWE}}{{\mathit{\sigma }}_{\mathrm{0}}}\right),\end{array}$

where SWE is in mm and σ0 (mm) is the standard deviation of SWE at the pre-melt maximum accumulation. The σ0 constrains the spatial variability of a snowpack and how it relates to SCA depletion. Snow cover with high spatial variability will have a longer duration of patchiness and therefore advection will contribution to more of the total snowmelt. Other parameterizations of SCA exist and this was selected for its simplicity, relative success in describing observed SCA curves, and derivation in similar environments with regard to what is being modelled.

### 2.1.4 Snow geometry

Perimeter–area relationships and patch area distributions of snow and snow-free patches show fractal characteristics that can be exploited to simplify the representation of snow-cover geometry needed to calculate advection. There are two commonly used scaling relationships. From application of Korcak's law by Shook et al. (1993a) the fraction of snow patches greater than a given area, F(Ap), is given as a power law distribution

$\begin{array}{}\text{(10)}& F\left({A}_{\mathrm{p}}\right)={\frac{{A}_{\mathrm{p}}}{{c}_{\mathrm{1}}}}^{-{D}_{k}/\mathrm{2}},\end{array}$

where c1 is a threshold value (given as the smallest patch size observed, and hereafter taken as 1 m2), Ap (m2) is patch area, and Dk (–) is the scaling dimension. The scaling dimension is the same between snow and snow-free patches, relatively invariant with time, and ranges between 1.2 and 1.6 (Shook et al., 1993b) and is not a fractal dimension (Imre and Novotn, 2016). A Hack's law relationship between linear dimension and area of landscape features was established by Rigon et al. (1996) and this was extended to Ap and L of snow patches by Granger et al. (2002) as

$\begin{array}{}\text{(11)}& L={c}_{\mathrm{2}}\cdot {A}_{\mathrm{p}}^{\frac{{D}^{\prime }}{\mathrm{2}}},\end{array}$

where c2 is a constant taken as 1 and D was fitted by Granger et al. (2002) to be 1.25.

The relationships of Eqs. (10) and (11) were exploited to develop a probability distribution of L. The exceedance fraction, Eq. (10), was converted to a probability distribution with calculation of probabilities for discrete intervals; this also entailed appropriate selection of intervals. The patch area probability (p(Ap)) is also equivalent to the probability associated with the probability of patch length (p(L)), therefore

$\begin{array}{}\text{(12)}& p\left({L}_{i}\right)=p\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right)=F\left({A}_{{\mathrm{p}}_{\mathrm{i}}-\mathrm{1}}\right)-F\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right),\end{array}$

where i is the index for intervals of Ap that span a range constrained as ${c}_{\mathrm{1}}\le {A}_{\mathrm{p}}<\mathrm{\infty }$. A discrete bin width of ≤1 m is advised to capture the large change in F(Ap) at the more frequent small values of Ap. To estimate an areal average advection exchange the normalized areal extent of each patch size was calculated. The limited number of the largest patches will dominate the exchange surface extent. Thus $p\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right)$ is transformed to give a normalized areal fraction of the unit area that is represented by each patch size $f\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right)$ as

$\begin{array}{}\text{(13)}& f\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right)=\frac{p\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right){A}_{{\mathrm{p}}_{\mathrm{i}}}}{\sum p\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right){A}_{{\mathrm{p}}_{\mathrm{i}}}}.\end{array}$

The transformation of the probability of occurrence to a fractional area of patch size is visualized in Fig. 4.

Figure 4Probability of patch size occurrence and its transformation to fractional area patch sizes for a range in patch sizes from 1 to 1000 m2.

Using the above-described parameterizations of $f\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right)$, L, SCA, Fwater and INF, as well as boundary layer integration HA and LEA parameterizations, the areal average advection, $\stackrel{\mathrm{‾}}{{Q}_{\mathrm{A}}}$ (W), can be calculated as

$\begin{array}{ll}\stackrel{\mathrm{‾}}{{Q}_{\mathrm{A}}}=& {f}_{\mathrm{s}}\left(\mathrm{1}-\mathrm{SCA}\right)\sum _{i=\mathrm{1}}^{i={A}_{max}}f\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right){H}_{\mathrm{A},\mathrm{sf}}\\ & +\left(\mathrm{1}-{f}_{\mathrm{s}}\right)\mathrm{SCA}\sum _{i=\mathrm{1}}^{i={A}_{max}}f\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right){H}_{\mathrm{A},\mathrm{sc}}\\ & +{f}_{\mathrm{s}}\left(\mathrm{1}-\mathrm{SCA}\right)\sum _{i=\mathrm{1}}^{i={A}_{max}}f\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right){\mathrm{LE}}_{\mathrm{A},\mathrm{sf}}\\ \text{(14)}& & +\left(\mathrm{1}-{f}_{\mathrm{s}}\right)\mathrm{SCA}\sum _{i=\mathrm{1}}^{i={A}_{max}}f\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right){\mathrm{LE}}_{\mathrm{A},\mathrm{sc}}.\end{array}$

The terms from left to right represent the HA from snow-free patches, HA to snow patches, LEA from snow-free patches, and LEA to snow patches. All summation terms constitute HA and LEA for the range of patch areas expected, from 1 m2 to an environment appropriate maximum expected patch size (Amax, m2). Calculation of HA and LEA use Eq. (2) with application of appropriate a and b parameterizations from Table 1 and L as calculated with Eq. (11) from the range of Ap. Advection fluxes for the range of patch sizes encountered are weighted by $f\left({A}_{{\mathrm{p}}_{\mathrm{i}}}\right)$, Eq. (13), to give an areal average maximum flux. The advection process must be constrained to snow-free or snow surfaces over which exchange takes place, hence the scaling of the maximum advection by (1−SCA) and SCA from snow-free patches and to snow patches respectively. The fs and (1−fs) terms quantify the relative contribution from snow-free patches and to snow patches over snowmelt and SCA depletion. The primary controls on the model behaviour are the horizontal gradients of humidity and temperature, as well as wind speed.

## 2.2 Re-evaluation of snow-geometry scaling relationships

The coefficients for the snow-cover geometry relationships are based on oblique terrestrial photography or aerial photography with coarse resolution and limited temporal sampling (Shook et al., 1993b). Recent advances in UAV technologies provide a tool to re-evaluate these relationships with georectified high-resolution imagery. During the 2015 and 2016 snowmelt seasons, 0.035 m × 0.035 m spatial resolution red–green–blue (RGB) imagery was collected daily during active melt. This imagery was classified into snow and non-snow areas with pixel-based supervised thresholding of blue band reflectance. Cells that share the same classification and were connected via any of the four mutually adjacent cell boundaries were grouped into snow and non-snow patches. The SDMTools R package (VanDerWal et al., 2014) was used to calculate patch areas. Patch length is a challenging to define and quantify. For this analysis a similar approach to Granger et al. (2002) was used in which the patch length was calculated as the mean of the height and width of the minimum rotated bounding box that contained the entire snow patch. Patches with areas less than 1 m2 were removed from the analysis as noise and classification artefacts are associated with such small patch sizes. The 1 m2 area threshold is consistent with the existing literature on advection and snow-cover geometry (Granger et al., 2002; Shook et al., 1993b, a). When SCA was less than 50 %, snow patch metrics were quantified, and when SCA was greater than 50 %, snow-free patch metrics were quantified. An example is provided in Fig. 5.

## 2.3 Model dynamics

The influence of the advection model upon snowmelt dynamics was explored with two approaches. The first approach is a scenario and sensitivity analysis where inputs are fixed and a selection of process parameterizations are employed to illustrate the relationship between HA and LEA and the snow-free surface humidity dynamics and snowmelt implications. The second approach coupled the SLHAM with an existing one-dimensional snowmelt model to estimate the influence of including, or not including, the advection process on snowmelt simulations.

Figure 5Example of snow-cover geometry scaling properties, exceedance faction versus patch area (b, Dk=21) and patch length versus patch area (c, ${D}^{\prime }=\mathrm{122}$), for snow-covered area classification at 1 m resolution from 29 March 2016 (a, axes are UTM 13N northing and eastings). Red lines are the best-fit scaling relationships, where slope provides the scaling constant.

### 2.3.1 Scenario analysis

To explore the dynamics of modelled advection contributions several scenarios were implemented with the model. The first scenario (No Advection) constitutes a baseline for a typical one-dimensional model that assumes no advection, the second (dry surface) includes advection from a warm dry surface, the third (wet surface) includes advection from a warm wet surface, and the fourth (dry to wet surface) includes advection from a warm surface that transitions from dry to wet as a function of the INF–SretFwater relationships. To understand the implications upon snowmelt for each scenario, input variables were held constant and the model was run until an assumed isothermal snowpack was fully depleted. A constant melt energy, Qnet (W m−2), was applied which represents the net snow surface energy balance as estimated via a typical one-dimensional model. The initialized SWE was ablated, leading to infiltration excess, detention-storage, runoff, or sublimation. The relative dynamics of the various scenarios are sensitive to the inputs and parameters used, as summarized in Table 2, and demonstrate the relationships between HA and LEA and the snow-free surface humidity conceptualization and snowmelt implications from a theoretical perspective.

### 2.3.2 Sensitivity analysis

The sensitivity of SLHAM to Twat, Ta, Tsoil, u, and RH variability is also explored to understand the implications upon SWE and SCA depletion, Fwater, HA, LEA, and net advection. The dry to wet surface scenario, holding the input variables constant and varying each variable in turn as detailed in Table 2, was employed to understand the dynamics of input variability. A common assumption is that Twat is 0 C as meltwater immediately after discharge from an isothermal snowpack is 0 C and underlying frozen soils are ≤0C. Unlike the snow surface the maximum temperature of ponded water is unconstrained by phase change so values ≥0C are expected because of possible low water surface albedos and high shortwave irradiance (SW${}_{\mathrm{atm}}^{↓}$) during the daytime. Analysis of available thermal images from a FLIR T650 thermal camera was used to correct for atmosphere conditions and water surface emissivity. This analysis showed that daytime Twat was generally >0 and <2C. This range in Twat was used to test the sensitivity of the Twat upon SLHAM dynamics. Intermittency of observations and inherent uncertainties in thermography prevented a more precise estimation of Twat. The ranges of all other variables were selected to represent conditions commonly experienced during snowmelt on the Canadian Prairies.

Table 2Input variables for scenario analysis of SLHAM dynamics.

### 2.3.3 Coupled advection and snow–stubble–atmosphere snowmelt model simulations

Conditions controlling advection processes are not constant over snowmelt; therefore SLHAM was coupled with a one-dimensional snowmelt model (SSAM) to estimate the role of advection contributions over a snowmelt season. Briefly, SSAM describes the relationships between shortwave, longwave and turbulent exchanges between a snow surface underlying exposed crop stubble and the atmosphere. The surface energy balance was coupled to a single-layer snow model to estimate snowmelt. A slight modification of SSAM, or any one-dimensional model that computes areal average snowmelt, is needed to include advection. The energy terms of one-dimensional energy balance models are represented as flux densities (W m−2) over an assumed continuous snow cover and therefore need to be weighted by an SCA parametrization (Eq. 9) to properly simulate the areal average melt energy available to the fraction of the surface comprised of snow. The SSAM was run with and without SLHAM to explore the impact of advection simulation on SWE. Simulation performance was quantified via root mean square error (RMSE) and model bias (MB) of the simulated SWE versus snow survey SWE observations. The relative contribution of advection was quantified through estimation of the energy contribution to total snowmelt. A commonly used snowmelt model, the Energy Balance Snowmelt Model (EBSM) of Gray and Landine (1988), was also run to benchmark performance. The EBSM has been widely applied in this region and simulation is deployed as an option within the Cold Region Hydrological Modelling (CRHM) platform (Pomeroy et al., 2007).

The SSAM, SSAM–SLHAM and EBSM simulations were driven by common observed meteorological data, parameters and initial conditions obtained from intensive field campaigns at a research site near Rosthern, Saskatchewan, Canada (52.69 N, 106.45 W). The data for the 2015 and 2016 snowmelt seasons reflect relatively flat agricultural fields characterized by standing wheat stubble, with 15 and 24 cm stubble heights for the respective years. Observations of Tsoil required for SLHAM come from infrared radiometers (Apogee SI-111) deployed on mobile tripods to snow-free patches. Unfortunately, no time series of Twat observations are available and values or models to describe Twat for shallow ponded meltwater in a prairie environment have not been discussed in the literature. Like snowpack refreezing, ponded meltwater can also refreeze at night as the heat capacity of this shallow water is limited. In this framework, as observations or models of Twat are unavailable, a simple physically guided representation of Twat takes the following form:

$\begin{array}{}\text{(15)}& {T}_{\mathrm{wat}}=\begin{array}{ll}{T}_{\mathrm{sc}}& {T}_{\mathrm{sc}}<\mathrm{0}\phantom{\rule{0.125em}{0ex}}{}^{\circ }\mathrm{C},\\ \mathrm{0.5}\phantom{\rule{0.125em}{0ex}}{}^{\circ }\mathrm{C}& {T}_{\mathrm{sc}}=\mathrm{0}\phantom{\rule{0.125em}{0ex}}{}^{\circ }\mathrm{C}.\end{array}\end{array}$

A description of the field site and data collection methodologies is detailed in Harder et al. (2018).

3 Results and discussion

## 3.1 Performance of extended GM2002

The extended GM2002 proposed here was tested using advection estimates from vertical air temperature and water vapour profiles as reported in Harder et al. (2017); the results are summarized in Table 3. The model slightly overestimated HA and LEA on 30 March 2015, likely due to the limiting assumptions of the GM2002 model. A key missing component of GM2002 is the influence of differences in surface roughness upon the growth of the internal boundary layer. A simple power law relationship with respect to distance from transition is employed in the model. Further work by Granger et al. (2006) demonstrated that boundary layer growth has a positive relationship with upwind surface roughness and that the parametrization employed in GM2002 overestimates the boundary layer depth, by up to a factor of 2 when upwind surface roughness is negligible. The GM2002 is based upon the integrated difference in temperature through the boundary layer depth; thus, a greater boundary layer depth will increase the estimated advection. This partly explains why the model overestimates values in the situation of a rough upwind surface. Other potential limiting assumptions include homogenous surface temperatures, uniform eddy diffusivities for different scalars, and no vertical advection. Despite the model limitations, the acceptable performance in simulating the 18 and 30 March observations gives confidence that this simple model is reasonable for some applications and provides guidance for future improvements.

Table 3Model parameters, estimates and observations for evaluation of the extended GM2002.

a Estimated from thermography. b Roughly estimated from application of a 1 : 100 sensor height to flux footprint ratio (Hsieh et al., 2000) as applied to concurrent UAV imagery.

## 3.2 Re-evaluation of snow-cover geometry

Differences exist between the originally reported parameters and those found from the analysis of UAV imagery (mean coefficients summarized in Table 4). Early work applying fractal geometry to natural phenomena (Mandelbrot, 1975, 1982) discusses the Korcak exponent as a fractal dimension. More recent work suggests that the Korcak law describing the area–frequency relationship is not a fractal relationship but rather a mathematically similar, but distinct, scaling law (Imre and Novotn, 2016). Therefore, the Dk value is not necessarily ≥1 or ≤2 and the identified exponent terms in Table 4 near or greater than 2 are plausible. The D terms are very similar to those previously reported (Granger et al., 2002). From this analysis, it is apparent that application of these parameters between sites must be done with caution as local topography and surface conditions may influence the snow patch size distribution. The lack of a temporal trend of these terms (time series of Dk in Fig. 6 and D in Fig. 7) over the course of snowmelt and equivalence in scaling of snow and snow-free patches implies that locally specific parameters may be applied as constants over the course of the melt and irrespective of patch type. The resolution of the underlying imagery, differences in classification methodologies and surface characteristics may contribute to some of the differences in terms observed and those previously reported. An illustrative comparison is that of a tall and short stubble surface. The tall stubble surface snow-cover geometry is heavily influenced by the early exposure (and hence classification as non-snow from nadir imagery) of stubble rows, which leads to very long and narrow patches even if snow is still present within the stubble. In contrast the oblique imagery of Shook et al. (1993b) and Granger et al. (2002) will not quantify the snow between stubble rows and larger and less complex snow patches would be represented by the previously reported coefficients. Further work is needed to calculate the scaling properties of patches over a more comprehensive variety of topography and vegetation types.

Table 4Updated mean snow-cover geometry parameters.

Figure 6Time series of fitted Dk parameter with respect to snow and soil patches for various land covers over the course of snowmelt.

Figure 7Time series of fitted D parameter with respect to snow and soil patches for various land covers over the course of snowmelt.

## 3.3 Implications of including advection in snowmelt models

### 3.3.1 Advection dynamics in scenario simulations

Figure 8Modelled snow water equivalent depletion for various advection scenarios.

Figure 9Latent heat (green), sensible heat (red) and net (blue) advection components for the SLHAM scenarios plotted with snow-covered area (black).

Figure 10Sensitivity of snow water equivalent and snow-covered area depletion, ponded water fraction, sensible heat advection, latent heat advection and net advection with respect to variation in water surface temperature.

### 3.3.2 Sensitivity to input variables

The influence of the input variables on the SLHAM model is evaluated through a sensitivity analysis (Fig. 10). It is apparent from the variability in SWE depletion that the Tsoil and u have the largest influence on advection contributions to snowmelt. This is expected as u and Tsoil variables quantify the first-order controls driving advection, the air mass movement and horizontal scalar gradients respectively. In contrast the Twat, Ta, and RH variables have considerably less variability for the ranges simulated as they have less influence upon the scalar profile differences between upwind and downwind locations. A critical model feedback relates to the influence dynamic upwind surface temperature and humidity and is articulated in this sensitivity analysis. If melt rates exceed the frozen soil infiltration capacity ponding occurs, Fwater>0, which forces the upwind surface to the assumed water surface temperature. The consequent sign of the surface humidity gradient will influence whether LEA induces condensation (increased melt rate) or sublimation (decreased melt rate), which influences the net advection and melt rate. This feedback is manifested in the sensitivity of all variables. The transition of the upwind surface from dry and warm to cooler and saturated tempers the advection contributions to melt. Generally, any change in a variable that increases the profile gradient or increases energy exchange will lead to increased SWE and SCA depletion rates and increased extent and duration of Fwater. Changes in HA and LEA tend to be compensatory, resulting in relatively small increases in net advection fluxes.

Sensitivity to any variable is only expressed towards the end of the snowmelt, when SWE <50 mm and SCA is depleting rapidly. Differences in melt rate are limited by the rapid reduction in the SCA exchange surface at the end of snowmelt. Whilst clearly important for simulating the dynamics of advection and sources of energy driving snowmelt, Twat, Ta, and RH have a relatively limited influence upon overall SWE depletion compared to Tsoil and u. In the absence of Twat models or observations, the assumptions outlined in Eq. (15) will have a relatively limited influence upon simulation of SWE with the fully coupled SSAM–SLHAM model.

The scenario analysis demonstrates the melt response to variations in surface wetness but actual snowmelt situations have forcings that vary diurnally and with meteorological conditions. Snowmelt simulations with three models of varying complexity provide insight into the implications of process representation. SSAM and SSAM–SLHAM show considerable improvement when compared to EBSM (Fig. 11 and Table 5). The SSAM simulation is by itself a significant improvement upon EBSM for SWE prediction during melt. The addition of SLHAM does not change the SWE simulation performance appreciably but does increase the physical realism of the model with its more complete surface energy balance. The SSAM–SLHAM simulations including advection, relative to SSAM simulations without advection, led to lower areal average melt rates in 2015 and higher rates in 2016. Lower wind speeds in 2015 led to lower advection contributions than 2016, which had relatively higher wind speeds. The comparison of the simulated melt with snow survey SWE observations showed that the differences are minimal (Fig. 11 and Table 5). While the SSAM–SLHAM simulations do not change melt rates or total amount of energy, the sources of energy driving snowmelt do change. Early melt displays no differences as SCA remains relatively homogenous. Differences appear due to decreases in the turbulent radiation fluxes, with a decrease in the SCA exchange surface due to increases in advection fluxes with increasing horizontal scalar gradients and surface heterogeneity. The cumulative net energy from advection for these two seasons contributed energy to melt 4 and 5 mm of SWE in 2015 and 2016, respectively (Fig. 12). The advection energy contribution represents 6.5 % and 10.6 % of total snowmelt in 2015 and 2016, respectively.

Figure 11Snow water equivalent simulation for EBSM (red line), SSAM (green line) and SSAM–SLHAM (blue line) with respect to snow survey mean (black points) and 95 % percentile sampling confidence interval (black lines).

Table 5Error metrics of snow water equivalent simulation versus snow survey observations for EBSM, SSAM and SSAM–SLHAM models.

## 3.4 Energy balance compensation

An unappreciated dynamic of local-scale advection during snowmelt is that LEA and HA may be of opposite sign and therefore will compensate for one another, leading to a lower net advection contribution. This occurs when the gradients of T and q between a snow-free and snow-covered surface are opposite in sign; a warm but dry snow-free surface upwind of a cool and wet snow-covered surface driving snow surface sublimation. This was evident in the reduction of the advection energy due to a negative LEA throughout the dry surface scenario and early melt of the dry to wet surface scenario (Fig. 9). In the 2015 and 2016 snowmelt simulations, the accumulated LEA was negative for much of the melt period, which compensated for the consistently positive HA term (Fig. 12). LEA only increased, enhancing the positive HA contribution, near the end of melt in 2015 when increased surface wetness led to a positive LEA term.

Figure 12Cumulative sensible (red), latent (green) and net (blue) advection terms in terms of energy (MJ: left) and equivalent melted snow water equivalent (mm SWE: right axis).

Table 6Cumulative energy from sensible, latent and net exchange for 2015 and 2016 snowmelt simulations with (SSAM–SLHAM) and without (SSAM) advection.

The simulation of snowmelt with, and without, advection gave minimal differences in the resulting SWE simulation. This demonstrates system insensitivity to processes that on their own appear to be important. This may explain why EBSM, like many other physically based snowmelt models (Jordan, 1991; Lehning et al., 1999; Marks et al., 1998), does not accommodate heterogeneous snow cover yet successfully simulates SWE depletion. In EBSM the simulation of an areal average albedo rather than a snow albedo performed relatively well in simulating SWE (Fig. 11) without considering SCA depletion or advection controls. The modelling challenges of asnow are not limited to EBSM as other asnow parameterizations, especially temperature-dependent ones, typically underestimate asnow during melt and therefore indirectly, and perhaps unintentionally, account for advected energy contributions (Pedersen and Winther, 2005; Raleigh et al., 2016). While modelled asnow values that underestimate actual asnow values are effective parameterizations for simulation of SWE, they cannot realistically incorporate the impacts of dust on snow or changes in snow albedo with grain size or wetness. Hence, SCA constraints and advection process conceptualizations are necessary to improve confidence in and applicability of snowmelt models. This is evident when comparing the more accurate and physically complete SSAM–SLHAM simulation of SWE to the EBSM simulation of SWE (Fig. 11).

Understanding the implications of land-use and climate changes on variables beyond SWE are needed to fully inform coupled modelling of land–atmosphere and radiation feedbacks between land surface and numerical weather or climate models. The framework presented explicitly considers advection and scales it with SCA, u and horizontal gradients, which are the primary controls of advection. A simple indication that a more appropriate model conceptualization is being used in this advection framework is that the minimum albedo value simulated is 0.75. This is consistent with that for clean, melting snow (Wiscombe and Warren, 1980), whilst the 0.2 in EBSM is not. While the SWE simulation differences are not particularly large, the new model is getting the “right” answer for the “right” reasons and without calibration. By including a more appropriate suite of physical processes, this model can produce realistic melt simulations in areas or years where the variables governing advection deviate from the conditions observed during model development.

## 3.6 Limitations and future research needs

The SLHAM framework replaces the large uncertainty deriving from physically unrealistic albedo parametrizations (Gray and Landine, 1987; Raleigh et al., 2016) and ignored SCA dynamics (Essery and Pomeroy, 2004) with a more physically realistic framework. The individual process parametrizations still have uncertainties that need to be constrained. The advection versus patch length parametrization of GM2002 lacks inclusion of surface roughness differences and the valid bounds of the parametrizations need clarification. Observations of stable atmospheric profiles over snow patches (Fujita et al., 2010; Mott et al., 2015, 2016; Shook and Gray, 1997) complicate energy exchange. The goal of this simple model was to develop an easy-to-implement advection framework with stability represented by the Weisman (1977) stability parameters. Future work will need to revaluate the stability assumptions of Granger et al. (2002) and Weisman (1977) or devise more appropriate schemes to account for the stability influence. The SCA model of Essery and Pomeroy (2004) is challenged by exposure of vegetation in shallow snow. The conceptual surface water ponding model developed in this work requires field observations or further parameterizations to accurately quantify the relevant variables. The transition of advection mechanism from snow-free sources to snow patch sources uses a conceptualized relationship to SCA. A targeted field campaign is needed to assess the validity of the conceptualized fs and its possible relation to the advection efficiency term of Marsh and Pomeroy (1996). An estimate of Tsf is needed to implement this framework and will limit application of SLHAM in its current form, as modelling Tsf is non-trivial and observations are often unavailable. Ideally a multisource land surface scheme with explicit representation of soils and ponded water is used to represent Tsoil and Twat. In the interim, the Twat assumptions in Eq. (15) may be used but need to be tested further. A regression of Tsoil to incoming shortwave radiation and Ta is presented in the Appendix to provide a simple and physically guided solution to remove this limitation when modelling snowmelt in agricultural regions on the Canadian Prairies. These uncertainties will be addressed in future work and will require additional field observations and model validation, testing, or development.

4 Conclusions

To date, the development of easily implementable and appropriate models to estimate the advection of HA and LEA to snow during melt has proved elusive. The formulation presented here is an initial framework that can be used to augment existing one-dimensional snowmelt models. When tested against observations the extended GM2002 model provides reasonable estimates of both HA and LEA and opportunities for improvement of the method are discussed. The scaling parameters necessary to describe the spatial heterogeneity of snow and snow-free patches were re-evaluated with UAV data. Coupling of the simple advection model with snow-cover geometry scaling laws, SCA depletion, frozen soil infiltration and a surface detention fractional water area parameterization resulted in a model that meets the objective of a formulation that can account for LEA and HA to snow as an areal average contribution. A scenario-based analysis of the model revealed the compensatory influence of LEA from a warm but dry surface; the LEA-driven sublimation offsets HA inputs. Coupling SLHAM with SSAM demonstrated that advection constitutes an important portion of melt energy: 11 % of the melt observed in the 2016 snowmelt season. The reduced radiation exchange to the snow surface fraction, due to decreasing SCA, is compensated for with an increase in net sensible and latent heat exchange that leads to minimal differences in the SWE depletion. This compensatory dynamic has sometimes allowed one-dimensional energy balance snowmelt models to provide adequate simulation of SWE despite using the “wrong” process conceptualizations. The advection model framework proposed here can be easily coupled to existing one-dimensional energy balance models and is expected to improve the prediction of snowmelt in areas dominated by heterogeneous snow cover during melt. Such adoption will permit successful use of more realistic albedo parameterizations. This work provides a guiding framework to address the long-identified need to develop “bulk methodologies” for calculating sensible and latent heat terms for patchy snow-cover conditions (Gray et al., 1986).

Code and data availability
Code and data availability.

The data and code discussed in this paper are available through the corresponding author, Phillip Harder (phillip.harder@usask.ca).

Appendix A

The SLHAM framework requires a Tsoil value which is a challenging variable to explicitly model during snowmelt. To provide an interim solution a multiple linear regression is developed to estimate Tsoil from SW${}_{\mathrm{atm}}^{↓}$ and Ta. This empirical parameterization is appropriate to the snowmelt situation on the Canadian Prairies when the surface is comprised of crop residues and should be treated with caution in other domains. The developed regression is physically guided as the main variables controlling Tsoil is the net radiation, whose variability is dominated by SW${}_{\mathrm{atm}}^{↓}$, and turbulent fluxes, which are dependent upon the Ta gradients. During nighttime Tsoil is very similar to Ta while during daytime the additional energy from SW${}_{\mathrm{atm}}^{↓}$ heats the surface to temperatures above Ta. A multiple regression that contains these parameters provides a simple but effective way to estimate Tsoil in a manner consistent with energy balance interactions. A full description of the observations used to parameterize this relationship can be found in Harder et al. (2018). Briefly the Ta is observed with a shielded Campbell Scientific HMP45C212 and SW${}_{\mathrm{atm}}^{↓}$ is observed with a Campbell Scientific CNR1 with both sensors 2 m above the ground surface. The Tsoil observations from Apogee SI-111 sensors, mounted on mobile tripods to ensure consistent representation snow-free surfaces, sampled surfaces of tall wheat stubble (0.35 m) and short wheat stubble (0.2 m) in 2015 and wheat stubble (0.24 m) and canola stubble (0.24 m) in 2016. Hereafter they are referred to as tall stubble, short stubble, wheat and canola, respectively. All observations were logged at 15 min intervals. The empirical representation of Tsoil (C) in relation to SW${}_{\mathrm{atm}}^{↓}$ (W m−2) and Ta (C) is

$\begin{array}{}\text{(A1)}& {T}_{\mathrm{soil}}=\mathrm{0.00339}{\mathrm{SW}}_{\mathrm{atm}}^{↓}+\mathrm{0.977}{T}_{\mathrm{a}}-\mathrm{1.22}.\end{array}$

Model performance was assessed with the root mean square error (RMSE) and model bias (MB). Each test provides a different perspective on model performance: RMSE is a weighted measure of the difference between the observation and model (Legates and McCabe, 2005), and MB indicates the mean over- or underprediction of the model versus observations (Fang and Pomeroy, 2007). The Tsoil regression provides good estimates of the diurnal variability and magnitudes with respect to observations (Fig. A1). The highest values during daytime are simulated well, which is critical for the appropriate simulation of advection processes. There is low bias for all simulations; MB < 1.09 C. The RMSEs between 1.39 and 1.94 C are negligible as most surface temperature models will simulate errors at a similar magnitude (Aiken et al., 1997). This parametrization provides a simple but effective workaround if Tsoil observations are unavailable or unmodelled. This empirical relation should be treated with caution if implemented outside of the conditions found during snowmelt in cropland areas of the Canadian Prairies. In such cases locally derived relationships should be developed or Tsoil should be explicitly modelled.

Figure A1Soil surface temperature observed versus modelled as scatter plots (a) and time series (b).

Author contributions
Author contributions.

PH collected all field data, conceptualized and coded the SLHAM model, performed simulations and analysis, and wrote the paper. JWP and WDH provided guidance and reviewed and revised model formulations and the paper.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Special issue statement
Special issue statement.

Acknowledgements
Acknowledgements.

Funding was provided by the Natural Sciences and Engineering Research Council of Canada through discovery grants, research tools and Instruments and the Changing Cold Regions Network and the Canada Research Chairs programme. Field and technical assistance from Bruce Johnson, Chris Marsh, Kevin Shook and Michael Schirmer is gratefully acknowledged. This work would not have been possible without the cooperation of Nathan Janzen and Robert Regehr, the farmers who accommodated the intensive field campaigns.

Edited by: Sean Carey
Reviewed by: Rebecca Mott and Richard L. H. Essery

References

Aiken, R. M., Flerchinger, G. N., Farahani, H. J., and Johnsen, K. E.: Energy Balance Simulation for Surface Soil and Residue Temperatures with Incomplete Cover, Agron. J., 89, 404–415, 1997.

Brun, E., Martin, E., Simon, V., Gendre, C., and Coleou, C.: An Energy and Mass Model of Snow Cover Suitable for Operational Avalanche Forecastiong, J. Glaciol., 35, 333–342, 1989.

Essery, R. and Pomeroy, J. W.: Implications of spatial distributions of snow mass and melt rate for snow-cover depletion: theoretical considerations, Ann. Glaciol., 38, 261–265, https://doi.org/10.3189/172756404781815275, 2004.

Essery, R., Granger, R. J., and Pomeroy, J. W.: Boundary-layer growth and advection of heat over snow and soil patches: modelling and parameterization, Hydrol. Process., 20, 953–967, 2006.

Fang, X. and Pomeroy, J. W.: Snowmelt runoff sensitivity analysis to drought on the Canadian prairies, Hydrol. Process., 21, 2594–2609, https://doi.org/10.1002/hyp.6796, 2007.

Fujita, K., Hiyama, K., Iida, H., and Ageta, Y.: Self-regulated fluctuations in the ablation of a snow patch over four decades, Water Resour. Res., 46, 1–9, https://doi.org/10.1029/2009WR008383, 2010.

Garratt, J. R.: The Internal Boundary Layer – A Review, Bound.-Lay. Meteorol., 50, 171–203, 1990.

Granger, R. J. and Male, D. H.: Melting of a prairie snowpack, J. Appl. Meteorol., 17, 1833–1842, 1978.

Granger, R. J., Male, D. H., and Gray, D. M.: Prairie Snowmelt, in: Symposium of the Water Studies Institute, 9–11 May 1978, Saskatoon, SK, Canada, 1978.

Granger, R. J., Pomeroy, J. W., and Parviainen, J.: Boundary-layer integration approach to advection of sensible heat to a patchy snow cover, Hydrol. Process., 16, 3559–3569, 2002.

Granger, R. J., Essery, R., and Pomeroy, J. W.: Boundary-layer growth over snow and soil patches: Field Observations, Hydrol. Process., 20, 943–951, 2006.

Gray, D. M. and Landine, P. G.: Albedo Model for Shallow Prairie Snow Covers, Can. J. Earth Sci., 24, 1760–1768, https://doi.org/10.1139/e87-168, 1987.

Gray, D. M. and Landine, P. G.: An Energy-Budget Snowmelt Model for the Canadian Prairies, Can. J. Earth Sci., 25, 1292–1303, 1988.

Gray, D. M., Pomeroy, J. W., and Granger, R. J.: Prairie Snowmelt Runoff, in: Water Research Themes, Conference Commemorating the Official Opening of the National Hydrology Research Centre, Canadian Water Resources Association, Saskatoon, Saskatchewan, Canada, 49–68, 1986.

Gray, D. M., Toth, B., Zhao, L., Pomeroy, J. W., and Granger, R. J.: Estimating areal snowmelt infiltration into frozen soils, Hydrol. Process., 15, 3095–3111, https://doi.org/10.1002/hyp.320, 2001.

Harder, P., Pomeroy, J. W., and Helgason, W. D.: Local scale advection of sensible and latent heat during snowmelt, Geophys. Res. Lett., 44, 9769–9777, https://doi.org/10.1002/2017GL074394, 2017.

Harder, P., Helgason, W. D., and Pomeroy, J. W.: Modelling the snow-surface energy balance during melt under exposed crop stubble, J. Hydrometeorol., 19, 1191–1214, https://doi.org/10.1175/JHM-D-18-0039.1, 2018.

Hsieh, C. I., Katul, G., and Chi, T. W.: An approximate analytical model for footprint estimation of scalar fluxes in thermally stratified atmospheric flows, Adv. Water Resour., 23, 765–772, https://doi.org/10.1016/S0309-1708(99)00042-1, 2000.

Imre, A. R. and Novotn, J.: Fractals and the Korcak-law: a history and a correction, Eur. Phys. J. H, 41, 69–91, https://doi.org/10.1140/epjh/e2016-60039-8, 2016.

Jordan, R. E.: A one-dimensional temperature model for a snowcover: Technical documentation for SNTHERM.89, Special Rep. 91-16, U.S. Army Cold Regions Research and Engineering Laboratory, Hanover, NH, USA, 49 pp., 1991.

Legates, D. R. and McCabe, G. J.: Evaluating the Use of “Goodness of Fit” Measures in Hydrologic and Hydroclimatic Model Validation, Water Resour. Res., 35, 233–241, https://doi.org/10.1029/1998WR900018, 2005.

Lehning, M., Bartelt, P., Brown, B., Russi, T., Stockli, U., and Zimmerli, M.: SNOWPACK model calculations for avalanche warning based upon a new network of weather and snow stations, Cold Reg. Sci. Technol., 30, 145–157, 1999.

Liston, G. E.: Local advection of momentum, heat and moisture during the melt of patchy snow covers, J. Appl. Meteorol., 34, 1705–1715, 1995.

Mandelbrot, B. B.: Stochastic models for the Earth's relief, shape and fractal dimension of coastlines, and number-area rule for islands, P. Natl. Acad. Sci. USA, 72, 3825–3838, 1975.

Mandelbrot, B. B.: The Fractal Geometry of Nature, Freeman, New York, USA, 1982.

Marks, D., Kimball, J., Tingey, D., and Link, T. E.: The sensitivity of snowmelt processes to climate conditions and forest cover during rain-on-snow: a case study of the 1996 Pacific Northwest Flood, Hydrol. Process., 12, 1569–1587, 1998.

Marks, D., Domingo, J., Susong, D., Link, T. E., and Garen, D.: A spatially distributed energy balance snowmelt model for application in mountain basins, Hydrol. Process., 13, 1935–1959, 1999.

Marsh, P. and Pomeroy, J. W.: Meltwater Fluxes At an Arctic Forest-Tundra Site, Hydrol. Process., 10, 1383–1400, 1996.

Marsh, P., Pomeroy, J. W., and Neumann, N.: Sensible heat flux and local advection over a heterogeneous landscape at an Arctic tundra site during snowtnelt, Ann. Glaciol., 25, 132–136, 1997.

Mott, R., Egli, L., Grünewald, T., Dawes, N., Manes, C., Bavay, M., and Lehning, M.: Micrometeorological processes driving snow ablation in an Alpine catchment, The Cryosphere, 5, 1083–1098, https://doi.org/10.5194/tc-5-1083-2011, 2011.

Mott, R., Gromke, C., Grünewald, T., and Lehning, M.: Relative importance of advective heat transport and boundary layer decoupling in the melt dynamics of a patchy snow cover, Adv. Water Resour., 55, 88–97, https://doi.org/10.1016/j.advwatres.2012.03.001, 2013.

Mott, R., Lehning, M., and Daniels, M.: Atmospheric Flow Development and Associated Changes in Turbulent Sensible Heat Flux over a Patchy Mountain Snow Cover, J. Hydrometeorol., 16, 1315–1340, https://doi.org/10.1175/JHM-D-14-0036.1, 2015.

Mott, R., Paterna, E., Horender, S., Crivelli, P., and Lehning, M.: Wind tunnel experiments: cold-air pooling and atmospheric decoupling above a melting snow patch, The Cryosphere, 10, 445–458, https://doi.org/10.5194/tc-10-445-2016, 2016.

Mott, R., Schlögl, S., Dirks, L., and Lehning, M.: Impact of Extreme Land Surface Heterogeneity on Micrometeorology over Spring Snow Cover, J. Hydrometeorol., 18, 2705–2722, https://doi.org/10.1175/JHM-D-17-0074.1, 2017.

Pedersen, C. A. and Winther, J. G.: Intercomparison and validation of snow albedo parameterization schemes in climate models, Clim. Dynam., 25, 351–362, https://doi.org/10.1007/s00382-005-0037-0, 2005.

Pomeroy, J. W., Gray, D. M., Brown, T., Hedstrom, N. R., Quinton, W. L., Granger, R. J., and Carey, S. K.: The cold regions hydrological model: a platform for basing process representation and model structure on physical evidence, Hydrol. Process., 21, 2650–2667, 2007.

Raleigh, M. S., Livneh, B., Lapo, K., and Lundquist, J. D.: How Does Availability of Meteorological Forcing Data Impact Physically Based Snowpack Simulations?, J. Hydrometeorol., 17, 99–120, https://doi.org/10.1175/JHM-D-14-0235.1, 2016.

Rigon, R., Rodriguez-Iturbe, I., Maritan, A., Giacometti, A., Tarboton, D. G., and Rinaldo, A.: On Hack's Law, Water Resour. Res., 32, 3367–3374, 1996.

Sauter, T. and Galos, S. P.: Effects of local advection on the spatial sensible heat flux variation on a mountain glacier, The Cryosphere, 10, 2887–2905, https://doi.org/10.5194/tc-10-2887-2016, 2016.

Shook, K. and Gray, D. M.: Snowmelt Resulting from Advection, Hydrol. Process., 11, 1725–1736, 1997.

Shook, K., Gray, D. M., and Pomeroy, J. W.: Geometry of patchy snowcovers, in: 50th Annual Eastern Snow Conference, 8–10 June 1993, Quebec City, Canada, 89–98, 1993a.

Shook, K., Gray, D. M., and Pomeroy, J. W.: Temporal Variation in Snowcover Area During Melt in Prairie and Alpine Environments, Nord. Hydrol., 24, 183–198, 1993b.

VanDerWal, J., Falconi, L., Januchowski, S., Storlie, L. S., and Storlie, C.: SDMTools: Species Distribution Modelling Tools: Tools for processing data associated with species distribution modelling exercises, available at: https://cran.r-project.org/package=SDMTools (last access: 20 December 2018), 2014.

Weisman, R. N.: Snowmelt: A Two-Dimensional Turbulent Diffusion Model, Water Resour. Res., 13, 337–342, 1977.

Wiscombe, W. J. and Warren, S. G.: A Model for the Spectral Albedo of Snow. 1: Pure Snow, J. Atmos. Sci., 37, 2712–2733, 1980.