Snow depth mapping in high-alpine catchments using digital photogrammetry

. Information on snow depth and its spatial distribution is crucial for numerous applications in snow and avalanche research as well as in hydrology and ecology. Today, snow depth distributions are usually estimated using point measurements performed by automated weather stations and observers in the ﬁeld combined with interpolation algorithms. However, these methodologies are not able to capture the high spatial variability of the snow depth distribution present in alpine terrain. Continuous and accurate snow depth mapping has been successfully performed using laser scanning but this method can only cover limited areas and is expensive. We use the airborne ADS80 optoelectronic scanner, acquiring stereo imagery with 0.25 m spatial resolution to derive digital surface models (DSMs) of winter and summer terrains in the neighborhood of Davos, Switzerland. The DSMs are generated using photogrammetric image correlation techniques based on the multispectral nadir and backward-looking sensor data. In order to assess the accuracy of the photogrammetric products, we compare these products with the following independent data sets acquired simultaneously: (a) manually measured snow depth plots; (b) differential Global Navigation Satellite System (dGNSS) points; (c) terrestrial laser scanning (TLS); and (d) ground-penetrating radar (GPR) data sets. We demonstrate that the method presented can be used to map snow depth at 2 m resolution with a vertical depth accuracy of ± 30 cm (root mean square error) in the complex topography of the Alps. The snow depth maps presented have an average accuracy that is better than 15 % compared to the average snow depth of 2.2 m over the entire test site.


Introduction
Snow is an important resource in alpine regions not only for tourism (e.g., Elsasser and Bürki, 2002;Nöthiger and Elsasser, 2004;Rixen et al., 2011) but also for hydropower generation, water supply (e.g., Marty, 2008;Farinotti et al., 2012) and ecological aspects of the local mountain flora and fauna (e.g., Wipf et al., 2009). Snow is also important in the context of natural hazards, such as avalanches prevention or flood forecasts in spring and early summer for the valleys downstream. For the latter it has been shown that the snow distribution at the winter maximum before the beginning of the melting period strongly determines the temporal evolution of the remaining snow resources and -if converted to snow water equivalent (Jonas et al., 2010) -the potential melt water runoff during the melting period . Several studies reported a very high spatial variability of snow depth and other snow pack parameters at different spatial scales in mountainous regions (e.g., Elder et al., 1991;Schweizer et al., 2008;Lehning et al., 2008;Grünewald et al., 2010;Egli, 2011). This high variation in snow cover distribution on very small scales requires a high spatial resolution of snow samples to measure different parameters of the snow pack such as, e.g., the areal mean snow depth on complex alpine topography and the temporal evolution of snow-covered areas during melt with high areal representativeness and low absolute uncertainty. In other words, snow pack monitoring in alpine terrain requires an area-wide observation with a large number of snow depth point measurements distributed over the area of interest.
Published by Copernicus Publications on behalf of the European Geosciences Union.

