Mapping snow depth in open alpine terrain from stereo satellite imagery

. To date, there is no deﬁnitive approach to map snow depth in mountainous areas from spaceborne sensors. Here, we examine the potential of very-high-resolution (VHR) optical stereo satellites to this purpose. Two triplets of 0.70 m resolution images were acquired by the Pléiades satellite over an open alpine catchment (14.5 km 2 ) under snow-free and snow-covered conditions. The open-source software Ame’s Stereo Pipeline (ASP) was used to match the stereo pairs without ground control points to generate raw photogrammetric clouds and to convert them into high-resolution digital elevation models (DEMs) at 1, 2, and 4 m resolutions. The DEM differences (dDEMs) were computed after 3-D coregistration, including a correction of a − 0.48 m vertical bias. The bias-corrected dDEM maps were compared to 451 snow-probe measurements. The results show a decimetric accuracy and precision in the Pléiades-derived snow depths. The median of the residuals is − 0.16 m, with a standard deviation (SD) of 0.58 m at a pixel size of 2 m. We compared the 2 m Pléiades dDEM to a 2 m dDEM that was based on a winged unmanned aircraft vehicle (UAV) photogrammetric survey that was performed on the same winter date over a portion of the catchment (3.1 km 2 ). The UAV-derived snow depth map exhibits the same patterns as the Pléiades-derived snow map, with a median of − 0.11 m and a SD of 0.62 m when compared to the snow-probe measurements. The Pléiades images beneﬁt from a very broad radiometric range (12 bits), allowing a high correlation success rate over the snow-covered areas. This study demonstrates the value of VHR stereo satellite imagery to map snow depth in remote mountainous areas even when no ﬁeld data are available.


