Spatial Snow Depth Assessment Using LiDAR Transect Samples and Public GIS Data Layers in the Elbow River Watershed, Alberta

Abstract

This study illustrates the potential to combine LiDAR remote sensing and GIS techniques for the purpose of estimating instantaneous winter snowpack volume within the mountainous Elbow River Watershed (ERW) upstream of Calgary, Alberta. Two LiDAR (Light Detection and Ranging) datasets, one during snow-free and the other during late winter were used to evaluate a procedure for snow depth sampling. These data were also used to classify terrain and canopy cover attributes to enable snow depth estimation in areas that were not directly sampled but for which equivalent land classifications could be derived via other means. The mean snow depth from 1675 field measurements collected coincident with the winter LiDAR survey (late March, 2008) in snow-covered areas only was 0.28 m ( = 0.27 m). The mean LiDAR-based snow depth in snow-covered areas was comparable with the field values at 0.26 m ( = 1.2 m), or 0.18 m when averaged across both snow-covered and snow-free areas. Using field measurements of snow density, a GIS routine was employed to estimate total watershed snow water equivalent (SWE) from ten snow accumulation units (SAUs) using elevation, aspect and canopy cover. The total watershed SWE estimate was 46.0 106 m3. This volume of water can also be expressed as 0.058 m of water depth across the entire basin, or approximately 18% of the total 2008 runoff yield. Further work is needed to improve LiDAR-based snow depth estimation in areas of shallow snowpack where the influence of noise in the data is highest and to optimize the methods of sampling and extrapolation. At the present level of airborne LiDAR sophistication, positional uncertainties in LiDAR data (though small) are such that high confidence in the watershed snowpack volume estimate, would only be achieved during deep snowpack years; which also tend to be the years where accurate data are least required. However, given the availability of LiDAR base maps is ever growing, and the accuracy and costs associated with the technology are constantly improving, this approach to snow depth sampling has the potential to become a useful tool to support headwater snowpack resource assessment in water-stressed regions of Canada.

