Modelled sensitivity of the snow regime to topography , shrub fraction and shrub height

Recent studies show that shrubs are colonizing higher latitudes and altitudes in the Arctic. Shrubs affect the wind transport, accumulation and melt of snow, but there have been few sensitivity studies of how shrub expansion might affect snowmelt rates and timing. Here, a three-source energy balance model (3SOM), which calculates vertical and horizontal energy fluxes – thus allowing within-cell advection – between the atmosphere, snow, snow-free ground and vegetation, is introduced. The three-source structure was specifically adopted to investigate shrub–tundra processes associated with patchy snow cover that singleor two-source models fail to address. The ability of the model to simulate the snow regime of an upland tundra valley is evaluated; a blowing snow transport and sublimation model is used to simulate premelt snow distributions and 3SOM is used to simulate melt. Some success at simulating turbulent fluxes in point simulations and broad spatial pattern in distributed runs is shown even if the lack of advection between cells causes melt rates to be underestimated. The models are then used to investigate the sensitivity of the snow regime in the valley to varying shrub cover and topography. Results show that, for domain average shrub fractional cover ≤ 0.4, topography dominates the preand early melt energy budget but has little influence for higher shrub cover. The increase in domain average sensible heat fluxes and net radiation with increasing shrub cover is more marked without topography where shrubs introduce wind-induced spatial variability of snow and snow-free patches. As snowmelt evolves, differences in the energy budget between simulations with and without topography remain relatively constant and are independent of shrub cover. These results suggest that, to avoid overestimating the effect of shrub expansion on the energy budget of the Arctic, future large-scale investigations should consider wind redistribution of snow, shrub bending and emergence, and sub-grid topography as they affect the variability of snow cover.


Introduction
The effects of shrub expansion or retreat on tundra surface energy balance have garnered much attention over the past decade.Increasing evidence from field observations, remote sensing and models suggests that warming in the Arctic has led to a "greening" that is partly due to the densification and expansion of shrub patches (Sturm et al., 2005a;Jia et al., 2006;Raynolds et al., 2006;Tape et al., 2006;Loranty et al., 2011;Tremblay et al., 2012).Chapin III et al. (2005) estimated that reductions in the duration of snow cover for Arctic Alaska has increased local atmospheric heating by 3 W m −2 per decade, which is similar in magnitude to local warming predictions spanning multiple decades from a doubling of atmospheric CO 2 , and continuation of current trends in shrub expansion could amplify this heating by two to seven times.
In addition to climate-induced changes, grazing can also control shrub coverage and height.In northern Finland, intensive year-round reindeer grazing prevents shrub growth such that albedo during the snow season is higher than in neighbouring Norway, where more moderate seasonal grazing management practices do not limit shrub height (Kitti et al., 2009;Cohen et al., 2013).For this reason, the contribution of grazing as one of many local solutions to control Published by Copernicus Publications on behalf of the European Geosciences Union.
Differences and changes in shrub cover also affect the snow mass balance of tundra because the height, density and location of shrubs affect snow distribution by reducing near-surface wind speeds within and downwind of shrub patches (Essery et al., 1999;Sturm et al., 2001), trapping wind-blown snow around isolated patches or at the edges of large patches (Essery and Pomeroy, 2004b).Pomeroy et al. (2004) found that the standard deviation of snow water equivalent (SWE) in a sub-Arctic basin increased with decreasing shrub height, decreasing cover and increasing exposure to wind.Liston et al. (2002) predicted that replacement of lowgrowing vegetation over a 4 km 2 domain in Arctic Alaska by shrubs with a 50 cm snow-holding capacity would decrease domain-averaged sublimation by 68 % and increase snow depth by 14 %.Essery and Pomeroy (2004b) argued that increased snow amount in shrubs is limited by the supply of wind-blown snow from open areas and predicted that the average SWE over an area in Arctic Canada with 30 % shrub cover would not increase further for increases beyond 1 m in shrub height.
Ground-based measurements, remote sensing and modelling studies in the Yukon (Pomeroy et al., 2006;Bewley et al., 2010), the Northwest Territories (Marsh et al., 2010), Alaska (Sturm et al., 2005a), Fennoscandia (Cohen et al., 2013) and the pan-Arctic (Loranty et al., 2011) have all shown that shrub branches exposed above snow decrease the surface albedo and increase the absorption of solar radiation.Shrubs can even absorb radiation while buried because short-wave radiation penetrates snow (Warren, 1982;Baker et al., 1991;Hardy et al., 1998).Sturm et al. (2005a) estimated that transitions from shrub-free to shrubby conditions could increase absorption of solar radiation over the snow-covered period by 69 to 75 %.Once exposed, branches can be 20 • C warmer than the surrounding snow (Pomeroy et al., 2006), increasing turbulent heat and long-wave radiation fluxes from the exposed shrub canopy.Tall shrubs reduce the short-wave radiation reaching the snow surface by shading but increased long-wave radiation and sensible heat fluxes from the canopy to the snow can give higher melt rates for snow beneath shrubs than for unshaded snow (Bewley et al., 2010).For example, Sturm et al. (2005b), Pomeroy et al. (2006) and Marsh et al. (2010) observed higher melt rates where shrubs were exposed above the snowpack than where shrubs were buried.
Improved understanding of shrub-tundra processes from field investigations has motivated recent model developments (e.g.Bewley et al., 2007Bewley et al., , 2010;;Marsh et al., 2010;Ménard et al., 2012).However, one of the remaining difficulties in modelling shrub-tundra is how to represent sparse canopies and horizontal interactions between shrub and non-shrub surfaces.Land surface models (LSMs) operate at large scales and represent the Earth surface as a series of grid boxes.Sub-grid heterogeneity is often addressed by classifying dif-ferent surfaces within a grid box as "tiles".The energy balance for each tile is calculated separately and area weighted to produce the grid-box-average energy balance.This model structure does not allow for horizontal advection of heat to the snowpack from shrub canopies or bare ground patches even though discontinuous snow cover is known to substantially enhance snowmelt (e.g.Liston, 1995;Neumann and Marsh, 1998;Granger et al., 2002Granger et al., , 2006;;Pomeroy et al., 2003).The solution proposed in many LSMs (see e.g.Sellers et al., 1986for SiB, Verseghy, 1991 for CLASS) is to calculate separate energy balances for vegetation and bare ground or snow surfaces within a grid box (so-called "two-source" models).However, this does not solve the problem for landscapes with exposed vegetation and patchy snow cover.Models that calculate a single energy balance for composite snowcovered and snow-free ground are known to fail to simulate late-lying snow patches (Marsh and Pomeroy, 1996;Bewley et al., 2010).This causes two issues in model performance: excessive melt rates and modelled sensible heat fluxes that are in the opposite direction to observations for patchy snow cover.In addition, these models do not consider horizontal advection of snow mass between grid boxes by wind.
A number of modelling studies (e.g.Kaplan and New, 2006;Lawrence and Swenson, 2011;Bonfils et al., 2012) have investigated the effects of shrub expansion on snow availability, distribution and/or energy fluxes, but none have considered the effects of topography at basin or smaller scales with a high-resolution model.This study aims to address this gap by setting two objectives.The first is to introduce and evaluate a new model (designated "3SOM" for 3source model) that addresses issues described above for energy fluxes over sparse shrub tundra canopies and heterogeneous snow cover.3SOM is a distributed model adapted from the single-point two-source (shrub and snow) model of Bewley et al. (2010), which was itself adapted for cold environments from Blyth et al. (1999), to solve three energy balances for bare ground, snow and shrubs simultaneously within each grid box.The second objective is to investigate the sensitivity of snow accumulation and ablation to shrubs and topography.This is achieved by initializing 3SOM with outputs from an existing distributed blowing snow model, DBSM (Essery and Pomeroy, 2004b).The combined models are used to estimate changes in net radiation, turbulent heat transfer and spatial distributions of SWE associated with variations in shrub cover with or without the influence of topography.

