the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

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

### Phillip Harder

### John W. Pomeroy

### Warren D. Helgason

Local-scale advection of energy from warm snow-free surfaces to cold snow-covered surfaces is an important component of the energy balance during snow-cover depletion. Unfortunately, this process is difficult to quantify in one-dimensional snowmelt models. This paper proposes a simple sensible and latent heat advection model for snowmelt situations that can be readily coupled to one-dimensional energy balance snowmelt models. An existing advection parameterization was coupled to a conceptual frozen soil infiltration surface water retention model to estimate the areal average sensible and latent heat advection contributions to snowmelt. The proposed model compared well with observations of latent and sensible heat advection, providing confidence in the process parameterizations and the assumptions applied. Snow-covered area observations from unmanned aerial vehicle imagery were used to update and evaluate the scaling properties of snow patch area distribution and lengths. Model dynamics and snowmelt implications were explored within an idealized modelling experiment, by coupling to a one-dimensional energy balance snowmelt model. Dry, snow-free surfaces were associated with advection of dry air that compensated for positive sensible heat advection fluxes and so limited the net influence of advection on snowmelt. Latent and sensible heat advection fluxes both contributed positive fluxes to snow when snow-free surfaces were wet and enhanced net advection contributions to snowmelt. The increased net advection fluxes from wet surfaces typically develop towards the end of snowmelt and offset decreases in the one-dimensional areal average melt energy that declines with snow-covered area. The new model can be readily incorporated into existing one-dimensional snowmelt hydrology and land surface scheme models and will foster improvements in snowmelt understanding and predictions.

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.

The differences in surface energetics between snow-covered and snow-free areas lead to a heterogeneous distribution of surface temperature and near-surface water vapour. These horizontal gradients drive a lateral exchange of heat (sensible heat advection) and water vapour (latent heat advection when considering the induced condensation or sublimation) over the leading edge of a snow patch. Advection contributions to snowmelt are not negligible as sensible heat advection has been estimated to account for up to 55 % of the snowmelt energy balance (Granger and Male, 1978), resulting in areal melt rates being greatest when snow cover is between 40 % and 60 % (Shook, 1995; Marsh et al., 1997). Advection is very challenging to directly observe due to the dynamic nature of snow-cover ablation. Direct observations of the melt implications of advection have utilized repeat terrestrial laser scanning to identify and quantify larger melt rates on the leading edge of snow patches (Mott et al., 2011). The development of internal boundary layers as air flows over heterogenous surfaces provides an alternate approach to measure the advection energy flux directly (Garratt, 1990). Measurements of these internal boundary layers across snow surface transitions reflect established power laws of boundary layer height (Shook, 1995; Granger et al., 2006) and can quantify the latent and sensible heat advection through boundary layer integration (Granger et al., 2002, 2006; Harder et al., 2017). In contrast to these findings the formation of boundary layers has also been identified as a potential cause of atmospheric decoupling of the atmosphere from the snow surface, leading to the suppression of sensible heat advection in the Alps (Mott et al., 2013, 2016, 2017). The reader is referred to Fujita et al. (2010), Granger et al.(2002, 2006), Mott et al. (2013, 2016, 2017) and Sauter and Galos (2016) for discussions on the complexities of boundary layer development over snow during advection situations and how this may influence overall energy exchange. Due to the complexity of the process and difficulties in observation, modelling has been the focus of much more work on this topic.

To model advection to snow patches, Weisman (1977) applied
mixing length theory, implicitly accounting for both latent heat advection
(LE_{A}) and sensible heat advection (*H*_{A}). 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.

