Edinburgh Research Explorer Measurements and modelling of snowmelt and turbulent heat fluxes over shrub tundra

. Measurements of snowmelt and turbulent heat ﬂuxes were made during the snowmelt periods of two years at two neighbouring tundra sites in the Yukon, one in a sheltered location with tall shrubs exposed above deep snow and the other in an exposed location with dwarf shrubs covered by shallow snow. The snow was about twice as deep in the valley as on the plateau at the end of each winter and melted out about 10 days later. The site with buried vegetation showed a transition from air-to-surface heat transfers to surface-to-air heat transfers as bare ground became exposed during snowmelt, but there were daytime transfers of heat from the surface to the air at the site with exposed vegetation even while snow remained on the ground. A model calculating separate energy balances for snow and exposed vegetation, driven with meteorological data from the sites, is found to be able to reproduce these behaviours. Averaged over 30-day periods the model gives about 8 Wm − 2 more sensible heat ﬂux to the atmosphere for the valley site than for the plateau site. Sensitivity of simulated ﬂuxes to model parameters describing vegetation cover and density is investigated.


Introduction
Shrub tundra occupies the latitudes and altitudes above the coniferous forest treeline. Many species are present in shrub tundra, although willow (Salix) and birch (Betula) dominate; 40% of the Alaska-Yukon Arctic region comprises upland low shrub canopies of Betula nana (dwarf birch) and Salix pulchra (Jorgenson and Heiner, 2003), and these species are Correspondence to: R. Essery (richard.essery@ed.ac.uk) also present in Western Alaska and Northern Scandinavia. Tall shrub canopies occupy less areal extent but are prevalent in both upland and sheltered riverine areas, consisting mainly of willow and alder (Alnus) species such as Alnus crispa (Jorgenson and Heiner, 2003). Epstein et al. (2004) suggest that shrubs are sensitive indicators of climate warming and appear to be the most responsive of a wide range of tundra and treeline biomes considered; this pattern is reflected in repeat photography and satellite imagery suggesting recent expansions in Arctic shrub cover (Tape et al., 2006).
Tundra may be snow-covered for 6 to 8 months per year in sub-arctic regions  and up to 9 months in the Arctic . Snowmelt is the major annual hydrological event in these areas, providing between 40% and 90% of the annual flow of freshwater to streams, lakes and oceans at high latitudes. However, the processes controlling the rates and magnitude of snow ablation under vegetation canopies remain one of the greatest uncertainties in snowmelt calculations (Link and Marks, 1999;Koivusalo and Kokkonen, 2002;Gelfan et al., 2004).
Previous field and modelling studies have identified the effect of shrubs on snow accumulation, snowmelt, hydrology and surface energy balance. Shrubs have a strong control on snow accumulation through suppression of wind transport; for instance, McFadden et al. (2001) showed from measurements in Alaska that shrubby tussock tundra and water track margin areas in which shrub birch was common had 27% deeper snow than tussock tundra, and Pomeroy et al. (1997) found that shrub tundra north of Inuvik, NWT accumulated four to five times more snow than sparsely vegetated tundra. Essery and Pomeroy (2004) and Liston et al. (2002) used blowing snow models to simulate erosion of snow from sparse tundra and trapping by shrubs. The amount of snow accumulation around shrubs was found to be sensitive to Published by Copernicus Publications on behalf of the European Geosciences Union. 1332 D. Bewley et al.: Heat fluxes over shrub tundra both the extent and height of shrubs (Essery and Pomeroy, 2004). Sturm et al. (2001) also found that the deeper snow in shrubs has different structural and thermal characteristics from snow in open areas, such as a greater proportion of weakly-bonded depth hoar of lower thermal conductivity, and proposed a positive feedback loop in which larger and more abundant shrubs trap more snow and reduce sublimation losses, leading to deeper snow cover; this promotes higher winter soil temperatures and microbial activity which ultimately increases the plant-available nitrogen required for further shrub growth.
The effect of shrubs on snowmelt was identified by Price and Dunne (1976). Emergence of low-albedo shrubs from a receding snowpack increases the absorption of solar radiation, increasing the energy available for snowmelt. Pomeroy et al. (2003) found little difference in net radiation on north and south facing shrub tundra slopes until the earlier shrub emergence on the south facing slope, after which higher values of net radiation led to earlier snowpack removal by over a month. Hydrologically, McCartney et al. (2006) showed that tall shrub canopies exert an inordinately large control on the timing and magnitude of streamflow discharge in tundra basins because of their enhancement of snow accumulation and rapid meltwater production which overwhelms the infiltration capacity of soils. An increase in abundance of shrub tundra during climate warming may therefore impact the timing and magnitude of freshwater runoff to the major drainage basins.
Accurate surface energy balance calculations are required as boundary conditions in atmospheric models. Melt rates also determine the degree of vegetation exposure and fractional snow covered area, which in turn strongly influence albedo. Both Strack et al. (2003) and Viterbo and Betts (1999) found that numerical weather predictions were too cold (by as much as 6 • C) when the role of vegetation over snow was not represented correctly. Estimates of the increase in heating associated with the transition from shrubfree tundra to shrubland range from ∼10 Wm −2 for the growing season (e.g., Chapin et al., 2000Chapin et al., , 2005 up to 37 Wm −2 for the non-growing season (Sturm et al., 2005). Although energy-balance snowmelt models have been used to simulate shrub effects on surface energy and moisture fluxes during snowmelt, representations of important natural processes in such models are tentative and require evaluation by field measurements.
This paper presents measurements of snow accumulation, melt and turbulent fluxes of heat and moisture at two nearby sites: one with short vegetation that is completely buried by snow prior to melt and another with tall shrubs that remain partially exposed throughout the winter. The measurements are then used to evaluate the performance of an energy balance model in simulating melt rates and energy fluxes to the atmosphere for the two sites, and the model is used to investigate how fluxes depend on shrub canopy structure.

Site descriptions
The 195 km 2 Wolf Creek Research Basin (60 • 35 N, 135 • 11 W) near Whitehorse, Yukon Territory, Canada, lies in the mountainous headwaters of the Yukon River and has a sub-arctic continental climate characterized by a large seasonal variation in temperature, low relative humidity and relatively low precipitation. The mean annual temperature is −3 • C, and annual precipitation is 300-400 mm with approximately 40% falling as snow. Wolf Creek is situated within the Boreal Cordillera ecozone (Environment Canada, 1995) and consists of three principle ecosystems: boreal forest, shrub tundra and alpine tundra cover 22%, 58% and 20% respectively of the total basin area (Fig. 1). The shrub tundra zone is snow-covered from early October to May with snowdrifts persisting until June.
Measurements were conducted in the melt periods of 2003 and 2004 at two sites within the Granger Creek sub-basin ( Fig. 1), with automatic weather stations (Fig. 2) providing half-hourly meteorological data. One station was located at 1363 m a.s.l. in a valley bottom with a discontinuous willow canopy up to 2.5 m in height overlying a ground cover of grasses and mosses. The shrubs in this area retain a number of marcescent leaves and become exposed during the melt season both gradually as snow ablates and rapidly by springup of buried branches when there is insufficient weight of snow to maintain their burial. Many of the shorter shrubs have been observed to be completely buried prior to melt  in high snowfall years, but part of the canopy remains exposed. Another station was located 330 m from the valley site at 1424 m a.s.l. on an exposed plateau. Shrub cover here is sparse and dominated by birches of limited height (average 0.32 m) that become completely buried by snowfall and blowing snow.

Meteorological measurements
Incoming shortwave and longwave radiation, air temperature and wind speed measured on the plateau are shown for 30day periods beginning on day 105 (mid-April to mid-May) in Fig. 3 for 2003 and Fig. 4 for 2004. Both sites are nearly level and, despite differences in sky view, Sicart et al. (2006) found that differences in incoming radiation between the sites were generally not significant. Air temperatures at the sites are also very similar, although the valley site is often a couple of degrees colder at night, perhaps due to cold air drainage. The major difference in the meteorology between the sites is

Snow measurements
Snow depths were measured automatically at each site with sonic ranging sensors, giving point measurements that may not represent good areal estimates. In addition, manual measurements were made on several dates to obtain areal-mean snow depth, snow water equivalent (SWE) and snow-covered fraction. Near the valley site, snow depths were measured on a 61-point transect with 1 m spacing in 2003 and on a 7×7 grid with 5 m spacing in 2004. There were no distributed snow measurements in the immediate vicinity of the plateau site in 2003, but snow depths were measured there on a 25point transect with 5 m spacing in 2004. Snow density was measured at a selection of points by weighing cores in an ESC-30 snow tube, and snow-covered fractions were calculated as the fraction of the total number of points that were snow-covered for each survey. The observed relationships between SWE and snowcover fraction snow shown in Fig. 5 can be approximated by a snowcover depletion curve of the form which has been used in previous studies (Yang et al., 1997;Essery and Pomeroy, 2004) for snow depletion around short vegetation. Adjusting the parameter X to minimize mean squared differences between Eq. (1) and the observations, a value of 44 mm gives a reasonable fit to snow distributions both in the valley and on the plateau. SWE and snow depths for both sites and both years are shown in Figs. 6 to 9. The results show consistent differences in pre-melt snow accumulation between the sites, with more snow accumulating in the valley than on the plateau. In 2003, rapid snow ablation began on day 113 and continued uninterrupted until day 120, corresponding with a period of clear skies, high solar irradiance and warm air temperatures that exceeded 10 • C during each afternoon and were above freezing when averaged diurnally. The shallow snow around the plateau station all melted by about day 117. Between days 121 and 127, the sky was overcast and air temperatures only rose above freezing in the afternoons. Thereafter, diurnal mean air temperatures again rose above freezing, and this led to the remaining snow at the valley station melting by day 131. Averaged over the snow-covered period, the melt rates were 4.3 mm d −1 on the plateau and 7.3 mm d −1 in the valley; the latter value incorporates the week-long period of near-zero melt at the beginning of May.
There was more snow in 2004 than in 2003. Rapid ablation started in the valley on day 118 during a period of clear skies and rising air temperatures that persisted until the end of April. The melt rate on the plateau was initially slower than in the valley but increased to a similar rate by day 121 when sufficient shrubs were exposed; the snow at this site remained until day 124, giving a mean diurnal melt rate of 8.6 mm d −1 . The remaining snow at the valley site melted steadily through May; patterns of cloud cover were less distinct than in 2003, with no extended periods of clear or overcast sky. Air temperatures were cooler between days 123 and 130 resulting in slightly lower melt rates, but the warmer period thereafter was sufficient to melt all remaining snow at the valley site by about day 134. Because more snow melted in a shorter period than in 2003, the overall valley melt rate was higher in 2004 at 10.9 mm d −1 .

Flux measurements
Eddy covariance systems (Campbell Scientific CSAT3 sonic anemometers and KH 2 O or LiCor-7500 absorption hygrometers) were installed 3 m above the ground surface in the valley and 2.3 m on the plateau to provide measurements of sensible and latent heat fluxes between the surface and the atmosphere, calculated by finding 30-min covariances between vertical wind speed fluctuations and temperature or vapour density fluctuations sampled at 10 Hz. Results for both sites and both years are shown in Figs. 6 to 9, with the convention that fluxes to the atmosphere are positive. Coordinate rotations have not been applied, but the double rotation method recommended by Finnigan (2004) only gave small adjustments for a sample of the flux data. Use of the eddy covariance technique for measurement of fluxes over snow has been reviewed recently by Reba et al. (2009).
On the plateau, sensible heat fluxes were small and mostly negative while the ground remained completely covered with snow but increased rapidly upon the emergence of shrubs and bare ground. Positive sensible heat fluxes were observed even during the period of continuous snowcover in the valley, where relatively dense shrubs were always exposed for long fetches in the predominant along-valley wind direcions. In contrast, latent heat fluxes were small and, assuming that transpiration was minimal from the leafless shrubs, were entirely due to snow sublimation and evaporation of soil moisture. The snow on the plateau was relatively well  exposed to turbulent transfer, with daily mean sublimation rates of 0.42 mm d −1 in 2003 while snow was still present, after which the surface dried and latent heat fluxes decreased. In the valley, sublimation was suppressed beneath the dense canopy, and the daily mean sublimation was 0.38 mm d −1 .

Model description
To simultaneously calculate energy fluxes from an exposed shrub canopy to the atmosphere and energy fluxes to underlying snow that determine the melt rate, it is necessary to use a model that represents separate energy balances for the canopy and the snow. Several land surface models, such as CLASS (Verseghy et al., 1993) and CLM (Oleson et al., 2004), already have this type of structure but have not been widely evaluated in comparison with field observations of energy fluxes and snow melt. The dual-source model of Blyth et al. (1999) was modified for comparison with observations here.
Sensible heat fluxes are parameterised as being proportional to temperature differences and inversely proportional to aerodynamic resistances. The sensible heat flux between the surface and the canopy air space is the flux between the vegetation and the canopy air space is and the flux above the canopy is where ρ and c p are the density and heat capacity of air, T s is the snow surface temperature, T v is the canopy temperature, f v is the vegetation fraction, T c is the canopy air temperature, and T a is the air temperature measured above the canopy at height z t . Assuming that there is no transpiration from the shrubs, the latent heat flux is simply where L is the latent heat of sublimation, Q a is the measured specific humidity of the air and Q sat (T s ) is the saturation humidity at the snow surface temperature. Considering the uncertainty in parameterisations of turbulent transport through vegetation canopies, Blyth et al. (1999) argued for the use of simple logarithmic forms for the aerodynamic resistances in Eqs.
(2) to (5), although this is not physically correct within dense canopies. The resistances are for turbulent transfers between the surface and the canopy air space, between the vegetation and the canopy air space, and above the canopy, where U is the wind speed measured at height z u , k is the von Kármán constant and z 0s and z 0v are roughness lengths for the surface and the vegetation. Following Mason (1988), the effective momentum roughness length z 0 is calculated here as With the original form of Eq. (8), increasing vegetation cover was found to give increasing sublimation rather than suppressing turbulent transport from underlying snow, so a modified form defined by is used here, where the dense canopy exchange coefficient C is given the value 0.004 (Zeng et al., 2005). Setting f v equal to 0 or 1, Eqs. (6-10) reduce to the usual logarithmic forms for resistances over homogeneous surfaces. Net shortwave and longwave radiation absorbed by the vegetation are calculated as and LW where SW ↓ and LW ↓ are the measured incoming shortwave and longwave radiation fluxes, α v is the canopy albedo and σ is the Stefan Boltzmann constant. Neglecting multiple reflections with the canopy, the net shortwave radiation absorbed at the surface is where α s is the surface albedo and τ is the canopy transmissivity (multiple reflections amount to a renormalization of τ ), and the net longwave radiation is where v f is the sky view through canopy gaps. For partial snowcover, the surface albedo is calculated as where α snow is the snow albedo, α ground is the albedo of snow-free ground and f snow is given by Eq. (1). Energy conservation requires at the surface, where G is the surface heat flux, for the canopy, and for the canopy air space. Substituting Eqs.
(2) to (5) into Eqs. (16) to (18) and linearizing gives a system of three equations that can be solved for the three unknowns T s , T v and T c ; Blyth et al. (1999) gives the solution. As the model was originally developed for savannah vegetation, it had to be modified to include snowmelt. If the initial solution gives a snow surface temperature in excess of 0 • C, the temperature is fixed at 0 • C and the solution is repeated to find the snowmelt heat flux as the residual in Eq. (16).

Parameter values
Roughness lengths are set to typical values of 1 mm for snow, 1 cm for snow-free ground and a tenth of the canopy height for vegetation. Incoming and reflected shortwave radiation measurements gave values of 0.85 for the albedo of snow and 0.15 for the albedo of snow-free ground at the plateau site, and 0.11 for the albedo of shrubs at the valley site. To make some allowance for the influence of radiation absorption by buried branches and leaf litter, the snow albedo at the valley site is set to a slightly lower value of 0.8. Considering the mature state of the snowpacks and the short simulation periods, no snow albedo decay is included in the model. The parameters f v , v f and τ are often parameterised as functions of leaf area index, but are difficult to measure and are poorly defined for discontinuous canopies. From hemispherical photography, an aerial photograph and an array of 9 radiometers deployed under the canopy near the valley site, Bewley et al. (2007) estimated that f v increased from 0.20 to 0.71 as snow melted and shrubs emerged between days 112 and 129 in 2004, v f decreased from 0.77 to 0.59 and τ decreased from 0.69 to 0.57. Simulations were initially made with these parameters set to median values of f v =0.46, v f =0.68 and τ =0.63 for the valley site. Assuming that there is no exposed vegetation at the plateau site while there is snowcover, the canopy parameters are simply set to f v =0 and v f =τ =1 for plateau simulations.

Model evaluation
Model results for both sites and both years are compared with observations in Figs. 6 to 9. The model was initialized with the measured pre-melt SWE and snow density in each case. For the plateau site ( Figs. 6 and 8), the model captures the transition from small downwards sensible heat fluxes before melt to large upwards sensible heat fluxes once the snow has disappeared, with smaller and predominantly upwards daytime latent heat fluxes throughout the simulation periods. The onset of melt is well-represented, although the model removes the snow a day or two early. For three days during melt in each year the model predicts sensible heat fluxes of opposite sign from the observations. These are periods when the snowcover was observed to have become patchy; a model calculating simultaneous energy balances for snow and snow-free ground may be better able to represent this situation.
For the valley site (Figs. 7 and 9), the model matches the observations in producing downwards sensible heat fluxes at night and upwards fluxes during daytime even while snow remains on the ground. The onset and initial rates of melt are captured well; where the automatic and manual depth measurements differ, the model is generally closer to the manual measurements. In 2004, the date at which the model becomes snow-free matches the manual observations closely, but in 2003 the model continues to melt snow during the cold period following day 120 during which the observations suggest that melting halted. By this time, about one-third of the measurement points were snow-free. Excessive melt rates for patchy snowcover are a known feature of models that calculate a single energy balance for the composite snow and snow-free surface (Liston, 2004;Essery et al., 2005); a model with separate energy balances may, again, produce better results. Although some of the quantitative error statistics for flux simulations shown in Table 1 are poor, the qualitative match to observed fluxes at both sites was obtained without any model calibration.

Sensitivity to canopy parameters
To explore how fluxes depend on the structure of the canopy, the model was run using the valley meteorology for 2004 but with the vegetation fraction and canopy transmissivity varied independently between the limits 0 and 1. Sky view was set equal to transmissivity, giving longwave radiation transmission equal to shortwave radiation transmission; this is likely to be a better approximation for overcast periods with purely diffuse shortwave radiation rather than clear periods when the shrubs cast shadows (Bewley et al., 2007). Figure 10 shows contour plots of radiative, turbulent and melt energy fluxes averaged over the period of snowcover for varying canopy parameters (these parameters are not fully independent in reality; the lower left corners of the plots correspond with the unrealistic situation of a canopy that covers a small area but efficiently shades the surface). As formulated in Eq. (13), the shortwave radiation at the snow surface depends solely on τ , so the contours are horizontal. Longwave radiation at the surface, given by Eq. (14), largely depends on v f , but the contours slope due to the dependence on canopy temperature; as the vegetation fraction increases, the canopy temperature increases due to less cooling of the canopy air space by transfer of sensible heat to the snow, and the longwave radiation to the snow correspondingly increases slightly. For the high albedo of snow, the increase in longwave radiation dominates over the decrease in shortwave radiation (the contours are tighter in the longwave plot than in the shortwave plot), and the net radiation at the snow surface increases as the vegetation fraction increases or the canopy becomes denser ; snowmelt increases correspondingly. As the canopy is the dominant source of heat transfer to the atmosphere, the sensible heat flux largely depends on vegetation fraction and increases as f v increases; however, the sensible heat flux also increases gradually as the canopy becomes denser and less heat is transported to the snow beneath. Finally the latent heat flux, which is derived entirely from the surface, increases as the canopy density decreases and the radiative energy available for driving sublimation increases. As the vegetation fraction increases, the latent heat flux first increases as the transport of moisture from the canopy air space to the atmosphere increases but then decreases as the transport of moisture from the surface into the canopy air space is suppressed. The simple model used here does not include stability corrections to aerodynamic resistances, but Lee and Mahrt (2004) found snowmelt predictions to be sensitive to sub-canopy mixing parameterisations at a more southerly site with shorter shrub vegetation.

Conclusions
Two field campaigns involved measurements of snowmelt, near-surface meteorology and energy fluxes at two neighbouring sites through the melt period. The sites had different vegetation structures and topographic settings: one was in a valley with tall shrubs and deep snow, and the other was on a plateau with shallow snow but no exposure of the dwarf shrubs there prior to melt. Topography and vegetation together determine the pre-melt snow accumulation by controlling the redistribution of snow by wind (Essery and Pomeroy, 2004).
Snow surface temperatures cannot exceed 0 • C, so there is a downward temperature gradient above a snow surface when air temperatures rise above 0 • C during spring melt, giving downward sensible heat fluxes from the atmosphere that contribute to melting. Due to their low albedo, shrubs can be warmed by absorption of solar radiation to temperatures well above those of the air or the snow surface, giving upward sensible heat fluxes to the atmosphere even while there is melting snow on the ground. Exposed shrubs can also greatly increase the roughness of the surface, increasing turbulent fluxes. Increases in incoming longwave radiation beneath shrubs can outweigh decreases in shortwave radiation due to shading, giving greater net radiation at snow surfaces below shrubs than for exposed snow. Heating of exposed shrubs at the valley site led to upward sensible heat fluxes approaching 200 Wm −2 while snow remained on the ground. Conversely, the plateau snow surface was well exposed to turbulent transfer, and sublimation was the main mechanism by which net radiation was dissipated while the snowcover was continuous, but sensible heat flux rapidly increased upon the emergence of shrubs over melting snow and bare ground between snow patches.
The meteorological observations were used to drive a model that calculates separate energy balances for a snow surface and an exposed shrub canopy. The model was initialized with pre-melt snow observations and then evaluated in comparison with the observations of snowmelt and sensible and latent heat fluxes. Differences between fluxes before and after melt and between the sites were represented well, although the model showed some deficiencies in calculating fluxes and melt rates for the few days in each year during which the snowcover was patchy. It is suggested that additionally calculating separate energy balances for snow and snow-free ground and considering horizontal advection of energy would help to address these deficiencies.
The model parameters include vegetation fraction, canopy transmissivity and sub-canopy sky view that were estimated from field observations. These parameters vary with the extent to which shrubs are buried by snow and would generally have to be parameterised. One complication is that shrubs are often observed to be bent over and buried by a depth of snow that is less than their height when erect. The mechanisms involved are uncertain, but are likely to be related to the temperature-related properties of branch stiffness, snow cohesion and strength that control the interception and unloading of snow in forest canopies (Schmidt and Pomeroy, 1990;Hedstrom and Pomeroy, 1998). Further observations are required of how shrubs of different species trap and bend under snow.