Sites and data
Model evaluation and sensitivity analyses were performed using data from a 1 km × 1 km area around a sub-Arctic tundra valley, the Granger Basin (GB) in the Yukon Territory, Canada.GB lies within the larger Wolf Creek Research Basin (WCRB), a 179 km 2 drainage basin 15 km south of Whitehorse that has been the subject of more than 100 peer-reviewed papers and technical reports to date on cold-region hydrology (Richard Janowicz, Yukon Environment, personal communication, 2013; see e.g.Pomeroy et al., 2003Pomeroy et al., , 2004Pomeroy et al., , 2006;;MacDonald et al., 2009;Bewley et al., 2010).A digital elevation model, canopy height map and fractional vegetation cover map of GB (Fig. 1) were derived from lidar data obtained during an airborne campaign in August 2007.Details of the campaign and of the lidar data processing are available in Chasmer et al. (2008) and Hopkinson and Chasmer (2009).
The valley site is situated in the sub-alpine ecozone of WCRB and is characterized by a southwest-facing slope (13 • ), a northeast-facing slope (17 • ) and a valley bottom with a stream flowing in a southeasterly direction.Tall (1.80 m on average) willows cover the riparian zones of the valley bottom; some of these become buried under snow in the winter and spring up during melt (Pomeroy et al., 2006;Ménard et al., 2012).Birch shrubs (0.3 to 1.5 m tall) are also widespread within the basin, generally being shorter on the exposed, well-drained plateau and taller near the wet valley bottom.
Heat fluxes and snowmelt simulated at points by 3SOM were evaluated using data from spring 2003 and 2004, and the combined DBSM and 3SOM models were driven on a high-resolution grid with data from winter and spring 2007-2008.Driving and evaluation data are used from seven automatic weather stations: 1.A station located amongst tall willow shrubs (> 2 m) in the riparian zone at 1363 m a.s.l., which has been measuring meteorological data since spring 1997.
2. Two stations located above short shrubs (birches of 0.32 m on average, Bewley et al., 2010)   in Fig. 1.Air temperature T a , wind speed u, wind direction u dir , relative humidity RH, incoming short-wave radiation SW in , incoming long-wave radiation LW in and snow depth S d were measured every half hour at the tall and short shrub stations.Data from the alpine station were used to fill gaps in wind speed and direction data required by the blowing snow model.As wind data were missing from both the second plateau and the alpine stations for 42 consecutive days from January to mid-February 2008, the missing data were infilled with data collected at WIA. Daily snowfall measurements were not made within the Wolf Creek Basin during the study period, so snowfall data were obtained from the WIA meteorological station.A correction factor for differences in snowfall between GB and WIA due to the difference in elevation was applied to the snowfall rate, as detailed in Ménard et al. (2012).Average air temperature during the 2007-2008 snow season, when the blowing snow model was run, was −10.9 • C with a total of 36 h (no more than 5 consecutive) above 0 • C; as 24 h were below 1 • C and all were below 2 • C it was deemed appropriate to assume no winter melt event.Visual field observations established that occasional overnight snowfall events occurred during the melt season in 2003, 2004 and 2008 but SR50 measurements did not show any difference in snow depth before and after the event, suggesting that they had no or very little effect on the snow mass balance at the site.Sensible and latent heat fluxes measured by eddy covariance towers at the tall and short shrub stations in spring 2003 and 2004 are compared with point simulations in Sect.3.4.A ground heat flux plate also provided measurements at the short shrub station in 2003.Data from the slope stations are used to evaluate the modelled distribution of incoming shortwave radiation (discussed in Sect.4.2).Snow depths were measured every 5 m and snow density every 25 m with an ESC-30 snow tube along three transects (Fig. 1) of approximately 400 m length.The transects cover the main topographic features in GB (i.e. the slopes and the valley) and a wide range of shrub height (0 to 2.60 m) and cover (0 to 0.63).In 2008, the transects were monitored every three to four days during the melt season until all the snow, except in drifts on the northeast-facing slope, had melted.The transect measurements are used to evaluate simulations of snow depth before and during melt in Sect.4.3.