La prsente tude dmontre la possibilit d'utiliser la tldtection par LiDAR en combinaison avec des technologies SIG pour estimer le volume instantan du manteau neigeux en hiver dans la zone montagneuse d'Elbow River Watershed (bassin versant d'Elbow River) situe en amont de Calgary, Alberta. Deux ensembles de donnes LiDAR (dtection et tlmtrie par ondes lumineuses), l'un collect en priode sans neige et l'autre vers la fin de l'hiver, ont servi valuer une procdure pour chantilloner lpaisseur de la couche de neige. Ces mmes donnes ont aussi servi la catgorisation des attributs du terrain et du couvert forstier et, partant, ont permis dstimer lpaisseur de la couche de neige dans des zones qui n'avaient pas t directement chantillones. La profondeur moyenne des 1675 mesures effectues sur le terrain en troite concordance avec le lev LiDAR d'hiver (fin Mars 2008) dans les seules zones enneiges tait de 0,28 m ( = 0,27 m). La profondeur moyenne de la couche de neige mesure base de LiDAR dans les zones enneiges tait comparable aux valeurs constates sur le terrain, savoir, 0,26 m ( = 1,2 m) soit 0,18 m toutes zones, enneiges et sans neige, confondues. Une routine SIG utilisant des mesures de la densit de la neige faites sur le terrain, a permis d'estimer lquivalent total en eau de la neige (EEN) avec altitude, orientation et couvert de la surface du sol. L'estimation du total EEN tait de 46,0 106 m3. Ce volume d'eau peut aussi bien tre exprim comme tant gal 0,058 m d'eau sur la surface entire du bassin, ou bien comme 18% du ruissellement en 2008. Un travail supplmentaire est indiqu pour amliorer les estimations base de LiDAR de lpaisseur de la neige dans des zones de manteau de neige peu pais o l'influence de bruit est la plus marque, et pour optimiser les mthodes dchantillonage et extrapolation. Au niveau actuel de sophistication du LiDAR aroport, des incertitudes de position dans les donnes LiDAR (quoique faibles) sont telles qu'un haut degr de confiance en l'estimation du volume de neige accumule dans le bassin ne serait ralis que pour les annes de neige profonde alors que, de faon gnrale, ce sont ces annes-l o le besoin de donnes prcises est moindre. Vu que les cartes ralises base du LiDAR sont de plus en plus disponibles et que les cots et la prcision des donnes associs cette technologie vont toujours s'amliorant, cette mthode dchantilloner lpaisseur de la neige est d'une utilit potentiellement trs souhaitable pour appuyer lvaluation des ressources en neige accumule du cours suprieur d'une rivire dans les zones de stress hydrique au Canada.

Introduction

In Southern Alberta's Bow River basin (BRB) (~26,000 km²), most of the runoff originates as snowpack in the mountainous headwaters of the Canadian Rockies. The importance of water, and therefore snow, in this region where supply appears insufficient for projected needs is illustrated in the recent cessation of new water allocations (Province of Alberta, 2007; Pentney and Ohrn, 2008). Snow accumulations in the Bow River headwaters are currently monitored by Alberta Environment at 20 snow course stations to provide information on summertime water availability and the potential for spring flooding. This number of sampling locations is insufficient to provide an accurate spatial map of snow depth but long term relationships between headwater snow course and downstream runoff observations can be generated such that contemporary snow course data provide an index of the following season's water availability. Such indices are reliable in so far as the climate or land surface characteristics are stable.

Global climate model (GCM) predictions for the mid 21st century suggest a more active precipitation regime, higher temperatures and a potential reduction in flow for the BRB (Martz et al., 2007). Furthermore, changes in land surface attributes, particularly forest canopy structure, due to insect and disease (Coops et al., 2010), wildfire and land use, have the potential to influence both snow accumulation (Boon, 2007) and melt timing (Musselman et al., 2008; Winkler and Boon, 2010), thus altering any long term relationship. Given the likelihood for continued environmental and developmental changes in parts of the BRB (and other headwater supply regions), there is some uncertainty regarding the adequacy of snow course networks to provide reliable indices of future water availability. Consequently, there is a need to explore alternative snowpack monitoring methods, such as using remote sensing techniques (e.g., Derksen et al., 2005), that can more directly quantify total headwater snow pack accumulations.

Previous studies have demonstrated that snowpack depth variation can be assessed at the meso-scale with airborne LiDAR (Light Detection and Ranging) (Hopkinson et al., 2004; Deems et al., 2006; Fassnacht and Deems, 2006; Minoru and Hiroshi, 2006; Trujillo et al., 2007). In a study conducted over the Marmot Creek Watershed in the headwaters of the BRB it was found that LiDAR estimates of snow depth in alpine, forested slopes and valley locations demonstrated mean depths within 0.13 m of corresponding field data (Hopkinson et al., 2011). Alpine slopes demonstrated the highest accuracy, presumably due to reduced system error propagation (Goulden and Hopkinson, 2010), while forest-covered slopes demonstrated the highest uncertainty, likely due to signal interference by the overlying canopy and understory vegetation. Furthermore, the LiDAR snow depth model (LSDM) clearly illustrated that the watershed hypsometric mean snow depth reached its maxima at treeline around 2250 m a.s.l. (Hopkinson et al., 2011). [For an in depth introduction to the basics of airborne LiDAR technology, the reader is referred to Baltsavias (1999) and Wehr and Lohr (1999).]

Given the growing database of high resolution and accuracy LiDAR topographic base map coverage within Alberta (e.g., Airborne Imaging, 2010), the City of Calgary and the provincial government funded a study to develop a LiDAR snow depth sampling strategy within a headwater region of the BRB. The criteria for the study were to devise a minimal cost method of sampling that did not require repeat mapping of the entire watershed, while allowing the use of readily available data layers for extrapolation of snow depth to the remaining area of the watershed. The approach used was to sample the watershed with systematic LiDAR transects and aggregate snow depth observations into discrete snow accumulation unit (SAU) classes based on dominant and controlling land surface attributes at the watershed scale. If successful, such an approach is seen as a potential supplement to snow course monitoring as it reduces the reliance on empirical models based on a few point measures and should thus facilitate a more explicit estimate of basin wide snow accumulation. The analysis is accompanied by a discussion of some of the uncertainties and limitations identified during the study.

Study Area

The ideal location to showcase LiDAR snow depth sampling would be an area of typically deep and widespread snow accumulation in the upper reaches of the BRB such as exist upstream of Banff or Lake Louise in Banff National Park. However, for this study the Elbow River Watershed (ERW) ( Figure 1 ) was chosen for a number of strategic reasons: 1) The ERW (1210 km²) drains into the Glenmore Reservoir (3.8 km²), which supplies the City of Calgary (~1.1 million people) with ~24% of its drinking water (Pernitsky and Guy, 2008). The reservoir also acts to buffer spring flood waters and provide an important recreational capacity to the people of Calgary; 2) Unlike the protected National Park setting of the Upper Bow, the ERW experiences forestry operations at intermediate elevations and agricultural land uses in the lower reaches; 3) The government of Alberta Sustainable Resource Development department (SRD) monitor and inventory land use and forest cover within the ERW; 4) Provincially owned LiDAR base map coverage from a snow free period in 2006 was already available for approximately 40% of the ERW, whereas only a fraction of the Upper Bow was available from previous research-based data collections.

Figure 1. Elbow River Watershed (ERW) study area in Alberta showing field and LiDAR sampling locations. ERW area background illustrated as terrain shaded relief.

Characteristic of the BRB in general, spring snowmelt from the ERW typically contributes the highest sustained period of inflow to the Glenmore Reservoir. Observed and modeled changes in the climatic regime over the ERW indicate increasing winter temperatures and a shift in the precipitation regime to more rainfall and less snow in the lower reaches (Valeo et al., 2007). The annually variable timing and magnitude of the ERW spring freshet requires that water levels in the Glenmore Reservoir be modified to ensure sufficient supply for prolonged periods while also holding back potential flood waters during spring melt and heavy rainfall runoff events (Valeo et al., 2007). Data quantifying the snowpack resource within the ERW is therefore valuable as it can inform water level management decisions to optimize surface water supplies at minimal hazard to the City of Calgary.

Snow precipitation and depth is monitored in and around the ERW by Environment Canada and the Government of Alberta. A search of the Environment Canada online climate data and information archive (climate.weatheroffice.gc.ca; last accessed February 2011) indicates that while there have been up to five meteorological stations simultaneously collecting data within ERW, there is currently no available data. However, the database does indicate that within 50 km of Bragg Creek, there are eight active stations, with at least four of these collecting snow data. Alberta Environment currently operate one snow pillow and collect monthly snow course measurements during winter at the Little Elbow snow monitoring station. This station is located in the westernmost upper reaches of the watershed at 2225 m a.s.l. in an area of relatively high snowpack accumulation. Given most of the snowpack in the ERW is found in the mountainous headwaters of the basin, the LiDAR sampling study was carried out across the 790 km² area upstream of Bragg Creek hydrometric gauging site 05BJ004 (Water Survey of Canada) ( Figure 1 ).