Subsequent models have had variable complexity.
Marsh and Pomeroy (1996) estimated areal average
*H*_{A} via a simple advection efficiency term related to SCA. Another
approach to estimate areal average estimates of advection applied internal
boundary layer integration approaches (Granger et
al., 2002, 2006) to tile models (Essery et al., 2006) and
accounted for the scaling properties of natural snow cover
(Shook et al., 1993a). Other investigations have employed
complex atmospheric boundary layer models (Liston, 1995) and
large eddy simulation (Mott et al.,
2015, 2017; Sauter and Galos, 2016) to quantify the non-linear relationships
between snow patch characteristics/geometry and advected energy. Numerical
models provide the most detailed description of the processes but are
constrained to idealized boundary conditions. The problem with all these
modelling approaches is that none have had validation with actual advection
observations and LE_{A} has largely been ignored. Only Liston (1995) considered LE_{A}, but in that study did not explicitly partition
advected energy into *H*_{A} or LE_{A} components and assumed a
constant saturated snow-free surface. This is in contrast to observations of
LE_{A} that show a dependency upon the dynamic extent of ponded meltwater
(and associated dynamic near-surface humidity content) which is prevalent in
areas of low topographic relief and limited snowmelt infiltration due to
frozen soils (Harder et al., 2017). The main challenges
in modelling these dynamics is to constrain the areal extent over which the
advection exchange takes place, quantify the gradients in scalar between
upwind and downwind surfaces, and relate the scalar gradients to advection
fluxes.

There remains a pressing need for an approach that can estimate areal
average *H*_{A} and LE_{A} contributions during snowmelt that
can easily integrate with existing one-dimensional snowmelt models. This
work seeks to understand the implications of including local-scale *H*_{A}
and LE_{A} with one-dimensional snowmelt models. To address this
objective, this paper presents a simple and easily implementable
*H*_{A} and LE_{A} 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.

The methodology to address the research objectives is briefly outlined here.
A conceptual and quantitative model framework extended the Granger et
al. (2002) advection model, hereafter referred to as the extended GM2002, to
also consider LE_{A}. The performance of the extended GM2002 was
evaluated with respect to *H*_{A} and LE_{A} observations as
reported in Harder et al. (2017). Snow-cover geometry scaling relationships
employed in the model framework (Granger et al., 2002; Shook et al., 1993b),
originally based on SCA classifications from coarse-resolution or oblique
imagery, were re-evaluated with high-resolution unmanned aerial vehicle (UAV)
imagery. The complete model framework, hereafter referred to as the Sensible
and Latent Heat Advection Model (SLHAM), was then used to explore the
dynamics of the extended GM2002 when coupled with parameterizations for frozen soil infiltration and the relationship between surface detention storage and fractional water area. Snowmelt simulation
performance and the implications of including *H*_{A} and
LE_{A} were explored by coupling SLHAM to the
Stubble–Snow–Atmosphere snowmelt Model (SSAM) (Harder et al., 2018). The
SSAM accounts for the dynamic influence of crop stubble emergence on the
sensible and latent heat and shortwave and longwave radiation terms of the
snow surface energy balance that is coupled to the mass balance of a
single-layer snowpack model to simulate snowmelt. Development and validation
for SSAM focused on representing snowmelt of shallow snowpacks in the
agricultural regions of the Canadian Prairies. SLHAM is coupled to SSAM here
as a demonstration of its ability to be coupled to existing snowmelt energy
balance models that assume continuous snow cover. Other snowmelt models could
similarly be easily coupled to SLHAM. The model performance of SSAM and
SSAM–SLHAM was also compared against the Energy Balance Snowmelt Model
(EBSM; Gray and Landine, 1988), a snowmelt model commonly implemented for
snowmelt prediction on the Canadian Prairies. In EBSM the contribution of
advection energy is indirectly addressed through simulation of an areal
average albedo that varies from a maximum of 0.8 pre-melt, a continuous snow
surface, to approach a low of 0.2 at the end of melt, which represents bare
soil rather than old snow (Gray and Landine, 1987). The areal average net
radiation, greater than typically received by a continuous snow surface, is
assumed to contribute to areal average snowmelt, thereby implicitly
accounting for advection. While a simple approach to include advection energy
for snowmelt, it is unconstrained by SCA dynamics and will overestimate melt
for low values of SCA. The implications of including advection were evaluated
with initial conditions and driving meteorology observed over two snowmelt
seasons from a research site located in the Canadian Prairies.