Three-source model description and evaluation
This section describes the surface, sub-surface and snow schemes of 3SOM, represented schematically in Fig. 2, and presents an evaluation of point simulations against heat flux and snowmelt measurements in spring 2003 and 2004.

Surface and sub-surface schemes
The surface scheme has separate energy balances for snow (F s ), bare ground (F g = 1 − F s ) and exposed vegetation fractions (F v , independent of F s and F g ) coupled through a network of resistances, allowing horizontal and vertical transfer of heat between sources within a grid box.Following the convention that radiative fluxes are positive towards the surface and turbulent and ground heat fluxes are positive away from the surface, the separate energy balance equations for vegetation, snow-free ground and snow, identified by subscripts v, g and s, respectively, are where R is the net radiation, H and LE are the sensible and latent heat fluxes respectively, M is the snowmelt heat flux, G g is the heat flux from the soil to the soil surface and G s from the snowpack to the snow surface.
The total moisture flux, ignoring evaporation from dormant vegetation, to the atmosphere is and, thus allowing advection from snow-free to snowcovered surfaces, the sensible heat flux is Neglecting multiple reflections between the vegetation and the ground or snow, net solar radiation absorbed at the surface is for the snow tile, for the vegetation tile and for the ground tile, where SW in is the incoming short-wave radiation, τ is the canopy transmissivity and α is the albedo.Unlike Bewley et al. (2010), who used a time-averaged τ calculated from hemispherical photographs for their point model evaluation, τ is calculated here as a function of exposed vegetation fraction following Bewley et al. (2007).Long-wave radiation is emitted both upwards and downwards from the vegetation canopy, so the net long-wave radiation for vegetation temperature Including absorption of long-wave radiation from the vegetation and the atmosphere, the net long-wave radiation for the snow and vegetation surfaces is and where T s and T g are the snow and ground surface temperatures.
The sensible heat fluxes at the snow, bare ground and vegetation surfaces are respectively and where ρ is air density, C p is the heat capacity of air, T c is the temperature of the canopy air space, and r as , r ag and r av are aerodynamic resistances for snow, ground and vegetation sources.Transpiration from the dormant vegetation is assumed to be negligible (E v = 0), but moisture fluxes over the snow and ground surfaces are calculated as and where Q sat (T ) is saturation humidity at temperature T and atmospheric pressure P , and Q c is the specific humidity of the canopy air space.
The resistance for evaporation of soil moisture from bare ground is parameterized as where θ u is the unfrozen soil moisture content and θ c is the critical volumetric soil moisture content, as in Best et al. (2011).
The aerodynamic resistance for energy exchange between the canopy air space and reference height for temperature, z t , is given by where k is the von Karman constant, u * is the friction velocity and d is the displacement height, taken to be 2/3 of the exposed vegetation height.The aerodynamic resistance for exchanges between the vegetation and the canopy air space is given by where z 0v is the roughness length of the vegetation, assumed to be 1/10 of the exposed vegetation height.Bewley et al. (2010) found that for understorey resistances, the form used in Eqs. ( 20) and ( 21) leads to an increase in sublimation with increasing vegetation cover rather than suppressing turbulent transport from the underlying snow.They proposed a modified form, which adapted for 3SOM gives 1 where C is a dense canopy exchange coefficient given the value of 0.004 (Zeng et al., 2002) and u * is the friction velocity such that where z u is the reference height for wind speed minus snow depth and z 0 is the momentum roughness length for which, as for the aerodynamic resistances, there exist a number of different formulations.For F v = 0 or 1, Eqs. ( 22) and ( 23) are reduced to the usual logarithmic forms for resistances over homogeneous surfaces.Following Mason (1988), an "effective" roughness length, intended to represent a spatial average in a heterogeneous terrain, is calculated by weighting the individual roughness values of the three sources such that Heat fluxes into the ground and snow surfaces are calculated as and where λ soil and λ snow are the thermal conductivities of surface soil and snow layers of temperature T g1 and T s1 and thickness z g1 and z s1 .Details of the soil thermodynamics can be found in Best et al. (2011).The three surface sources share a single soil column of four layers extending to 1.5 m depth per grid box.3SOM uses a diffusive heat conduction approach from the ground and snow sources to the top soil layer to calculate the top soil temperature.Transfers of heat and water between soil layers are neglected; this makes the model unsuitable for distributed hydrological modelling over complex topography but has little impact on the simulations of turbulent fluxes and snowmelt that are the focus of this paper.
An implicit solution is used to find increments in surface temperatures over each model time step.The flux parameterizations are linearized in temperature and humidity increments (using the Clausius-Clapeyron equation to expand Q sat ) and substituted into the balance equations ( 1), ( 2), ( 3), ( 6) and ( 5).This gives five equations for five unknowns (increments in vegetation, snow and ground surface temperatures, and temperature and humidity increments in the canopy air space) that are solved using LU decomposition (Press et al., 1996).A first solution is found assuming no snowmelt; if this gives a snow surface temperature greater than 0 • C, the solution is repeated assuming T s = 0 • C and the residual of the energy balance is used to melt snow at rate M s .

