Mapping snow depth in open alpine terrain from stereo satellite imagery

Abstract

To date, there is no definitive 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 $\mathrm{m}$ resolution images were acquired by the Pléiades satellite over an open alpine catchment (14.5 ${\mathrm{km}}^{\mathrm{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 $\mathrm{m}$ resolutions. The DEM differences (dDEMs) were computed after 3-D coregistration, including a correction of a $-$0.48 $\mathrm{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 $\mathrm{m}$, with a standard deviation (SD) of 0.58 $\mathrm{m}$ at a pixel size of 2 $\mathrm{m}$. We compared the 2 $\mathrm{m}$ Pléiades dDEM to a 2 $\mathrm{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 $\mathrm{km}$${}^{\mathrm{2}}$). The UAV-derived snow depth map exhibits the same patterns as the Pléiades-derived snow map, with a median of $-$0.11 $\mathrm{m}$ and a SD of 0.62 $\mathrm{m}$ when compared to the snow-probe measurements. The Pléiades images benefit 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 field 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 . Snow cover is important for hydropower production, irrigation, urban supply, risk assessment, and recreation . 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 ${\mathit{\rho }}_{\text{s}}$, and the snow water equivalent (SWE) . 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 . 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 .

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 . 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 before 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 . In addition, many studies showed that the snow density is much less variable in space than the snow depth . 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 . GPS interferometry has been recently used to measure the HS at decimetre resolution 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..

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 . 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 . The vertical accuracy ranges from centimetres to a few decimetres . This technique is being extended for operational purposes in the USA (; 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 $\mathrm{m}$) . However, holes in the dataset caused by convex landforms such as hills or moraines may limit the spatial covering of the TLS acquisition . 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.

Airborne and terrestrial photogrammetry has been investigated on snow surfaces since the 1960s . Nevertheless, their successful assessment has been achieved only recently . Airborne photogrammetry represents a relatively inexpensive alternative to lidar to generate accurate and precise HS maps. However, its use implies the presence of an operator to drive an unmanned aircraft vehicle (UAV) or a pilot to fly an airplane .

Satellite snow cover observations, including operational applications, have been performed for many decades e.g.. Numerous satellite-derived products exist at the global scale . Snow-covered area maps are routinely produced from visible or near-visible bands (e.g. MODIS products; ). 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.. 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;. 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 $\mathrm{km}$ , and the presence of liquid water in the snowpack. Another limitation is the SWE threshold, which impedes SWE retrieval for deep snowpacks (; $>$ 0.15 $\mathrm{m}$–0.20 $\mathrm{m}\mathrm{w}.\mathrm{eq}.$). Several attempts have been made to retrieve spatially distributed HS or SWE data from space by radar imagery . 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.

Data sources and description. ADS means Airbus Defence and Space. GEODE and CESBIO are both laboratories of the Toulouse University (France). GeoFalco is a French start-up specialized in UAV data acquisition and processing.

Data sources	Acquisition	Institution	Ground sampling	Photogrammetric	Products
	date	(acquired by)	distance ($\mathrm{m}$)	information	(resolution)
Pléiades-1B triplet	26 Oct 2014	ADS	0.70–0.73 $\mathrm{m}$	$B$ $/$ $H$ $=$ 0.22; 0.23; 0.45	snow-free DEM (1 $\mathrm{m}$;2 $\mathrm{m}$; 4 $\mathrm{m}$)
Snow probe	10 Mar 2015	GEODE	10–30 $\mathrm{m}$	-	validation dataset
measurements		CESBIO
UAV photographs	10 Mar 2015	GeoFalco	0.10–0.40 $\mathrm{m}$	70 % end lap, 70 % side lap	winter DEM (0.1 $\mathrm{m}$;2 $\mathrm{m}$)
Pléiades-1A triplet	11 Mar 2015	CNES	0.70–0.73 $\mathrm{m}$	$B$ $/$ $H$ $=$ 0.22; 0.26; 0.48	winter DEM (1 $\mathrm{m}$;2 $\mathrm{m}$; 4 $\mathrm{m}$)
UAV photographs	13 Jul 2015	GeoFalco	0.10–0.40 $\mathrm{m}$	70 % end lap, 70 % side lap	snow-free DEM (0.1 $\mathrm{m}$;2 $\mathrm{m}$)

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 .

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-of-summer glacier surfaces , lake deposits and dunes , or landslide areas . Pléiades-derived DEMs exhibited sub-metre accuracy 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 $\mathrm{km}$${}^{\mathrm{2}}$), which is an open alpine terrain in the north-eastern Pyrenees (Fig. ). Bassiès is one of the main sub-basins of the Upper Vicdessos Valley, which has a long history of hydropower production . The elevation ranges between 1156 and 2676 $\mathrm{m}$ a.s.l. (median elevation 1659 $\mathrm{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${}^{\circ }$C and the mean annual precipitation is 1640 $\mathrm{mm}$, of which at least 30 % falls as snow . The snow season generally starts in November–December and ends in May–June (Fig. ). 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).

Datasets

Pléiades images

The satellites Pléiades-1A and 1B fly on the same near-polar sun-synchronous orbit at an altitude of 694 $\mathrm{km}$ with a 180${}^{\circ }$ phase and descending node at 10:30 a.m. The CCD optical sensors acquire images in push broom mode by using 5 $\times$ 6000 pixel arrays and a maximum of 20 integration lines (TDI) for the panchromatic band (480–830 $\mathrm{nm}$) . 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 $\mathrm{m}$ but are oversampled at 0.50 $\mathrm{m}$ before image delivery by a post-processing algorithm that was implemented by the French Space Agency (CNES).

Top: localization of study site relative to Europe and to the Pyrenees mountain. Bottom: Bassiès catchment (14.5 ${\mathrm{km}}^{\mathrm{2}}$) on a TOP25 topographic map ^©IGN (Institut national de l'information géographique et forestière).

[Figure omitted. See PDF]

Comparison of terrestrial oblique pictures taken by an automatic camera (see Fig.  for localization). On the left, the pictures were taken on 26 October 2014. On the right, the pictures are from 10 March 2015. These photographs were taken at the same time of day as the Pléiades images taken on 11 March 2015 (11:00 LT).

[Figure omitted. See PDF]

Two Pléiades triplets were acquired over the Bassiès catchment, which is the area of interest in this study (Table ). 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 ${\mathrm{km}}^{\mathrm{2}}$, which was centred on the Bassiès catchment. The images were acquired with viewing angles of 11.9, 0.7, and $-$10.9${}^{\circ }$ in the along-track direction with respect to the nadir and $-$4.8, $-$4.3, and $-$3.7${}^{\circ }$ 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${}^{\circ }$, azimuth 167${}^{\circ }$). 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 ${\mathrm{km}}^{\mathrm{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${}^{\circ }$ in the along-track direction with respect to the nadir and with viewing angles of 0.4, $-$2.7, and $-$6.4${}^{\circ }$ 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${}^{\circ }$, azimuth 157${}^{\circ }$).

UAV images

Two winged-UAV photogrammetric surveys were performed over a central subset of the Bassiès catchment (3.15 $\mathrm{km}$${}^{\mathrm{2}}$) to determine the snow depth by DEM differencing (Table ). The UAV was a real-time kinematic (RTK) ^© eBee that was equipped with a 12 $\mathrm{MP}$ camera:

• in winter, on 10 March 2015, 785 images during four parallels flights with 70 % lateral and longitudinal overlaps were acquired by a Canon IXUS 127 HS camera (4608 $\times$ 3465 pixels; sensor dimension: 6.170 $\mathrm{mm}$ $\times$ 4.628 $\mathrm{mm}$; focal length: 4.380 $\mathrm{mm}$);

• in summer, on 13 July 2015, 964 images during four parallels flights with 70 % lateral and longitudinal overlaps were acquired by a Sony DSC-WX220 camera (4896 $\times$ 3672 pixels; sensor dimension: 6.170 $\mathrm{mm}$ $\times$ 4.628 $\mathrm{mm}$; focal length: 4.572 $\mathrm{mm}$).

The flight altitude was maintained at approximately 150 $\mathrm{m}$, which provided a mean ground sampling distance (GSD) from 0.10 to 0.40 $\mathrm{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 $\mathrm{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 ). 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 $\mathrm{m}$. We used two types of snow probes with lengths of 2.2 and 3.2 $\mathrm{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 $\mathrm{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 (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.

Methods

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 (Fig. ). The ASP was primarily designed to create DEMs of ice and bare-rock 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) . Spatio-triangulation 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 $\mathrm{m}$ cell sizes from the merged point cloud through the point2dem command. Resolutions lower than 1 $\mathrm{m}$ are not relevant given the original image resolution and resolutions higher than 4 $\mathrm{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) .

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. .

Input data	Reference data	Type of	Values	Comments
		coregistration	of adjustments
4 $\mathrm{m}$ Pléiades	4 $\mathrm{m}$ Pléiades	$xy$ relative	$-$5.2 $\mathrm{m}$ north	Workflow metrics
winter DEM	summer DEM	coregistration	$+$2.8 $\mathrm{m}$ east	(same shifts applied to the 1 $\mathrm{m}$
		$\mathrm{\Delta }X\mathrm{rel}.\mathrm{\Delta }Y\mathrm{rel}.$		and 2 $\mathrm{m}$ Pléiades winter DEMs )
1 $\mathrm{m}$ Pléiades	1 $\mathrm{m}$ Pléiades	$xy$ relative	$-$5.2 $\mathrm{m}$ north (SD $=$ 0.70 $\mathrm{m}$)	Verification
winter orthoimage	summer orthoimage	coregistration	+3.2 $\mathrm{m}$ east (SD $=$ 0.50 $\mathrm{m}$)	metrics
1, 2, 4 $\mathrm{m}$ Pléiades	dDEM snow-free	$z$ relative	$b_{\mathrm{1}\mathrm{m}}=-0.46\mathrm{m}$ (SD $=$ 0.25 $\mathrm{m}$)	Workflow
dDEMs	football field	coregistration	$b_{\mathrm{2}\mathrm{m}}=-0.48\mathrm{m}$ (SD $=$ 0.20 $\mathrm{m}$)	metrics
		$b$	$b_{\mathrm{4}\mathrm{m}}=-0.44\mathrm{m}$ (SD $=$ 0.15 $\mathrm{m}$)
2 $\mathrm{m}$ Pléiades dDEMs	78 widespread points	$z$ relative	$\mathrm{Median}b=-0.70\mathrm{m}$	Verification
	over snow-free areas	coregistration	$\mathrm{Mean}b=-0.74\mathrm{m}$	metrics
		$b$	$\mathrm{SD}b=0.26\mathrm{m}$
1 $\mathrm{m}$ Pléiades	6 widespread points on	$xy$ absolute	+3 $\mathrm{m}$ north (SD $=$ 0.38 $\mathrm{m}$)	Workflow metrics
summer orthoimage	the 0.50 $\mathrm{m}$ IGN orthoimage	coregistration	$-$0.8 $\mathrm{m}$ east (SD $=$ 0.35 $\mathrm{m}$)	(same shifts applied to all
		$\mathrm{\Delta }X\mathrm{abs}.\mathrm{\Delta }Y\mathrm{abs}.$		the Pléiades dDEMs)
0.1 $\mathrm{m}$ UAV dDEM	353 widespread points	$\mathrm{\Delta }Z$ correction based	RMSE: 0.34 $\mathrm{m}$	Post-treatment correction.
	over snow-free areas	on a trend surface		Same correction applied
		of order 3		on the 1 $\mathrm{m}$ and 2 $\mathrm{m}$ UAV dDEMs

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”).

Data source	Pixel				Percentage of
	size	No data	$\mathrm{\Delta }Z<\mathrm{0}\mathrm{m}$			$\mathrm{\Delta }Z>15\mathrm{m}$
			in the catchment	on snow	on sunny snow	in the catchment	on snow
Pléiades tri-stereo	1 $\mathrm{m}$	2.4 %	22.4 %	14.7 %	9.4 %	0.14 %	0.09 %
	2 $\mathrm{m}$	1.7 %	24.5 %	17 %	11.3 %	0.15 %	0.1 %
	4 $\mathrm{m}$	1.2 %	22 %	14.5 %	9.8 %	0.17 %	0.1 %

Successive phases in Pléiades triplets processing and associated metrics of assessment.

[Figure omitted. See PDF]

Snow-free and winter DEMs at 4 m were horizontally co-registered by iteratively shifting the winter DEM with respect to the summer DEM (reference) by minimizing the standard deviation (SD) of the elevation difference distribution . The final horizontal shifts were applied to the winter DEM were $-$5.2 m in northing and $+$2.8 m in easting (Table ). We obtained similar results by computing the optimal shift at 1 and 2 $\mathrm{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 $\mathrm{m}$ for a circular error at a confidence level of 90 % (CE90) for Pléiades-1A and 4.5 $\mathrm{m}$ for Pléiades-1B . The same shift was applied to the 2 and 1 $\mathrm{m}$ winter DEMs.

Winter and snow-free nadir images were rectified at 1 $\mathrm{m}$ resolution from their respective DEMs, before co-registration. By picking six widespread corresponding points on the snow-free and winter images, the mean shifts were $-$5.2 $\mathrm{m}$ in northing (SD $=$ 0.7 $\mathrm{m}$) and $+$3.2 $\mathrm{m}$ in easting (SD $=$ 0.5 $\mathrm{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 ). 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 $\mathrm{m}$ spatial resolution by subtracting the snow-free DEM from the winter DEM on a pixel by pixel basis: $$\mathrm{\Delta }{Z}_{\mathrm{0}}=Z_{\text{w}}-Z_{\text{s}},$$ where $Z_{\text{w}}$ is the pixel value in the winter DEM, and $Z_{\text{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 $\mathrm{m}$. The shift between the snow-free Pléiades orthoimage and the IGN orthophoto was $+$3 $\mathrm{m}$ (SD $=$ 0.38 $\mathrm{m}$) in northing and $-$0.8 $\mathrm{m}$ (SD $=$ 0.35 $\mathrm{m}$) in easting (Table ). 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 $\mathrm{\Delta }{Z}_{\mathrm{0}}$ (Eq. ) to obtain the final dDEMs: $$\mathrm{\Delta }Z=\mathrm{\Delta }{Z}_{\mathrm{0}}-b,$$ 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 $\mathrm{km}$ from the mountain refuge (Fig. ). The value of $b$ was assumed to be equal to the median of the dDEM distribution on the football field (Table ). After this bias correction, dDEM pixels with negative values were classified as “no data”, which include 8 to 10 pixels that correspond to a snow-probe measurement (Table ). We classified the percentage of negative dDEM pixel values over the Bassiès catchment according 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 $\mathrm{m}$ Pléiades dDEM before bias correction (Eq. ). We photo-interpreted these points on snow-free rock areas, roads, or bare soil in the absolute georeferenced winter orthoimage by avoiding the steepest slopes ($<$ 30${}^{\circ }$) and by covering a large elevation range (790–2510 $\mathrm{m}$). We did not use this information to remove the bias because we aimed to evaluate a simple workflow that could become operational (Table ).

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 structure-from-motion algorithm . The focal 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 $\mathrm{m}$ cell sizes for both the snow-free and winter DEMs. Subsequently, 0.1, 1, and 2 $\mathrm{m}$ dDEMs were obtained by differencing the corresponding snow-free DEM from the winter DEM. The UAV-image acquisition, the UAV-image processing, and the UAV DEM generation were performed by a private company (Table ).

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 $\mathrm{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 $\mathrm{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 $\mathrm{m}$; for order 1, RMSE is 0.44 $\mathrm{m}$; for order 2, RMSE is 0.39 $\mathrm{m}$; for order 3, RMSE is 0.34 $\mathrm{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 $\mathrm{m}$: $+$0.33 $\mathrm{m}$; 1 $\mathrm{m}$: $+$0.43; 2 $\mathrm{m}$: $+$0.41 $\mathrm{m}$).

Pléiades and UAV dDEM assessments and comparison

We compared the Pléiades dDEM at 1, 2, and 4 $\mathrm{m}$ resolutions and the UAV dDEM at 0.1, 1, and 2 $\mathrm{m}$ resolution to the snow-probe measurements. We calculated the values of the residual vector $R_{\mathrm{\Delta }Z}$ as follows: $$R_{\mathrm{\Delta }Z}=\mathrm{\Delta }{Z}_i-\mathrm{HS}.$$ $\mathrm{\Delta }{Z}_i$ is the subset of the dDEM values, where $\mathrm{\Delta }Z\ge \mathrm{0}$ after bias correction (Eq. ), 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 ${\mathit{\sigma }}_{\mathrm{probe}}=0.15\mathrm{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_{\mathrm{\Delta }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 : $$\mathrm{NMAD}=1.4826\mathrm{median}\left(\right|(R_{\mathrm{\Delta }Z}-m_{R_{\mathrm{\Delta }Z}}|),$$ where $m_{R_{\mathrm{\Delta }Z}}$ is the median of the residuals.

We also assessed how the $\mathrm{\Delta }Z$ and HS values correlate through a rank correlation method. We used the Spearman correlation factor, called cor${}_{\text{s}}$, which is sensitive to neither the presence of outliers nor the existence of nonlinear correlations .

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 $\mathrm{m}$ UAV dDEM from the 2 $\mathrm{m}$ Pléiades dDEM. We visually compared both dDEM maps and the dDEM differences. We performed a SW–NE transect (1.6 $\mathrm{km}$ long) and compared the dDEM values along that transect.

Residual analysis on the Pléiades data

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. ). The Pléiades panchromatic images had a pixel size of 0.5 $\mathrm{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 $\mathrm{m}$ Pléiades dDEM as follows: $$\begin{array}{}& & R_{Z_{\text{w}}}=Z_{\text{w}}-Z_{w,\mathrm{DGPS}},& & R_{Z_{\text{s}}}=Z_{\text{s}}-(Z_{w,\mathrm{DGPS}}-\mathrm{HS}).\end{array}$$ We assume all the random errors to be uncorrelated. The random error on the DGPS measurements was ${\mathit{\sigma }}_{\mathrm{DGPS}}=0.1\mathrm{m}$ for all the points. The second term of the Eq. () $(Z_{w,\mathrm{DGPS}}-\mathrm{HS})$ provides snow-free reference elevation values at the 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 snow-probe measurements (${\mathit{\sigma }}_{\mathrm{DGPS}+\mathrm{probe}}=0.18\mathrm{m}$; see the Supplement for more details on the random error calculation).

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 .

Data source	dDEM	Number of	Median	Standard	NMAD	Spearman correlation
	pixel size	snow-probe sampling	($\mathrm{m}$)	deviation (m)	($\mathrm{m}$)	cor${}_{\text{s}}(\mathrm{\Delta }Z,\mathrm{HS})$
Pléiades tri-stereo	1 $\mathrm{m}$	443	$-$0.15	0.62	0.47	0.71${}^{*}$
	2 $\mathrm{m}$	442	$-$0.16	0.58	0.45	0.72${}^{*}$
	4 $\mathrm{m}$	441	$-$0.12	0.69	0.51	0.67${}^{*}$
Pléiades front/nadir stereo pair	4 m	411	$-$0.54	0.64	0.53	0.62${}^{*}$
Pléiades nadir/back stereo pair	4 $\mathrm{m}$	450	0.13	0.61	0.47	0.73${}^{*}$
UAV photographs	0.1 $\mathrm{m}$	343	$-$0.07	0.63	0.38	0.8${}^{*}$
	1 $\mathrm{m}$	336	$-$0.15	0.62	0.36	0.79${}^{*}$
	2 $\mathrm{m}$	339	$-$0.11	0.62	0.35	0.79${}^{*}$

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.

[Figure omitted. See PDF]

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 $\mathrm{m}$ (Fig. ). 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. ). Here, we summarize the different ranges of the main topographical variables that are associated with the snow-probe data:

• The sampling-elevation range is 1645–2000 $\mathrm{m}$ a.s.l. However, 70 % of the snow-probe values are between 1645 and 1700 $\mathrm{m}$ a.s.l. The elevation range of the catchment is 1156–2676 $\mathrm{m}$ a.s.l. (median elevation of 1930 $\mathrm{m}$). Therefore, we did not assess the residuals' distribution (Eq. 1) according to the elevation.

• The slope, which was derived from the 2 $\mathrm{m}$ snow-free Pléiades DEM, was associated with the snow-probe measurements and ranged continuously from almost 0${}^{\circ }$ (flat areas) to 25${}^{\circ }$. 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${}^{\circ }$, 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 land-cover 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 snow-probe measurements.

Hillshade of snow-free 2 $\mathrm{m}$ Pléiades (top left) and 1 $\mathrm{m}$ UAV (bottom left) DEMs and winter 2 $\mathrm{m}$ Pléiades (top right) and 1 $\mathrm{m}$ UAV (bottom right) DEMs. The outlines of the UAV extents are indicated in the Pléiades extents.

[Figure omitted. See PDF]

2 $\mathrm{m}$ Pléiades dDEM map interpreted as snow height (metres) in Bassiès catchment.

[Figure omitted. See PDF]

Results

Pléiades dDEM assessments

The snow-free DEM and winter DEM are shown in Fig. . 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. ) has the typical gamma or log-normal distribution shape that is reported in the literature e.g.. Considering the whole Bassiès catchment, the mean of $\mathrm{\Delta }Z$ is 2.15 $\mathrm{m}$ and its SD is 1.72 $\mathrm{m}$.

The Pléiades 2 $\mathrm{m}$ dDEM is composed of 1.7 % of no-data entries in the Bassiès catchment (2.4 and 1.2 % for the 1 and 4 $\mathrm{m}$ dDEMs respectively) (Table ). These no-data entries originate from data gaps in the raw points clouds, which are produced by the ASP before rasterization. Considering the Bassiès catchment area, 25 % of the pixels of the 2 $\mathrm{m}$ dDEM exhibit negative values (23 and 22 % for the 1 $\mathrm{m}$ and 4 $\mathrm{m}$ dDEMs respectively). The percentage of negative 2 $\mathrm{m}$ dDEM pixel values on the SCA is 17 % (14.7 and 14.5 % for the 1 and 4 $\mathrm{m}$ dDEMs respectively). This fraction is less important if we do not consider the snow pixels in the shaded areas (direct shadow from the surrounding cliffs): 11.3 % for the Pléiades 2 $\mathrm{m}$ dDEM (9.4 and 9.8 % for the 1 and 4 $\mathrm{m}$ dDEMs respectively). The Pléiades 2 $\mathrm{m}$ dDEM pixels with values above 15 $\mathrm{m}$ represent a very limited fraction, which is negligible on snow (less than 0.1 %). These values should most probably be interpreted as inconsistent and classified as no data.

We calculated a constant vertical bias $b$ from a snow-free football field (see section ). The value of $b$ for each dDEM resolution is $b_{\mathrm{1}\mathrm{m}}=-0.46\mathrm{m}$, $b_{\mathrm{2}\mathrm{m}}=-0.48\mathrm{m}$, and $b_{\mathrm{4}\mathrm{m}}=-0.44\mathrm{m}$. The bias value distribution of the football field has a mean value that is close to the median (1 $\mathrm{m}$: $-$0.43 $\mathrm{m}$; 2 $\mathrm{m}$: $-$0.45 $\mathrm{m}$; and 4 $\mathrm{m}$: $-$0.42 $\mathrm{m}$) and a low SD (1 $\mathrm{m}$: 0.25 $\mathrm{m}$; 2 $\mathrm{m}$: 0.20 $\mathrm{m}$; 4 $\mathrm{m}$: 0.15 $\mathrm{m}$). The bias assessment which was performed over the entire Pléiades dDEM (110 $\mathrm{km}$${}^{\mathrm{2}}$) and was based on 78 widespread values (see section ) indicates a median of $-$0.70 $\mathrm{m}$, a mean of $-$0.74 $\mathrm{m}$, and an SD of 0.26 $\mathrm{m}$. The low SD value and the median difference confirm the possibility to remove a constant bias from a unique area, with small random and systematic errors: $$\mathrm{median}\left(\mathrm{football}\mathrm{field}\right)-\mathrm{median}\left(\mathrm{entire}\mathrm{dDEM}\right)=-0.22\mathrm{m}.$$ The comparison with the snow-probe data indicates that the Pléiades dDEMs are consistent with the snow depth measurements (Table ). The median values of the residuals distribution $R_{\mathrm{\Delta }Z}$ are relatively low (between $-$0.12 and $-$0.16 $\mathrm{m}$) and close to the mean of the distribution at each pixel resolution ($\pm$0.05 $\mathrm{m}$ between the median and mean). A slight influence from the pixel size is present (Table ). For our validation dataset, the 2 $\mathrm{m}$ Pléiades dDEM exhibits slightly better precision and accuracy. For this dDEM, the SD is 0.58 $\mathrm{m}$ and the NMAD is 0.45 $\mathrm{m}$. The $\mathrm{\Delta }Z$ and HS datasets are significantly correlated at each pixel resolution (cor${}_{\text{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 ($\mathrm{\Delta }{Z}_{\mathrm{2}\mathrm{m}}=0.90\cdot \mathrm{HS}$) (Fig. ).

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 .

Variable	Interval	Number of	Median	Standard	NMAD	Spearman correlation
	bins	snow-probe sampling	(m)	deviation ($\mathrm{m}$)	($\mathrm{m}$)	cor${}_{\text{s}}\left(\right|R_{\mathrm{\Delta }Z}|,\mathrm{HS})$
Snow depth	[0 m; 0.5 m]	25	0.24	0.22	0.31	0.3*
	[0.5 $\mathrm{m}$; 1 m]	65	$-$0.01	0.46	0.33
	[1 $\mathrm{m}$; 1.5 m]	94	$-$0.07	0.44	0.39
	[1.5 $\mathrm{m}$; 2 $\mathrm{m}$]	114	$-$0.24	0.60	0.34
	[2 $\mathrm{m}$;2.5 m]	72	$-$0.32	0.68	0.54
	[2.5 $\mathrm{m}$; 3.2 m]	46	$-$0.63	0.56	0.39
Slope	[0${}^{\circ }$ ; 5${}^{\circ }$]	150	$-$0.10	0.42	0.32	0.26${}^{*}$
	[5${}^{\circ }$; 10${}^{\circ }$]	117	$-$0.19	0.53	0.41
	[10${}^{\circ }$; 15${}^{\circ }$]	81	$-$0.30	0.53	0.59
	[15${}^{\circ }$; 20${}^{\circ }$]	63	$-$0.30	0.79	0.7
	$>$ 20${}^{\circ }$	31	$-$0.18	0.93	0.75
Aspect	North	159	$-$0.20	0.6	0.43	-
	East	113	$-$0.15	0.63	0.48
	South	134	$-$0.16	0.55	0.46
	West	43	$-$0.12	0.55	0.39
Land cover	All classes	442	$-$0.16	0.58	0.47	0.72
	Mineral	56	$-$0.2	0.79	0.60	0.74
	Water	21	$-$0.32	0.55	0.50	0.67
	Low grass	140	$-$0.16	0.49	0.35	0.74
	Shrub	140	$-$0.15	0.63	0.51	0.68
	Peatland	84	$-$0.15	0.51	0.42	0.69

Figure  shows the spatial distribution of the snow depth measurements and the residuals of the 2 $\mathrm{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. ).

Overall, the snow-probe dataset exhibits a low systematic error and is spatially homogeneously distributed.

Top left: histogram of the distribution of the residuals between the 2 $\mathrm{m}$ Pléiades dDEM and the snow-probe measurements. Bottom left: histogram of the distribution of the residuals between the 2 $\mathrm{m}$ UAV dDEM and the snow-probe measurements. Top right: scatter plot between the 2 $\mathrm{m}$ Pléiades dDEM and the snow-probe measurements. Bottom left: scatter plot between the 2 $\mathrm{m}$ UAV dDEM and the snow-probe measurements. Blue line is the 1 : 1 line. Red line is the least-square best fit of a linear function with a zero intercept $y=ax$.

[Figure omitted. See PDF]

Left: map of snow depth sampling over the Pléiades winter orthoimage. Middle: map of 2 $\mathrm{m}$ Pléiades residuals $R_{\mathrm{\Delta }Z}$. Right: map of 2 $\mathrm{m}$ UAV residuals $R_{\mathrm{\Delta }Z}$.

[Figure omitted. See PDF]

UAV dDEM assessments

The SD and the NMAD indicate a decimetric random error, at each pixel resolution (SD${}_{\mathrm{2}\text{m}}=0.62$ $\mathrm{m}$, NMAD${}_{\mathrm{2}\text{m}}=0.35$ $\mathrm{m}$). The $\mathrm{\Delta }Z$ and HS values are significantly correlated (mean cor${}_{\text{s}}$ 0.79). The median value of the residual distribution $R_{\mathrm{\Delta }Z}$ is slightly negative and ranges from $-$0.07 to $-$0.15 $\mathrm{m}$ according to the pixel size.

Figure  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${}_{\text{s}}$ between the Pléiades and UAV dDEM values is 0.62 and significant at 95 % confidence ($N={527.10}^{\mathrm{3}}$, number of values of the sample size). The datasets were not co-registered. The comparison between the 2 $\mathrm{m}$ Pléiades and the 2 $\mathrm{m}$ UAV dDEMs is characterized by a residual distribution with a median of $-$0.14 $\mathrm{m}$ (mean $-$0.06 $\mathrm{m}$), an SD of 1.47 $\mathrm{m}$, and an NMAD of 0.78 $\mathrm{m}$.

The 2 $\mathrm{m}$ Pléiades and UAV dDEM maps exhibit very similar patterns (Fig. ). 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. ). A transect over a common area that is covered by both the 2 $\mathrm{m}$ Pléiades and 2 $\mathrm{m}$ UAV dDEMs highlights the consistency in both $\mathrm{\Delta }Z$ variations. Over this transect, the SD of the residuals between the Pléiades and UAV dDEMs is 0.78 $\mathrm{m}$ and the median is $-$0.16 $\mathrm{m}$.

Pléiades and UAV dDEM and dDEM differences. dDEMs and $\mathrm{\Delta }$dDEM (dDEM differencing) values along the transect.

[Figure omitted. See PDF]

Top: residuals of the comparison between the 2 $\mathrm{m}$ Pléiades winter DEM and the winter DGPS measurements (see Eq. (5), Sect. 4.4.2), after removal of the bias (median of the residuals). Middle: residuals of the comparison between the snow-free 2 $\mathrm{m}$ Pléiades DEM and the estimated summer surface elevation (see Eq. (6), Sect. 4.4.2), after removal of the bias (median of the residuals). Bottom: residuals of the comparison between the 2 $\mathrm{m}$ Pléiades dDEM (black bars) and the snow-probe measurements according to the probe Id ranked in the ascending HS (red line) order (see Eq. (3), Sect. 4.3), and after removal of the bias (median of the residuals).

[Figure omitted. See PDF]

Residual analysis on the Pléiades data

Photogrammetric processes

The density values of the raw point clouds (pts. $\mathrm{m}$${}^{-\mathrm{2}}$) from the correlation process based on the front nadir/stereo pair are close to the maximum achievable value (4 pts. $\mathrm{m}$${}^{-\mathrm{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 snow-free and winter DEMs to the dDEM residuals (Eqs.  and , Fig. ). The medians of $R_{Z_{\text{w}}}$ and $R_{Z_{\text{s}}}$ distributions are $-0.91$ and $-0.25$ $\mathrm{m}$ respectively, leading to a difference of $$\mathrm{median}\left(R_{Z_{\text{w}}}\right)-\mathrm{median}\left(R_{Z_{\text{s}}}\right)=-0.66\mathrm{m}.$$ This value is consistent with the median of $-$0.64 $\mathrm{m}$ for the $R_{\mathrm{\Delta }{Z}_{\mathrm{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 $\mathrm{m}$). The $R_{Z_{\text{w}}}$ and $R_{Z_{\text{s}}}$ values in Fig.  are corrected from the bias by removing the median. The SDs of $R_{Z_{\text{w}}}$ and $R_{Z_{\text{s}}}$ are 0.32 $\mathrm{m}$ and 0.66 $\mathrm{m}$ respectively. These estimations are consistent with the SD of the 2 $\mathrm{m}$ residual distribution $R_{\mathrm{\Delta }Z}$ (0.58 $\mathrm{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_{\mathrm{\Delta }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 ).

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 $\mathrm{m}$). The most important dispersions are associated to the mineral and the shrub classes (SD $=$ 0.79 and 0.63 $\mathrm{m}$).

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 . 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 snow-free 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 .

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. ). This observation could be due to the difficulty in treating micro-topography with the native GSD of Pléiades (0.7 $\mathrm{m}$ at nadir). Snow-covered areas offer a smoother surface compared to vegetated or stony snow-free 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 ). In both the snow-free 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 $\mathrm{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 ). 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 multi-view stereo possibility of the ASP (only available since version 2.5.0), which limited our correlation to successive pair matching. 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 4 $\mathrm{m}$ dDEM; SD $=$ 0.64 and 0.61 $\mathrm{m}$ for the bi-stereo). The accuracy was slightly better for the tri-stereo dDEM (median $=$ $-$0.12 $\mathrm{m}$ for tri-stereo; median $=$ $-$0.54 and +0.13 $\mathrm{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 assess the entire snow depth variability at a plot scale, typically 10 $\mathrm{m}$ $\times$ 10 $\mathrm{m}$ . 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${}_{\text{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, sampling 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 snow-free 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 “non-convergent”, which could also result in marked deformation . Winged UAVs are potentially less stable than UAVs with rotors , although recent works have highlighted their great potential for snow mapping in high-alpine catchments even in relative windy conditions . 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 .

Recent works based on UAV systems to map snow depth highlight much better performance than the results reported in this study (2 $\mathrm{m}$ UAV dDEM: SD $=$ 0.62 $\mathrm{m}$, NMAD $=$ 0.35 $\mathrm{m}$, median $=$ $-$0.11 $\mathrm{m}$; see Table ). 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}^{-\mathrm{3}}$ $\mathrm{m}$) over a small mountainous terrain (0.07 ${\mathrm{km}}^{\mathrm{2}}$) with thick vegetation cover. A comparison with a reduced sample of snow-probe measurements ($N=20$) highlighted an RMSE of 0.096 $\mathrm{m}$ using GCPs and 0.184 $\mathrm{m}$ without (0.084 $\mathrm{m}$ with one point of co-registration). In , an UAV-octocopter was used to collect imagery at two alpine sites of the region of Davos in the swiss Alps (1940 $\mathrm{m}$ and 2500 $\mathrm{m}$ a.s.l.). The images were acquired with a customized Sony NEX-7 camera with an overlap of 70 % along and across track. Reference data were constituted by plots of 1 m${}^{\mathrm{2}}$ with five manual snow depth measurements. Four snow depth maps were produced and assessed with the manual plots (between 12 and 22 plots according to the map). Accuracies of 0.07 to 0.15 $\mathrm{m}$ RMSE are reported in a detailed analysis, according to the study sites and the land-cover classes. Considering all the reference plots in the valley bottom site, the HS RMSE is 0.25 $\mathrm{m}$ and there is an average systematic underestimation of HS by 0.20 $\mathrm{m}$. In , a Sensefly Ebee RTK UAV was used to collect imagery at a cultivated agricultural Canadian Prairie and a sparsely vegetated Rocky Mountain alpine ridgetop site (2 300 $\mathrm{m}$ a.s.l.). In the alpine site, the images were acquired with a Canon IXUS, with a lateral overlap of 85 %, a longitudinal overlap of 75 %, and a flight altitude of 100 $\mathrm{m}$. Multiple acquisitions (43) were performed with careful flight plans. The snow depth was measured with five snow depth measurements in a 0.4 $\mathrm{m}$ $\times$ 0.4 $\mathrm{m}$ square at the locations of the GNSS survey locations. The average snow depth of the five values was then compared to the snow depth determined by the UAV, with a number of snow depth measurements between 3 and 20 measurements per flight. The reported snow depth accuracy is characterized by a RMSE of 0.085 m.

In the case of our study, the DEM of the SCA was generated from a unique flight plan. Some problematic flights were reported by (5 from 43 flights for all sites, or 11.6 %) with DEMs showing an RMSE of up to 0.32 $\mathrm{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 , 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$\mathit{\sigma }$ 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 . 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 airborne-derived 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 state-of-the-art of digital aerial photography. The performance highlighted by the UAV system mentioned in the previous section are very satisfactory . Nevertheless their spatial coverage is limited to several hectares. have created a DEM over a SCA of a mountainous terrain in Australia from a high-spatial-resolution camera (GSD up to 0.05 $\mathrm{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 $\mathrm{m}$ with an SD of 0.08 $\mathrm{m}$. employed an optoelectronic line scanner (ADS 80) that was mounted on an aeroplane to map the snow depth at 2 $\mathrm{m}$ resolution (GSD of 0.25 $\mathrm{m}$) in the Swiss Alps. A comparison between the ADS and different individual HS measurements revealed both RMSE and NMAD of approximately 0.3 $\mathrm{m}$, which is equivalent to GSD of 1 of the input images. Over the polar snow of Alaska, generated dDEMs over rather flat areas from a consumer-grade 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 $\mathrm{m}$ (GSD of 0.06 to 0.2 $\mathrm{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 $\mathrm{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 $\mathrm{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 $\mathrm{km}$ for a stereo acquisition and 25 $\mathrm{km}$ for a tri-stereo acquisition (195 $\mathrm{km}$ and 80 $\mathrm{km}$ respectively for a $B$ $/$ $H$ of 0.4) . Considering a $B$ $/$ $H$ of 0.2, areas of up to 1600 ${\mathrm{km}}^{\mathrm{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 $\mathrm{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 $\mathrm{m}$). In glaciology, DEMs that are generated from optical stereos are often considered to be inaccurate in accumulation areas . 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 . 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 ${\mathrm{km}}^{\mathrm{2}}$) was 2.15 $\mathrm{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 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 $\mathrm{mm}$ over the hydrological year 2011–2012 and 1150 $\mathrm{mm}$ in 2014–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 $\mathrm{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 $\mathrm{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 snow-probe 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 ${\mathrm{km}}^{\mathrm{2}}$). The distribution of the residuals between the 2 $\mathrm{m}$ Pléiades dDEM values and the corresponding snow-probe measurements presents a median of $-$0.16 $\mathrm{m}$ and an SD of 0.58 $\mathrm{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 snow-probe measurements (median of the residuals is $-$0.11 $\mathrm{m}$, SD is 0.62 $\mathrm{m}$). The comparison between the 2 $\mathrm{m}$ Pléiades and the 2 $\mathrm{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 $\mathrm{m}$ (mean is $-$0.06 $\mathrm{m}$), SD is 1.47 $\mathrm{m}$, and NMAD is 0.78 $\mathrm{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éiades-derived 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 $\mathrm{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.

Acknowledgements

This manuscript was greatly improved thanks to the constructive comments of Guillaume Chambon (editor), Yves Bühler (referee), and Emmanuel Thibert (referee). The Pléiades images acquisition was supported by the CNES through the ISIS Pléiades programme (project CRYOPYR–4500107601). This work was supported by the Région Midi-Pyrénées and the University of Toulouse through the CRYOPYR project. We thank the OHM Vicdessos Human-Environnemental Observatory and the Labex DRIIHM (Dispositif de Recherche Interdisciplinaire sur les Interactions Hommes-Milieux) for logistic and financial support. We also warmly thank Audrey Chone, Yoann Malbeteau, Yoann Moreau, and Vincent Rivalland who helped us to collect the snow-probe measurements. The authors deeply acknowledge Michael J. Willis for his help in ASP parameterization to produce the DEMs. We also thank Joaquín María Muñoz-Cobo Belart for his useful comments that improved the manuscript.Edited by: G. Chambon