## 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 ≤0 ^{∘}C 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).

To scale any estimate of fetch length advection to an areal average
representation the geometric properties and extent of exchange are needed.
Over the course of melt, SCA declines from completely snow-covered to
snow-free conditions, with the intermediate periods defined by a
heterogeneous blend of both. Conceptually the advection of energy to snow
therefore is bounded by the areas of snow-free and snow-covered surfaces
that constrain energy transfer. Initial advection contributions to melt are
dominated by energy advecting from emerging snow-free patches to the
surrounding snow (Fig. 2a). The total amount of energy advected will be
limited by the smaller snow-free surface source area available to exchange
energy; all energy entrained by air movement across isolated snow-free
patches will be completely advected to the surrounding snow surfaces. At the
end of snowmelt, snow patches remain in a snow-free domain, and some energy
is advected from the warm surrounding snow-free surface to isolated snow
patches (Fig. 2b). The amount of energy advected is limited by the smaller
snow surface area available to exchange energy. When the snow surface is the
most heterogeneous, with a complex mixture of snow and snow-free patches,
advection occurs between isolated snow-free patches, surrounding snow cover,
snow-free surfaces, and isolated snow patches at the same time. Conceptually
there will be a gradual transition from isolated snow-free patch to isolated
snow patch advection constraints. Marsh and
Pomeroy (1996) and Shook et al. (1993b) found that magnitude of the snowmelt
advection flux will be greatest when SCA is 40 %–60 % and this range was
used to bound the transition of advection constraints. The advection
mechanism transitions over the course of the melt and was conceptually
related to SCA by a fractional source (*f*_{s}) term that assumes a linear
weighting between 60 % and 40 % SCA as

An *f*_{s} of 1 implies the exchange of advection energy is limited by the
snow-free patch areas and an *f*_{s} 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.

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,
LE_{A} 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 *H*_{A} and
LE_{A} 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.

*e*_{sc}: snow surface vapour pressure (kPa);
*T*_{sc}: snow surface temperature (K); *e*_{sf}: snow-free
surface vapour pressure (kPa); *T*_{sf}: 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}); *q*_{sc}: snow surface specific
humidity (kg kg^{−1}); *z*_{0s}: snow surface roughness
(0.005 m); *q*_{sf}: 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, *Q*_{A} (W m^{−2}), was presented as a power
function of patch length, *L* (m), downwind of a surface transition as

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 (*H*_{A}) and latent (LE_{A}) 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 *H*_{A} contributions to snow. To extend this
approach to LE_{A} the *a* and *b* parameterizations of GM2002 were
assumed to remain valid. The parameterization for coefficient *a* in the
case of LE_{A} 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 LE_{A} 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 LE_{A} term. The *e*_{sc} was estimated by assuming
saturation at the *T*_{sc}. The *e*_{sf} is more challenging as it varies
with the surface fraction of ponded water (*F*_{water}, no unit) as

The surface water vapour pressure for water surfaces (*e*_{wat}, kPa) was
estimated by assuming saturation at the surface temperature of the ponded
water (*T*_{wat}, K). Assuming negligible evaporation from dry soil
surfaces during snowmelt, the surface water vapour of soil (*e*_{soil}, kPa)
can be taken to be the same as actual vapour pressure observed above the
surface. The *T*_{sf} was also weighted by *F*_{water} as