Vegetation and snow cover fractions
The exposed vegetation fraction, allowing for bending and burial of shrubs by snow, is taken from Liston and Hiemstra (2011) such that where F v0 is the snow-free vegetation fraction, h c is canopy height and B = 0.85 is a bending parameter estimated by Ménard et al. (2012) for the valley site in the 2007-2008 season.
Although the lognormal distribution has been used in many studies to describe snow cover fraction (e.g.Pomeroy et al., 1998b;Liston, 2004;DeBeer and Pomeroy, 2009), Essery and Pomeroy (2004b) argue that a closed-form solution is preferable to the parametric form of the lognormal snow cover distribution curve.Instead they propose a parameterization, adapted from Yang et al. (1997) and used in this study, which closely fits the lognormal distribution such that where S d is snow depth (m) and σ 0 is the premelt standard deviation of snow depth.DBSM was used to calculate the standard deviation of snow depth for the points and the distributed simulations; σ 0 was found to equal 0.20 at both the short and tall shrub sites.

Snow scheme
Snow albedo decay during melt takes an exponential form following Verseghy (1991) such that while cold albedo decay is linear following Douville et al.
(1995) such that where where α fresh is fresh snow albedo and Sf thresh is the minimum amount of snowfall (10 kg m −2 ) required to refresh albedo to α fresh .
The snow compaction and thermodynamic schemes used in 3SOM are taken from the Joint UK Land Environment Simulator (JULES).A simple description of the scheme is provided below for convenience but the reader is referred to Best et al. (2011) for a full description.The number of snow layers is adjustable and up to three snow layers are used here.Separate temperatures, densities, ice and water contents are calculated for each layer and updated at each time step.The thickness of the snow layers varies as snow compacts or fresh snow accumulates.Less-dense freshly fallen snow (100 kg m −2 ) is accumulated in the surface layer while mechanical compaction leads to higher snow densities in deeper layers.Upon growth of the snowpack, the first snow layer increases in thickness until it reaches a prescribed depth, at which point the layer splits in two.This process is repeated upon accrual of each layer until the maximum number of prescribed layers is reached.At this point, the thickness of the top snow layer stays fixed and subsequent increases in snow depth are accommodated by increasing the depth of the bottom layer.During snowmelt, snow depth is accommodated by decreasing the depth of the deepest snow layer first.When liquid water in excess of liquid capacity is reached, water is moved to the layer below.Refreezing of liquid water within the snowpack and the resulting latent heat release occurs when water flows to a layer where snow temperature < 0. The temperature of the layer is then recalculated accordingly.Fresh snow and frost are added as an interim layer 0 and the snow layers and snow pack properties are recalculated accordingly following the conservation of mass and energy.

Evaluation of the model at points
Considering the single-point two-source-model presented by Bewley et al. (2010) as a benchmark against which 3SOM is compared, 3SOM is first evaluated against observations at the tall and short shrub stations throughout the melt periods of 2003 and 2004.The model was initialized using manual measurements of vegetation height, premelt snow depth and snow water equivalent at both stations.Results are presented in Fig. 3. Melt started earlier in 2003 and the snow at the short shrub site had disappeared by 30 April, but melt at the tall shrub site stagnated because of a drop in air temperature between 30 April and 6 May.The small increase in measured snow depth (2 cm) and SWE (6 kg m −3 ) on 3 May 2003, which persists for the next two days, suggests that a snowfall or wind-induced snow redistribution event occurred during this 7-day cold period.Although this may have a small effect on end-of-melt modelled snow, the snow pack is mostly shallower than measurements because melt rates are overestimated at both stations in 2003.Despite this, 3SOM is able to reproduce late-lying snow at the tall shrub site, which was one of the features that the two-source model presented by Bewley et al. (2010) could not reproduce.Modelled melt rates are 49 % larger at the tall than at the short shrub site in 2003 and 56 % larger in 2004 on average due to advection of heat from the shrubs to the snow.
In both years, the model is able to reproduce the direction and magnitude of energy exchanges with the atmosphere over a dynamic surface.Contributions of the individual sources to energy fluxes are shown in Fig. 4. At the short shrub site, H g from a small snow-free fraction is positive during daytime from the beginning of the run but H v is negligible because little vegetation is exposed and total sensible heat fluxes (measured and modelled) are predominantly negative until the snow has melted.At the tall shrub site, H is positive during daytime owing to large upwards heat fluxes from exposed vegetation and negative at night, with heat transferred from the atmosphere to both snow and vegetation.Downwards sensible heat fluxes at night are slightly overestimated by the model, particularly at the tall shrub site.Latent heat fluxes from snow and bare ground are upwards during the day, when they are overestimated by the model, and small at night for both sites.From 23 to 25 April 2003 and 1 to 3 May 2004 at the short shrub site, the two-source model proposed by Bewley et al. (2010) simulated sensible heat fluxes of opposite directions, but of similar magnitude, to measurements leading to model errors of up 250 W m −2 .Such errors do not occur with 3SOM because, as Bewley et al. (2010) suggested, calculating separate energy balances for the snow and snow-free ground produces a more accurate representation of the turbulent fluxes during snowmelt over a heterogeneous surface.
Table 2 shows quantitative assessments of the modelled fluxes.Although just above half of the performance metrics are improved compared to the two-source model, the addition of a third source has solved some of the issues identified with the two-source model, namely the simulation of latelying snow patches and the direction of sensible heat fluxes during snowmelt.Furthermore, although some of the errors are large, turbulent fluxes modelled at points have been compared here with measurements that are influenced by heterogeneous upwind surface conditions.Ideally, a footprint model should be used to weight the modelled fluxes for comparison with measurements.However, accurate determination of the measurement footprint and flux sources would constitute an entirely separate study.Given these improvements in point scale process representations, the model was distributed to perform simulations over a landscape.