This paper reports on the field and LiDAR sampling strategy and the GIS methodology adopted to estimate snow depths within areas of the watershed that were outside the LiDAR sampling transects. It is not the intent of this paper to discuss in detail the process or results of LiDAR snow depth measurement and validation in a mountainous environment, as such analyses have been presented elsewhere (Deems et al., 2006; Fassnacht and Deems, 2006; Trujillo et al., 2007; Hopkinson et al., 2011). To support this project, however, complimentary LiDAR snow depth and land surface type SAU analyses conducted immediately west of the ERW along the slopes of the Kananaskis River Valley (KRV), which includes Marmot Creek watershed, are presented.

LiDAR Snow Depth Mapping

The Challenges

The common approach used to create a LiDAR-based snow depth model (LSDM) is to generate a digital elevation model (DEM) of the snow free ground surface and a digital surface model (DSM) of the snow covered ground and assess the difference (Hopkinson et al., 2004):

Clearly, two LiDAR data collections must take place to survey each of these surfaces. This leads to two self evident requirements: 1) Either both data collections need to be planned ahead of time, or some baseline dataset describing the snow free terrain surface must already be available at the time of LiDAR snow depth sampling; 2) The two acquisitions must be at least a few months apart to capture a full winter of snow accumulation but as the intervening time period increases, so too does the possibility for terrain surface and landcover modification. As of the time of the study, the LiDAR base map coverage in Alberta did not extend far into the eastern front ranges of the Canadian Rockies. For the ERW above Bragg Creek there was approximately 40% base LiDAR coverage, which was limited to the northern half of the watershed and the lower 1000 m of relief ( Figure 1 ). The southern half and the upper 400 m of relief (or the upper 11% of basin hypsometry) are not represented. Within meso-scale mountainous watersheds, snow depth observations and simulations can vary widely as a result of different landcover and terrain features exerting variable levels of control (e.g., Elder et al., 1998). Therefore, if our sample set has no representation for the upper 400 m of the watershed, this constitutes a serious limitation. While it was not possible to directly represent this part of the ERW using publicly available LiDAR data, we were fortunate to have access to a research-based LiDAR dataset (DeBeer and Pomeroy, 2010; Hopkinson et al., 2011) collected over the Kananaskis River Valley (KRV) and adjacent slopes immediately to the west ( Figure 1 ). The range in landcover and elevation of the KRV survey encompasses that of the upper westernmost headwaters of the ERW and thus provides a useful proxy.

In addition to spatial sampling challenges, collecting, processing and then comparing high resolution LiDAR surface models has many opportunities for error propagation (Hodgson et al., 2005; Deems and Painter, 2006; Goulden and Hopkinson, 2010). Therefore for each LiDAR surface there is a need to check for and reduce systematic positional bias and ensure comparable data resolutions prior to subtraction. While contemporary airborne LiDAR data accuracies are frequently quoted to be <15 cm (Optech Incorporated, 2005), internal variability within individual datasets can occur and so the process of spatial and vertical co-registration over immobile features can be necessary. Moreover, even a highly accurate LiDAR DEM or DSM can contain some spatial bias within the noise level of the data such that when compared to another dataset, this bias becomes a component of the change. It is important that DSMs are created and compared at equivalent resolutions, so that feature scaling influences do not propagate into the change surface. A challenge of particular relevance in mountainous terrain is that of horizontal positional uncertainty propagating a vertical error over steep terrain surfaces (Hodgson and Bresnehan, 2004; Hodgson et al., 2005). In a study of glacier surface downwasting in the Canadian Rockies it was reported that vertical errors in the change surface of <100 m were observed at cliff edges due to a horizontal shift of less than one DEM grid node (Hopkinson and Demuth, 2006). Completely removing such effects over steep terrain and cliff edges is challenging and can require manual manipulation of the DEM or DSM. Therefore, to accurately reproduce snow depth in mountainous environments, the above sources of error must be anticipated and addressed. If they are not, then the resultant LSDM change surface is likely to contain areas of systematic error that reflect properties of the underlying terrain and other uncertainties in the data.

LiDAR Data Preparation

Recognizing that airborne LiDAR monitoring has the potential to be costly over large areas, a sampling strategy was devised to minimize air time while representing a range of terrain and landcover attributes within the ERW. Two LiDAR datasets were required to perform the analysis; the first was collected as part of a Provincial base mapping initiative during snow-free and green foliage conditions in September 2006; while the second was commissioned specifically for this study during anticipated deep watershed snow accumulation. Both surveys were flown at an altitude of 3500 m a.s.l. using the same Airborne Laser Terrain Mapper (ALTM) 3100 sensor (Optech; Toronto, Ontario) owned and operated by Airborne Imaging Inc. based in Calgary, AB. In both cases, the pulse repetition frequency used was 33 kHz, and the average point spacing at ground level was between 1 m and 2 m, with actual point density increasing and swath width decreasing with terrain elevation. Sensor calibration and validation was performed before and after each flight at the Springbank Airport runway 20 km north-east of Bragg Creek and resulted in a vertical R.M.S. error less than 0.1 m.