where *T*_{soil} (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 *F*_{water} and SCA.

### 2.1.2 Fractional coverage of ponded water

To estimate *F*_{water}, 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

where *C* (2.1, no unit) is a coefficient representing prairie soils, *S*_{0} (–)
is a surface saturation (generally assumed to be 1), *S*_{i} (–) is the
antecedent soil saturation, *T*_{si} (K) is the initial soil temperature,
and *t*_{0} (h) is the infiltration opportunity time. The *t*_{0} term
is estimated as the cumulative hours of active snowmelt over the course of
the snowmelt period. Excess meltwater (*M*_{excess}, mm) is calculated as

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

To relate *M*_{excess} to *F*_{water}, 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 *F*_{water} is
given by the solution of

where the ratio of filled detention storage (*S*_{ret}, no unit) is determined
from

where a user-defined *S*_{max} (mm) is the maximum detention storage of the
surface. Any *M*_{excess} that is greater than *S*_{max} is removed as
runoff and thereafter unavailable to future infiltration.

### 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,

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*(*A*_{p}),
is given as a power law distribution

where *c*_{1} is a threshold value (given as the smallest patch size
observed, and hereafter taken as 1 m^{2}), *A*_{p} (m^{2}) is patch
area, and *D*_{k} (–) 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
*A*_{p} and *L* of snow patches by Granger et al. (2002) as

where *c*_{2} 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*(*A*_{p})) is also equivalent to the
probability associated with the probability of patch length (*p*(*L*)),
therefore

where i is the index for intervals of *A*_{p} 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*(*A*_{p}) at the more frequent small values of *A*_{p}. 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

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

### 2.1.5 Areal average advection

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

The terms from left to right represent the *H*_{A} from snow-free
patches, *H*_{A} to snow patches, LE_{A} from snow-free
patches, and LE_{A} to snow patches. All summation terms
constitute *H*_{A} and LE_{A} for the range of patch areas
expected, from 1 m^{2} to an environment appropriate maximum expected
patch size (*A*_{max}, m^{2}). Calculation of *H*_{A} and
LE_{A} 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 *A*_{p}. 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
*f*_{s} and (1−*f*_{s}) 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 m^{2} were removed from the analysis
as noise and classification artefacts are associated with such small patch
sizes. The 1 m^{2} 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 *H*_{A} and
LE_{A} 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.

### 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–*S*_{ret}–*F*_{water} 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, *Q*_{net} (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 *H*_{A} and LE_{A}
and the snow-free surface humidity conceptualization and snowmelt
implications from a theoretical perspective.

### 2.3.2 Sensitivity analysis

The sensitivity of SLHAM to *T*_{wat}, *T*_{a}, *T*_{soil}, *u*, and RH
variability is also explored to understand the implications upon SWE and
SCA depletion, *F*_{water}, *H*_{A}, LE_{A}, 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 *T*_{wat} is 0 ^{∘}C as meltwater immediately after discharge from an isothermal
snowpack is 0 ^{∘}C and underlying frozen soils are ≤0 ^{∘}C. Unlike the snow surface the maximum temperature of ponded
water is unconstrained by phase change so values ≥0 ^{∘}C are
expected because of possible low water surface albedos and high shortwave
irradiance (SW${}_{\mathrm{atm}}^{\downarrow}$) 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 *T*_{wat} was generally >0 and
<2 ^{∘}C. This range in *T*_{wat} was used to test the
sensitivity of the *T*_{wat} upon SLHAM dynamics. Intermittency of
observations and inherent uncertainties in thermography prevented a more
precise estimation of *T*_{wat}. The ranges of all other variables were
selected to represent conditions commonly experienced during snowmelt on the
Canadian Prairies.

### 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 *T*_{soil} required for SLHAM come
from infrared radiometers (Apogee SI-111) deployed on mobile tripods to
snow-free patches. Unfortunately, no time series of *T*_{wat} observations
are available and values or models to describe *T*_{wat} 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 *T*_{wat} are unavailable, a simple physically
guided representation of *T*_{wat} takes the following form:

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

## 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 *H*_{A} and LE_{A} 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.

## 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 *D*_{k} 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 *D*_{k} 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.

## 3.3 Implications of including advection in snowmelt models

### 3.3.1 Advection dynamics in scenario simulations