Description and evaluation of distributed simulations
This section describes how meteorological data are adjusted for driving distributed simulations and presents results obtained using DBSM to simulate the evolution of accumulating snow from September 2007 to April 2008 and then using 3SOM to simulate snowmelt through April and May 2008.The sequential coupling of the models captures the dominant snow processes in each season but neglects any winter melt or spring redistribution events; a full coupling of the models is in development.
DBSM calculates changes in SWE within a grid box over time as where q s is the rate of sublimation from blowing snow and q t is the rate of mass transport by blowing snow.Parameterizations of the blowing snow fluxes depending on local wind speed, vegetation height, erodability of the snow surface, air temperature and relative humidity are described by Essery et al. (1999), Pomeroy and Li (2000) and Essery and Pomeroy (2004a).

Effect of grid-box resolution on process representation
DBSM was run at 4, 8, 12, 16 and 20 m resolution in order to determine the grid-box size at which snow, vegetation and slope scale heterogeneity at GB could be best represented.Investigations into the representation of snow heterogeneity showed that σ 0 decreases as a power function of increasing grid-box size (Fig. 5).At higher resolution, the surface features -such as slope, aspect, shrub patches -which are responsible for the heterogeneity of the snow distribution are explicitly represented; for example, DBSM redistributes more snow on a grid box downwind of a grid box with high shrub cover and height than on one upwind.On the other hand, at lower resolution, the grid-box scale exceeds the scale at which surface features affect snow redistribution, thus having the net effect of homogenizing the landscape and, by extension, the snow variability dependent upon it.
With regard to the representation of shrub cover heterogeneity, Hopkinson et al. (2005) showed that there could be large errors in assessing short-vegetation parameters from single laser pulse airborne lidar sensors because the open canopy structure within and between shrubs allows numerous penetrations of the laser pulse into the foliage.As a result, the number of laser pulse returns per grid box used to distinguish between ground and shrubs needs to be high enough to minimize the noise from "false" within-shrub returns but low enough to allow for the sparse between-shrub canopy structure to be identified.As a consequence, both DBSM and 3SOM were run at 8 m resolution as this was found to be most representative of both snow heterogeneity and of canopy structure.

Distribution of the driving data and initial conditions
Some of the slopes in GB are in excess of 20 • but the altitude range is less than 260 m, so only the incoming shortwave radiation and the wind speed are distributed for model driving; Pomeroy et al. (2003) measured all of the meteorological forcing variables on different slopes and at different elevations in the basin and found this to be a reasonable assumption.
Measurements of total incoming short-wave radiation made with a levelled radiometer at the tall shrub site were partitioned into diffuse and direct components following the empirical method proposed by Erbs et al. (1982) and adjusted for the slope and aspect of each grid box according to the solar elevation and azimuth for each time step as described by Oke (1987) and Liston and Elder (2006).Comparisons between predictions and measurements made parallel to the slopes at the slope stations in 2003 are shown in Fig. 6.
Wind speeds were distributed using the Mason and Sykes (1978) model of flow over topography, which requires wind speed and direction at an exposed location as inputs.These data were obtained from measurements at the plateau station when available and gaps were filled with data from the alpine or WIA stations.Average wind directions for periods of overlapping data at these stations were within the 45 • discretization used by DBSM.Average wind speed at WIA was lower than the plateau's by only 0.16 m s −1 but was 1.18 m s −1 higher at the alpine station.Wind speeds at the alpine station were multiplied by a factor of 0.77 found by minimizing the difference between measured and DBSM modelled snow depths used as initial conditions for 3SOM.
DBSM was run from 1 September 2007 to 19 April 2008 to generate SWE and snow depth grids as initial conditions for melt simulations by 3SOM, which was then run from 19 April to 28 May.The standard deviation of premelt snow depth was calculated, using the snow depth grids produced by DBSM, from each grid box and its adjoining boxes (nine in total).Lacking distributed measurements or an adequate model of below-ground heat and water transport, the soil temperature was initialized to the measured temperature at the tall shrub site (the soil temperature at the plateau station on 19 April only differed by 2 • C despite different snow conditions).As a simple method of making the valley bottom wetter than the plateau, initial soil moisture content θ was scaled according to a topographic index (Kirkby and Weyman, 1974;Beven and Kirkby, 1979) such that θ = θ sat TI/TI max , where θ sat is the volumetric soil moisture content at saturation and TI max is the maximum value of topographic index found within the domain.

