Simulation of the snowmelt runoff contributing area in a small alpine basin

Simulation of areal snowmelt and snowcover depletion over time can be carried out by applying point-scale melt rate computations to distributions of snow water equivalent (SWE). In alpine basins, this can be done by considering these processes separately on individual slope units. However, differences in melt timing and rates arise at smaller spatial scales due to the variability in SWE and snowpack cold content, which affects the timing of melt initiation, depletion of the snowcover and spatial extent of the snowmelt runoff contributing area (SRCA). This study examined the effects of variability in SWE, internal energy and applied melt energy on melt rates and timing, and snowcover depletion in a small cold regions alpine basin over various scales ranging from point to basin. Melt rate computations were performed using a physically based energy balance snowmelt routine (Snobal) in the Cold Regions Hydrological Model (CRHM) and compared with measurements at 3 meteorological stations over a ridge within the basin. At the point scale, a negative association between daily melt rates and SWE was observed in the early melt period, with deeper snow requiring greater energy inputs to initiate melt. SWE distributions over the basin (stratified by slope) were measured using snow surveys and repeat LiDAR depth estimates, and used together with computed melt rates to simulate the areal snowcover depletion. Comparison with observations from georeferenced oblique photographs showed an improvement in simulated areal snowcover depletion curves when accounting for the variability in melt rate with depth of SWE in the early melt period. Finally, the SRCA was characterized as the product of the snowcovered area and the fraction of the SWE distribution undergoing active melt and producing an appreciable runoff quantity on each slope unit. Results for each slope Correspondence to: C. M. DeBeer (cmd225@mail.usask.ca) were then aggregated to give the basin scale SRCA. The SRCA is controlled by the variability of melt amongst slope units and over individual SWE distributions, the variability of SWE, and the resulting snowcover depletion patterns over the basin.