The dynamics of the various scenarios are expressed through visualizations
of SWE depletion (Fig. 8) and magnitudes of the *H*_{A}, LE_{A}
and net advection terms (Fig. 9). A critical consequence of including
SCA in snowmelt calculations is that areal average melt rates will vary
between a continuous and heterogeneous snow surface. The *Q*_{net} driving
melt in a one-dimensional melt model is in terms of a flux density; an
energy flux with a unit area dimension (W m^{−2}) where exchange is
limited to the SCA. As the SCA decreases the corresponding areal average
energy to melt snow will also decrease, which will decrease the areal average
melt rate. This is evident in the melt rate of the No Advection scenario,
which decreases with time as the SCA decreases. Including energy from
advection, for the dry surface, wet surface, and dry to wet surface
advection scenarios, causes the SWE to deplete faster as there is now an
additional energy component that increases as SCA depletes. In these
scenarios the additional energy gained from advection is greater than the
reduction of areal average *Q*_{net} as SCA decreases. LE_{A} from a
constant wet surface is greater than any other advection scenario. Despite a
reduction in *H*_{A} from the cooler surface and therefore an overall slower
melt, the consistently positive LE_{A} towards the snow leads to a large
net advection flux. In contrast, a consistently warm dry surface has a much
higher *H*_{A} flux, and faster melt rate, than the wet surface that is
partly compensated by a negative LE_{A} due to sublimation and a decrease
in the overall energy for melt from advection. When the surface wetness is
parameterized by detention storage and frozen soil infiltration capacity,
dry to wet surface, the snow-free surface is dry and warm in the early
stages of melt and LE_{A} is negative and limits melt, as in the dry surface scenario. As melt proceeds and *F*_{water} begins to increase, the
upwind *T*_{sf} cools and the humidity gradient switches, resulting in
positive LE_{A} and a decrease in *H*_{A}, which compound to slow melt
relative to the dry surface scenario. There are clear implications for the
timing of melt and thus snow hydrology depending upon the upwind condition.
It is evident that SLHAM can quantify the key advection behaviours in
relation to the upwind surface dynamics.

### 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 *T*_{soil} and *u* have the largest influence
on advection contributions to snowmelt. This is expected as *u* and
*T*_{soil} variables quantify the first-order controls driving advection,
the air mass movement and horizontal scalar gradients respectively. In
contrast the *T*_{wat}, *T*_{a}, 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, *F*_{water}>0, which forces the upwind surface to the assumed water surface
temperature. The consequent sign of the surface humidity gradient will
influence whether LE_{A} 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 *F*_{water}. Changes in *H*_{A} and
LE_{A} 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,
*T*_{wat}, *T*_{a}, and RH have a relatively limited influence upon overall
SWE depletion compared to *T*_{soil} and *u*. In the absence of *T*_{wat}
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.

### 3.3.3 Advection dynamics in coupled advection and snowmelt models

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.

## 3.4 Energy balance compensation

An unappreciated dynamic of local-scale advection during snowmelt is that
LE_{A} and *H*_{A} 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 LE_{A} 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 LE_{A} was negative for much of the
melt period, which compensated for the consistently positive *H*_{A}
term (Fig. 12). LE_{A} only increased, enhancing the positive
*H*_{A} contribution, near the end of melt in 2015 when increased
surface wetness led to a positive LE_{A} term.