The LiDAR DSM acquisition occurred on March 28^th, 2008. Flight lines were flown in east to west transects across the watershed and separated by approximately 5 km for optimal sample coverage with minimal flying time ( Figure 1 ). Given the snow-free LiDAR coverage in the ERW was incomplete and the highest terrain elevations unrepresented, snow-free and snow covered data collected from the KRV, located immediately to the west, were also used in this study. These additional datasets were flown using equivalent survey settings and ALTM sensors. The first snow-free dataset was flown in summer 2007 by the Canadian Consortium for LiDAR Environmental Research (C-CLEAR), while the second was flown by Airborne Imaging on March 29^th, 2008. The only difference in survey configurations was that the KRV surveys had complete spatial coverage and did not follow a transect sampling pattern ( Figure 1 ).

Before snow depth could be mapped, the LiDAR point cloud data underwent the following pre-processing steps to reduce uncertainty in the final LSDM:

1.	Ground classification of the LiDAR point clouds to remove vegetation and buildings from the terrain (Axelsson, 1999). This was performed using the industry-standard TerraScan software (Terrasolid, Finland).
2.	Horizontal and vertical co-registration of the snow-covered and snow-free LiDAR datasets. The two LiDAR point clouds were visually checked for spatial alignment in areas of no snow cover. Several profiles across the two datasets were extracted and compared throughout each transect to ensure spatial correspondence. Highway, building and cliff edge features were used to assist with the fine alignment of data.
3.	The classified and corrected point cloud data were gridded to a 1 m resolution raster surface using a Triangular Irregular Network (TIN) interpolation procedure. This gridding technique was used as it maintains point position integrity and is less susceptible to artificial smoothing of break line features such as cliffs and gorges (Keckler, 1995).

Following creation of the bare ground DEM and snow covered DSM, the LSDM was generated in ArcMap (ESRI; Redlands, CA). A histogram of the LSDM grid node values was generated so that systematic biases could be identified and to enable subsequent snow depth summaries for certain land surface classes. Bin widths of 0.1 m were chosen for the snow depth increments, as this was close to the observed precision in the LSDM (Hopkinson et al., 2011). A snow free threshold of 0.05 m was used to minimize the number of cells erroneously classified as bare snow. This depth threshold resulted in a snow covered area of 70% for the transect areas sampled. The remaining 30% of the depths were clustered around 0.0 m with some outliers (<1%) illustrating negative depths of up to several meters. Upon further scrutiny, it was found that all negative outliers occurred in areas of steep slope and were the result of small horizontal offsets between the DEM and DSM, as observed in Hopkinson and Demuth (2006) and explained by Hodgson et al. (2005). It was reasoned that the same outlying behaviour would also lead to erroneous positive outliers, so to mitigate this all data associated with vertical break lines in the DEM (slope > 80) were removed from further analysis. The final stage of data cleaning was to remove any negative depths from the LSDM and assign a zero value, as negative snow depth is logically impossible. In practice, the areas displaying negative values represented <0.1% of the total area sampled and had a negligible influence on the depth statistics.

Snow Depth Validation, Stratification and Extrapolation

Field Sampling

Field snow depths were sampled both in the ERW and in the KRV study areas over a five day period starting one day prior to the airborne LiDAR surveys and ending two days after. Within both the ERW and KRV the intent of the field campaign was to sample snow depths that were coincident with airborne LiDAR estimates while representing the range of elevation and canopy conditions experienced in ERW. Data from the more easily accessible KRV sites were a valuable supplement to the ERW analyses, as the terrain and land covers are similar to ERW so the SAU controls on relative accumulations are comparable. In practice, the KRV data were used to validate the LSDM approach, as flight lines over the ERW were offset from the field data due to a real time malfunction in the LiDAR navigation system. Consequently, the ERW field data were used to evaluate land surface type SAU influences on snow depth and to provide a comparative sample estimate of snow depth, instead of the intended correlative analysis.

Field data were collected at 25 spatially distributed sites (12 ERW and 13 KRV) at elevations ranging from <1300 m a.s.l. to >2300 m a.s.l. (see Figure 1 ) using either ground or helicopter transportation. At each of the sites, at least two profiles of nested snow depth measurements were made ( Figure 2 ). Profile lengths varied from 25 m up to 100 m in length and five measurements of snow depth were made at every 5 m increment along the profile (one centre and four in each cardinal direction, 1 m out from centre). In total, 865 individual depth measurements were made within the ERW and 1394 within the KRV. Positions of all samples were noted with hand held GPS to give approximate location to with approximately ten m. At 894 of the measurement locations in the KRV, the sampling positions were surveyed with differential GPS to provide cm-level positioning accuracy and to allow direct comparison with LSDM estimates. At all 12 ERW and three of the KRV field sites, depth-integrated bulk snowpack density measurements were made on site so that snow water equivalent (SWE) could be calculated. For profiles under forest canopy within the ERW, 34 digital hemispherical photographs (DHP) were collected so that canopy cover at each location could be evaluated (LeBlanc et al., 2005) and related to local snow depth. Each field snow depth and density measurement was entered into a GIS database and attributed with spatially corresponding DEM and canopy characteristics, so that controls on depth and density could be evaluated, and to support the definition of snow accumulation unit (SAU) properties.