Introduction
The seasonal snow cover in mountainous areas sustains mountain glaciers, alters frozen ground through its insulating effect, and plays a major role in mountainous ecosystems and plant survival (Keller et al., 2005). Snow cover is important for hydropower production, irrigation, urban supply, risk assessment, and recreation (Barnett et al., 2005). The seasonal snow on the ground can be characterized by various metrics, including the snow-covered area (SCA), the snow height (HS), the snow density ρ s , and the snow water equivalent (SWE) (Fierz et al., 2009). A key moment to evaluate the snow cover as a water resource in an alpine catchment is the accumulation peak, when the SWE reaches its maximum value. In the Pyrenees, the accumulation peak associated to the persistent snow pack is generally between March and April (López-Moreno and García-Ruiz, 2004;López-Moreno et al., 2013). Even for small mountain catchments with areas of a few square kilometres, the spatial variability of the snow height and water equivalent is high due to a number of different processes: elevation gradient of snowfall, preferential deposition of precipitation induced by topographical effects, redistribution of snow by wind, sloughing, and avalanching (Grünewald et al., 2014).
Various techniques exist to monitor the HS and SWE at specific locations. The snow course is a standard protocol that is used to measure the SWE in the catchment areas of dams in many countries (DeWalle and Rango, 2008). An operator measures the HS with a snow probe at a number of predefined waypoints. The survey is repeated a few times during winter to obtain the amount of accumulated snow be-fore spring freshets. The snow density is also estimated during a snow course, but this measurement is not conducted at every point because coring and weighing the snowpack takes a longer time than snow depth measurements (Sturm et al., 2010). In addition, many studies showed that the snow density is much less variable in space than the snow depth (Pomeroy and Gray, 1995;Marchand and Killingtveit, 2005;Jonas et al., 2009;López-Moreno et al., 2013). The snow course remains a time-consuming task which can be dangerous because of the risk of avalanches. Even in small catchments, this approach does not enable field operators to routinely sample the entire catchment area. Automatic measurements that are based on snow pillows, sonic rangers, and nuclear snow gauges are widely used in addition to manual measurements (Egli et al., 2009). GPS interferometry has been recently used to measure the HS at decimetre resolution (Larson et al., 2009;Gutmann et al., 2012) and could represent an alternative in snow-dominated regions, where geodetic GPS receivers are already operating for various purposes (e.g. plate deformation or weather monitoring). All these point-scale observations must be extrapolated by using statistical models and/or remotely sensed data (e.g. Martinec and Rango, 1981;Luce et al., 1999;Molotch et al., 2005;López-Moreno and Nogués-Bravo, 2006;Grünewald et al., 2013).
Remote sensing techniques are particularly suitable for monitoring snowpacks at the catchment scale under satisfactory safety conditions. Recent advances in the fundamental understanding of the distribution of mountain snow depth have been achieved through airborne lidar (light detection and ranging) campaigns (Deems and Painter, 2006;Deems et al., 2013). Lidar provides an accurate measurement of the snow depth with a very high spatial resolution, which is perfectly suited for monitoring snowpacks in mountainous areas, including forested areas (Hopkinson and Sitar, 2004;Grünewald et al., 2013). The vertical accuracy ranges from centimetres to a few decimetres (Grunewald and Scheithauer, 2010;Deems et al., 2013). This technique is being extended for operational purposes in the USA (Painter and Berisford, 2014; Airborne Snow Observatory, http://aso.jpl.nasa.gov/). However, airplane surveys are costly and do not allow global coverage. Terrestrial laser scanners (TLSs) are relatively less expensive than an airborne laser scanner (ALS) and offer comparable resolution and accuracy at mid-range distances (up to 300-500 m) (Prokop, 2008;Grünewald et al., 2010). However, holes in the dataset caused by convex landforms such as hills or moraines may limit the spatial covering of the TLS acquisition (Bühler et al., 2016). The beam divergence of TLS is generally lower over steep terrain but coarser over flat areas, which highlights the complementary nature of both ALS and TLS techniques in mountainous terrain.
Satellite snow cover observations, including operational applications, have been performed for many decades (e.g. Rango A, 1976Rango A, , 1994Dietz et al., 2012). Numerous satellitederived products exist at the global scale (Frei et al., 2012). Snow-covered area maps are routinely produced from visible or near-visible bands (e.g. MODIS products; Hall et al., 2002). When combined with a distributed snowmelt model, the SWE can be reconstructed from the monitoring of the SCA, provided that the last day of snow on the ground is known (e.g. Molotch and Margulis, 2008). An important limitation of this method for operational purposes is that it requires the user to wait until the end of the snow season.
Microwave remote sensing techniques have been demonstrated to be effective for monitoring snowpack-related metrics (SCA, HS, SWE, wet/dry state; Sokol et al., 2003). Numerous spaceborne radiometers with appropriate frequency channels have been in orbit since the 1960s (e.g. SMMR 1978;SSM/I 1987;AMSR-E 2002). However, the application of passive microwaves to snowpack monitoring in alpine regions is limited by the coarse resolution of spaceborne sensors, which are typically 10-25 km (Clifford, 2010), and the presence of liquid water in the snowpack. Another limitation is the SWE threshold, which impedes SWE retrieval for deep snowpacks (Dozier et al., 2016; > 0.15 m-0.20 m w.eq.). Several attempts have been made to retrieve spatially distributed HS or SWE data from space by radar imagery (Papa et al., 2002;Leinss et al., 2014;Rott et al., 2014;Dedieu et al., 2014). However, the optimal frequency channels (Ku, Ka) are still absent from current SAR satellites. Radar can operate even under cloudy conditions, but snow penetration from band X or band C complicates these measurements, and large areas may remain masked because of the oblique view of the imager.
Satellite altimetry (e.g. ICESat) could potentially accurately determine the snow depth, but the large footprint is not optimal for small alpine catchments. Errors may arise from signal saturation and beam penetration. To date, there is no definitive approach to map snow depth in mountainous areas from spaceborne sensors (Lettenmaier et al., 2015).
The objective of this paper is to assess the potential of stereo images from a very-high-resolution (VHR) satellite to retrieve the snow depth. Recently, digital elevation models (DEMs) that were derived from Pléiades satellites have been assessed over various types of surfaces, such as end-ofsummer glacier surfaces (Marti et al., 2014;Berthier et al., 2014), lake deposits and dunes (Schuster et al., 2014;Lucas et al., 2015), or landslide areas (Stumpf et al., 2014;Lacroix et al., 2015). Pléiades-derived DEMs exhibited sub-metre ac- curacy in the elevation of these rugged topographies, which opens the possibility to sense the snow depth from space by subtracting a DEM that was obtained under snow-free conditions from a DEM that was obtained near the peak of snow accumulation. This study's goals are as follows.
-Generate, co-register, and differentiate two Pléiades DEMs in a small mountainous catchment without ground control points: a snow-free DEM and a DEM that was acquired near the snow accumulation peak.
-Assess the quality and accuracy of the difference in the Pléiades DEMs (dDEMs) based on two datasets: (i) snow-probe measurements and (ii) another dDEM generated from two UAV surveys.
-Discuss the influence of the topography and land cover on the residuals between the Pléiades dDEM and the snow-probe measurements.

Study site
The study area is the Bassiès catchment (14.5 km 2 ), which is an open alpine terrain in the north-eastern Pyrenees (Fig. 1). Bassiès is one of the main sub-basins of the Upper Vicdessos Valley, which has a long history of hydropower production (Taillefer, 1939;Antoine et al., 2012). The elevation ranges between 1156 and 2676 m a.s.l. (median elevation 1659 m) with a contrasted relief: while steep slopes delimit the watershed, the valley bottom is rather flat and exhibits gentle slopes in its central part. The catchment is ungauged, but the streamflow at the outlet is diverted toward a hydropower plant operated by "Electricité de France" (EDF). The average annual temperature in the area is 6.6 • C and the mean annual precipitation is 1640 mm, of which at least 30 % falls as snow . The snow season generally starts in November-December and ends in May-June (Fig. 2). The catchment is 65 % covered by subalpine meadow and 25 % by vegetation-free rock and bare soils. The last 10 % is composed of intermediate vegetation (scattered short-conifer, 5 %), forest (2 %), and water surfaces (lakes and rivers, 3 %) (see the Supplement for the land-cover map).

Pléiades images
The satellites Pléiades-1A and 1B fly on the same near-polar sun-synchronous orbit at an altitude of 694 km with a 180 • phase and descending node at 10:30 a.m. The CCD optical sensors acquire images in push broom mode by using 5 × 6000 pixel arrays and a maximum of 20 integration lines (TDI) for the panchromatic band (480-830 nm) (Poli et al., 2015). The system can achieve stereoscopic imaging with an additional quasi-vertical image (tri-stereoscopy), which is particularly suited for dense urban and mountainous areas. The tri-stereo mode can combine three stereo pairs to generate multiple DEMs, namely front/nadir, nadir/back, and front/back stereo pairs. The Pléiades's pixel depth at acquisition is 12 bits, and the panchromatic images have an initial resolution of 0.70 m but are oversampled at 0.50 m before image delivery by a post-processing algorithm that was implemented by the French Space Agency (CNES). Two Pléiades triplets were acquired over the Bassiès catchment, which is the area of interest in this study (Table 1). The snow-free acquisitions were programmed on 26 October 2014 (10:53:10, 10:53:31, and 10:53:52 LT). Each snow-free image covered a surface area of approximately 117 km 2 , which was centred on the Bassiès catchment. The images were acquired with viewing angles of 11.9, 0.7, and −10.9 • in the along-track direction with respect to the nadir and −4.8, −4.3, and −3.7 • in the across-track direction. Consequently, the base to height (B / H ) ratios were 0.22 (front/nadir pair), 0.23 (nadir/back pair), and 0.45 (front/back pair). The northern slopes were exposed to large shadows (approximately 10 % of the catchment area) and exhibit poor image contrast because of the sun's position during autumn (sun elevation 34 • , azimuth 167 • ). No saturation or cloudiness were observed in the snow-free images.
The second triplet was acquired on 11 March 2015 (10:56:42, 10:57:03, and 10:57:27 LT), when the snow accumulation was presumably close to its maximum peak. Each winter image covered a surface area of approximately 115 km 2 , centred on the Bassiès catchment, as achieved for snow-free images. The images were acquired with viewing angles of 10.5, −0.7, and −14 • in the along-track direction with respect to the nadir and with viewing angles of 0.4, −2.7, and −6.4 • in the across-track direction. Consequently, the estimated B / H were 0.22 (front/nadir pair), 0.26 (nadir/back pair), and 0.48 (front/back pair). The images had a very low cloudiness (<2 %). Saturated zones represented less than 3 % of the images and were located almost exclusively along the southern-exposed slopes. The northern slopes also exhibited abundant shadows (approximately 5 % of the catchment area), but these shaded areas with low contrast were less extensive than those in the snow-free acquisitions (sun elevation 41 • , azimuth 157 • ).

UAV images
Two winged-UAV photogrammetric surveys were performed over a central subset of the Bassiès catchment (3.15 km 2 ) to determine the snow depth by DEM differencing ( The flight altitude was maintained at approximately 150 m, which provided a mean ground sampling distance (GSD) from 0.10 to 0.40 m. Both the winter and snow-free acquisitions were achieved under very clear sky conditions. Onboard RTK corrections were performed at a frequency of 20 hz. UAV orientation was improved during the winter survey through the use of a GPS base, which was installed on the flat dropping zone of the mountain refuge during the survey. Five georeferenced ground targets were placed in the valley bottom during the summer and identified on the UAV images to improve the absolute positioning accuracy.

Snow probing
We collected up to 501 hand-probed depth measurements on 10 March 2015, at the time of the UAV survey, and 1 day before the Pléiades acquisition (Table 1). Because of the limited available time on the field, we attempted to cover an area that could represent of a large part of the catchment topography. The distance between each sample ranged from 10 to 30 m. We used two types of snow probes with lengths of 2.2 and 3.2 m. The snow probing coordinates were recorded by using a differential GPS (DGPS) with a mean of 15 acquisitions (one per second) per probe location. We used the Trimble Geo XH 2008 (GPS) and Geo XH 6000 (GPS and Glonass). Post-treatment corrections were collected from a base that was 21 km away, specifically the French "Réseau Géodésique Permanent" (RGP) network (base "Mercus-Garrabet"). This process enabled us to achieve 0.1 m accuracy in the horizontal and vertical directions of the snow probing locations.

Land-cover map
A 2008 land-cover map, which was updated by a field survey in July 2015, was generated through an object-based approach and expert interpretation of aerial photographs (Sheeren et al., 2012;Houet et al., 2015) (see the Supplement for the land-cover map). The vegetation types were aggregated into seven classes to reflect the type of land cover that may influence the comparison between the Pléiades dDEM and the snow-probe measurements: mineral surfaces (bare soil and rocks), water surfaces (rivers and lakes), peatland, low grass (rangeland, grassland, and subalpine meadows), shrubs, trees (conifer and deciduous), and unknown.

Production of DEMs, orthoimages, and dDEMs from Pléiades images
A tri-stereoscopic acquisition was considered to (i) limit the areas potentially masked by the rugged topography of the studied catchment, (ii) improve the correlation by providing different B / H ratios, and (iii) obtain a nearly nadir image to improve the orthorectification process and accuracy of the absolute co-registration offset. Snow-free and winter Pléiades DEMs were generated from the image triplets through the Ames Stereo Pipeline (ASP,version 2.4.8.), an open-source automated stereogrammetry software by NASA (Broxton and Edwards, 2008;Moratto and Broxton, 2010;Willis et al., 2015) (Fig. 3). The ASP was primarily designed to create DEMs of ice and barerock surfaces. The ASP supports any Earth imagery that uses the rational polynomial coefficient (RPC) camera model format. The RPC model is an analytical model, provided here as meta-data by Airbus Defense and Space (ADS), which gives a relationship between the image coordinates and the ground coordinates with z as the height above an ellipsoid and includes both a direct model (image to ground) and an indirect model (ground to image) (ASTRIUM, 2012). Spatiotriangulation was based on the RPC model, which was refined from an automated tie points generation without including ground control points (GCPs). We parameterized the ASP to project the images into an epipolar geometry to reduce the search range before the correlation (normalized cross correlation) and triangulation steps. We generated three point clouds from the three stereoscopic pairs from the stereo command and merged them. The DEMs were rasterized at 1, 2, and 4 m cell sizes from the merged point cloud through the point2dem command. Resolutions lower than 1 m are not relevant given the original image resolution and resolutions higher than 4 m will smooth out most of the interesting snow depth features. The elevation values at a given grid point were obtained as a weighted average of the elevations of all points in the cloud within the search radius of the grid point, with the Gaussian curve as weighting function (see the Supplement for the ASP's parameters) (NASA, 2015).
Snow-free and winter DEMs at 4 m were horizontally coregistered by iteratively shifting the winter DEM with respect to the summer DEM (reference) by minimizing the standard deviation (SD) of the elevation difference distribution (Berthier et al., 2007). The final horizontal shifts were applied to the winter DEM were −5.2 m in northing and +2.8 m in easting (Table 2). We obtained similar results by computing the optimal shift at 1 and 2 m resolution. This result is consistent with the expected localization precision that was provided by the RPCs from the Pléiades images. Without GCPs, the horizontal location accuracy of the images was estimated at 8.5 m for a circular error at a confidence level of 90 % (CE90) for Pléiades-1A and 4.5 m for Pléiades-1B (Lebegue et al., 2010;Gleyzes et al., 2013). The same shift was applied to the 2 and 1 m winter DEMs.
Winter and snow-free nadir images were rectified at 1 m resolution from their respective DEMs, before coregistration. By picking six widespread corresponding points on the snow-free and winter images, the mean shifts were −5.2 m in northing (SD = 0.7 m) and +3.2 m in easting (SD = 0.5 m), which are consistent with shifts from the DEM co-registration technique. The low SD values indicate that  the horizontal shift was almost constant in the image. A classification of the image pixels into snow and snow-free classes based on intensity thresholds was performed on the winter orthoimage (Table 3). Two intensity thresholds were visually adjusted in order to treat specifically the case of the shaded snow surfaces from the general case. dDEMs were produced at 1, 2, and 4 m spatial resolution by subtracting the snow-free DEM from the winter DEM on a pixel by pixel basis: where Z w is the pixel value in the winter DEM, and Z s is the pixel value in the snow-free DEM. An absolute horizontal shift in the Pléiades DEMs was estimated from six widespread points that were identified on an aerial orthophoto from IGN (Institut national de l'information géographique et forestière), which presents an absolute accuracy of approximately 2 m. The shift between the snow-free Pléiades orthoimage and the IGN orthophoto was +3 m (SD = 0.38 m) in northing and −0.8 m Table 2. Summary of the different co-registrations and bias corrections performed to produce the Pléiades and the UAV DEM and dDEM maps. SD means standard deviation. The term workflow metrics refer to the data presented in Fig. 3 Table 3. Percentage of potential outliers and no data in the dDEM values, considering the catchment area, the snow-covered area of the catchment, and the snow-covered area of the catchment located out of the shadows due to the high cliffs (here indicated as "sunny snow").  Table 2). The dDEMs were then shifted based on this absolute horizontal offset to be consistent with the DGPS and the georeferenced snow-probe measurements. Then, we removed a constant vertical bias from Z 0 (Eq. 1) to obtain the final dDEMs: where b is a constant vertical bias, which is determined from a unique, stable, and flat area of the satellite winter and autumn images that is easy to interpret. We chose to evaluate b from a snow-free football field in the image that was 5 km from the mountain refuge (Fig. 1). The value of b was assumed to be equal to the median of the dDEM distribution on the football field (Table 2). After this bias correction, dDEM pixels with negative values were classified as "no data", which include 8 to 10 pixels that correspond to a snowprobe measurement (Table 4). We classified the percentage of negative dDEM pixel values over the Bassiès catchment ac-cording to the presence of snow and excluded shadow areas from steep rocks or cliffs. Verifying whether a vertical bias that is measured over a small portion of a dDEM at low elevations (football field) can be used to correct an entire dDEM is very important. To test this assumption, we extracted 78 widespread values from the 2 m Pléiades dDEM before bias correction (Eq. 1). We photointerpreted these points on snow-free rock areas, roads, or bare soil in the absolute georeferenced winter orthoimage by avoiding the steepest slopes (< 30 • ) and by covering a large elevation range (790-2510 m). We did not use this information to remove the bias because we aimed to evaluate a simple workflow that could become operational (Table 2).

Production of UAV DEMs and dDEMs from the UAV images
UAV DEMs were generated from the overlapping drone images by using the PIX4D software, which uses a structurefrom-motion algorithm (Westoby et al., 2012). The focal www.the-cryosphere.net/10/1361/2016/ The Cryosphere, 10, 1361-1380, 2016 length as well as the lens distortion modelling parameters of the cameras were adjusted for each flight during the automatic PIX4D workflow. Five GCPs were available in summer to improve the snow-free images orientation. Except the position of the GPS-base, no GCPs were available during the winter survey; thus the winter images were co-registered to the summer images to improve their orientation. Generated point clouds were rasterized at 0.1, 1, and 2 m cell sizes for both the snow-free and winter DEMs. Subsequently, 0.1, 1, and 2 m dDEMs were obtained by differencing the corresponding snow-free DEM from the winter DEM. The UAVimage acquisition, the UAV-image processing, and the UAV DEM generation were performed by a private company (Table 1).
After an initial comparison with the snow-probe measurements, a marked planar bias-oriented SW-NE was identified on the dDEMs. Comparing the winter UAV DEM values to the winter DGPS measurements (N = 343) showed that the bias resulted from a bad stereo orientation, which led to some deformations in the winter DEM. To correct that bias, we extracted 353 widespread values from the 0.1 m UAV dDEM at locations where the snow depth was supposed to be 0 based on the winter orthoimage (emerging bare rock). We generated trend surfaces of order 1, 2, and 3 based on these values and subtracted them from the 0.1 m dDEM. The trend surfaces are defined by fitting a polynomial function to the sample points. Here we tried polynomial functions of order 1, 2, and 3. This processing was done using ArcGIS Spatial Analyst toolbox. The results improved significantly at each polynomial order, so we corrected the dDEM with the order 3 trend, which best fit the dDEM values from the emerging bare rocks (root mean square error (RMSE) before trend removal is 0.96 m; for order 1, RMSE is 0.44 m; for order 2, RMSE is 0.39 m; for order 3, RMSE is 0.34 m). The results presented below are based on the de-trended dDEM values at each pixel resolution. An extra point that was located on the flat dropping zone of the mountain refuge was used to correct a constant bias after the trend removal (0.1 m: +0.33 m; 1 m: +0.43; 2 m: +0.41 m).

Pléiades and UAV dDEM assessments and comparison
We compared the Pléiades dDEM at 1, 2, and 4 m resolutions and the UAV dDEM at 0.1, 1, and 2 m resolution to the snowprobe measurements. We calculated the values of the residual vector R Z as follows: Z i is the subset of the dDEM values, where Z ≥ 0 after bias correction (Eq. 2), which were sampled by snow probing. HS are the snow-probe measurements. We considered that the measurements from the snow probes had a random error of σ probe = 0.15 m but did not introduce a systematic error term.
The metrics that were used to describe the quality of the dDEMs were the percentage of no-data values after the stereo processing and the statistics of R Z : (i) the mean and the median, which were used to evaluate the vertical accuracy of the dDEMs, and (ii) the SD and the normalized median absolute deviation (NMAD), which were used to characterize its vertical precision. The NMAD is a metric for the dispersion of data that are not as sensitive to outliers as the SD (Höhle and Höhle, 2009): where m R Z is the median of the residuals. We also assessed how the Z and HS values correlate through a rank correlation method. We used the Spearman correlation factor, called cor s , which is sensitive to neither the presence of outliers nor the existence of nonlinear correlations (Chueca et al., 2007;Borradaile, 2013).
The snow depth was greater than the snow-probe length for 50 occurrences. These cases where the operator did not reach the ground were excluded from these statistics and were only exploited as binary information to assess the dDEMs (see the Supplement).
We snapped and subtracted the 2 m UAV dDEM from the 2 m Pléiades dDEM. We visually compared both dDEM maps and the dDEM differences. We performed a SW-NE transect (1.6 km long) and compared the dDEM values along that transect.

Photogrammetric processing
We calculated the density of the summer and winter raw point clouds that were generated during the correlation process based on the front nadir/stereo pair (Fig. 3). The Pléiades panchromatic images had a pixel size of 0.5 m, so a mean density of 4 points per square metre would indicate a correlation success at the minimum achievable interval. Areas with lower density values require a higher search range in the interpolation of the raster cell value from the point cloud.

DEM contributions
To identify whether the final systematic and random errors were due to the snow-free DEM or the winter DEM, we computed two distinct residuals terms for the 2 m Pléiades dDEM as follows: We assume all the random errors to be uncorrelated.  Figure 4. Distribution of the snow-probe sampling, according to the snow depth, the elevation, the slope, the aspect, the curvature, and the land-cover classes. On the whole snow-probe sampling (N = 501), the snow probe did not reach the ground for 50 occurrences.
snow-probe locations. This term has a random error from both uncertainties in the snow probe and the DGPS. Hence, two random error terms exist because of the DGPS and snowprobe measurements (σ DGPS+probe = 0.18 m; see the Supplement for more details on the random error calculation).

Snow height, topography, and land-cover influences
Various factors limit the acquisition of snow-probe measurements, such as exposure to avalanches, human mobility on this challenging terrain, and the available time on the field. The snow depths that were obtained with the snow probes ranged from 0 to 3.2 m (Fig. 4). We assessed the influence of HS on the residuals between the dDEM and HS (Eq. 1). The snow heights from the snow probes do not represent the entire topographic variability of the catchment (Fig. 4). Here, we summarize the different ranges of the main topographical variables that are associated with the snow-probe data: . Therefore, we did not assess the residuals' distribution (Eq. 1) according to the elevation.
-The slope, which was derived from the 2 m snow-free Pléiades DEM, was associated with the snow-probe measurements and ranged continuously from almost 0 • (flat areas) to 25 • . A dozen snow-probe values were recorded in steeper zones but were not considered as statistically representative. The median slope of the catchment was 26 • , and a variety of slope values are present in the catchment, from flat area to cliff.
-The different aspect classes were well sampled during the snow-probe survey.
-The snow depth sampling range according to the curvature was quite limited because of the difficulty in performing snow probing in marked convex or concave areas. Therefore, we did not assess the residuals' distribution (Eq. 1) according to the curvature.
The distribution of the residuals between the dDEM and HS values was analysed according to the different landcover classes. The land-cover classes in the snow-probe data were minerals (12 %), water surfaces (4 %), low grass (32 %), shrubs (33 %), and peatland (19 %). The peatland class is overrepresented and mineral surfaces are underrepresented in the probe dataset with respect to the Bassiès catchment area.

Contribution of the tri-stereoscopy
To our knowledge, the added value of tri-stereoscopy relative to bi-stereoscopy has not been clearly established for an open alpine terrain. To provide a preliminary assessment of this contribution, we generated two seasonal DEMs from two individualized stereo pairs in both snow-free and winter cases. The first Pléiades pair consists of backwards and almost nadir images, and the second pair consists of forward and almost nadir images. Consequently, we generated a dDEM map for each stereo pair. We compared these dDEMs to the snowprobe measurements.

Pléiades dDEM assessments
The snow-free DEM and winter DEM are shown in Fig. 5. The small-scale topographic features are well captured by the high-spatial resolution of the DEMs. The winter DEM is characterized by a smoother texture. The distribution of the dDEM values (inset in Fig. 6) has the typical gamma or log-normal distribution shape that is reported in the literature (e.g. Winstral and Marks, 2014). Considering the whole Bassiès catchment, the mean of Z is 2.15 m and its SD is 1.72 m.
The Pléiades 2 m dDEM is composed of 1.7 % of no-data entries in the Bassiès catchment (2.4 and 1.2 % for the 1 Table 4. Statistics relative to the comparison between the Pléiades and the UAV dDEMs to the snow-probe measurements, according to the pixel resolution. Significant correlations (p values < 0.05) are marked with asterisks. NMAD means normalized median absolute deviation (Höhle and Höhle, 2009  and 4 m dDEMs respectively) ( Table 3). These no-data entries originate from data gaps in the raw points clouds, which are produced by the ASP before rasterization. The comparison with the snow-probe data indicates that the Pléiades dDEMs are consistent with the snow depth measurements (Table 4). The median values of the residuals distribution R Z are relatively low (between −0.12 and −0.16 m) and close to the mean of the distribution at each pixel resolution (±0.05 m between the median and mean). A slight influence from the pixel size is present (Table 4). For our validation dataset, the 2 m Pléiades dDEM exhibits slightly better precision and accuracy. For this dDEM, the SD is 0.58 m and the NMAD is 0.45 m. The Z and HS datasets are significantly correlated at each pixel resolution (cor s ranges between 0.67 and 0.72). The linear regression between the dDEM values and the snow-probe measurements is close to the 1 : 1 line ( Z 2m = 0.90 · HS) (Fig. 7). Figure 8 shows the spatial distribution of the snow depth measurements and the residuals of the 2 m dDEM. No obvious pattern is present in the residuals, although the absolute residuals are higher in the southern part of the surveyed area, where the slopes are the steepest (see Sect. 5.4.3).
Overall, the snow-probe dataset exhibits a low systematic error and is spatially homogeneously distributed.

UAV dDEM assessments
The SD and the NMAD indicate a decimetric random error, at each pixel resolution (SD 2 m = 0.62 m, NMAD 2 m = 0.35 m). The Z and HS values are significantly correlated (mean cor s 0.79). The median value of the residual distribution R Z is slightly negative and ranges from −0.07 to −0.15 m according to the pixel size. Figure 8 shows the spatial distribution of the snow depth measurements and the residuals of the 2 m UAV dDEM. From this map, no obvious pattern is present in the residuals.

Comparison of the Pléiades and UAV dDEMs
The Spearman correlation factor cor s between the Pléiades and UAV dDEM values is 0.62 and significant at 95 % confidence (N = 527.10 3 , number of values of the sample size). The datasets were not co-registered. The comparison between the 2 m Pléiades and the 2 m UAV dDEMs is characterized by a residual distribution with a median of −0.14 m (mean −0.06 m), an SD of 1.47 m, and an NMAD of 0.78 m.
The 2 m Pléiades and UAV dDEM maps exhibit very similar patterns (Fig. 9). Similar snow features are identifiable in both dDEM maps, such as a marked over-accumulation of snow along a topographic ridge that stretches from the refuge to the lake, snow traps for wind-blown snow and snow cornices. These features are also observable in the terrestrial photography (Fig. 2). A transect over a common area that is covered by both the 2 m Pléiades and 2 m UAV dDEMs www.the-cryosphere.net/10/1361/2016/ The Cryosphere, 10, 1361-1380, 2016 Table 5. Statistics relative to the comparison between the 2m Pléiades dDEM (tri-stereo) and the snow-probe measurements, according to the snow depth, slope and aspect, and the land-cover classes. Significant correlations (p values < 0.05) are marked with asterisks. NMAD means normalized median absolute deviation (Höhle and Höhle, 2009 highlights the consistency in both Z variations. Over this transect, the SD of the residuals between the Pléiades and UAV dDEMs is 0.78 m and the median is −0.16 m.

Photogrammetric processes
The density values of the raw point clouds (pts. m −2 ) from the correlation process based on the front nadir/stereo pair are close to the maximum achievable value (4 pts. m −2 ) in both the winter and snow-free DEMs at the snow-probe locations (see the Supplement for the density maps, Figs. 2  and 3). Therefore, the dDEM assessment should not be influenced by the interpolation process that creates the raster DEMs at the first order.

Pléiades DEM assessment
We decompose the respective contributions from the snowfree and winter DEMs to the dDEM residuals (Eqs. 5 and 6, Fig. 10). The medians of R Z w and R Z s distributions are −0.91 and −0.25 m respectively, leading to a difference of This value is consistent with the median of −0.64 m for the R Z 0 distribution that was identified with the HS probe measurements before the bias correction (the bias that was identified on the football field was −0.48 m). The R Z w and R Z s values in Fig. 10 are corrected from the bias by removing the median. The SDs of R Z w and R Z s are 0.32 m and 0.66 m respectively. These estimations are consistent with the SD of the 2 m residual distribution R Z (0.58 m).

Influences of topography and land cover
The correlations between the residuals distribution and the snow depth or the terrain slope are weak but significant (0.3 and 0.26). The deviation of the residuals distribution R Z increases slightly with the slope. However, the number of snow-probe measurements varies by interval and thus limits the interpretation of the statistics (Table 5).
The snow-probe measurements associated to the low grass and peatland classes present the lower deviation in the residuals' distribution (SD = 0.49 and 0.51 m). The most important dispersions are associated to the mineral and the shrub classes (SD = 0.79 and 0.63 m). 6 Discussion

Production of DEMs and dDEMs from Pléiades images
The method that was proposed here is based on VHR satellite stereo imagery. The agility of the Pléiades satellites provides a wide range of B / H ratios, including small values, which are necessary for alpine topography. We programmed a B / H of 0.2 between two consecutive stereo pairs to improve the correlation success rate and limit the shading effect of topography. The snow-free and winter front/back pairs (B / H = 0.4) created less dense photogrammetric clouds. Thus, the number of no-data pixels would have increased in the final DEM for a bi-stereo acquisition that was based on a B / H of 0.4 instead of 0.2. The stereo-orientation from the RPC ancillary data was sufficient to adjust the relative orientation of the images prior to their projections in the epipolar geometry. The affine epipolar transform of both the left and right images is based on automated tie-point measurements, whose effect is equivalent to rotating the original cameras which took the pictures (NASA, 2015). The command bundle adjust could probably improve the relative stereo orientation. We did not intentionally use GCPs to avoid the need for a field survey in the workflow. We did not remove outliers from the 3-D triangulated point cloud, which could be done by parameterizing the ASP ("near and far universe-radius parameters"; see the Supplement). The images that show which pixels were matched by the stereo correlation, which are called "good pixel maps" in the ASP, highlight a significant correlation in both the snowfree and winter DEMs. For steep slopes and/or a limited density of raw photogrammetric clouds, the map projection of the images through the mapproject tool on a coarse DEM before the Stereo pre-processing stage of ASP could improve the correlation success. Another option could be the direct calculation of the distance between the snow-free and winter raw point clouds instead of a raster representation (Westoby et al., 2015;Passalacqua et al., 2015).
The statistics, which were calculated separately for both DEMs, highlight the better performance in the elevation determination of the snow-covered images compared to the snow-free images (Fig. 10). This observation could be due to the difficulty in treating micro-topography with the native GSD of Pléiades (0.7 m at nadir). Snow-covered areas offer a smoother surface compared to vegetated or stony snowfree surfaces. The results on bare rock may be directly connected to the slope influence because most of this type of surface is located on steep slopes (Table 5). In both the snowfree and winter acquisitions, the shadow areas were the most challenging for the correlation process and appeared as very noisy surfaces with more no-data entries because of the correlation failures and outliers, such as negative dDEM values after vertical bias removal. The resolution of 2 m presents the most favourable statistics according to our validation dataset and potentially highlights a good compromise between the horizontal accuracy and the smoothing of the snow height.
Snow areas under shadows from high cliffs constitute a large erroneous fraction of negative dDEM pixel values (Table 3). Together with emerging steep rock, these areas should be treated as no-data entries with a sufficient buffer to limit the uncertainties on the mean HS retrieval.
No snowfall occurred in the Bassiès catchment during the 20 h between the field survey and the Pléiades acquisition. Fresh snow probably may have complicated the correlation stage and increased the number of saturated pixels. During the triangulation stage, we did not exploit the multiview stereo possibility of the ASP (only available since version 2.5.0), which limited our correlation to successive pair matching. Berthier et al. (2014) showed for the Mont Blanc area that a simple combination of the different DEMs derived from the three images of a tri-stereo can reduce the percentage of data voids and slightly improve the precision of the merged DEM. In our case, we did not notice an improvement in the dDEM precision through the comparison with the snow-probe measurements (SD = 0.69 m for the tri-stereo

Residuals (dDEM-HS)
- dian = −0.12 m for tri-stereo; median = −0.54 and +0.13 m for bi-stereo). The medians were of opposite signs for the front/nadir and nadir/back stereo pairs, which may explain the median values for the tri-stereo. The density maps from the point clouds exhibited similar patterns, because the correlation failed for both stereo pairs in the shadow areas.

Comparison to the snow-probe measurements
The validation dataset was strongly limited by the measurement protocol. To cover the largest extent in a limited time, we did not apply an optimal sampling strategy to as-sess the entire snow depth variability at a plot scale, typically 10 m × 10 m (López-Moreno et al., 2011;Bühler et al., 2015). The dDEM pixel values were therefore assessed by a unique snow depth measurement, which could explain the modest correlation between the dDEM values and the snow-probe measurements (mean cor s = 0.7 for Pléiades). The snow probes were too short to measure the highest snow depth, and we only provided binary information in these cases (see the Supplement). We did not survey the highest crest where drifted snow accumulates, which led to the highest snow accumulations. Even with longer snow probes, sam- The pling the snow depth in these areas would not have been safe. Increases in slope have a clear influence on the magnitude of the dispersion of the residuals between the dDEMs and the snow-probe measurements. However, the snow-probe dataset was not sufficiently representative to determine the influence of the slope.

Comparison to the UAV dDEM
A bias was identified in the winter UAV DEM. We could remove this bias in the final UAV dDEM thanks to the snowfree bare rock areas, which provided a valuable opportunity to generate widespread vertical offsets. However, this strategy for bias correction has obvious limits, and identifying and correcting the sources of this bias would have been better. The RTK signal was repeatedly lost during the survey, which negatively affected the photographs' orientation. The acquisition mode of UAV photographs is largely "nonconvergent", which could also result in marked deformation (Westoby et al., 2015). Winged UAVs are potentially less sta-ble than UAVs with rotors (Bühler et al., 2016), although recent works have highlighted their great potential for snow mapping in high-alpine catchments even in relative windy conditions (Harder et al., 2016). We noted large mismatches between the Pléiades and UAV dDEM maps for steep slopes, which could be due to incorrect flight plans or lens calibration co-registration errors (James and Robson, 2014). Recent works based on UAV systems to map snow depth highlight much better performance than the results reported in this study (2 m UAV dDEM: SD = 0.62 m, NMAD = 0.35 m, median = −0.11 m; see Table 4). Jagt et al. (2015) used a DSLR camera mounted on a multi-rotor UAV platform to map the snow depth at a very high spatial resolution (GSD of 6.10 −3 m) over a small mountainous terrain In the case of our study, the DEM of the SCA was generated from a unique flight plan. Some problematic flights were reported by Harder et al. (2016) (5 from 43 flights for all sites, or 11.6 %) with DEMs showing an RMSE of up to 0.32 m. The results mentioned above were extracted from multiple surveys with well spread GCPs and more dedicated surveys. We did not use GCPs during the winter survey and only five GCPs in summer, not well spread (bottom of the valley only). According to Harder et al. (2016), GCPs are needed to achieve the sub-decimetric accuracy, and a bias correction may also be necessary. Furthermore, residuals of the comparison between the UAV dDEM and the HS manual snow measurements were not filtered (e.g. a statistic criteria like 1σ threshold, the land-cover classes, or the slope). Therefore, despite the discrepancies observed in this study, we consider that the UAV dDEM map was a valuable independent source to evaluate the Pléiades snow depth map because the comparison revealed similar snow depth patterns, while the random and systematic errors of both dDEMs are comparable.

Limitations and perspectives
The digital photogrammetric determination of snow depth in mountainous areas has been a longstanding issue (Cline, 1993(Cline, , 1994Ledwith and Lundén, 2001). Until recently, terrestrial and aerial photographs and optical satellites images have been used almost exclusively to determine the spatial distribution of SCAs. Identifying conjugate, ground control points, contrast, and lighting issues were the main factors that have impeded the production of DEMs of SCAs.
Recent works have highlighted the potential of airbornederived techniques to produce centimetric and decimetric vertical accuracy and precision in DEM generation over SCAs and in dDEM generation from snow-free and winter DEM differencing. Pléiades-derived snow heights do not have the same accuracy and precision compared to this stateof-the-art of digital aerial photography. The performance highlighted by the UAV system mentioned in the previous section are very satisfactory (Jagt et al., 2015;Bühler et al., 2016;Harder et al., 2016;De Michele et al., 2016). Nevertheless their spatial coverage is limited to several hectares. Lee and Jones (2008) have created a DEM over a SCA of a mountainous terrain in Australia from a high-spatial-resolution camera (GSD up to 0.05 m) and an enhanced radiometric dynamic (12 bits) on a GPS/inertial motion unit (IMU) airplane system. An assessment by 183 GPS measurements revealed a mean of the residuals of +0.14 m with an SD of 0.08 m. Bühler et al. (2015) employed an optoelectronic line scanner (ADS 80) that was mounted on an aeroplane to map the snow depth at 2 m resolution (GSD of 0.25 m) in the Swiss Alps. A comparison between the ADS and different individual HS measurements revealed both RMSE and NMAD of approximately 0.3 m, which is equivalent to GSD of 1 of the input images. Over the polar snow of Alaska, Nolan et al. (2015) generated dDEMs over rather flat areas from a consumergrade camera that was coupled to a dual-frequency GPS on a manned aircraft without the use of an IMU. The comparison of the dDEMs to 6000 snow-probe measurements highlighted an SD of the residuals of 0.1 m (GSD of 0.06 to 0.2 m). These techniques that are based on airborne platforms remain suitable if clouds are present above the flight altitude. However, these approaches present serious constraints absent from satellite acquisition: the need for a pilot, a ground operator, or the use of a specific sensor and an ad hoc installation.
In remote areas such as high mountain catchment, these requirements could seriously compromise the acquisition process.
Pléiades, along with GeoEye-1, WorldView-1, WorldView-2, and QuickBird, belongs to class 6 satellites (GSD of 0.40-0.75 m). The main limitation of the images that are derived from these satellites could be the surveying of large areas because of the relatively limited swath (20 km for Pléiades). The maximum length of Pléiades stereoscopic coverage from the same orbit with a B / H of 0.2 is 80 km for a stereo acquisition and 25 km for a tri-stereo acquisition (195 km and 80 km respectively for a B / H of 0.4) (Gleyzes et al., 2012). Considering a B / H of 0.2, areas of up to 1600 km 2 may be imaged repeatedly in any part of the world that is covered by the Pléiades satellite constellation. Pléiades images do not exhibit the best spatial resolution of this class. However, its main advantage is its pixel depth at acquisition of 12 bits, while other VHR sensors have a pixel depth at acquisition of 11 bits. With 4096 shades of grey by pixel instead of 2048, subtle nuances, especially at the beginning or end of the spectrum, can be distinguished. As for all optical sensors, the main drawback of the Pléiades constellation is the need for clear-sky or with limited cloud cover conditions to obtain suitable images. Snow-free images can be acquired over a large temporal window, and repeating these acquisitions each time a dDEM must be processed is unnecessary. Winter images are more constrained because the key moment to evaluate the snow cover height is the vicinity of the accumulation peak, which may span several weeks. However, the daily revisit interval of the Pléiades satellite constellation increases the possibility of obtaining cloud-free and valuable images. Winter datasets can also be acquired at the end of various winters for interannual comparisons of snow depth.
The method that was proposed here does not provide any information on the snow thickness under trees. The ALS remains the only technique to extract high-resolution HS information in forested terrain. In the study area, this point is not critical because most of the catchment is open terrain. In general, most of the snow in the Pyrenees accumulates above the tree line near 1600 m a.s.l. .
Despite the above-mentioned limitations and given the results of this first study, we believe that satellite photogrammetry is a promising alternative to recently developed techniques that are based on lidar or aerial digital photogrammetry to retrieve snow depth. This conclusion is especially true in areas where field or airborne campaigns are not feasible or too expensive and where the snow accumulation is significant (above 2 m). In glaciology, DEMs that are generated from optical stereos are often considered to be inaccurate in accumulation areas (Schiefer et al., 2007;Racoviteanu et al., 2010). However, Pléiades DEMs that are acquired at the beginning and end of accumulation seasons could be used to evaluate the seasonal components of the glacier mass balance (Berthier et al., 2014). In hydrology and water resource applications, there remains a substantial uncertainty on the final snow volume at the watershed scale that need to be better assessed. In our study site, the mean dDEM value in the Bassiès catchment area (14.45 km 2 ) was 2.15 m. The corresponding coefficient of variation (CV) value was 0.80 (CV is the ratio of the SD to the mean snow depth). This CV agrees with the classification that was proposed by Liston (2004) since it falls in the category 9 "mid-latitude, treeless mountain (e.g. Rocky Mountains, alpine)". In terms of accumulation, the 2011-2012 winter was very comparable to the 2014-2015 winter in the Bassiès catchment. According to a Météo-France meteorological reanalysis, the precipitation was 1130 mm over the hydrological year 2011-2012 and 1150 mm in 2014-2015. Szczypta et al. (2015) used a distributed snowpack model to simulate the snowpack and its temporal evolution on a regular grid over the Bassiès catchment at a spatial resolution of 25 m during the 2011-2012 snow season. At the accumulation peak, the mean monthly snow depth that was simulated over the entire catchment in April was 2.2 m. Although both mean values cannot be readily compared, the order of magnitude appears to be consistent with the mean dDEM value that was found for the 2014-2015 winter and was based on Pléiades data.

Conclusions
We generated a DEM difference map that was based on winter and snow-free tri-stereoscopic Pléiades satellite images. The comparison of this Pléiades dDEM map to 451 snowprobe measurements, which were collected simultaneously, shows that the snow height can be retrieved from space with decimetric systematic and random errors and a metric horizontal resolution at the scale of a small mountain watershed (14.5 km 2 ). The distribution of the residuals between the 2 m Pléiades dDEM values and the corresponding snow-probe measurements presents a median of −0.16 m and an SD of 0.58 m. An independent dDEM map was generated through a winged UAV photogrammetric survey on the same date based on a similar workflow. Despite some outliers, the UAV dDEM map was also successfully validated by the snowprobe measurements (median of the residuals is −0.11 m, SD is 0.62 m). The comparison between the 2 m Pléiades and the 2 m UAV dDEMs is characterized by a relatively scattered distribution of the residuals mainly due to some outliers in the UAV dDEM: median is −0.14 m (mean is −0.06 m), SD is 1.47 m, and NMAD is 0.78 m. The snow cover features that were obtained by Pléiades DEM differencing were consistent with those that were derived from the UAV acquisition. The correlation between the snow heights from both techniques is statistically significant, even though some discrepancies were present on the steepest slopes.
The accuracy might be insufficient in areas where the snowpack remains thin even at peak accumulation (North American prairies, semiarid mountains) and for the study of small-scales snow features like sastrugi or penitents. Further studies should focus on influences of the snow height, the topography, and the land cover on the accuracy of Pléiadesderived snow heights based on lidar-derived snow height maps. Our validation dataset limited the analysis to gentle slopes or relatively flat areas and snow heights up to 3.2 m. The shadows that are projected onto slopes create a lack of radiometric contrast in both snow-free and winter images and constitute an inherent limitation to optical sensors. Other limitations include obstructions by the forest canopy and cloud cover.
These results are promising because they introduce the possibility of retrieving the snow height at a metric horizontal resolution in remote mountainous areas that are difficult to access. Indeed, the processing of the Pléiades data does not require mandatory field data like ground control points, although such reference measurements are always highly desirable. An adjustment on a snow-free flat surface, which can be located kilometres apart and at lower elevations, is needed to correct a vertical bias in the Pléiades DEM difference. The size of the study area could vary from several square kilometres to several hundreds of square kilometres.
The Supplement related to this article is available online at doi:10.5194/tc-10-1361-2016-supplement.