Y. Bühler et al.: Snow depth mapping in high-alpine catchments
Currently, in the Swiss alpine region snow depth is measured at specific locations by automated weather stations or observers in the field, while both observations are restricted to flat sites exhibiting a rather homogeneous snow cover (Bründl et al., 2004;Egli, 2008). These flat field point measurements are assumed to represent snow cover characteristics for a larger area around the stations and are therefore interpolated over large distances and are combined with snow cover information from optical satellites (Foppa et al., 2007). This method is unable to capture the small-scale variability of snow depth. Investigations into the representativeness of point snow depth measurements on snow depth for entire catchments are sparse (Grünewald and Lehning, 2014).
Remote-sensing instruments have been used for snowrelated studies since such data became available (e.g., Rango and Itten, 1976;Dozier, 1984;Hall and Martinec, 1985). A very common parameter measured by remote-sensing instruments is snow-covered area (SCA). Operational products on a global scale, such as MODIS snow-cover products (Hall et al., 2002) or GlobSnow (Koetz et al., 2008), are widely used today (Frei et al., 2012). For example, Dozier (1989), Nolin and Dozier (1993), Fily et al. (1997) and Dozier et al. (2009) published investigations into snow grain size with finer spatial resolution on a regional scale. Snow depth and snow water equivalent (SWE) has been assessed using passive microwave sensors (e.g., Ulaby and Stiles, 1980;Chang et al., 1982;Pulliainen, 2006). However, due to the coarse spatial resolution of these sensors (25 km), the results do not display small-scale snow cover characteristics of alpine catchments. Active microwave sensors use much smaller wavelengths (millimeters to centimeters) and achieve finer spatial resolutions of up to 20 m (e.g., Schanda et al., 1983;Shi and Dozier, 2000). However, this method is limited to dry snowpacks and faces problems in steep high-alpine terrain (Buchroithner, 1995). Nolin (2011) and Dietz et al. (2012) give an overview of recent advances in the remote sensing of snow.
Terrestrial laser scanning (TLS) was previously used to derive spatially continuous snow depth (Prokop, 2008;Grünewald et al., 2010). Even though the accuracy of such measurements is very good (usually better than 0.1 m, depending on laser footprint and distance from sensor), largescale catchments such as the Dischma Valley ( Fig. 1) cannot be covered completely. Data acquisition with TLS is time and manpower consuming and only possible at easily accessible spots under fair conditions (depending on the avalanche situation, weather) for areas within the line of sight of the measurement location. This results in limited coverage and many data gaps, e.g., behind bumps. Airborne laser scanning (ALS) from helicopters or airplanes can cover larger areas in a shorter time also under difficult avalanche danger situations. Recent studies demonstrate that the accurate mapping of snow depth is possible (Deems et al., 2013;Melvold and Skaugen, 2013). However, the costs to cover larger areas are still high  and overflights are, as with digital photogrammetry, restricted to fair weather conditions. Previous attempts to map snow depth using scanned aerial imagery were already made 50 years ago (Smith et al., 1967) and the topic was investigated in detail by Cline (1993Cline ( , 1994. However, their results suffer from image saturation and insufficient reference data leading them to the conclusion that photogrammetry has much potential but is not yet accurate enough for large-scale snow depth mapping. Ledwith and Lunden (2010) used scanned aerial imagery to derive digital elevation models over glaciated and snow-covered areas in Norway. They report a mean accuracy of 2.8 m in comparison with differential Global Navigation Satellite System (dGNSS) transects, which is clearly too low for meaningful snow depth mapping in alpine regions. Lee et al. (2008) used a digital mapping camera (DMC) digital frame camera to cover an area of approximately 2.3 km 2 with a very high mean ground sampling distance (GSD) of 0.08 m. The reported mean differences compared to dGNSS measurements are approximately 0.15 m, stressing the great potential of digital photogrammetry for accurate snow depth mapping. However, no snow depth mapping has been performed and compared to different reference data sets, covering larger areas.
In this investigation we apply digital photogrammetry based on high spatial resolution aerial imagery (0.25 m) to calculate digital surface models (DSM) of winter and summer terrain. Traditional photogrammetry using analogue aerial imagery and 8 bit digital sensors faced problems over snow-covered areas mainly due to saturation and the homogeneous surface (Kraus, 2004). Modern digital sensors can acquire data with 12 bit radiometric resolution to overcome these limitations. We calculate spatially continuous snow depth maps using the summer and winter DSMs for two test sites near Davos, Switzerland (145 km 2 in total). This technology is much more economical for covering large areas than ALS or TLS but still has an acceptable spatial resolution to map small-scale spatial variability. To assess the accuracy of our results, we compare the calculated snow depths to hand measurements, dGNSS points, TLS measurements and ground-penetrating radar (GPR) transects acquired simultaneously with the aerial imagery.

Test sites Wannengrat and Dischma, Davos, Switzerland
The two areas covered by the ADS80 sensor on a Pilatus Porter airplane are located close to the winter sport resort Davos in the eastern part of Switzerland (Fig. 1). The Wannengrat test site is located to the north of Davos and covers an area of approximately 3.5 × 7.5 km (26.25 km 2 ). The valley bottom is about 1500 m a.s.l.; the highest peaks reach up to 2780 m a.s.l. (Amselflue at the southwestern part of the test site). The large ski resort Davos-Parsenn is located at the northeastern edge of the test site. The covered mountain chain is characterized by high-alpine meadows, rock faces and scree-covered areas. The area below 2000 m a.s.l. is covered by sparse forest which becomes dense as of ca. 1800 m a.s.l. The Wannengrat area is used as test site for various research projects at the WSL Institute for Snow and Avalanche research SLF, mainly because of the very good accessibility from Davos even if the avalanche danger level is considerable. We collected hand-measured snow depth plots, dGNSS points and TLS data sets close to the Wannengrat peak as reference data sets (see Sect. 3.2) on the day of the ADS80 data acquisition. The Dischma test site is a high-alpine valley branching off from the main valley of Davos (1500 m a.s.l.) in a southeasterly direction and rising up to 2000 m a.s.l. at the end of the valley. It covers an area of ca. 7 × 17 km (119 km 2 ) and contains the complete catchment of the Dischma Creek where several hydrological studies have been performed (Bavay et al., 2009). The peaks surrounding this catchment reach up to 3130 m a.s.l. (Piz Grialetsch). Forest covers the lower part of the valley up to 2000 m a.s.l. The southeastern two thirds of the valley are completely forest free. We collected GPR snow depth measurements at the valley bottom in the northwestern part of the test site as reference data on the day of the ADS80 data acquisition. Because the central flight strip at the valley bottom was corrupted in the summer 2010 data set, resulting in a low-quality summer DSM, we repeated the flight in summer 2013.

Sensors and data sets
To measure spatially continuous snow depth and to validate these measurements, we used independent state-of-theart technologies. It is a difficult task to measure multiple, spatially widely distributed snow depths in high-alpine areas within a short time span. Several teams were deployed in the field on the day of the ADS80 data acquisition, guaranteeing a small temporal offset to the ADS80 imagery because snow depth can change very quickly under spring conditions.

Airborne optoelectronic scanner ADS80
Two optoelectronic line scanner data sets were acquired with the ADS80-SH52 sensor. The acquisition of the summer images was carried out on 12 August 2010 (Wannengrat) and 3 September 2013 (Dischma). Winter imagery of the snowcovered sites was acquired on 20 March 2012 (close to the maximum snow cover, peak of winter). The covered area consists of 12 overlapping image strips (approx. 70 % overlap across track) flown during approximately 90 minutes at an elevation of approximately 4000 m a.s.l. (1500 m above mean ground elevation). The mean ground sampling distance (GSD) of the imagery is 0.25 m, limited by the minimal flying height for high-alpine terrain . The ADS80 scanner simultaneously acquires four spectral bands (red: 604-664 nm; green: 553-587 nm; blue: 420-492; near infrared: 833-920 nm) and a panchromatic band (465-676 nm) with a radiometric resolution of 12 bits and two viewing angles (nadir and 16 • backward). The nadir and forward-looking panchromatic bands were not used due 232 Y. Bühler et al.: Snow depth mapping in high-alpine catchments to saturation issues caused by the broader sensitivity of these bands. GNSS/IMU (Global Navigation Satellite System/inertial measurement unit)-supported orientation of the image strips supplemented by the use of ground control points achieve a horizontal accuracy (x, y) of 1-2 GSD (0.25-0.5 m). The sources of the ground control points used are a combination of GNSS ground surveys and already existing oriented stereo images (with unknown absolute accuracy). We tried to distribute the ground control points (GCPs) regularly; however, they are denser at the lower altitudes. We applied between 11 and 33 ground control points per acquisition date showing residuals of 3 to 21 cm in x, 4 to 17 cm in y and 10 to 33 cm in z direction. The ADS sensor was successfully used to detect avalanche deposits in the area of Davos (Bühler et al., 2009). Sandau (2009) gives more detailed information on the Leica ADS optoelectronic scanner.

Manual snow depth measurements
Simultaneously with the ADS80 data acquisition, a field team acquired manual snow depth measurements using a 3.2 m avalanche probe at 15 different plot locations within the test site Wannengrat. A plot consists of 5 × 5 probe measurements with a distance of 2 m between points (Fig. 2a), resulting in 375 single-probe measurements localized using dGNSS data of the corner points. Because snow depth can vary substantially within the distance of some decimeters if there is, e.g., a rock at the surface (Lopez-Moreno and Nogues-Bravo, 2006), we use the average snow depth and the standard deviation to compare it to the corresponding ADS80 snow depth values within this 10 × 10 m area (Fig. 2b). The acquisition of field measurements is very challenging because the terrain is steep and human mobility is limited. The avalanche danger for wet snow avalanches rises quickly during the day due to sunny spring weather conditions, limiting the time the field team can move within the test sites. Therefore, the number of performed field measurements at 15 plots distributed over an area of 1 × 1.5 km is close to the possible maximum that can be obtained with the number of those participating in the experiment. Because this number is, in our opinion, not sufficient to assess the potential of the method proposed, we apply further reference data sets.

Differential Global Navigation Satellite System (dGNSS) measurements
During ADS80 data acquisition on 20 March 2012, 137 dGNSS points were measured with the Leica GPS 1200 device in the test site Wannengrat (Fig. 2a). The points were measured with real-time correction using the virtual reference station of the swisstopo Automated GNSS Network for Switzerland (AGNES) in Davos. The surveyed points show a horizontal accuracy better than 1 cm (1 standard deviation) and a vertical accuracy better than 2 cm (1 standard deviation). Measured points represent the top of the snow cover in meters above sea level.

Terrestrial laser scanning (TLS)
In the last decade, terrestrial laser scanning has been increasingly applied for continuous snow depth mapping (e.g., Deems, 2013;Schirmer et al., 2011;Prokop, 2008;Prokop et al., 2008). To calculate snow depth, an elevation model of the bare ground and another of the snow-covered winter surface is produced. Snow depth is then obtained by subtracting the two surfaces from each other. In this study, we use the Riegl LPM-321 device operating at 905 nm. This device has been proven to accurately measure snow depth in alpine terrain (Prokop, 2008;Prokop et al., 2008). Grünewald et al. (2010) compared TLS measurements to tachymeter measurements and found a mean vertical deviation of 4 cm with a standard deviation of 5 cm at a distance of 250 m using the LPM-321. To assure the quality of the laser scans, we additionally performed reproducibility tests. A laser scan acquired in a coarse resolution (three points per square meters at a distance of 300 m) was compared with the full-resolution acquisition (eight points per square meters at a distance of 300 m). This allows detecting misalignments between the two data sets due to an instable scan setup (unstable tripod, wind influence, etc). Scans which showed a mean difference larger than 10 cm were excluded. The upper end of the Steintaelli was scanned once in summer 2011 and a second time on 20 March 2012 during the ADS80 data acquisition (Fig. 2c). Fixed installed reflector points were used to match the summer and winter TLS data sets.

Ground-penetrating radar (GPR)
GPR data were collected using a MALÅ ProEx system configured for synchronous measurements with four pairs of separable shielded 400 MHz antennas. The antennas were set up as a common-midpoint (CMP) array with separation distances of 0.31, 0.95, 1.6 and 2.8 m. The GPR antennas were mounted on two pulkas, which were rigidly connected to one another to guarantee fixed relative antenna positions throughout the measurements. This assembly was pulled along a transect of 4.8 km length. After the initial stacking of four individual traces, data were recorded every 0.5 s, which resulted on average in one record every 30 cm along the transect. GPS coordinates were taken every second along the transect using an onboard GPS receiver as well as an external Trimble GeoExplorer 6000 dGNSS system. GPS data were slightly smoothed before associating them with the GPR data records. Snow depth data were obtained using standard CMP analysis procedures, partly involving the commercial software package ReflexW 7.0 (Sandmeier, 2013). Along the GPR transect, we obtained 130 manual snow depth readings. These data were used for cross-validation of the GPR data. Concurrent GPR and manual snow depth ranged from 0.76 to 2.70 m. Correlation between both data sets resulted in an R 2 of 0.96 and an RMSE of 0.07 m.

Generation of summer and winter digital surface models
For DSM generation we use the "automatic terrain extraction" (ATE) as part of the SOCET SET software version 5.4.1 from BAE Systems. The software implements an areabased algorithm calculating similarity measures with a twodimensional cross-correlation approach. ATE has no need for user input on specific image-matching strategies and parameters as a function of terrain type. ATE uses an "inference engine" which adaptively generates image-matching parameters depending on facts such as terrain type, signal power, flying height or X and Y parallax. A user-given post spacing distance is used to control image correlation spacing (e.g., 2 m); hence, cross-correlation is not calculated for every image pixel (Zhang and Miller, 1997). We use the green, red and near-infrared bands of the sensor as input. The nearinfrared band absorbs a larger part of the incoming radiation over snow, and the reflected signal is sensitive to grain size variation within short distances (Bühler et al., 2015). This improves the performance of the ATE point-matching algorithm, in particular over old snow covers not recently covered by new snow. ATE SOCET SET gave the best results regarding blunders and completeness. We also tested NGATE from SOCET SET from Inpho. XPro and MatchT use semi global matching techniques (SGM) for image correlation. Although this is the state-of-the-art method for dense image matching (especially in urban areas with a very high image overlap), the results on snow surface were comparable to or even worse than ATE SocetSet. MatchT gave similar results to ATE but was much slower regarding calculation time. The stereo blocks of each year were orientated separately. Although jointly adjusted image blocks would increase the relative accuracy between the blocks, it was not possible to adjust blocks in this way due to different visibilities of ground control points in different years. We want to demonstrate the workflow for future campaigns where a reorientation of all existing blocks together is not feasible.
As this study focuses on snow depth mapping for widescale applications, we set the spatial resolution of the derived DSM to 8 × GSD (8 × 0.25 m), which results in significantly lower demand in CPU usage compared to a resolution at pixel level. Additionally we apply a 3 × 3 low-pass In our research setup, single buildings and forest or scrub cannot be modeled with sufficient horizontal accuracy due to the limited spatial resolution of the input imagery. Slight differences in x, y positions of such objects in the summer and winter DSM would lead to big outliers in the snow depth product. Therefore, all buildings and forest or scrub areas were masked out. For the detection of forest or scrub areas, a combination of NDVI (normalized differenced vegetation index) and a canopy height layer was applied. With this approach, all visible vegetation in the winter images and vegetation higher than 1.5 m in the summer images was masked out. The detection of buildings (settlements) only from spectral or elevation information is not feasible since rock-covered areas return an identical spectral signature as settlements and are prone to big outliers. Therefore, we use the building layer from the topographic landscape model (TLM) of the Swiss Federal Office of Topography. This step might not be necessary if the input imagery had a higher spatial resolution (15 cm or better).
Large-scale imagery of mountainous, snow-covered landscapes show a maximal range of radiometric-image information over a short distance, which is highly demanding for image correlation processes. For this reason, generating a complete DSM from one entire image strip is not expected to give optimal results for snow-covered areas. As a response to this challenge, we divided the test site into 809 tiles, for which DSMs were calculated separately. Another well-known difficulty in steep mountain areas is a suboptimal viewing angle or even occlusion in an image strip. Considering this difficulty, we calculated two DSMs for each tile, using the "most nadir" and the "second-most nadir" color infrared (CIR) image strips (near infrared, red, green) to increase the chance of a good image match for a given point on the ground. For the generation of the final DSM, we calculated the mean slope for every processed DSM tile. By selecting the DSM with the smaller mean slope for every given tile, big blunders caused by a nonoptimal viewing angle or occlusion could mostly be automatically eliminated. We used the approach described to process all DSMs.
The orientation of ADS80 image strips has to be considered as a critical point especially for winter images. All processing and evaluation efforts are worthless if there is a lack of accuracy in image orientation. Due to a small number of highly accurate reference points in remote areas and sometimes almost unrecognizable ground control points in snow-covered, high-alpine regions (e.g., east part of Dischma Valley without any anthropogenic features), orientation quality shows certain limitations. For the areas mentioned, orientation during the post-processing of image strips (software Leica XPro) could not be substantially improved, resulting in a final orientation accuracy of about 1 GSD. Welldistributed artificial reference points measured at the ground with dGNSS could improve the orientation quality substantially but were not available for the winter 2012 imagery.

Results and validation
To quantify the accuracy of the digital photogrammetry products, we use the following measures recommended by Höhle and Höhle (2009) to compare elevation data sets from different sources: -(a) the root mean square error: ; this measure is often used and simple to calculate but very prone to outliers.
-(b) normalized median absolute deviation: , where h j denotes the individual errors and m h is the median of the errors.
-(c) additionally, we use the empirical correlation coefficient to assess how well two snow depth measurements from different sources correlate.
To make the comparison of elevations in digital elevation model (DEM) products possible, it is crucial that a coherent coordinate system is applied for all data sets. We use the swisstopo LV03 LN02 (CH1903) system with the elevation reference point at Repère Pierre du Niton H(RPN) = 373.6 m a.s.l. in Geneva, Switzerland (swisstopo, 2008).

Photogrammetric summer DSMs (DSM ADS )
Three DSM ADS (winter 2012, summer 2010 and 2013) were processed for this study. For a quantification of the quality of the derived DSM ADS , we perform an accuracy assessment using a digital terrain model ( in summer 2009 by an airborne laser scanner (Riegl LMS-Q240i) mounted on a helicopter as a reference and assuming the changes in terrain to be negligible (which might not be true for areas prone to erosion and deposition). The average point density acquired was two to three points per square meter from an average flight height of 300 m above ground. Airborne laser scanning is reported as a very accurate method for DTM generation in various studies (e.g., Aguilar and Mills, 2008;Höhle and Höhle, 2009), also on snow (Deems et al., 2013) and in high-alpine terrain (Bühler and Graf, 2013). The quantification of the accuracy is described by the distributions of vertical deviations between the two data sets (886 000 points). Vegetation and buildings were excluded for the analysis. The statistical measures in Table 1 show a good correspondence between the DTM ALS and DSM ADS . The RMSE value without outlier removal indicates the presence of big outliers. Since the mean values of the deviations with and without outlier removal differ only by 3 cm, these big outliers are both negative and positive. A detailed quality assessment of DSMs derived by ADS80 image strips in very steep and complex alpine terrain showed that the accuracy of photogrammetric DSMs decreases significantly in terrain steeper than 50 • , explaining the occurrence of the above-mentioned outliers .
In Fig. 3, on the right, image correlation completeness in terms of correlated and interpolated points is shown for a section of the test site at Wannengrat for winter 2012. Imagematching completeness for the whole test site is given in Table 2 (Wannengrat and Dischma without buildings and vegetation). These results show the high matching success with the 12 bit imagery, in particular over snow-covered areas.

Snow depth maps
The snow depth maps are calculated by subtracting the photogrammetric winter DSM from the summer DSM. The spatial resolution is 2 m as for the input DSMs. Because negative snow depths cannot occur, values smaller than 0 are set to "no data". Consulting the input orthophotos of the winter data acquisitions allows identifying whether a certain area is snow free or not. Overall, 19.42 % of all pixels are classified as trees and scrubs and 1.65 % as buildings. Of the remaining pixels, 4.83 % were classified as "no data".
The snow depth maps generated (Figs. 4 and 5) reveal a very high spatial variability of snow depth even within small distances. Snow depth can vary by more than 5 m within a few meters. Snow traps for wind-blown snow and deposits from past avalanche events are clearly visible. We identify the same snow trap features in the Wannengrat area, reported by Schirmer et al. (2011), that were measured in winter 2008. This indicates that snow traps and cornices are persistent over different winters due to dominant main wind directions. High snow depths due to avalanche deposits are persistent in tracks where avalanches occur several times each winter but are not where avalanches occur with return periods of more than 1 year.
The large area at the northern edge of the Dischma test site (Fig. 5) classified as "no data" is Lake Davos. This natural lake is used for power generation during winter, and the surface level is lowered by up to 50 m. By subtracting the winter DSM from the summer DSM we get clearly negative values in this area, which are classified as outliers. The large outlier areas at the southern edge of the investigation area are the glaciers of the Grialetsch range. These small glaciers lost a significant part of their volume between winter 2012 (winter DSM) and summer 2013 (summer DSM), and their surface elevations were lowered (Zemp et al., 2006). Therefore, highly positive values occur and are classified as outliers. Further outliers occur in very steep terrain (> 50 • ) because the footprint of the sensors is very small in such areas , demonstrating the limitation of the proposed method for snow in rock faces. These areas are less relevant for most snow depth applications because little snow usually accumulates in very steep terrain (e.g., Fischer et al., 2011).

Differential Global Navigation Satellite System (dGNSS) measurements
A comparison of the ADS-derived Winter 2012 DSM with 137 dGNSS points, describing elevations in meters above sea level (top of the snow cover), results in an RMSE of 0.37 m and a NMAD of 0.28 m. With a mean of 0.21 m, the ADS DSM models the surface of the snow cover systematically higher than dGNSS measurements. For the area Wannengrat in Fig. 1a, it can therefore be assumed that snow cover thickness is overestimated using photogrammetric methods, mainly because of orientation inaccuracies. A bias introduced during the dGNSS survey could be caused by the penetration of the dGNSS device into the soft snow cover by a few centimeters, which could explain some of the mean differences in elevation values between photogrammetry and dGNSS measurements.

Terrestrial laser scanning (TLS)
We compare the independently acquired TLS-derived snow depth (TLS winter minus TLS summer) with the ADSderived snow depth (Fig. 6a, b). In total we look at 55 272 pixels of 2 m resolution. It is hard to detect differences between the two snow depth products at first sight. All prominent snow features such as filled channels, cornices or areas where the snow has been blown away are clearly visible in both products. In the difference image between the two snow depth products, four regions with large deviations of up to 2 m stand out (marked by black circles in Fig. 6c). Three areas with significantly negative deviations (red, TLS higher than ADS) are located in small depressions. In these areas the incident angle of the laser beam is very flat, resulting in lower accuracies. The ADS sensor looks at these spots from a nadir view, producing more reliable snow depth values. On the ridge at the southern edge of the subset a large cornice was formed by wind during the winter (see Fig. 2c in the background). Because of the nadir-viewing angle, the ADS data set maps this cornice with snow depth values that are too large. The TLS sensor sees the overhanging cornice from below, producing better snow depth measurements than the ADS. However, the correlation analysis for the two snow depth measurement methods results in cor e = 0.94, the RMSE is 0.33 m and the NMAD 0.26 m. This proves the quality of the ADS snow depth measurements, especially concerning the complex, representative terrain of this subset (mean slope angle of 27 • , ranging from 0 to 81 • , elevations ranging from 2332 m to 2639 m a.s.l.).

Hand-measured plots
The comparison of the snow depth values derived from the ADS80 DSMs to the manual plot measurements is given in Table 3. In 3 out of the 15 plots, the snow depth exceeds the length of the avalanche probe (3.2 m), and the correct values could not be measured at all 25 points (measurements deeper than 3.2 m: plot1, 14; plot11, 5; plot 13, 5). The hand measurements could also be distorted by the penetration of the snow cover not being exactly vertical (especially in deep snow packs); by thick ice layers in the snowpack, which cannot be penetrated by the avalanche probe; by rough bedrock or by inaccuracies of the positioning of the dGNSS. Therefore, we average the 25 single measurements and compare the mean and standard deviation of an entire plot to the ADS80 DSM-based snow depth values (mean of all cells within the plot area). The RMSE is 0.35 m for the mean snow depth, and the standard deviation is 0.13 m over all plots. The NMAD is 0.22 (mean) and 0.06 m (SD). The correlation coefficient cor e for the mean snow depth is 0.92 and 0.81 for the standard deviation. If we eliminate the three plots (1, 11 and 13), which contain unreliable measurements, the RMSE is reduced to 0.19 (mean) and 0.11 (SD) and the NMAD to 0.18 m (mean) and 0.06 m (SD). The correlation coefficients shift to 0.95 (mean) and 0.76 (SD). The standard deviation is underestimated by the DSM ADS -derived snow depth values due to the smoothing effect of the 2 m pixel size. However, these results indicate the feasibility of the proposed method for snow depth mapping.

Ground-penetrating radar (GPR)
To allow a comparison between GPR snow depth measurements and the ADS measurements, we assigned all individ- ual 18 136 GPR point measurements to the 2 × 2 m ADS raster and calculated the mean of all GPR values within each cell, resulting in 1522 cells with GPR-based comparison data. The variability of the GPR snow depth within these cells amounted to between 0.1 and 0.3 m. Parts of the GPR data have been obtained close to taller vegetation such as trees and bushes. However, heavily affected measurements have been masked out before comparison, as ADS data cannot represent snow depth under forest canopy.
Comparing GPR to ADS data results in an overall RMSE of 0.43 m and an NMAD of 0.36 m. This is approx. 0.1 m worse compared to the reference data sets acquired at the Wannengrat area. The overall correlation coefficient between both data sets is 0.45 (Fig. 7a) only; note, however, that the GPR data set features a significantly lower range in snow depth when compared to the TLS data set (Fig. 6), mainly because it was acquired at the valley bottom. When analyzing different segments of the GPR data set we find considerable differences. While the correlation is acceptable for individual GPR segments that feature large snow depth variability (Fig. 7b), it appears less favorable for GPR segments with a small variability in snow depth (Fig. 7c). By comparing the profiles of the snow depth values along the two segments numbers 1 and 5 (Fig. 7d, e), we find the ADS values to be too low over large parts of the transects.  ADS snow depth values are too low. In the profile number 5 (Fig. 7d), the first 200 m of the segment covers meadow. The second part is on a road, running along a slope. While the GPR snow depth values remain quite constant, the ADS snow depth values show a large variability. While all GPR measurements are made strictly on the road, the 2 × 2 m ADS pixels include adjacent areas on both sides of the road which could be nearly snow-free or covered by deep snow at the edge of the road. Another explanation for the worse accordance between GPR and ADS snow depth values might be the greater distance of the ADS sensor to the ground. While the Wannengrat reference data sets were collected at an altitude of approximately 2400 m a.s.l., the valley bottom of the Dischma, where the GPR data were collected, has an elevation of approximately 1600 m a.s.l. This results in a coarser effective ground sampling distance (GSD) and therefore in a lower accuracy of the corresponding ADS data set. This finding indicates that the spatial resolution of input imagery matters for the accuracy of the resulting snow depth estimates. Table 5. Price ranges in thousands of Swiss Franks (kCHF) and relative differences derived from quotations of three independent data providers. We asked to cover the investigation area of this paper (145 km 2 ) with airborne laser scanning (ALS) and digital photogrammetry with a spatial resolution of 2 m and a vertical accuracy of approx. 30 cm.

Discussion
Compared to airborne laser scanning the proposed method is expected to be slightly less accurate but more economical if large areas (> 100 km 2 ) have to be covered repeatedly. To assess the economic advantage of digital photogrammetry, we requested quotations from three independent data providers offering digital surface models generated by airborne laser scanning and digital photogrammetry to cover the investigation area of this study (145 km 2 ). We asked for a GSD of 2 m for the final DSM and a vertical accuracy of approx. 30 cm (RMSE). Table 5 presents an overview of the answers we received. Digital photogrammetry is 40-50 % more economical than ALS in data acquisition, mainly because of the more efficient flight pattern resulting in reduced flight time for a given area. Data processing is 10 to 40 % more economical, resulting in a significant total price reduction of 25 to 37 %. Now, the successor sensor Leica ADS100 is available, incorporating almost twice as many detectors than the ADS80 sensor, resulting in a better spatial resolution for the same flying height above ground. Digital photogrammetric DSMs can be generated using unmanned aerial vehicles (UAVs) flying close to the ground and producing higher spatial-resolution imagery (Mancini et al., 2013) on the order of centimeters, resulting in more accurate (better than 10 cm in vertical direction) and much more economical snow depth maps. However, the feasibility of UAVs in high-alpine terrain has to be investigated further. Winged UAVs might not be stable enough under windy conditions, which usually prevail in alpine terrain. Furthermore, it might be difficult to find suitable starting and landing spots due to the rough terrain. UAVs with rotors are much more stable and can acquire data under windy conditions if the wind is not gusty. However, they have very limited flight times due to high energy consumption and the batteries have to be changed very often (approx. every 5 min). UAVs with rotors are not yet able to efficiently cover areas larger then a few square kilometers in alpine conditions and the risk of crashing the UAV in rocky terrain is high.
Challenging for image correlation on snow-covered terrain are the big spectral differences of surface cover properties between bright snow-covered slopes and rocky ter- -High precision -Very limited coverage -Point measurements -Extreme terrain inaccessible -Need for being in the field -Expensive device rain in shadow. If terrain properties change within short distances, the probability of big outliers or even complete failures of image-matching rises. We modeled only 0.25 km 2 per step to decrease these differences within the correlated images. With this approach massive failure of image matching could mostly be averted. For some tiles, issues with big outliers remained, showing a certain limitation to the modeling of snow-covered areas with the image correlation software used. For future investigations the choice of more advanced image correlation algorithms, such as methods of the semi-global matching family, has potential to solve part of this limitation. The modeling of steep slopes (> 50 • ) using image-matching techniques is not accurate mainly due to the small footprint of the sensor . But because snow accumulation is reduced on such steep slopes (Schweizer et al., 2008;Fischer et al., 2011), these areas are less important for applications in hydrology and avalanche science. The methodology proposed does not work in forested terrain or in regions covered by scrubs. Therefore, these areas were masked out prior to the snow map calcula-tion. This is not possible for areas with high grass in summer; therefore, we clearly underestimate the snow depth with the ADS data in such areas (see Fig. 7d, e). In forested terrain ALS has a strong advantage compared to photogrammetry because the terrain surface can be measured between the trees if the forest cover is not too dense. The accuracy of final DSM products depends heavily on the image strip orientation quality. Here we faced two major limitations: (a) we could gather only a small number of reference points, measured with high accuracy in x, y and z, and b) in areas covered by deep snow without anthropogenic signs visible, the recognition of clearly identifiable reference points is sometimes almost impossible. Therefore, we see great potential to increase the quality of final products by collecting more accurately measured reference points and by marking reference points in remote parts of the covered area for upcoming data acquisition campaigns. But such fieldwork can be costly if several people have to be deployed in the field to cover large areas and different elevation levels in difficult terrain, reducing the economic advantage of photogrammetry. The results presented demonstrate the potential of digital photogrammetry for catchment-wide snow depth mapping. The extensive validation using independent data sets acquired simultaneously reveals an accuracy of approximately 30 cm (RMSE, NMAD), equivalent to ∼ 1 GSD of the input images (Table 4). Due to the high radiometric resolution of the images (12 bit) and the use of the near infrared band, the images were not saturated over bright, snow-covered areas and information could be acquired even in shadow. The image correlations work even over very homogeneous areas. Table 2 reveals almost the same correlation success with winter images compared to summer images. The resulting snow depth maps visualize the high spatial variability of snow depth even within short distances of a few meters. Snow traps for wind-blown snow, cornices and deposits from past avalanche events can be identified easily by high snow depth values of up to 15 m. In this paper we applied six different methodologies to map snow depth in high-alpine terrain. Table 6 lists the major strengths and weaknesses of these methods based on the experience of the authors. However, which method should be applied in a specific case depends on many different factors and should be evaluated with care.

The
We plan to acquire similar data sets at the end of upcoming winters for an interannual comparison of snow depth. This would also open the door for investigations into the representativeness of snow depth measurements at given points, for example at automated weather stations. Future comparisons between snow depth maps generated by lidar and digital photogrammetry will provide more detailed information on the specific strengths and weaknesses of the two methods.