Evaluation of distributed simulations
Modelled snow depths and standard deviation of snow depths are compared with manual measurements and discretized per melt period and three topographic features: the northeastfacing slope (NFS), the valley and the southwest-facing slope (SFS) (Fig. 7).The modelled standard deviation was extrapolated to 4 m resolution following the relationship between grid-box size and σ 0 shown in Fig. 5 in order to provide a more direct comparison between measurements and model results.Measured and modelled snow depths along the transect with the longest persisting snow cover are also shown (Fig. 8) to provide a more spatially explicit comparison.Figures 7 and 8 show that the broad features of snow distribution and disappearance are captured; DBSM reproduces the drifts at each side of the transect and 3SOM melts the southwestfacing slope faster than the northeast-facing slope.Nevertheless, errors in both space and time remain; errors in snow depths increase during melt but decrease at the end of melt.
Figure 7 shows that DBSM and 3SOM overestimate the standard deviation of snow depth on the northeast-facing slope but Fig. 8 shows that this is due to DBSM underestimating the width of the snow drift at the top of the slope; this error then persists throughout the 3SOM simulation.The simulated disappearance of snow from the southwest-facing slope later than observed is consistent with the overestimation of snow depth before the start of melt.This may occur because the model only accounts for atmospheric advection of heat within and not between grid boxes.Warming of air by upwards sensible heat fluxes over snow-free patches and warming of snow by downwards heat fluxes as the air then passes over snow patches is a process that has been well documented at GB (Granger et al., 2002(Granger et al., , 2006;;Essery et al., 2006).Figure 9 shows the southwest-facing slope on three dates during the first two weeks of the melt season; the red arrow points to a snow-free patch close to one of the transects that was small enough on 20 April to be contained within a single model grid box but had grown much larger by 27 April.Tests simulations were performed to assess whether increasing the resolution of the grid box to a size more representative of larger shrub patches would improve model results.These simulations revealed that the more homogeneous snow cover causes the turbulent exchanges regime to be become dominated by latent heat (28 W m −2 climatological area average at 20 m against 14 W m −2 at 8 m) instead of sensible heat fluxes (−6 W m −2 at 20 m against 28 W m −2 at 8 m) and thus did not solve the advection.Section 3.4 shows that 3SOM performs well in point simulations when initial snow depth, SWE, vegetation fraction and vegetation height are specified from direct measurements.A number of factors could contribute to errors and some of the poor quantitative statistics in Table 3 in the distributed simulations.Firstly, scales between model (8 m × 8 m) and observations (snow probe at 5 m interval) differ; the sensitivity of σ 0 to grid-box size shows that scaling errors are inevitable.Secondly, manual snow depth measurements along the transects are also prone to sampling errors, as can be seen in Fig. 7 where mean snow depth at the end of melt on the northeast-facing slope is higher than during melt; although snow surveys were conducted as close as possible to the previous surveys, individual points generally differed by a few centimetres (between 5 cm to 1 m) because of the  destructive nature of sampling snow density and depth.Finally, given the high resolution of the model, small errors in the lidar mapping of topography and vegetation will have some influence on model results.

Effects of shrub expansion
Despite issues discussed above, the DBSM and 3SOM models are able to capture the evolution of broad spatial snow patterns along transects representative of the topographic features within GB and diurnal and seasonal variations in energy fluxes at two points representative of short and tall shrub cover.The models are now used in a sensitivity study investigating the effects of shrub expansion on snow distribution and energy balance, with and without topography.Shrub expansion can proceed by infilling and lateral growth of existing shrub patches, increase in the height of shrubs and colonization of areas beyond the shrubline (Tape et al., 2006;Myers-Smith et al., 2011).Without going to the complexity of introducing an ecological model, shrub expansion by the first two mechanisms was simulated by iteratively increasing the area and height of existing shrubs in the Granger valley for perturbed simulations.In each iteration, the vegetation fraction and height in each model grid box were increased by a random amount up to the maxima in any of the eight neighbouring boxes.This process was repeated 20 times, saving vegetation fraction and height maps after 1, 3, 5, 10 and 20 iterations; the vegetation fraction increases from 8 % in the lidar-derived map to 52 % after 20 iterations.DBSM and 3SOM runs were performed without vegetation, with existing vegetation, and with each of the five increased vegetation scenarios.Two runs were performed in each case: one with the existing topography and one on a flat domain.Model outputs were averaged over the central 1 km 2 of the domain.
Premelt conditions on 22 April are shown in Fig. 10.About one third of the vegetation is buried by snow in each case, with little impact of topography.With no vegetation and no topography, DBSM diagnoses sublimation of blowing snow but no net redistribution within the domain.As the vegetation fraction is increased for the flat domain, the domain-average premelt SWE increases because the reduction of near-ground wind speed by shrubs decreases blowing snow sublimation.In runs with topography, the increase in SWE with increasing F v is much less marked because deposition of snow also occurs in areas of decreased wind speed, such as hillslopes and depressions, and reduces the possible deposition in the tall shrubs of the valley bottom; Quinton et al. (2004)  premelt SWE is similar for runs with or without topography and is controlled by the shrubs.Spatial variability in premelt SWE can be characterized by the coefficient of variation (CV = standard deviation divided by mean).In simulations without topography, the snow cover is uniform when there is no vegetation (CV = 0 when F v = 0) but CV first increases when shrubs are introduced because they cause some spatial variability in snow accumulation; as shrub cover increases CV decreases because the shrubs suppress wind-induced redistribution of snow.Variability reaches a maximum without vegetation in simulations with topography and drops to CV = 0.29 at F v = 0.52.In comparison, Pomeroy et al. (2004) found CV = 0.27 from snow surveys for well-vegetated sites in WCRB.
Both the spatial variability in premelt SWE and average premelt SWE are important to capture: the former determines how much snow-free ground is exposed after a certain amount of snow has melted (Pomeroy et al., 1998a;Essery and Pomeroy, 2004a) and hence influences the surface energy partitioning, the latter determines how much energy is required to melt all of the snow.Time series of snow cover fraction and snow cover depletion curves (F s plotted against SWE) are shown for selected simulations in Fig. 11.With no topography and no vegetation, the premelt SWE and spatial distribution of melt energy are uniform and the transition from complete snow cover occurs in one time step.The simulation with topography but without vegetation has the highest premelt CV, giving a flattened snow cover depletion curve and a much more gradual decrease in snow cover; some ground is exposed early as shallow snow melts but some snow cover persists late into May in deep drifts.Increasing vegetation fraction increases premelt SWE and decreases CV, delaying the onset of snow cover depletion but increasing its rate once it has begun because of advected energy from exposed shrubs.
Figure 12 shows net radiation, sensible heat and latent heat fluxes averaged over early melt (22 April-4 May), main melt (5-16 May) and late melt (17-28 May) periods in simulations with and without topography.During the early melt period, the snow cover in simulations without topography is nearly complete (F s ≈ 1 in Fig. 10) and the sensible heat is dominated by downward fluxes from the atmosphere to the snow; this is offset by increasing upwards sensible heat fluxes from the vegetation with increasing vegetation fractions.In simulations with topography and low vegetation fractions, the snow cover is already incomplete (F s < 1) on 22 April, giving upwards sensible heat fluxes from snow-free ground that offset the downwards fluxes to snow and reduce the sensitivity of overall sensible heat to vegetation fraction.In simulations with and without topography, the increase in net radiation with increasing shrub cover is largely balanced by less-negative sensible heat fluxes.Latent heat fluxes increase slightly with increasing shrub cover due to advection of heat from exposed vegetation to snow within the same grid box.
The relationship between simulations with and without topography changes during the main melt period; differences in R n are constant independently of F v and sensible heat fluxes in simulations without topography are now almost independent of vegetation fractions less than 0.4.All simulations estimate less than 50 % snow cover by the end of this period, and those without topography have smaller snow cover fractions than those with topography and the same vegetation fractions because of lower initial SWE and more rapid snow depletion.Advection between grid boxes, if introduced in the model, would likely have the greatest influence during this period when there are significant fluxes from vegetation and snow-free ground but significant snow cover still remains.There are only small differences in fluxes averaged over the late melt period between simulations with and without topography because the surface is dominated by snow-free ground and vegetation.The difference in sensible heat fluxes between simulations with and without topography are due to late-lying topographical drifts, as indicated by longer persisting snow cover in Fig. 11.Sensible heat fluxes increase due to increasing surface roughness with increasing vegetation fractions but latent heat fluxes decrease slightly due to a decrease in the fraction of late-lying snow patches.Because the model lacks a hydrology module, partitioning of available energy between latent and sensible heat fluxes for snow-free ground is uncertain.

Discussion and conclusions
Corroborating previous findings (e.g.Chapin III et al., 2005;Liston and Hiemstra, 2011;Bonfils et al., 2012) this study suggests that expansion and densification of tundra shrub patches in a warming climate will have a positive feedback on warming through decreases in surface albedo and increases in sensible heat fluxes to the atmosphere.This change in surface energetics with warming is predicted despite the inclusion of a shrub bending parameterization which reduces the exposed vegetation fraction and increases the albedo at the beginning of the snowmelt season.However, topography was found to moderate the magnitude of the effects of shrub expansion on premelt energy budgets and snow accumulations; for the domain studied here, wind-blown snow from the exposed plateau can be trapped in a topographic drift on the northeast-facing slope before it reaches shrubs in the sheltered valley bottom.This suggests that the positive feedback identified in studies of level Arctic plains may in fact be dampened in Arctic hills and mountain location, such as in upland Alaska, the Yukon Territory and the Northwest Territories in Canada.
The high resolution of the grid used here allowed influences of topography and vegetation on snow accumulation and melt to be explicitly resolved and shows dampening feedbacks due to small-scale topography that are likely to upscale and affect biome-scale energy fluxes.These findings have implications for studies investigating shrub expansion over larger scales because of the number of processes which are poorly or not represented in land surface or climate models but which were found to affect the energy budget, namely (i) the effects of sub-grid topography on snow distribution, (ii) the effects of sub-grid redistribution of snow by wind, (iii) the simplistic sub-grid snow cover fraction parameterization as a function of SWE or snow depth that many models adopt, and (iv) the effects of shrubs on redistribution of snow which are neglected (except in CLM4 in Lawrence and Swenson, 2011).
The crucial role of spatial variability in snow on snowmelt is well known (Pomeroy et al., 1998a(Pomeroy et al., , 2004;;Clark et al., 2011;Egli et al., 2012) but its implications for shrub expansion studies have not been assessed until now.This study suggests that modelling shrub expansion without considering the effects of topography on wind redistribution of snow and on snowmelt rates may lead to overestimation of winter and early spring energy budget responses.Although there have been fewer investigations of high-latitude end-of-winter energy budget, many have shown large differences in net radiation and sensible heat fluxes between low albedo (generally trees or shrubs) and snow-covered surfaces despite limited solar radiation (e.g.Harding and Pomeroy, 1996;Chapin et al., 2000;Sturm et al., 2005a).Given the contrasting effect of topography between premelt and melt processes, further research is needed to understand the significance of shrub expansion in both complex and flat terrains during the whole snow season.
Partial snow cover was explicitly represented by the models at an 8 m resolution, a scale that is orders of magnitude smaller than climate model grid boxes; analysing model results in terms of 1 km areal averages identified the spatial variability of snow water equivalent as the most important factor affecting snowmelt energetics at the kilometre scale.Not only do the results presented in this study show basinscale behaviour in snow depletion and energetics due to very fine-scale processes but previous studies have argued for inclusion of more realistic snow cover depletion parameterizations (e.g.Pomeroy et al., 1998a;Roesch et al., 2001;Liston, 2004;Essery and Pomeroy, 2004b;Dornes et al., 2008;Clark et al., 2011), confirming that our findings should also be considered in larger-scale studies.
The model proposed here was specifically developed to investigate snow redistribution and snowmelt energetics in sparse canopies at high latitudes.3SOM addressed the need, expressed by Sturm et al. (2005b), Pomeroy et al. (2006) and Bewley et al. (2010), to account for the bending of shrubs under the snowpack in energy balance calculations by incorporating an exposed vegetation fraction parameterized from a shrub-bending model.The model also addressed the known limitation of dual-source models in reproducing snow melt rates for discontinuous shrub and snow cover (Liston, 2004;Essery et al., 2005;Bewley et al., 2010) by calculating separate energy balances for snow, bare ground and vegetation.The study was conducted at a single location in a mountain valley which is well understood because of the numerous research campaigns conducted at the site over the years; further studies are required to apply and confirm the relevance of these findings in other sub-Arctic and Arctic environments.In addition, further work should focus on year-round changes to the energy budget associated with shrub cover and topography.Further model developments, such as adding a hydrology module, accounting for heat advection between grid boxes and fully coupling 3SOM and DBSM, will be required to further improve our understanding of the surface and soil processes associated with shrub expansion.

Figure 1 .
Figure 1. 8 m resolution lidar-based digital elevation map of the 1 km × 1 km area around the valley in the Granger basin, overlaid by a vegetation fraction map and a contour map (1360 to 1480 m at 20 m interval).Four meteorological stations (from north to south: south-facing, tall shrubs, north-facing and short shrub stations) and the three snow depth transects are shown in black.

Figure 2 .
Figure 2. Structure of 3SOM with reference to the heat exchanges and the resistance network.

Figure 3 .
Figure 3. Measured (dots) and modelled (lines) (a) snow depth, (b) snow water equivalent, (c) sensible heat flux, (d) latent heat flux and (e) ground heat flux during snowmelt in 2003 and 2004.Snow depth and water equivalent measurements were manual (automatic snow depth at the short shrub site in 2003) and fluxes were measured using eddy correlation systems.

Figure 4 .
Figure 4. Modelled (a) sensible heat fluxes, (b) latent heat fluxes and (c) ground heat fluxes over individual sources (black = ground, red = snow, green = shrub) at the tall and short shrub sites in 2003 and 2004.

Figure 6 .
Figure 6.Modelled (black) and measured (red) incoming short-wave radiation at the two slope stations in spring 2003 (top panel: southwestfacing slope; bottom panel: northeast-facing slope).All the data above 710 W m −2 are unavailable on the southwest-facing slope because the programmed maximum data-logger voltage range was exceeded.

Figure 7 .
Figure 7. Modelled vs. measured average snow depths and standard deviation of snow depth on the northeast facing slope (NFS), the valley bottom and the southwest-facing slope (SFS) before, during and at the end of the melt season.Modelled standard deviation is extrapolated to 4 m resolution.

Figure 9 .
Figure 9. Area surrounding the two easternmost transects on the southwest-facing slope on 20 April, 27 April and 3 May 2008 from left to right.The visible portion of the transects is approximately 200 m long.The red arrow is pointing to a snow-free patch.

Figure 10 .
Figure 10.Exposed vegetation fraction, vegetation height, mean SWE and coefficient of variation in SWE on 22 April as functions of snow-free vegetation fraction for simulations with (dots) topography and without (crosses).

Figure 12 .
Figure 12.Area-averaged latent heat flux, sensible heat flux and net radiation as functions of snow-free vegetation fraction for simulations with topography (dots) and without (crosses) over three time periods.

Table 1 .
Specifications of the instruments most relevant to this study whose measurements were used to compile forcing meteorological data for DBSM and 3SOM.
melting and cold snow albedo decay time.Albedo is refreshed by snowfall, S f , such that, with increment α s such that δt is the length of the model time step in seconds, α s(min) is the minimum aged albedo and τ m and τ c are the Hydrol.Earth Syst.Sci., 18, 2375-2392, 2014 www.hydrol-earth-syst-sci.net/18/2375/2014/

Table 2 .
Bewley et al. (2010)ean square errors, bias and correlation coefficient of sensible, latent and ground heat fluxes at the tall and short shrub sites in 2003 and 2004.The ME, RMSE and r 2 of the two-source model described inBewley et al. (2010)are reproduced here in parentheses for convenience; the bias was calculated for this study.The bold numbers indicate higher quantitative values for 3SOM than the two-source model.ME and RMSE are in W m −2 .

Table 3 .
Coupled model bias, mean error, root mean square (rms) and rms errors normalized by the standard deviation of measured snow depths, for snow depths on the three dominant topographical units and at different stages of the melt season (Premelt = 21 April; Melt = 25 April, 28 April, 3 May and 7 May; End of Melt = 10 May, 16 May and 19 May.Only one of the three transects was surveyed on the last two dates).ME and RMSE are in metres.