Figure 2. Optimal field snow depth sampling design. Four radial depth measurements were made one metre out from at each sampling location along each profile. Due to local terrain and land cover constraints, most field sample profiles at each site did not intersect at the midpoint and many were limited to a length of 50 m. Only the locations in the KRV underwent differential GPS positioning.

Land Surface Classes and Snow Accumulation Units

In addition to the spatially coincident LSDM validation that was carried out for most of the KRV, a less direct and more subjective method of LSDM verification was conducted in the ERW and remaining KRV sites by stratifying field and LSDM data into common land surface classes and comparing the results. However, for the sake of being able to collect field snow depth data, the profiles were located in areas where snow cover was present. Occasional zero depths were noted along individual profiles but given the inability to visit the entire watershed and accurately assess snow covered area (SCA), average field depths within a given land surface class will a priori be greater than aerially aggregated LiDAR depths. The primary value in these data, however, is in assisting with the establishment of SAUs based on the variation of snow depths between land surface classes.

In mountain environments, there are many controls on snow depth, with some being more universally applicable than others. In this study, we chose aspect, elevation and canopy cover to classify and use as the basis for distinct SAUs. Even at a local or hill slope scale, the controls on snow depth distribution are complex and numerous (Pomeroy and Gray, 1995). However, the intent in this study is to identify general SAU properties that apply at the watershed scale. Greater snow accumulations tend to occur on north-facing slopes due to decreased levels of incoming solar radiation (Pomeroy and Gray, 1995; Anderton et al., 2004; Sicart et al., 2006), while a higher frequency and intensity of snowfall combined with decreased evaporation and melting generally lead to increasing snow depth with elevation (Pomeroy and Gray, 1995; Anderton et al., 2004). Increased canopy cover tends to reduce snow accumulation primarily due to canopy interception and sublimation (Hedstrom and Pomeroy, 1998; Pomeroy et al., 2002; Essery et al., 2003; Lpez-Moreno and Latron, 2008). Indeed, strong negative correlations between SWE and LiDAR-based forest canopy cover at forested sites in British Columbia have been demonstrated (Varhola et al., 2010).

Each of the three land surface types were stratified into appropriate classes based on observations in the LSDM and field data. These classes were then combined to derive unique SAU class properties, which could be used as a proxy to estimate snow depths in areas of the watershed not directly sampled. For the LSDM area, terrain aspect and elevation were directly classified from the bare earth LiDAR DEMs, while canopy cover was computed by dividing the sum of all canopy returns by the sum of all returns (e.g., Barilotti et al., 2006), and stratifying into open and closed canopy classes. While the SAU properties could be derived from the LiDAR data, they also needed to be capable of derivation from publicly available data sets to facilitate depth estimation in areas where LiDAR data were unavailable.

Watershed SWE Estimation

Given the operational need for a watershed estimate of SWE, and given available bare earth LiDAR covered ~40% of the watershed, LSDM data were extracted and summarized according to their SAU properties. Observed SAU depth properties could then be applied to corresponding SAUs across the remainder of the watershed. To facilitate snow depth extrapolation using the terrain attributes of aspect and elevation, the 1:50,000 Canadian Digital Elevation Data (CDED) topographic map sheet DEMs were downloaded from the Natural Resources Canada Geobase web portal (www.geobase.ca). The spatial resolution of the CDED DEM was ten metres. LiDAR DEM and derivative properties are generally found to be superior to publicly available DEM data sources (e.g., CDED or USGS DEMs) but at coarse resolutions the elevation and basic terrain properties are comparable even over complex land covers and alpine terrain surfaces (Rayburg et al., 2009; Hopkinson et al., 2009).

Canopy cover was mapped from the canopy closure data layer within the Alberta Vegetation Inventory (AVI) produced by the Department of Sustainable Resource Development (SRD) Alberta. The AVI canopy closure attribute is derived by photo interpretation from medium resolution aerial photography (1:40,000 or 1:60,000) and stratified into four classes: 030%, 3150%, 51-70% and 71100% (Alberta Government, 2005). These class divisions are subjectively derived from data collected at different times to the LiDAR and summarized for forest stand polygons. Therefore, while the AVI canopy closure metric is analogous to that derived from LiDAR or digital hemispheric photography, it cannot be reliably compared due to spatial and temporal inconsistencies. Also LiDAR and DHP cover estimates are floating percentages from 0100%, while the AVI has already been aggregated to four discrete classes (Alberta Government, 2005). To facilitate a practical utilization of the AVI, while account for the dominant canopy cover influence, it was decided to classify the data into open and closed canopy, and for the LiDAR, DHP and AVI the threshold was set to 30%. This value was initially chosen as it was the upper limit of the open canopy class in the AVI but was subsequently found to correspond closely to the median LSDM value in the ERW and KRV.

The approach of extrapolating LSDM data based on SAU properties assumed aspect, elevation and canopy cover from publicly available datasets captured the dominant sources of snow depth variability at the watershed scale. Using a GIS, the mean SAU snow depths observed in the LSDM transects, were applied to the entire watershed using SAUs derived from the CDED and AVI datasets. No attempt was made to reconstruct local snowpack distribution patterns as the level of detail in the public data layers and physical modeling sophistication required was not available for implementation at the entire ERW scale. Once the SAU-based estimates of depth were generated, field observations of bulk snow pack density were used to calculate watershed SWE. To place the results into context, comparisons were made to the SWE monitored at Little Elbow snow station, and with the runoff yield during the following spring melt period for the watershed outlet at Bragg Creek. Finally, A ten square kilometre section of one of the LSDM transects near the centre of the study area was reserved for comparison with the GIS extrapolated snow depths.