Introduction
The end-of-winter snowcover over open alpine terrain is characterized by considerable spatial variability in its snow water equivalent (SWE) depth, primarily as a result of snow redistribution by wind (Anderton et al., 2004;Elder et al., 1991;Winstral et al., 2002).Topographic depressions, areas of exposed alpine shrubs and vegetation, and leeward slopes tend to accumulate snow drifts up to several metres or more in depth over the winter, whilst adjacent and more exposed terrain is scoured by wind and may accumulate little or no snow.During the spring snowmelt period, a complex patchwork of snowcovered and snow-free areas emerges as a result of this heterogeneity and the spatial variation in applied melt energetics over the landscape.Both sources of variation control the timing and rate of areal snowcover depletion (SCD), which is important hydrologically since the fraction of snowcovered area defines the maximum extent of the snowmelt runoff contributing area (SRCA) within the basin.Thus, studies dealing with alpine runoff require sufficient understanding and representation of this variability and the SCD to realistically derive the snowmelt hydrograph.
Previous work has shown that snowmelt dynamics in complex terrain can be represented at intermediate spatial scales in models reasonably well by disaggregating the landscape into terrain units defined by features such as slope and aspect, and treating coupled mass and energy balances separately over these units (DeBeer and Pomeroy, 2009;Dornes et al., 2008a, b;Pomeroy et al., 2003).This avoids having to rely on fully distributed melt representations, whilst at the same time utilizing an objective means of landscape stratification rather than an arbitrary grid system that does not conform to the natural variability of the terrain and the processes involved.In this regard, the approach is useful for explicitly representing the variability in applied melt energetics between slopes as a result of differences in irradiation.Further, areal SCD can be represented within individual slope units by stochastically accounting for the pre-melt distribution of SWE and applying slope-based melt rate computations to the SWE distribution.DeBeer and Pomeroy (2009) demonstrated how this approach yields improved results in simulated SCD curves from the case of uniform applied melt energetics and a spatially lumped distribution of pre-melt SWE over the landscape.
DeBeer and Pomeroy (2009) noted, however, that problems with this approach may arise in situations where melt rates and timing are non-uniform within individual terrain units.For example, small scale statistical associations between SWE and melt rates may accelerate or prolong areal SCD.Faria et al. (2000) and Pomeroy et al. (2004) attributed small scale negative associations to the effects of vegetation.In areas with a cold, redistributed and highly non-uniform snowcover, such associations would also result from differences in the internal energy content of the snow.This may be particularly problematic in the early melt period as shallow snow becomes isothermal and begins to melt while adjacent deeper snow is still warming to 0 • C (Gray and O'Neill, 1974;Male and Gray, 1975).Essery and Pomeroy (2004) suggested integrating over the joint SWE and melt rate frequency distributions in the case of inhomogeneous melt to derive SCD curves, but noted that these integrals will generally be intractable.This is because both the frequency distributions of these variables and the association between them are not constant over time, and thus an analytical solution is not generally possible.
Other than studies using a fully distributed approach, there has been no attention given to the effects of small scale spatial variation of initial SWE and internal snowpack energy storage on the computation of melt rates over a SWE distribution.Physically-based (i.e.energy/mass balance) pointscale melt rate computations applied to these distributions have typically been based on a single initial snowpack state (depth, density, internal energy, etc.) and have tended to neglect differences in internal energy changes over the SWE distribution (Liston, 1999;Luce et al., 1999).Anderton et al. (2002) applied an energy balance snowmelt model in spatially distributed simulations of snowmelt over an alpine catchment, but considered the initial snowcover to be isothermal at the melting point throughout the catchment.Kuchment and Gelfan (1996) accounted for stochastic variations of SWE within sub-areas of a forest-steppe basin for computing basin averaged meltwater outflow.Their calculations were also based on a spatially uniform snowpack temperature of 0 • C. Similarly, most of the recent work dealing with spatial simulation of the snowcover in alpine terrain has considered the SRCA to be equal to the snowcovered area for estimating meltwater volume (e.g.Martinec et al., 1998).This is an oversimplification of reality for cold, redistributed snowcovers since not all areas will generate and release meltwater in the early stages of melt (Gray and O'Neill, 1974;Marsh and Pomeroy, 1996).Often, melt is first produced from the shallow snow, which may completely disappear before the areas with a deeper snowpack begin to produce melt (Male and Gray, 1975).
The purpose of this study is therefore to investigate the problems described above more thoroughly and to suggest a potential approach to resolve them in modelling applications.We first examine, at the point scale, how variability in initial SWE depth affects the computation of snowmelt rates and timing on variously orientated slopes over an alpine ridge within the Canadian Rocky Mountains.The effects of pointscale melt rate variation due to differences in internal energetics over a non-uniform snowcover are then examined with respect to the areal SCD and the temporal evolution of the SRCA for 2 different slopes in an adjacent cirque basin.We finally explore these features at the scale of the small headwater basin, which encompasses multiple slope units of different aspect, gradient, and elevation.The approach we develop and use in this study is meant to be applied over open, sparsely vegetated alpine slopes.This analysis provides an initial basis for predicting the magnitude and timing of the snowmelt hydrograph here and in similar environments.

Study area and field methods
This study was conducted at Fisera Ridge (∼2315 m) and the adjacent sparsely vegetated cirque basin within the Marmot Creek Research Basin, Kananaskis, Alberta, Canada (Fig. 1).Marmot Creek is located within the Front Ranges of the Canadian Rocky Mountains, where climatic conditions are dominated by continental air masses.Winters are long and cold, with an average temperature of -15 • C for the months of January through March for Fisera Ridge.The spring season is generally cool and wet, often producing late season snowfall events at high elevations.Average temperatures at Fisera Ridge for May and June (during the primary snowmelt period) are +2 and +8 • C, respectively.Historically, annual precipitation at Marmot Creek has been observed to average about 900 mm, increasing to over 1140 mm at treeline near the Mt.Allan cirque (Storr, 1967); recent observations corroborate these values.Roughly 60-75% of the annual precipitation falls as snow at Marmot Creek with the percentage increasing with elevation.
Continuous meteorological and snowpack condition observations were made for 3 consecutive years (2007, 08, 09) at a permanent meteorological station located on Fisera Ridge, and at 2 outlier stations on south-east (S-E) and and the incoming longwave component was measured with a CG3 pyrgeometer (expected accuracy of ±10% for daily sums).Air temperature and relative humidity were measured at each of the sites with a Vaisala HMP45C212 hygrothermometer housed inside a Gill radiation shield.Wind speed and direction were measured at the main station using an R.M. Young propeller anemometer (Model 05103-10), and wind speed at both of the outlier stations was measured with Met-One 014A 3-cup anemometers.Rainfall was measured at Fisera Ridge using a Campbell Scientific TB4-L tipping bucket rain gauge, and in August 2008 an Alter-shielded Geonor T-200B weighing precipitation gauge was installed in a sheltered site on the ridgetop near the main station to measure total precipitation.Measurements of total precipitation for the entire study period were also available from a similar Geonor gauge located at a mid-elevation (1843 m) forest clearing site within the Marmot Creek Basin roughly 2 km away.Snow depths were recorded using Campbell Scientific SR50 sonic ranging sensors at each of the stations (at the 2 outlier stations, the sensors were oriented perpendicu-lar to the sloping ground surface).To monitor the internal snowpack temperature, a series of fine-wire thermocouples was installed vertically at the 2 outlier stations, beginning at 2 cm above the ground surface and continuing up at 20 cm intervals.
A digital single lens reflex camera (Pentax model K110D) mounted to the main station on Fisera Ridge was used to take several pictures daily of the cirque.Photos taken at 12:00 p.m. local time each day were selected for analysis of the snowcovered area over the cirque, except in situations when low cloud cover or snowfall obscured the terrain.The images were projected orthogonally onto a 1 m resolution digital elevation model (DEM) of the cirque derived from airborne LiDAR (Light Detection And Ranging) using the Georeferencing Terrestrial Photography software tool (Corripio, 2004;see DeBeer and Pomeroy (2009) for a description of this procedure).From the corrected images, the snowcovered area fraction was determined for the cirque and for various slopes within it using a threshold procedure in ESRI® ArcMap TM 9.1 GIS.Snowcovered area was determined as the ratio of the number of 'snow' pixels to the total number of pixels over the relevant area (excluding pixels containing "no-data" values which were hidden from view of the camera).Snow surveys were repeated a number of times throughout the winter and spring along fixed transects over the ridge and within the cirque (Fig. 1).These surveys were done to characterize the spatial variability of SWE on different slopes over time.Snow depth was measured using a graduated aluminium rod every 3 m along each transect, and density measurements were taken roughly every fifth depth measurement where possible using an ESC-30 snow tube and scale.Several snow pits were located at various positions across the ridge and the within the cirque to examine the vertical snowpack structure and measure the density of the individual snow strata.Bulk snowpack density measured in the pits near the stations was used to convert SR-50 snow depth measurements into SWE measurements (Fig. 3).
In late March, 2008 a second LiDAR DEM was obtained over Marmot Creek for the purpose of examining spatial patterns of snow accumulation over the basin (Hopkinson et al., 2010).The DEM representing snow-free conditions was subtracted from the snowcovered DEM to generate a 1 m spatial resolution raster image of snow depth values.Comparison of point values of snow depth extracted from this raster image with values measured along the Fisera Ridge survey transect on 28 March 2009 showed good correspondence (R 2 =0.94; n=137; RMS error=0.17m) over the range of measured snow depths from 0 to 2.88 m, and the LiDAR depths reproduced the observed spatial pattern over the ridge very well.There were no apparent biases in error due to slope or aspect, but LiDAR-derived snow depth had greater error (up to ∼0.6 m) in locations with exposed alpine shrubs.Thus, the LiDAR snow depth raster appears to provide a useful tool for examining the spatial pattern and distribution of snow depth over most of the sparsely vegetated and open alpine cirque.

Snowmelt modelling and validation
Simulations of snowmelt rate and timing were carried out using the Cold Regions Hydrological Model (CRHM) platform (Pomeroy et al., 2007).CRHM is a flexible, object-oriented modelling system that can be used to develop, support, and apply dynamic hydrological process algorithms.These algorithms are applied over hydrological response units (HRU; i.e. homogeneous landscape units characterized by their terrain and surface properties), within which conditions and processes are represented by single sets of parameters, state variables, and energy and mass fluxes.Various component modules representing basin characteristics, observations, and hydrological processes are combined within CRHM to form an operational model of the system that has a level of complexity specified by the needs of the user.
Internal snowpack temperature and melt rates were computed using Snobal (snowmelt energy balance model; Marks et al., 1998Marks et al., , 1999Marks et al., , 2008)), which was incorporated as a module within CRHM.This module approximates the snowcover as being composed of 2 layers: a surface active layer of fixed maximum thickness, and a lower layer representing the remaining snowpack.The module solves for the temperature ( • C) and the specific mass (kg/m 2 ) or water equivalent depth per unit area (mm) of each layer for each timestep.The energy balance of the snowpack at a point is expressed as: where Q m is the energy available for snowmelt, Q * is the net radiation comprised of both short-and longwave components, Q H , Q E , and Q G are the sensible, latent, and ground heat fluxes, respectively, Q P is the energy added to the snowpack by precipitation, and U is the internal energy of the snowpack.The melt energy can be expressed as a depth of melt, m, by: where ρ is the density of the snow, h f is the latent heat of fusion (0.334 MJ kg −1 ), and β is the fraction of ice in snow (taken as 0.97).Melt is computed in either layer when the accumulated energy exceeds that required to bring the snowpack to 0 • C, at which point positive values of Q m result in snowmelt.A comprehensive description of the Snobal model is given by Marks et al. (1998Marks et al. ( , 1999Marks et al. ( , 2008)).
The model was run at the point scale at 15 min time intervals in CRHM using meteorological observations from the various stations as the external forcing variables.Observations of short and long-wave radiation components were only available from the main station, so these had to be corrected for slope, aspect, and sky exposure to be applied at the 2 outlier stations.For shortwave radiation this was done using the modules Global and Slope Qsi within CRHM, which use expressions from Garnier and Ohmura (1970) and List (1968) to determine clear sky direct and diffuse radiation on slopes.The ratio of measured shortwave radiation to the calculated theoretical clear sky direct and diffuse radiation on a horizontal plane is used to adjust the calculated clear sky shortwave  radiation value to the slope.Incoming longwave radiation was adjusted for differences in sky exposure between the main station and the outlier stations using a modified version of the parameterization given by Sicart et al. (2006).More detailed descriptions of these procedures and the formulae used are given by DeBeer and Pomeroy (2009).
The model was run for the late winter and spring periods of 2007 at the main station, and the same periods in 2008 and 09 at all stations.The first 2 years of observations were used to identify and select the most suitable parameter values, while observations from the final year were used to evaluate the model performance.This assessment was based the correspondence between simulated and measured snow depth, SWE, and internal snowpack temperature.
Parameter selection focused on 3 parameters related to snowpack properties represented by Snobal, and 4 parameters related to the decay function for snow albedo, α.The first 3 parameters are the roughness height of the snow surface, z 0 , the maximum active layer thickness of the snowpack, max z,s0 , and the maximum liquid water holding capacity of the snowpack, w c,max , (i.e. the irreducible water content).Roughness height was set as z 0 =2.0×10 −2 m.This value is not physically realistic for the actual roughness of the pure snow surface as measurements of z 0 over seasonal snowcovers generally indicate small values on the order of 10 −4 to 10 −3 m (Harding, 1986).However, this serves as an effective value and accounts for the fact that turbulent fluxes in this environment are affected by shear stresses created by local surface roughness as well as non-local processes (Helgason, 2009).Helgason (2009) found effective values of z 0 in the range of 2-5×10 −2 m using eddy covariance measurements in Marmot Creek Basin.The default value of the maximum active layer thickness in Snobal is 0.25 m, but a value of 0.1 m was used here.Marks et al. (2008) found that reducing the value of max z,s0 from 0.25 to 0.1 m improved simulations of sensible and latent heat fluxes under a pine canopy within the Fraser Experimental Forest.Further, the value of 0.1 m is more physically representative of the upper exchange layer of the snowpack within which turbulent energy exchange and short-wave radiation penetration occur.Maximum water holding capacity, which is defined as the ratio of the volume of liquid water to the volume of pore space within the snowpack, was set as 0.01 following Marks et al. (1998Marks et al. ( , 2008)).In the absence of ice layering, evidence suggests that it rarely exceeds 1% of the snowcover void space, although with ice layering in a wet, melting snowpack it can be as much as 5%.The value of 0.01 yielded good results, while higher values were found to delay simulated snowpack ablation by retaining too much liquid water.
The snow albedo was parameterized as an exponential decay during the melt period using the following expressions from Essery and Etchevers (2004).For each timestep with snowmelt, the albedo is updated according to: where t is the timestep length, τ is a time constant applied to melting snow, and α min is the minimum albedo that the decay declines asymptotically toward.For time steps with snowfall, the albedo is increased by: where S f is the snowfall amount during the timestep, and S min,α is the minimum snowfall amount required to refresh the albedo to α max .Values of τ =10 6 s, α min =0.3, S min,α =10 mm, and α max =0.85 were used in Eqs. ( 3) and ( 4).These values produced a reasonable correspondence between predicted and observed albedo at the main station for times when there was snowcover at this point.The observed albedo was determined on a daily basis using accumulated shortwave measurements from the upward, as well as a downward facing CM3 pyranometer.The value α min =0.3 allows the albedo to decline to unrealistically low values for pure snow, but this effectively represents an areal albedo that represents the broader surface supplying energy to the melting snow.
Figure 2 shows the observed air temperature and precipitation series at Fisera Ridge for the spring of 2008, and Fig. 3 provides an example of the model performance at the point scale for the S-E and N facing stations over this period.The model performs very well in terms of the timing and rate of snowmelt at both sites, although it has problems with simulation of the snow accumulation earlier in spring because it does not account for blowing snow redistribution.The model was capable of producing similarly reasonable results at these sites and the main ridgetop site for the other observation years.Simulations of internal snowpack temperature were somewhat problematic at times prior to the onset of melt, but this may be partly attributable to problems with the thermocouples and accurate measurements of the corresponding snowpack layers due to their position.However, the timing of the snowpack reaching an isothermal state of 0 • C at each site corresponded very well with the observations.Table 1 provides Nash-Sutcliffe model efficiency values (Nash and Sutcliffe, 1970) and RMS errors between measured and simulated snow depth.These provide a quantitative measure of the model performance and indicate that it is capable of producing reasonable results over multiple spring seasons and on slopes of different orientation.The success of the model for predicting snowmelt rates and timing under these various conditions at Fisera Ridge indicates that it can be used towards examining potential effects of differences in initial SWE here and elsewhere within the Mt.Allan cirque.

Simulated melt rate -SWE associations
To examine the sensitivity of snowmelt timing and melt rates to initial SWE and the associated differences in internal energy content, the model was run at the point scale for various initial depths of SWE.These simulations were carried out using the parameter values listed in the previous section and the observed 15 min meteorological conditions at Fisera Ridge and the 2 outlier stations for 2008.Outputs of SWE from the model each day at midnight were differenced to derive the simulated daily melt rates, which were then compared with SWE values at the beginning of each day.In each case, the The meteorological conditions during the spring in 2008 led to several distinct periods of melt interrupted by snowfall events and cold weather that refreshed the snowcover and temporarily delayed snowmelt (Figs. 2, 3).There were some early, but short-duration melting events in the month of April, when average air temperatures were above freezing for several days.Peak snow accumulation was around 12 May, while the main snowmelt period was during the remainder of May and June.This pattern during the melt period is characteristic of most springs in the alpine zone of Marmot Creek and common in the Rocky Mountain Front Ranges.
Figure 4 shows plots of simulated melt rates vs. SWE at the S-E and N facing sites for different times during the melt period.These results show that there are considerable differences in the rates and timing of snowmelt between the two locations on opposite sides of the ridge.For example, in the early melt period (26 April to 4 May), predicted SWE values in the range of ∼200-500 mm at the S-E facing site began melting in appreciable amounts nearly 1 week before the corresponding SWE depths at the N facing site.At all times in this period the simulated melt rates were greater in magnitude at the S-E facing location in comparison to the N facing site.Similar variation is evident at other times in the melt period, with the onset of melting conditions following snowfall events occurring earlier, and melt rates generally being of greater magnitude at the S-E facing site than at the N facing site.Differences in melt rates amongst individual slope units have been previously reported (DeBeer and Pomeroy, 2009;Dornes et al., 2008a;Pomeroy et al., 2003) and are explained primarily by differences in short-wave radiation receipt.
The results in Fig. 4 also show that differences in initial SWE have a significant effect on the rate and timing of simulated melt under the same applied melt energetics at each location.The effect is especially important early in the melt period as energy inputs to shallow snow are expended on melt while inputs to deeper snow go towards warming and ripening the pack.At the S-E facing location in the early melt phase of 2008, simulated shallow snow began melting nearly a week prior to deeper snow (i.e.>500 mm SWE), while minimal melt was predicted for the deepest snow.The effect was similar but even more pronounced at the N facing site during this time.A strong negative association between simulated melt rates and SWE is clear at both sites during this time.
The results for this early melt period were moderately sensitive to the maximum active layer thickness parameter in the model.We performed the same simulations using the default value of 0.25 m in Snobal to compare with these results.With other factors equal, differences of up to ∼2 mm/day were found.Greater melt rates were computed for the simulations using the 0.1 m maximum thickness, and the results were most sensitive for SWE values in the approximate range of 200-800 mm.Also, the results were more sensitive to max z,s0 at the S-E facing site than the N-facing site.As mentioned previously, the value of 0.1 m is more physically representative of the surface exchange layer of the snowpack, and yields the optimum results together with the other chosen parameter values in this study.
At later times in the snowmelt period, Fig. 4 shows the patterns to change from those observed earlier in the spring.A more pronounced phase of the snowmelt period began in mid-May following the peak accumulation (Figs. 2, 3), and within one or two days, all snow began actively melting and the negative melt rate-SWE association disappeared at both sites.This was due to the fact that all snow depths were either at, or near to isothermal conditions.A negative association became evident again for a brief time in the late melt period in June at both sites following a snowfall event.Again, the deeper snow would take greater energy inputs to bring it back to 0 • C following the short period of colder weather.Interestingly, at other times a slight positive association between melt rates and SWE was predicted by the model (e.g.18 May).This may be explained by the refreezing of liquid water within shallow snow overnight, which requires greater energy inputs than deeper snow (due to the greater thermal mass) to satisfy the latent heat requirement and restore melting conditions.
The variation in melt rates indicated by the model due to variable SWE depth alone is significant.Briefly, at times such as the early snowmelt period or following cold weather and snowfall events, it is of equal or greater magnitude as the variation between slopes resulting from differences in radiation receipt.The association between melt rates and SWE at these times tends to be non-linear and changes over time, making the effects of the association difficult to parameterize in terms of a constant linear covariance term (e.g.Pomeroy et al., 2004).This is an important result that has not been previously considered in most snowmelt modelling applications, but which has significant implications for simulating melt rates over a cold, highly redistributed snowcover.These implications are considered in the following section.

Areal SCD and Snowmelt Runoff Contributing Area
Simulated snowmelt rates were used to model the decline in snowcovered area and the temporal evolution of the snowmelt runoff contributing area (SRCA) over 2 opposing slopes within the adjacent cirque.Our approach is based on the theoretical lognormal frequency distribution, which can be expressed for SWE values in the following linear form: where SWE and CV are the mean and the coefficient of variation (standard deviation/mean) respectively, and K is the frequency factor for each particular SWE value calculated from the rank of observation (Chow, 1954).If the underlying distribution is lognormal, values of SWE plotted against K should approximate a straight line with a slope equal to the standard deviation of SWE and an intercept at K = 0 equal to SWE.The theoretical framework for determining K values from observed data and for deriving SCD curves is described by DeBeer and Pomeroy ( 2009), Faria et al. (2000), Pomeroy et al. (1998) and Shook (1995).Essentially, the value of K corresponding to SWE=0 (i.e.K min ) provides an index of the fraction of snowcovered area over an individual terrain unit, and by applying slope-based computed melt rates to the distribution, K min values can be continuously recalculated to generate the SCD curve.
Areal SCD was simulated using this framework on both a north-and a south-facing slope within the cirque basin (following DeBeer and Pomeroy (2009)).Measurements from the snow surveys and the LiDAR snow depth raster were used to define the parameters of the SWE distributions on these slopes (Table 2).Snowmelt modelling was carried out as described previously in this paper, with the short-wave and long-wave radiation inputs adjusted for the gradient, aspect, and sky view of the 2 slope units.To investigate the effects of the variation in computed melt rates due to differences in initial SWE on simulated SCD curves, daily melt rates were applied uniformly from simulations using a single initial SWE value in the mid-range of the distribution, and were also applied across the distribution using melt rates computed for different initial values of SWE depth, as in the previous section.In the latter case, SCD was simulated by discretizing the linear form of the SWE distribution (Eq.5) on each slope into a number of individual segments, which were defined by the initial SWE values used in the melt computations.As each of these points melted over time, the segments were readjusted between the points, and that representing the shallowest part of the distribution was used to determine K min , and thus, snowcovered area.The effects of snowfall events during the melt period, which refresh the snowcover until the new snow has melted, were handled using a rescaled depletion curve following Moore et al. (1999).
The SRCA on each slope was characterized as the product of the snowcovered area and the fraction of the remaining SWE distribution with an appreciable predicted quantity Hydrol.Earth Syst.Sci., 14, 1205Sci., 14, -1219Sci., 14, , 2010 www.hydrol-earth-syst-sci.net/14/1205/2010/ of daily runoff out of the base of the snowpack (i.e.nonzero SWE with runoff rates >∼5 mm/day for illustrative purposes).Simulated snowpack outflow accumulated over the day for each of the model runs were used to define runoff over the SWE distribution.To determine the fraction generating significant runoff, the relative proportion of the total SWE distribution within each of the initial segments of the linear representation of the distribution must be known.For the lognormal distribution, this can be found knowing the initial SWE and CV as follows: 1.For each value of initial SWE used to define the limits of these segments, determine the corresponding K value from Equation 5.The value of the frequency factor K provides an index of the exceedence probability of that particular value of SWE in the distribution.
2. To determine the exceedence probability, calculate the value of the transformed frequency factor K y (i.e.Equation (3) in DeBeer and Pomeroy ( 2009)) corresponding to each value of K. Theoretically, K y is calculated as: where s y is the standard deviation of the log-transformed data, which can be estimated from (Chow, 1954): The probability of K y being exceeded, P (K y ), is equivalent to the exceedence probability of SWE, and can be determined as: 3. The fraction of the total initial distribution represented by each segment can then be determined as the difference in exceedence probability between the lower and upper limits of each class.
At any time, the fraction of the remaining distribution represented by the individual segments in the plot of K vs. SWE can be found by multiplying the fraction of the total initial distribution represented between K values by the reciprocal of the snowcovered area.By this approach, the range of SRCA is limited in value between zero (in the case of no significant runoff computed for any SWE depth) and the fraction of snowcovered area (in the case of all SWE actively producing runoff).
Figure 5 shows the simulated and observed SCD curves, as well as the evolution of the SRCA on both the north and south facing slopes in the cirque during the early melt period of 2008.Table 3 provides RMS errors between the simulation results and observations on both slopes for this period.The results show an improvement in the simulated SCD curves when accounting for the variability of melt rates over different initial SWE depths, particularly on the south facing slope.This initial phase of SCD on the south facing slope was due mainly to the melt and disappearance of large areas of shallow snow, while adjacent deeper snow was still warming to 0 • C. Thus, the simulation using melt rates based on a single melt rate for average SWE poorly represented the decline of snowcovered area in this early snowmelt phase.On the north facing slope at this time, SCD was initially due to wind scouring of the snowcover (which is not accounted for in this approach).Melt of the shallow snow did not begin until 3 May, at which time the model based on SWE classes predicted a small decline in snowcovered area, corresponding to the observations.The simulation based on melt rates computed for a single initial SWE predicted a very minor decline in snowcovered area at this time.
The temporal evolution of the SRCA differs considerably between the north and south facing slopes during this early melt period as a result of differences in the SWE distributions and SCD curves, as well as the predicted runoff over these distributions.Simulated snowmelt runoff generation began on 26 April over the south facing slope, but the areal extent of the SRCA was limited since not all of the SWE distribution had began actively melting.This extent quickly grew larger as melt and runoff occurred over an increasing portion of the SWE distribution, which was characterized by a relatively shallow mean SWE.However, at the same time, the depletion of snowcovered area partially counteracted the effect of melt over an increasing fraction of the remaining distribution.The maximum extent of the SRCA occurred on 28 April (0.67) after which it declined due to further areal SCD.On the north facing slope, simulated snowmelt runoff generation did not begin until 3 May and reached its maximum extent on 5 May (0.74).At this time, the areal extent of the SRCA was larger than that over the south facing slope because of the limited areal SCD occurring on this slope.Subsequent snowfall occurring on 5 and 6 May (Figs. 2, 3) refreshed the snowcover and temporarily ended snowmelt runoff generation on both slopes.
Figure 6 shows the results of this modelling on both slopes carried out over the entire spring melt period of 2008.The simulations for the remainder of the spring, following the initial period of melt and SCD, were based on the peak SWE distributions around 12 May.The distribution parameters at this time were estimated from the previous spatial snow distributions measured with the repeat LiDAR and from surveys done on 16 and 17 May, just prior to any significant depletion of the snowcovered area.Simulated melting conditions resumed on 14 May (south facing slope) and 15 May (north facing slope), and the value of the SRCA rapidly approached the snowcovered area on both slopes as significant runoff generation was computed for all SWE depths.Over the course of the spring there were several more snowfall events, with similar effects on areal snowcover, melt rates, and extent of the SRCA.Both approaches (uniform applied melt and melt computed for different SWE depths) yielded similar SCD curves that corresponded reasonably well with the observations, with the exception of some errors that likely resulted from the rescaled depletion curve approach of Moore et al. (1999) and uncertainty in the SWE distributions later in spring.The similarity in these curves between the two approaches was due to the fact that the variation in melt rates across the distributions of SWE were less pronounced later in the snowmelt period.However, at several times this variation following cooler periods and new snowfall events did have a minor effect on the rescaled depletion curves by initially accelerating the areal SCD.Despite this, the improvements from including simulations of inhomogeneous melt over the entire snowmelt period in the spring were negligible (Table 3).This partly explains why the simulated SCD curves of DeBeer and Pomeroy (2009) were found in good agreement with observations later in the melt period using uniform slope-based melt rates.
This modelling framework was extended to all slopes within the upper Middle Creek basin to simulate areal SCD and SRCA at the basin scale.The basin was segregated into four slope-based HRUs, which include the north and south facing slopes.In addition, a broad east facing slope on Mt.Allan in the upper part of the basin was defined, and the small cirque floor area was distinguished from the other slope units.The floor is characterized by relatively level terrain with several small knolls and depressions and incised by several gulleys.Table 2 gives the statistical parameters of the SWE distributions on each HRU used to approximate the lognormal distribution in the late winter.The basin contains some forested areas in the lower region, where the tundra snow modelling framework in this study is not expected to apply.These areas were therefore excluded from the analysis.Further, there are a few local areas of steep cliffs in the upper part of the basin, which are permanently snow-free.To distinguish these areas and exclude them from the analysis, a terrain classification similar to that of Blöschl et al. (1991) was carried out to identify areas with a gradient steeper than 50 • .This value of the critical slope gave results that corresponded well with permanent snow-free areas identified in the corrected oblique photographs.
The simulation results on each slope unit were weighted according to the area of the slope and combined to give the basin scale SCD curve and SRCA, which are shown in Fig. 7.This basin scale SRCA is initially zero and grows in extent beginning 26 April due to runoff generation over the south facing slope.As melt and runoff generation begin to occur over an increasing portion of the distributions on the other slopes, the SRCA approaches the value of the basin scale rel-ative snowcovered area.Periods of new snowfall refresh the snowcover and reduce the snowmelt runoff generation area to zero until significant melt and runoff resumes following these events.As with the simulations on the individual slope units, improvements in the basin scale SCD curve resulting from including simulations of inhomogeneous melt are limited primarily to the early melt period in late April and early May (Fig. 7; Table 3).These results show that the basin scale SCD curve and SRCA are controlled by the relative extent of the various slope units and the changing snowcovered area and snowmelt dynamics over them with time.
As an example of the spatial variability of the SCD and SRCA within the basin, Fig. 8 shows the extent of snowcovered area and melting snow over several dates in the early melt period of 2008.These patterns were mapped from the LiDAR-derived snow depth raster image by accounting for the melt and depletion of different SWE depths over time on each of the slopes.Initially, simulated runoff generation began over limited parts of the south facing slope and the cirque floor.These areas then expanded and subsequently became depleted of snow as more parts of the other slopes began to produce melt.By the time the SRCA had expanded to include most of the area on the north facing and upper east facing slopes, the snowcovered area on the other slopes had declined considerably, thus limiting the runoff generating source area there.Thus, the differences in melt rates and SWE distributions amongst slope units controls not only the overall timing of the development of the SRCA, but the location of contributing areas as well.

Discussion and conclusions
This study has shown that differences in the warming and melt rates over a cold and highly redistributed snowcover have an important effect on areal SCD as well as the timing and development of the SRCA.The effects are most pronounced early in the melt period when only the shallow snow is actively melting.This leads to an initial acceleration of the SCD from the case of uniform applied melt rates based on a single initial SWE value in the mid-range of the distribution, due to earlier and faster melt of areas with a thin snowcover.As the melt period progresses and the effects of variable SWE and internal energy differences become reduced, melt rates over a given slope (and SWE distribution over the slope) become more uniform and consideration of melt differences amongst SWE classes becomes less important.
It is noted, however, that this variability can be important again at other times such as following spring snowfall events (which are common in this environment and refresh the snowcover) or periods of cold weather.Even night-time cooling and longwave radiation losses (particularly on clear nights) can have a considerable effect on internal snowpack energetics and become manifested as differences in melt rates amongst different SWE depth classes the following day.Male and Gray (1975) showed the importance of the internal energy term in applying the energy balance principle to shallow snowcovers.Shallow snow tends to undergo large diurnal variations in internal energy due to overnight cooling and refreezing, in contrast to deeper snow, which takes longer to warm in spring, but exhibits relatively damped diurnal variation of internal energy content due to its greater thermal mass.Fierz et al. (1997Fierz et al. ( , 2003) ) also showed the importance of explicitly representing internal processes within the snowcover in alpine terrain to simulate the evolution of initially subfreezing snowcovers and the effects of refreezing overnight.Thus, it is important to consider differences in melt rates due to variation in SWE throughout the entire melt period in such environments.
Realistic simulation of the SRCA depends on proper representation of the areal SCD, which limits the maximum contributing area, as well as the fraction of the SWE distribution that is actively producing runoff.Consideration of differences in melt timing and rates over a non-isothermal SWE distribution is therefore important in this regard.As described previously, earlier and initially faster melt of the shallow snow depletes the snowcover and exposes the ground more rapidly in these sites.Thus, the effect of the earlier melt of shallow areas of snow, which initially expands the SRCA, is somewhat counteracted by the fact that these areas become more rapidly snow-free and the runoff contributing area is then limited by the snowcovered area.At the basin scale, which generally encompasses multiple slope units of different aspect, gradient, and elevation, the pattern is even more complex.Snowmelt and areal SCD begin first over areas with shallow snow on south facing slopes or level lower elevation sites, and then proceed to deeper snow within these areas and shallow snow on other slopes.By the time that most of the snowcover on the remaining slopes begins to melt, much of the snowcover on the south facing slopes may have disappeared.Subsequent spring snowfalls add even more complexity to these patterns as the snowcover is reestablished over the landscape and the SRCA must again evolve.
The approach used here was to apply slope-based melt energetics to different SWE classes from a spatial frequency distribution of SWE over the different slope units in the basin.This yielded reliable prediction of areal SCD that corresponded well with observations, and it is reasonable to infer that the simulated SRCA was in good agreement with the actual areal extent and location of the snowcovered area producing appreciable melt quantities, based on the success of the model at the point scale and the known spatial distributions of SWE.The approach required far less data input and parameter information than a fully distributed approach, and yet was capable of simulating the complex spatial evolution of the snowcovered area and SRCA during the spring.This is because the processes and SWE distributions were represented at the appropriate spatial scales, areas and locations.This approach could easily be applied in other similar alpine basins with meteorological forcing adjusted for the slope/aspect and elevation, and information on the basic parameters describing the initial SWE distributions on the slopes (i.e.SWE and CV for the lognormal distribution).
The approach is limited and may not be suitable in some situations.For example, it is not expected to work where other factors have a significant influence on the association between melt rate and SWE, unless those factors can be accounted for by the model and related to SWE.In this study, we neglected densely vegetated parts of the basin where a more complex relation between melt rates and SWE may exist (Faria et al., 2000;Pomeroy et al., 2001).When applying this framework over mountainous terrain it is necessary to account for cliffs, which are subject to snow avalanches and do not retain snowcover (Blöschl et al., 1991).The lognormal distribution generally provides a useful approximation to the SWE distribution over the remainder of the slope.However, if the terrain is highly complex with extensive cliffs comprising much of the slope unit area, then this approximation may fail.Further measurements are needed to determine whether and how the presence of cliffs over larger parts of the landscape might affect SWE distributions.Similarly, where substantial redistribution of snow by avalanching occurs, it may also be difficult to describe the frequency distribution of SWE in a convenient parametric form such as the lognormal distribution.
The findings and developments in this work can be used towards an improved framework for simulating the spatial variability of meltwater generation and release from the snowpack in open alpine terrain.Realistic, physically-based representation of these processes must account for the spatial variability in melt rates and SWE at multiple scales (i.e.amongst individual slopes and within them).This leads to a further issue that should then be considered -the retention of meltwater within the snowpack and the transmission of meltwater fluxes through the pack.The rate of propagation of meltwater fluxes through snow is related to the magnitude of the flux (Colbeck, 1972;Colbeck and Davidson, 1973), while the travel time through the snowpack is a function of this rate as well as the depth of snow.This travel time can become important for some of the deepest snow found in this study site.Thus, the framework described in this study for simulating the areal SCD and SRCA should eventually be coupled with a variable meltwater routing scheme through snow that accounts for these principles.This will allow for a more realistic simulation of the meltwater release over the landscape, which can then be used as the input to hydrological models for streamflow simulation.

Fig. 1 .
Fig. 1.Shaded topographic map of Fisera Ridge and the Upper Middle Creek Basin within the Marmot Creek Research Basin, Alberta.Snow survey transects in 2008 are represented by the white dashed lines and location of the detailed Fisera Ridge inset map is indicated.Other insets include an aerial photograph of the basin and a map showing the location of Marmot Creek within western Canada.

Fig. 2 .
Fig. 2. Air temperature and precipitation at Fisera Ridge during the late winter and spring periods of 2008.

Fig. 3 .
Fig. 3. Comparison of measured and simulated snow depth, SWE, and internal snowpack temperature at (a) south-east facing station, and (b) north facing station during the spring period in 2008.

Table 1 .
Nash-Sutcliffe model efficiency values and RMS errors for snow depth at each site for the period 1 April until final snow disappearance during each of the three simulation years.Nash-Sutcliffe efficiency RMS error S-E facing N facing ridge-top S-E facing N facing ridge-

Fig. 4 .
Fig. 4. Plots of simulated melt rates vs. SWE on consecutive dates at various times throughout the melt period in 2008 at (a) the south-east facing station and (b) the north facing station on Fisera Ridge.

Fig. 5 .
Fig. 5. Simulated vs. observed SCD curves and the SRCA during the early melt period within the Mt.Allan cirque in 2008: (a) Southfacing slope unit; (b) North-facing slope unit.SCD curves were derived based on inhomogeneous melt (In) from computations based on different initial SWE values, and from uniform slope-based melt (U) based on a single value of SWE and applied across the SWE distribution.Observed SCD curves were derived from the orthogonally corrected photographs of the cirque.

Fig. 6 .
Fig. 6.Simulated vs. observed SCD curves and the SRCA over the entire spring period within the Mt.Allan cirque in 2008: (a) South-facing slope unit; (b) North-facing slope unit.SCD curves were derived based on inhomogeneous melt (In), and from uniform slope-based melt (U) as in Fig. 5.

Fig. 7 .
Fig. 7. Simulated basin-scale SCD curve and SRCA over the Upper Middle Creek Basin for the spring period of 2008.SCD curves were derived based on inhomogeneous melt (In), and from uniform slope-based melt (U) as in Figs. 5 and 6.

Fig. 8 .
Fig. 8. Spatial snowcover and SRCA patterns over the Upper Middle Creek Basin (showing the four different slope-based HRUs) in the early snowmelt period of 2008.Patterns were mapped based on the spatial distribution of snow depth and SWE derived from the repeat LiDAR data.Snow melt patterns were not simulated in the densely forested or cliff areas.

Table 2 .
Statistical parameters of pre-melt SWE distributions over individual slope units in 2008 based on snow surveys and the Li-DAR snow depth raster.March and the simulation ended once the snow disappeared.Thus, melt rates associated with shallow SWE later in the melt period were based on the remaining snowpack from simulations with greater initial SWE values, rather than initializing the model with shallow SWE at later times in the melt period.

Table 3 .
RMS errors between simulation results and measurements of fractional snowcovered area on different slopes for both early melt phase (25 April-5 May) and the entire spring melt period.Simulations were carried out using both the uniform and inhomogeneous melt approaches.