The advection fluxes may also be of opposite sign to the sensible
(*H*_{snow}) and latent (LE_{snow}) turbulent fluxes between the snow
surface and the atmosphere. Inclusion of the advection process therefore
influences the overall sensible and latent heat exchange at the snow surface
(net exchange). This interaction is further complicated by the varying
SCA of the SSAM–SLHAM model versus the complete snow cover assumption of
SSAM. Including advection decreased cumulative LE by 1.4 MJ in 2015 and by
3.9 MJ in 2016 (Table 6). Cumulative *H*, when including advection,
increased by 0.2 MJ in 2015 and by 5.7 MJ in 2016. The net exchange when
including advection shows that the inclusion of LE_{A} decreases the
influence of *H*_{A}; the change in net exchange is lower than the change in
*H* exchange (Table 6). The role of advection in modifying net exchange is
clearly complex and varies by season. Despite differences in magnitude, the
opposite signs of LE_{A} and *H*_{A} demonstrate that these energy
contributions partially compensate for one another, therefore reducing the
net influence of advection on snowmelt. This compensatory relationship has
been missed by the focus on *H*_{A} in snowmelt advection research, which
has therefore overemphasized the contribution of *H*_{A} to snowmelt. This
compensatory mechanism also helps to explain why observed latent heat fluxes
are often much smaller than model predictions in the meltwater-ponded
Canadian Prairies during melt (Granger et al., 1978). The
compensation of *H*_{A} by LE_{A} will be a more important interaction on
the Canadian Prairies, or similar level environments, but perhaps less so in
mountain regions where complex terrain leads to rapid meltwater runoff.

## 3.5 To advect or not to advect?

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 *a*_{snow} are not limited to EBSM as other
*a*_{snow} parameterizations, especially temperature-dependent ones,
typically underestimate *a*_{snow} during melt and therefore indirectly, and
perhaps unintentionally, account for advected energy contributions (Pedersen and Winther,
2005; Raleigh et al., 2016). While modelled *a*_{snow} values that
underestimate actual *a*_{snow} 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 *f*_{s} and its possible relation to the advection
efficiency term of Marsh and Pomeroy (1996). An estimate of *T*_{sf}
is needed to implement this framework and will limit application of SLHAM in
its current form, as modelling *T*_{sf} 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
*T*_{soil} and *T*_{wat}. In the interim, the *T*_{wat}
assumptions in Eq. (15) may be used but need to be tested further. A
regression of *T*_{soil} to incoming shortwave radiation and
*T*_{a} 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.

To date, the development of easily implementable and appropriate models to
estimate the advection of *H*_{A} and LE_{A} 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 *H*_{A} and LE_{A} 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 LE_{A} and
*H*_{A} to snow as an areal average contribution. A scenario-based
analysis of the model revealed the compensatory influence of LE_{A}
from a warm but dry surface; the LE_{A}-driven sublimation offsets
*H*_{A} 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).

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

The SLHAM framework requires a *T*_{soil} value which is a challenging
variable to explicitly model during snowmelt. To provide an interim solution
a multiple linear regression is developed to estimate *T*_{soil} from
SW${}_{\mathrm{atm}}^{\downarrow}$ and *T*_{a}. 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 *T*_{soil} is the net radiation,
whose variability is dominated by SW${}_{\mathrm{atm}}^{\downarrow}$, and
turbulent fluxes, which are dependent upon the *T*_{a} gradients.
During nighttime *T*_{soil} is very similar to *T*_{a} while
during daytime the additional energy from SW${}_{\mathrm{atm}}^{\downarrow}$
heats the surface to temperatures above *T*_{a}. A multiple regression
that contains these parameters provides a simple but effective way to
estimate *T*_{soil} 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
*T*_{a} is observed with a shielded Campbell Scientific HMP45C212 and
SW${}_{\mathrm{atm}}^{\downarrow}$ is observed with a Campbell Scientific CNR1
with both sensors 2 m above the ground surface. The *T*_{soil}
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 *T*_{soil} (^{∘}C) in relation to SW${}_{\mathrm{atm}}^{\downarrow}$ (W m^{−2}) and *T*_{a} (^{∘}C) is

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 *T*_{soil} 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
*T*_{soil} 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
*T*_{soil} should be explicitly modelled.

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.

The authors declare that they have no conflict of interest.

This article is part of the special issue “Understanding and predicting Earth system and hydrological change in cold regions”. It is not associated with a conference.

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

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.: Simulation of the Ablation of Prairie Snowcovers, PhD Thesis, University of Saskatchewan, Saskatoon, Saskatchewan, Canada, 1995.

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.