Results

Snow Depth

A summary of the good correspondence between LSDM and field data collected within KRV (including Marmot Creek) is illustrated in Figure 3 . The average field depth was found to be 0.54 m ( = 0.44 m), while the corresponding average LSDM value was 0.60 m ( = 0.44 m). The bias and uncertainty varied by site and landcover but overall, the correlation between field and LSDM estimates demonstrates that average snow depths can be mapped within mountainous environments to within about a decimeter as long as care is taken to ensure alignment of the ground DEM and snow surface DSM. However, the field and LSDM data collected within the KRV area were expected to display increased depth values relative to comparable land surface classes in ERW due to the precipitation shadow effect as one travels east towards the foothills and prairie lands. The two closest and comparable active provincial snow course stations are Little Elbow (2225 m a.s.l.) in the headwaters of the ERW and Three Isle Lake (2170 m a.s.l.) immediately west of Little Elbow and south of the KRV. The 30 year average winter SWE for Little Elbow is 0.299 m or 54% of the 0.550 m observed at three Isle Lake (Alberta Government, 2011).

Figure 3. LSDM plotted against corresponding field measured depth within the KRV sampling locations: (a) alpine sites; (b) valley sites; (c) intermediate elevation forested slopes; (d) grouped snow depth standard deviations for the forested sites. (e) Combined alpine, forest slope and valley profiles are. (adapted from Hopkinson et al. 2011) [Note: forest slope data thinned by averaging (4:1) to reduce sample heteroscedacity].

Stratified LSDM and field data for the ERW are illustrated in Figure 4 a and 4b, respectively, and stratified LSDM data for the KRV are illustrated in Figure 4 c. Overall, depths in the more westerly KRV are approximately double those observed at ERW, illustrating a similar trend to that in the Provincial snow course records noted above. The mean field depth for the ERW samples was found to be 0.28 m ( = 0.27 m). The mean LSDM value in snow-covered areas (defined as all areas where LSDM > 0.05 m) was comparable at 0.26 m ( = 1.2 m), or 0.18 m when zero depths in snow-free areas were considered. The higher standard deviation in the LSDM is due to outlying depth estimates typically found in areas of steep slope or associated with deep gullies. Indeed, it was observed in the field data that depressions in alpine areas can fill with snow up to several meters and the ability for the LSDM to accurately capture these deep accumulations in alpine areas is illustrated in Figure 3 a.

Figure 4. Snow depth sample data stratified by terrain aspect, canopy fractional cover, and terrain elevation. (a) LSDM results for sampling transects collected in the ERW. (b) Field sampling snow depth results collected in the ERW. (c) Temporally coincident LSDM results collected in the KRV area immediately west of the ERW.

Apart from the already documented reduction in snow depth magnitude in eastern areas of the front ranges, similarity in LSDM behaviour at KRV and ERW is apparent when the data are stratified by north vs. south aspect, canopy cover and elevation ( Figures 4 a and 4c). North slopes possess deeper snow than south; open canopies illustrate deeper snow than closed canopies; and snow depth increases with elevation up to treeline. These observations are consistent with documented observations for northern hemisphere meso-scale watershed environments. The observation for snow depth to peak at tree line requires further explanation. Snowpack above treeline is redistributed and ablated by blowing snow (Greene et al., 1999) and at lower elevations by canopy interception in evergreens (Pomeroy and Gray, 1995). However, the treeline in this region is characterised by a high proportion of deciduous larch trees which trap blowing snow and do not have significant interception losses (Fisera, 1977). Above treeline there is some uncertainty as to the average snow depth properties in the ERW, as the alpine zone was minimally sampled. Indeed, the little sample coverage there is suggests snow depth might reduce above treeline ( Figure 4 a). However, while a decrease might occur locally on some slopes, this is not the pattern that was observed over the larger more continuous sampling area of the KRV. Based on the immediate adjacency of the KRV and ERW areas, the terrain and land cover diversity in the upper ERW is similar to that of the KRV and thus it is likely that snow depth above treeline in the ERW is similarly distributed and also displays no elevation gradient.

As expected, stratified field snow depths were systematically greater than aerially integrated LSDM values due to field data not representing snow free areas, whilst the LSDM represented all areas with and without snow cover ( Figures 4 a and 4b). As with the bulk field vs. LSDM estimate, the differences when the SCA was factored in reduced significantly and were well within the 0.27 m standard deviation observed in all field results. Both canopy cover and elevation stratifications of the field data illustrate the same general tendencies as observed in the LSDM. However, the north vs. south aspect stratification of the field data did not. This is due to field samples being collected in areas where snow accumulated and were readily accessible for measurement, and these areas tended to be nearer to the base of slopes, in forested or sheltered areas where aspect exerts less control. Therefore, in this case, the slope aspect stratification observed in the ERW and KRV LSDMs could be more reliable indicator of snowpack behaviour than the field data, as these results are not influenced by field access limitations.

Watershed SWE

Based on the above observations, ten unique SAUs were created from all plausible permutations of: north (27090) and south (90270) aspect; closed (>30%) and open (<30%) canopy cover; and low (<1700 m a.s.l.), medium (1700 >2200 m a.s.l.) and high (above treeline or >2200 m a.s.l.) elevation classes. [Note: above treeline, there is no canopy cover, so the two closed canopy classes for north and south facing slopes are redundant.] The LSDM depth data were stratified into these ten SAUs and used to train the CDED- and AVI-based SAUs for the remainder of watershed ( Figure 5 ).

Figure 5. Classification scheme used to create 10 unique snow accumulation units (SAUs) based on elevation, aspect and canopy cover. Rectangular box near centre of watershed illustrates LSDM area used to test the spatial extrapolation approach.

Depth integrated field snow density measurements ranged from 18 g cm³ to 50 g cm³ with a mean of 26 g cm³ (= 8 g cm³). There were no discernible elevation, aspect or canopy cover trends in the density observations so the mean value was used to convert depth to SWE. By applying the trained snow depths and mean snowpack density to each of the ten watershed SAUs and then totaling the results, the total SWE estimate for the watershed was 46.0 10⁶ m³ (a depth of 0.058 m). The basin SWE depth distribution using this approach is illustrated in Figure 6 . Placing this into a water resource context, we note that this quantity of snowpack water storage was equivalent to 18% of the total annual yield monitored at the Bragg Creek gauging station for 2008 (Environment Canada, 2009), or 25% of the spring runoff from April to June for 2008 ( Figure 7 ). Furthermore, this estimate of watershed SWE is approximately 20% that of the Alberta Environment SWE measurement made at the end of March at the Little Elbow snow monitoring station in the headwaters of the ERW ( Figure 7 ). It is self-evident but worth reiterating that while snow monitoring data provides a useful long term proxy for water resource availability (assuming no significant climatic or watershed changes) the magnitude of SWE recorded at such sites can differ greatly from the mean snowpack conditions in the entire watershed.

Figure 6. Spatial distribution of estimated SWE across watershed.

Figure 7. Estimated watershed SWE at end of March relative to the increasing winter SWE at Little Elbow snow course station (2200 m a.s.l.) and the specific yield of the Elbow River Watershed at Bragg Creek.

The comparison of LSDM test data and GIS extrapolated snow depth across the full range of ten SAUs, indicated that while the correlation was reasonable (r² = 0.78), the GIS results tended to over-estimate the LSDM data by approximately ten percent near the centre of the watershed ( Figures 5 and 8 ). This demonstrates that there are noticeable spatial variations in terrain and landcover controls on snow depth within the watershed. This is further illustrated, if we recalculate watershed SWE using the average LiDAR snow depth from all sampling transect areas (0.18 m) and simply apply the observed snow density (0.26 g cm³). Using this simplified approach the total watershed SWE volume dropped to 37.4 10⁶ m³ (depth = 0.047 m). This 17% reduction relative to the land surface class SAU approach illustrates the need for a more sophisticated method than simply using the average LSDM sample depth; i.e., land surface classes that exert some control over snowpack depth were not equally represented within the sample data or the watershed as a whole.

Figure 8. Comparison of GIS extrapolated snow depth with LSDM sampled snow depth over all 10 land surface classes with the test area.

Discussion

Field and airborne snow depth sampling in the ERW took place during a relatively dry period when snow cover at the lower elevations was patchy and thin. The average LiDAR snow depth was estimated to be 0.18 m overall or 0.26 m in snow covered areas only. By stratifying the LiDAR snow depth data into land surface class SAUs and applying these depths to the same SAU classes derived from analogue data sources, an estimate of basin SWE (46.0 10⁶ m³) was generated. However, steep heterogeneous terrain and dense canopy covers are not the ideal conditions for point-level LiDAR snow depth assessment, with uncertainties at individual locations being potentially high. We contend that when multiple data points are aggregated across common land surface classes, some of the random point-level error is mitigated. Overall, it was found that LiDAR estimates, when stratified using the same criteria as field measurements, tended to slightly under-estimate snow depth in the ERW, whereas direct comparison at the KRV sites illustrated a slight over-estimation. In both cases, however, the LSDM precision was on the order of 0.1 m.

Based on similar LSDM observations elsewhere, it has been reported that ground-level vegetation tends to systematically elevate true ground surface by up to ~0.1 m (Hopkinson et al., 2005), whilst snowpack surfaces, being highly reflective and smooth, tend to be more accurately represented in LiDAR data. The net effect is an under-estimation of snow depth in areas of dense ground level foliage. Steep terrain is known to introduce random errors into the surface elevation due to the propagation of horizontal uncertainty (Hodgson et al., 2005; Hollaus et al., 2006). A slope raster created from the 2006 LiDAR DEM indicated that only 1% of the surface exceeded 45. The proportional effect of these depth uncertainties, therefore, would be limited, and most likely there would be some compensation of under- and over-estimated depths. A cautionary note, however, is that steeper slopes tend to occur higher in the watershed on the western side, where snow depths are expected to be higher. Therefore, it might be reasonable to expect that random errors in depth would increase in those areas of alpine watersheds that typically experience deeper snowpack.

While individual LSDM grid-level values of several metres were observed in some areas and zero depths occurred over approximately 30% of the watershed, the mean depth of 0.26 m was approximately two times the manufacturer quoted 0.15 m accuracy for a single LiDAR data set (Optech Incorporated, 2005). A certain magnitude of error is to be expected even over perfectly flat and unambiguous ground or snow surfaces. In an extreme example, then, if both LiDAR ground and snowpack surfaces possessed equal but opposite systematic errors of this magnitude, this could result in areas with snow depths of 0.3 m appearing devoid of snow and areas of no snow displaying depths of 0.3 m. Add to this small systematic depth under-estimation due to foliage cover and random terrain slope effects, and it is clear that an observed mean depth of 0.26 m could conceivably have 100% statistical uncertainty! Given the high correspondence displayed within the KRV and the mean field and LiDAR depth estimates were within 0.02 m, however, such a high level of error is not expected.

To calculate a meaningful snow depth error for the total area sampled by LiDAR would require an investigation of all the individual components of random error and systematic bias. These components would then need to be applied to each grid node of the LSDM, and not to the overall mean depth. It is possible that a mean depth of 0.26 m could contain a sufficiently high proportion of uncertainty that the estimate might not be statistically meaningful. However, this does not mean the approach of LiDAR-based snowpack monitoring has no value, it simply provides limits on when such approaches are viable. Without the type of uncertainty assessment suggested above, it would be impossible to quantify at what average snow depth the method provides results at a pre-determined level of confidence. However, given the types of bias and random behaviour discussed are likely to introduce uncertainties at the decimeter scale over most areas and potentially at the metre scale for small proportions (~1%) of the watershed, it is reasonable to assume that mean snow depths of approximately 1 m would produce reliable and useable results at a high level of confidence.

It was demonstrated that a relatively small spatial variation of <5 km in LiDAR depth observations led to the GIS results over-estimating the LSDM class-summaries by approximately 10% in the test area. The available base LiDAR data for LSDM creation was limited to the northern 40% of the watershed and had limited representation above tree line. While proxy data were available from the nearby Kananaskis River Valley area to provide some insight as to the expected snow depth patterns, it is known that snow depth can vary significantly at meso-scales (e.g., Elder et al., 1998). The controls on depth at the watershed scale are not always localized and can vary due to synoptic meteorological variations, orographic and precipitation shadow effects. Therefore, by having no sample representation in the southern part of the watershed, this created an unquantifiable level of uncertainty in the SWE estimate generated. If this method of snowpack water resource assessment were to be used in an operational setting, a more spatially complete base LiDAR coverage would be required to enable sampling over appropriately spaced and representative land surface classes covering the full geographic and elevation range of the watershed. LiDAR snow sampling transects would need to be spaced close enough to capture significant spatial variations but far enough apart to minimize sample redundancy and cost, while maximizing aerial coverage. Further research is needed to optimize the transect configuration (altitude, swath width, spacing and sampling density) for a given environment.

Conclusion

Despite the challenges and uncertainties discussed above, this study has demonstrated a remote sensing approach for mountainous headwater snowpack water resource assessment. The use of LiDAR snowpack sampling combined with a GIS land surface classification of snow accumulation units (SAUs) to fill the gaps between transects shows potential as a supplement to traditional in-situ helicopter-based snow depth sampling. Benefits of a LiDAR sampling and GIS extrapolation approach include a reduced need for field personnel deployment combined with spatially explicit and rapid (within two to three days of acquisition) estimations of snowpack depth and volume at the meso-scale. This relatively simple approach to snowpack monitoring could be of significant value to water resource managers when accurate and repeatable estimates of the spring runoff volume are needed for seasonal water supply, irrigation and power generation forecasts, as well as flood risk assessment.

At the current time it is not thought that the approach to snowpack monitoring described in this study is sufficiently cost-effective or accurate for operational water resource monitoring. However, the speed, aerial coverage, accuracy and costs of airborne LiDAR data have all improved greatly in recent years and continue to do so. Therefore, given growing water scarcity and potential flood risk challenges in parts of the Bow River basin, it is possible that the need for more accurate spatial estimates of headwater snowpack volumes will make the investment in LiDAR snow sampling worthwhile at sometime in the next decade. Indeed, active snow course monitoring already requires helicopters to transport snow monitoring crews and the Alberta LiDAR base coverage is gradually heading towards completion (Airborne Imaging, 2010). Therefore, it will soon be feasible to mount a small and cost effective LiDAR profiling sensor on the helicopter so that snow depth transects could be automatically collected en route between snow course sites. While there would be a cost associated with the hardware and the post-processing, the actual operational costs would be little to no more than they are at present.

Acknowledgements

The City of Calgary and Government of Alberta are gratefully acknowledged for funding this project. Funding for the 2007 LiDAR acquisition over Marmot Creek was provided by the Canadian Foundation for Climate and Atmospheric Sciences, and the Natural Sciences and Engineering Research Council. Sustainable Resource Development and Alberta Environment are acknowledged for providing ATV field transportation, helicopter time, and loan of sampling equipment. The authors would like to thank Wolf Survey and Mapping (a division of Destiny Resources) for the generous loan of survey equipment and the use of vehicles to perform field snow sampling. Dr John Diiwu is acknowledged for supporting this project at the conceptual stage. Chris DeBeer and Chad Ellis are acknowledged for snow depth field collection and analysis. Hubert Howe and James Churchill provided valuable assistance in the field. We are also grateful to Airborne Imaging Inc. for flying the 2008 snowpack LiDAR survey on a cost-recovery basis. Mike Demuth, Geological Survey of Canada, and the Canadian Consortium for LiDAR Environmental Applications Research provided aviation support in 2007.