Glacier topography and elevation changes derived from Pléiades sub-meter stereo images

In response to climate change, most glaciers are losing mass and hence contribute to sea-level rise. Repeated and accurate mapping of their surface topography is required to estimate their mass balance and to extrapolate/calibrate sparse field glaciological measurements. In this study we evaluate the potential of sub-meter stereo imagery from the recently launched Pléiades satellites to derive digital elevation models (DEMs) of glaciers and their elevation changes. Our five evaluation sites, where nearly simultaneous field measurements were collected, are located in Iceland, the European Alps, the central Andes, Nepal and Antarctica. For Iceland, the Pléiades DEM is also compared to a lidar DEM. The vertical biases of the Pléiades DEMs are less than 1 m if ground control points (GCPs) are used, but reach up to 7 m without GCPs. Even without GCPs, vertical biases can be reduced to a few decimetres by horizontal and vertical co-registration of the DEMs to reference altimetric data on ice-free terrain. Around these biases, the vertical precision of the Pléiades DEMs is ±1 m and even ±0.5 m on the flat glacier tongues (1σ confidence level). Similar precision levels are obtained in the accumulation areas of glaciers and in Antarctica. We also demonstrate the high potential of Pléiades DEMs for measuring seasonal, annual and multiannual elevation changes with an accuracy of 1 m or better if cloud-free images are available. The negative regionwide mass balances of glaciers in the Mont-Blanc area (−1.04± 0.23 m a water equivalent, w.e.) are revealed by differencing Satellite pour l’Observation de la Terre 5 (SPOT 5) and Pléiades DEMs acquired in August 2003 and 2012, confirming the accelerated glacial wastage in the European Alps.


Introduction
In a context of nearly global glacier wastage, new means to retrieve accurate and comprehensive measurements of glacier topography and elevation changes are welcome.Digital elevation models (DEMs) are needed to properly orthorectify satellite images and to extrapolate pointwise glaciological mass balance measurements to entire ice bodies (Kääb et al., 2005;Zemp et al., 2013).The geodetic method, based on the differencing of multi-temporal DEMs, has been used for decades to retrieve glacier-wide and region-wide glacier mass balances (e.g.Bamber and Rivera, 2007;Finsterwalder, 1954).This method reveals the spatial patterns of elevation changes over individual glaciers or entire regions.Geodetically derived mass balances are now included in global assessments of glacier mass loss (Cogley, 2009;Vaughan et al., 2013).Furthermore, the differences between multi-temporal DEMs derived from aerial photos and airborne lidar can be used to check and, if necessary, correct cumulative mass balances measured using the field-based glaciological method over periods of typically 5-10 years (e.g.Abermann et al., 2010;Jóhannesson et al., 2013;Soruco et al., 2009b;Zemp et al., 2013).However, airborne sensors cannot acquire data everywhere on Earth because of the logistical difficulties involved in flying an airplane over some remote regions (e.g.high-mountain Asia, polar regions).Lidar data from the Ice, Cloud and land Elevation Satellite (ICESat) mission (and from the future ICESat-2) remain too sparse to provide a comprehensive coverage of individual glaciers; hence, mass balances can be retrieved reliably only for sufficiently large regions (Arendt et al., 2013;Gardner et al., 2013;Kääb et al., 2012).The geodetic method has also been applied extensively to DEMs derived from spaceborne optical or radar sensors such as the Shuttle Radar Topographic Mission (SRTM) and SPOT DEMs (e.g.Gardelle et al., 2013).However, the vertical errors of these DEMs (5 to 10 m) and their resolution (40 to 90 m) limits the possibility of retrieving accurate glacier-wide mass balances on individual small to medium size glaciers covering typically less than 10 km 2 for time periods of a few years.DEMs derived from sub-meter stereo images have the potential to fill this gap between coarse spaceborne DEMs and very high resolution DEMs from aerial surveys.
After the launch of the first non-military sub-meter resolution satellite (IKONOS) in September 1999 and until recently, relatively little work was carried out on deriving DEMs from these images over glaciers, probably due to the cost of the images.However, over the last 2-3 years, interest in these data sets has grown due to easier accessibility by researchers (e.g.Haemmig et al., 2014;Marti et al., 2014;Sirguey and Cullen, 2014).In the US, WorldView-1 and WorldView-2 images are distributed by the Polar Geospatial Center (PGC) to US-funded researchers (Noh and Howat, 2014a, b).Since the launch of the Pléiades 1A and 1B satellites in December 2011 and 2012, their high-resolution images have become available for researchers from the European Union, Iceland, Norway and Switzerland through the ISIS program of the French Space Agency, CNES (http: //www.isis-cnes.fr/).In this context, the goal of the present study is to assess the accuracy of the DEMs retrieved from Pléiades stereo images and to test their potential to estimate glacier elevation changes at seasonal, annual and multiannual timescales.

Pléiades stereo images
Pléiades stereo pairs acquired in five different regions are used in this study (Fig. 1, Table 1).The study sites were selected to represent a variety of glacial settings ranging from the small (1 km 2 ) Agua Negra Glacier in the high (> 5000 m a.s.l.), steep and arid Andes of Argentina to the flat and 7 km wide Astrolabe Glacier, an outlet glacier of the East Antarctic ice sheet in Adélie Land.The main reason for selecting these glaciers was that they were all targets of ongoing field programs (Björnsson et al., 2013;Jóhannesson et al., 2013;Le Meur et al., 2014;Vincent et al., 2009Vincent et al., , 2014;;Wagnon et al., 2013) so that accurate reference data were available.A Pléiades image is shown for each study site in Fig. 2, with an enlargement of a small area in the upper part of Astrolabe Glacier.
The Pléiades 1A and 1B twin satellites were launched 17 December 2011 and 2 December 2012, respectively.Images are delivered at a ground sampling distance (pixel size) of 0.5 m for the panchromatic channel (wavelength in Table 1.Characteristics of the study sites and Pléiades images.For each site, the approximate geographic coordinates and maximum altitude (Z max , m a.s.l.) are indicated.All images were acquired by Pléiades-1A, except a Pléiades-1B stereo pair over Tungnafellsjökull.The base-toheight ratio (B/H), the ratio of the distance between two successive positions of the satellite to its height above ground, is an indicator of the sensitivity to topography.A single B/H is indicated for stereo pairs, whereas three values for the front/nadir, nadir/back and front/back pairs are provided for tri-stereos.The percentage of saturation in the image is given for the front/back images for stereo pairs and front/nadir/back for tri-stereos.The reference (noted Ref.) altimetric data used to evaluate the Pléiades DEMs are kinematic global navigation satellite system (kGNSS), stop and go GNSS measurements and, for the Tungnafellsjökull Ice Cap only, a lidar DEM.Identification codes of the Pléiades images are not listed for the sake of concision.the 480-830 nm range) and 2 m for the multi-spectral blue, green, red and near-infrared channels (http://smsc.cnes.fr/PLEIADES/).However, the actual resolution of the sensor is slightly coarser (0.7 and 2.8 m) and the images are therefore oversampled by the ground segment (Gleyzes et al., 2012).One advantage of the Pléiades satellites (compared to SPOT1-5 and ASTER, Advanced Spaceborne Thermal Emission and Reflection Radiometer) for glaciological studies is the fact that the panchromatic band images are coded on 12 bits (digital number range is 2 n , where n is the number of bits).Such a wide radiometric range gives much better image contrast over flat and textureless regions (such as snowcovered accumulation areas, Fig. 2) and the risk of saturation is reduced.As a direct result of the 12 bit encoding, the percentage of our Pléiades images with saturation, i.e.where digital numbers equal 4095, is low (Table 1).The maximum percentage is observed over the Mont Blanc images (5 %) and is due to several snowfalls during the days preceding the 20 September 2013 acquisition.This date excluded, the percentage is always lower than 1 %.However, new Pléiades images acquired in north-west Himalaya in August 2014 contained a higher percentage of saturated pixels, sometimes up to 10 %.This is probably due to a high solar elevation angle in August at this relatively low latitude (∼ 33 • N).In such cases, specific acquisition parameters (i.e.lower number of time delay integration, TDI, stages) may help to reduce the saturated areas.
A positive consequence of this 12 bit encoding is that nearly all images from the archive are now useful for the study of ice and snow surfaces, whereas with SPOT1-5 or ASTER stereo imagers (8 bit encoding), special acquisition plans with a low gain had to be defined to avoid saturation and to enhance image contrast over snow and ice (Korona et al., 2009;Raup et al., 2000).
The images of a Pléiades stereo pair are acquired along the orbit (along track) within a few tens of seconds due to the agility (pitch) of the platform.Triplets of images (referred to as tri-stereos) are available for two of our study sites, in the Andes and in the Alps.A tri-stereo is made of three images (front, nadir and back images) that can be combined in three stereo pairs for multiple DEM generation: front/nadir, nadir/back and front/back.The front/nadir and nadir/back pairs are acquired from closer points of view than the front/back pair and hence exhibit less distortion, facilitating the recognition of identical features in the images by automatic matching algorithms.However, the sensitivity to topography is reduced by a factor of about 2, so matching errors will lead to doubled altimetric errors (Toutin, 2008).As for all optical sensors, the main drawback of the Pléiades constellation is the need for clear-sky conditions to obtain suitable cloud-free images.

SPOT 5 DEMs
Two SPOT 5 DEMs were used in this study.
First, a DEM covering 30 km × 60 km and including the entire Mont-Blanc area was computed in a previous study from two SPOT 5-HRG images acquired 19 and 23 August 2003 (Berthier et al., 2004).The ground sampling distance of a SPOT 5-HRG image is 2.5 m and the DEM pixel size is 10 m.This DEM was subtracted from the Pléiades 19 August 2012 Pléiades DEM to measure the geodetic mass balance for the entire Mont-Blanc area.
Second, a DEM including the Astrolabe Glacier was extracted from the SPIRIT (SPOT 5 stereoscopic survey of polar ice: reference images and topographies) database (Korona et al., 2009).The ground sampling distance of a SPOT 5-HRS (high-resolution sensor) image is 10 m across track and 5 m along track (Berthier and Toutin, 2008;Bouillon et al., 2006).This SPOT 5 DEM has a pixel size of 40 m and covers an area of 120 km × 230 km.
For these two study sites (Mont-Blanc area and Astrolabe Glacier), the Pléiades and SPOT 5 DEMs are also compared qualitatively to determine the added value of the Pléiades broader radiometric range and higher resolution for DEM generation, in particular over the snow-covered accumulation areas.

Validation data: GNSS and lidar
Table 1 includes the type, date and precision of the reference data available to evaluate the Pléiades DEMs.Most of the reference data are differential GNSS (global navigation satellite system) measurements including the US GPS and Russian GLONASS constellations.They are processed relative to the closest available base station.A decimetric precision can be reached in the best cases (the Mont-Blanc area), whereas, in the worst case, the precision is closer to ±0.5 m in the Andes due to the fact that the base station (TOLO) is located 100 km away from Agua Negra Glacier.For Iceland, the Pléiades DEM was also compared to an airborne lidar DEM with a 2 m pixel size.The lidar DEM was acquired 7 and 8 August 2011, slightly more than 2 years before the acquisition of Pléiades images (Jóhannesson et al., 2013).The vertical accuracy of the lidar DEM and its horizontal positioning accuracy has been estimated to be well within 0.5 m (Jóhannesson et al., 2011).

Pléiades DEM generation
All Pléiades DEMs presented in this paper were calculated using the OrthoEngine module of PCI Geomatica 2013.The original orientation of each image was read in the ancillary data provided with the images in the form of rational polynomial coefficients.Without ground control points (GCPs), the horizontal location accuracy of the images was estimated at 8.5 m (CE90, Circular Error at a confidence level of 90 %) for Pléiades-1A and 4.5 m for Pléiades-1B (Lebègue et al., 2013).This initial orientation was then refined using GCPs when available.DEMs were generated with a pixel size of 4 m, a compromise offering relatively fast processing and sufficient resolution.In addition to pixel size, two other main parameters can be tuned during DEM generation with Or-thoEngine: the level of detail of the DEM (from low to very high) and the type of relief in the scene (from plain to mountainous).Unless otherwise mentioned, all DEMs evaluated in the present paper were obtained with the low and mountainous settings.The choice of these parameters is justified in Sect.4.1.Data voids were generally not filled by interpolation during DEM generation because our aim was not to obtain complete coverage but rather to determine where elevations were extracted and what their quality was.Thus, the statistics on elevation differences are only provided outside of Pléiades DEM data voids, except when mentioned explicitly.
The use of a commercial software is one of the limitations of our study.For example, it is not possible to know the exact DEM extraction parameters (e.g.correlation window size, correlation threshold) hidden behind the PCI processing parameter settings.Future analyses are planned to evaluate and inter-compare some open-source software able to generate DEMs from Pléiades stereo images such as, among others surface extraction through triangular irregular network-based minimisation (SETSM) (Noh and Howat, 2014b), the Satellite Stereo Pipeline (s2p) (De Franchis et al., 2014) and the Ames Stereo Pipeline (ASP) (Moratto et al., 2010).

Pléiades DEM evaluation
The metrics used to describe the quality of the DEMs are the percentage of data voids and various statistics on the elevation differences between the Pléiades DEMs and the reference (GNSS or lidar) data: (i) their mean, µ, and median to evaluate the vertical accuracy of the DEMs and (ii) their standard deviation (SD) and normalised median absolute deviation (NMAD) to characterise their vertical precision.The NMAD is a metric for the dispersion of the data (also at the 1σ confidence level) that is not as sensitive to outliers as the SD and is recommended to evaluate DEM precision (Höhle and Höhle, 2009).
where h j denotes the individual errors and m h is the median of the errors.All statistics were computed after horizontal coregistration of the DEMs to the reference data.These horizontal shifts were obtained (i) by minimising the SD of the elevation differences when two DEMs were compared (Berthier et al., 2007;Rodriguez et al., 2006); (ii) in other cases by estimating the shift between a Pléiades orthoimage and the GNSS measurements acquired along certain trails and roads clearly visible in the images (e.g.Wagnon et al., 2013).When GCPs were used to compute the DEMs, we always verified visually that no horizontal shift of more than one to two pixels (i.e.0.5 to 1 m) remained between the Pléiades orthoimages and the GNSS tracks acquired along roads and trails.Thus, no planimetric correction was required.For Tungnafellsjökull, the GCPs were derived from a shaded relief image of the lidar DEM and a small shift remained between the Pléiades and the lidar DEM (see Sect. 4.2).
All elevation differences between the Pléiades DEMs and the reference data (GNSS surveys or lidar DEM) are assumed to be due to errors in the Pléiades DEMs.In fact the total error budget also includes (i) on glaciers only, real and spatially varying elevation changes over the time span of a few days or weeks separating the acquisition of the Pléiades images and the reference data and (ii) everywhere, errors in the reference data themselves.Error (i) can also matter off glacier if snow was present during the acquisition of the Pléiades images or the reference data.Therefore, the present study provides an upper limit for the errors of the Pléiades DEMs.When the time between satellite acquisition and field measurements exceeds a few weeks (e.g. the lidar DEM in Iceland), only reference data off glaciers are considered for the evaluation of the DEM, and data on glaciers are used only for the study of glacier elevation changes.

Glacier outlines
Glaciers outlines are needed for the Mont-Blanc area and the Icelandic study site to perform separate statistics on/off glaciers and to compute the glacier-wide mean elevation changes.For Mont Blanc, glacier outlines were drawn manually on a SPOT 5 2.5 m orthoimage acquired 23 August 2003.Very little seasonal snow remained when this image was acquired due to the heatwave that baked Europe in early August 2003.For Tungnafellsjökull, the margin of the ice cap was digitised manually from a shaded relief image of the August 2011 lidar DEM.The surrounding snow patches were obtained through semi-automatic classification applied to the corresponding lidar intensity image.For others study sites, GNSS data points were overlaid on the Pléiades orthoimages and visually classified as on/off glaciers so that no glacier outline was necessary.

Accuracy and precision of the Pléiades DEMs
The stable terrain around the Tungnafellsjökull ice cap (Iceland), which was snow-free on 7-8 August 2011 when the lidar DEM was acquired, is our primary site to evaluate the Pléiades DEM because of the extensive coverage and high accuracy provided by the lidar DEM.For this site, we explore the influence of the different processing parameters in PCI Geomatica (Sect.4.1), of the number of GCPs (Sect.4.2) and check the occurrence of spatially varying biases in the Pléiades DEM (Sect.4.3).The Pléiades DEMs for the other study sites were then evaluated (Sect.4.4).

Processing parameter settings
Our criteria for selecting processing parameter settings for the DEM generation were (i) the percentage of coverage with valid data (covered area column in Table 2) and (ii) the dispersion of the elevation differences around the mean and median (quantified using the SD and NMAD, respectively).
With the parameters "Type of relief" set to mountainous, "DEM detail" set to low and without filling data voids, the dispersion is just slightly larger and the area covered is greatly improved (nearly 99 % versus less than 93 %).All DEMs examined in the rest of this study were generated using these parameter settings, which seems to be the best compromise between a low percentage of data voids and a low dispersion of the elevation differences.We acknowledge that the differences obtained for different processing parameter settings are not very large and hence that other settings may be more appropriate in some cases.For instance, to map crevasses using a Pléiades DEM, a high or extra high level of detail would probably be a better setting.Filling data voids by interpolation leads to much larger errors (the SD is multiplied by a factor of 3, Table 2) and is not recommended except if a complete DEM is really needed (e.g. for orthorectification of an image).

Influence of the numbers of GCPs and tie points
Our five study sites are the target of glacier field measurements; therefore, GCPs are already available from static GNSS measurements of prominent features such as large boulders or crossing roads or could be measured in the future during dedicated campaigns.However, GCPs are not available in all ice-covered regions and it is important to assess their influence on DEM quality and determine whether useful DEMs can be retrieved in remote regions where no GCPs are available.At the time of DEM processing, no GCPs were available for the Astrolabe (Antarctica) and Mera (Nepal) study sites.The best GCP coverage was around Tungnafellsjökull where 19 evenly distributed GCPs were identified manually on a shaded relief image (pixel size of 2 m) derived from the lidar DEM.Uncertainties in the latter GCPs results from (i) errors in pointing manually identical features in the Pléiades images and the lidar shaded relief image and (ii) the horizontal and vertical errors in the lidar data.Tungnafellsjökull was thus the site chosen to test the influence of the number of GCPs.
For the Tungnafellsjökull site, the Pléiades DEM derived without any GCPs exhibits a horizontal shift of about 3.3 m and is also about 3 m too high on the average (Table 3).The Table 3. Influence of the collection of ground control points (GCPs) and tie point (TPs) on the accuracy of the Pléiades DEMs.Statistics are computed for the elevation differences (dh) between the Pléiades and lidar DEMs around and on the Tungnafellsjökull ice cap (Iceland).The Pléiades DEMs are computed using different numbers of GCPs and TPs.The parameter settings used to generate all the DEMs are as follows: DEM detail is low, type of terrain is mountainous, pixel size is 4 m, data gaps are not filled.The number of pixels used in these statistics is over 3 000 000. a after horizontal co-registration; b after horizontal and vertical co-registration, i.e. correction of the mean horizontal and vertical shifts between the Pléiades and lidar DEMs estimated on the ice-free terrain.

Nb of GCPs
addition of GCPs reduces the horizontal shift to about 2 m and the vertical shift nearly vanishes.In fact, a single accurate GCP appears to be sufficient to eliminate most of the vertical bias.In contrast, the horizontal shift is never entirely removed, even with 19 GCPs, possibly because a systematic shift may arise from GCP identification in the shaded relief image of the lidar DEM.The last column of Table 3 corresponds to the mean elevation difference between August 2011 and October 2013 on the ice cap after horizontal and vertical co-registration (referred to as 3-D co-registration hereafter) of the Pléiades and lidar DEMs.The similarity of these values (within 0.03 m) illustrates the effectiveness of 3-D co-registration.It implies that, if the subtracted DEMs include a sufficient proportion of stable areas (i.e.ice-free terrain), accurate elevation change measurements can be retrieved even without GCPs.
In the case of the Icelandic study site, the collection of tie points (TPs, i.e. homologous points identified in both images of the stereo pair but with unknown geographic coordinates) had no clear influence on the quality of the DEM: the vertical bias is increased by 0.5 m and the dispersion is slightly lower.In Sect.4.4, we will show that this conclusion does not hold for the other study sites.

Spatial pattern of the off-glacier elevation differences
In the previous sections, we assessed the overall accuracy of the Pléiades DEM on the whole ice-free terrain surrounding Tungnafellsjökull.However, past studies have highlighted the occurrence of spatially varying elevation errors in ASTER and SPOT 5 DEMs due notably to unrecorded variations in satellite attitude (Berthier et al., 2012;Nuth and Kääb, 2011).To quantify these errors, we split the map of elevation differences between the Pléiades DEM (computed using 5 GCPs) and the lidar DEM into X × X tiles (with X the number of tiles in each direction, varying from 2 to 5) and computed the median elevation difference on the ice-free terrain of each tile.The median was preferred here because it is a metric of centrality less affected by outliers (Höhle and Höhle, 2009).The results are shown in Fig. 3 for X = 3.The maximum absolute departure from 0 is observed for the north-west tile where the median elevation difference between the Pléiades and lidar DEMs reaches  −0.15 m (µ = −0.20,SD = 0.79, N = 141 754), followed by the south-east tile (median = 0.10 m, µ = 0.12, SD = 0.37, N = 12 246).These two tiles are also the ones with the most limited data coverage.This maximum absolute median elevation difference is 0.08, 0.26, and 0.24 m when X = 2, 4, and 5, respectively.These statistics show very limited spatially varying elevation differences between the Pléiades and lidar DEMs, implying that, within each quadrant, elevation changes of a sufficiently large ice body could be retrieved with an uncertainty of about ±0.3 m or less.

Evaluation of the Pléiades DEMs from all study sites
Table 4 summarises the results of the evaluation of the Pléiades DEMs with GNSS measurements for all study sites.In the Andes and for Mont Blanc, 5-13 GCPs were available to compute the DEMs.
For Mera Glacier in Himalaya no accurate GCPs, i.e. measured in the field using static GNSS positioning, were available at the time of processing.Instead, a set of 22 GCPs was derived from a coarser resolution SPOT 5 data set (2.5 m orthoimage and 40 m DEM).These SPOT 5 data were previously co-registered to GNSS data acquired along the trails of the Everest base camp, outside of the Pléiades images (see Wagnon et al., 2013, for a complete description).Because of the coarse resolution of the SPOT 5 DEM we tried as much as possible to select GCPs over flat terrain.The horizontal precision of these GCPs is limited by the SPOT 5 pixel size (2.5 m) and their vertical precision is about ±5 m, the precision of the SPOT 5 DEM.For Astrolabe Glacier, no GCPs were available at the time of processing.

Precision of the DEMs
The vertical precision (quantified using the SD and the NMAD) is relatively homogeneous for all sites and generally between ±0.5 and ±1 m (NMAD values ranging from 0.36 to 1.1 m and SDs from 0.51 to 1.26 m).These precision values are consistent with a recent study for a small glacier in the Pyrénées (Marti et al., 2014) and slightly better than those obtained using Pléiades stereo pairs in two other studies (Poli et al., 2014;Stumpf et al., 2014).For the relatively steep and vegetated terrain around landslides in the southern French Alps (Stumpf et al., 2014), the precision of the Pléiades DEM is around ±3 m.For the urban landscape around the city of Trento in Italy (Poli et al., 2014), it is ±6 m or more.The seemingly lower precision in other studies is probably not due to differences in the processing of the Pléiades images but more likely due to the influence of the landscape on the DEM precision.A quasi-linear relationship has been found between DEM precision and terrain slope (e.g.Toutin, 2002).It is also problematic to extract precise DEMs in urban areas (e.g.Poli et al., 2014).We would therefore expect to obtain a better precision on smooth glacier topography.This is confirmed by the results for the two study sites (Agua Negra and Tungnafellsjökull), where GNSS data have been collected on and off glaciers (Table 4).The precisions are always higher on glaciers.The improvement is particularly spectacular on the Tungnafellsjökull study site where the SD of the elevation difference is 0.53 m on the ice cap and 1.33 m elsewhere.The decrease of the SD is not as strong for the Agua Negra study site, possibly because the glacier presented a rough topography (0.5 to 1 m high penitents) at the end of the ablation season when the Pléiades images and GNSS measurements were acquired.
It is notoriously difficult to extract reliable elevation measurements in the flat and textureless accumulation areas of glaciers using stereo photogrammetry.In our Pléiades DEMs, only a limited percentage of data voids is present in these accumulation areas which suggests a good correlation between the images.This is confirmed by the vertical precision of the DEMs which reaches the same level as in the ablation areas.For example, in the Mont-Blanc area, 16 GNSS points were measured in 2012 above 4000 m a.s.l. in the Col du Dôme area (Fig. 2, southernmost points on the Mont Blanc Pléiades image), well above the regional equilibrium line altitude of ∼ 3100 m a.s.l.during the last 10 years (Rabatel et al., 2013).For these 16 points, the mean bias of the August 2012 Pléiades DEM is 0.19 m (median = 0.14 m) and the SD of the elevation difference is 0.40 m.
These quantitative assessments are confirmed qualitatively when examining elevation contours and shaded relief images generated from the DEMs (Figs. 4 and 5).flects the smoothness of the Pléiades DEMs for the upper accumulation basin of the Mer de Glace (the so-called Glacier du Géant) and the Astrolabe Glacier in Antarctica.The noise level is much larger in the SPOT 5 DEMs of these areas.
Together with the larger dynamic range (12 bits for Pléiades vs. 8 bits for SPOT 5), we suggest that the higher resolution of Pléiades allows capturing some fine-scale surface features that facilitate the matching of the images.

Accuracy of the DEMs
When GCPs are available, the vertical bias (i.e. the mean or median of the elevation differences) is generally lower than 1 m.The GNSS surveys on glaciers used to evaluate the DEMs are not exactly temporally coincident with the acquisition dates of the Pléiades images (Table 1).Hence, part of these vertical biases may be explained by real (but unknown) glacier elevation changes during this time interval.For example, in the case of the Mont Blanc 2012 DEM, the ca. 1 m elevation difference could be easily explained by the thinning that likely occurred between 19 August 2012 (Pléiades DEM) and 5-8 September 2012 (GNSS survey) on the rapidly thinning tongues of Argentière Glacier and Mer de Glace.Ablation field measurements performed on stakes on Argentière Glacier between 16 August and 5 September 2012 gave values of −0.98 m water equivalent (w.e.) at 2400 m a.s.l., −0.75 m w.e. at 2550 m a.s.l. and −0.60 m w.e. at 2700 m a.s.l.These ablation values, measured during a similar time period, approach the 1 m elevation difference observed between the Pléiades DEM and the GNSS survey but this agreement can only be considered as a general indication because glacier elevation changes are the combination of surface mass balance and ice dynamics processes.
Without GCPs, the vertical biases of the DEMs are larger: about 1 m for the Agua Negra study site, 2 m for the Astrolabe Glacier (Antarctica) and as much as 7 m for the Au-gust 2012 Mont Blanc DEM.These results clearly demonstrate the necessity to vertically adjust the Pléiades DEMs built without GCPs on stable terrain before using them to retrieve elevation changes (Paul et al., 2014).
For the Agua Negra study site, we obtained similar vertical biases between the three different versions of the Pléiades DEMs (front/nadir, nadir/back, front/back) computed without GCPs, with mean vertical biases ranging from 0.99 to 1.33 m.Conversely, with our 5 GCPs, the mean vertical biases range from −0.33 to 1 m.There are two possible explanations for this.First, the coordinates of the GCPs were calculated using a GNSS base station located as far as 100 km away and are more uncertain than for other study sites.Second, the identification on the Pléiades images of the features (GCPs) measured in the field (e.g.large boulders) was sometimes ambiguous.

Necessity of tie points
We already mentioned above that TPs had no influence on the coverage and the precision of the DEM of the Tungnafellsjökull ice cap (Table 3).However, this does not hold for the Mont-Blanc area and the Astrolabe Glacier, two other sites where DEMs were generated without GCPs.In both cases, the automatic collection of about 20 TPs provided far better coverage probably by improving the relative orientation of the two images of the stereo pairs.The necessity of collecting TPs will depend on the accuracy of the navigation data (position and attitude of the satellite) provided with the images.Given that the latter accuracy is not known a priori, we recommend collecting TPs between the images.

Added value of a tri-stereo
There is a moderate added value of a tri-stereo instead of a simple stereo pair.In general, among the three possible pair combinations of the three images, the largest percentage The Cryosphere , 8, 2275-2291, 2014 www.the-cryosphere.net/8/2275/2014/ of data gaps is observed for the front/back pair, due to the stronger distortion between these images.This is especially true when the base-to-height ratio is high (> 0.5), for example in the case of the Mont Blanc 2013 tri-stereo where the data voids represent as much as 30 % for the nadir/back pair and only 15 % for the front/nadir and nadir/back pairs.For the latter two pairs, we combined the two DEMs in a merged DEM as follows: (i) for pixels where both DEMs were available, the mean of the two values was calculated; (ii) for pixels where only one DEM was available, this single value was retained; (iii) for pixels corresponding to data voids in both DEMs, a data void was kept in the merged DEM (i.e.no gap filling was used).The percentage of data voids in the Mont Blanc 2013 merged DEM dropped from 15 to 6 %, showing that data voids were not at the same location in the front/nadir and nadir/back DEMs.For Agua Negra Glacier (Andes), the same merging reduced the percentage of data voids by only 1 % (but the initial coverage in the DEMs was higher, over 96 %) with no significant improvements in vertical precision.

Seasonal elevation changes
A comprehensive GNSS survey of Tungnafellsjökull ice cap was performed using a snowmobile on 2 May 2013, 5 months before the acquisition of the Pléiades stereo pair.The Pléiades DEM, first co-registered horizontally and vertically to the lidar DEM, was compared to GNSS elevations to reveal the surface elevation changes during the 2013 melt season between 2 May and 9 October.As expected, the surface was lower in October due to snow, and to a lesser extent, firn compaction and surface melt during summer (Fig. 6).On the average, the surface lowering between May and October 2013 was 3.1 m (SD = 0.89 m, N = 4800).Part of this lowering is also due to ice dynamics in the accumulation area, whereas ice motion (i.e.emergent velocity) only partly counteracts the strong lowering due to melting in the ablation zone in summer.A clear pattern with elevation is observed, with greater thinning at lower elevations close to the margins of the ice cap (inset of Fig. 6).Precise elevation changes are available at cross-over points between the extensive GNSS survey in May 2013 and the more limited GNSS coverage in September 2013 (Fig. 2).At those locations, there is good agreement (within 0.4 m) between the Pléiades-GNSS and repeated GNSS elevation change measurements.

Annual elevation changes
For the Mont-Blanc area, two Pléiades DEMs were available with a time difference of slightly more than a year (19 August 2012 and 20 September 2013).These two DEMs were first co-registered by minimising the SD of their elevation differences on the ice-free terrain (e.g.Berthier et al., 2007).The corrected horizontal shifts were 1 m in easting and northing.The remaining vertical shift on the ice-free terrain after horizontal co-registration was only 0.1 m (median of the elevation differences) and the dispersion (NMAD) was 1.79 m.These very low horizontal and vertical shifts could be expected given that the 2012 and 2013 DEMs were built using the same set of GCPs.This negligible median elevation difference off glaciers tends to confirm that the 19 August 2012 DEM is not biased and that the ca. 1 m average elevation difference between the 19 August 2012 DEM and the 5-8 September 2012 GNSS measurements (Table 4) is due to real elevation changes, not errors in the 2012 Pléiades DEM.
Once co-registered in 3-D, the DEMs were differentiated to map the glacier elevation changes that occurred over the 13 months between 19 August 2012 and 20 September 2013 (Fig. 7).Due to the difference in seasonality, the glaciological interpretation of these changes, as well as their comparison to field measurements (performed annually around 10 September), is not straightforward.However, the map reveals a clear thinning for all glacier tongues, whereas thickening is observed on most glaciers above 3000 m a.s.l., with some localised elevation gains of over 5 m, probably due to avalanches.This high-elevation thickening cannot be confirmed by field measurements but is in line with the slightly above-average accumulation during 2012-2013 (unpublished GLACIOCLIM-LGGE, data).We did not attempt to compute a glacier-wide or region-wide annual mass balance over such a short time span (13 months) since it would likely be skewed by seasonal variations and because of the high uncertainties that would arise from the poorly constrained density of the material gained or lost for such a short period of time (Huss, 2013).Despite these issues, this result highlights the very strong potential of repeat Pléiades DEMs for accurate mapping of glacier elevation changes, even over short time periods.

Multi-annual elevation changes
To fully explore the potential of repeated high-resolution satellite DEMs for measuring glacier elevation changes and glacier-wide mass balances, the 19 August 2012 Pléiades DEM was next compared to a 10 m DEM derived from SPOT 5 images acquired 19 and 23 August 2003 over the Mont-Blanc area.The 19 August 2012 Pléiades DEM was preferred to the 20 September 2013 DEM because it was acquired at the same time of year as the SPOT 5 DEM and contains less data voids.To entirely cover the Mont-Blanc area, a Pléiades DEM was derived in its southern part from a second Pléiades stereo-pair also acquired 19 August 2012 (Fig. 2).This southern DEM was computed using only 2 GCPs, shared with the northern stereo-pair.The mean elevation difference between the southern and the northern Pléiades DEMs off glacier was 0.27 m (median = 0.30, SD = 0.98 m) and the horizontal shift was 0.2 m.After 3-D co-registration off glacier, the mean elevation difference on glacier between those two synchronous DEMs is 0.07 m (median = 0.06, SD = 1.20 m).The full 2012 Pléiades DEM is used as reference DEM for 3-D co-registration of the two DEMs on the stable terrain.The 2003 DEM exhibits a 4.6 m total horizontal shift (2.5 m in easting and 3.9 m in northing) and is, on the average, 2.3 m lower than the 2012 DEM.Finally, a plane is fitted to the map of elevation differences in order to remove a spatially varying bias between the DEMs, with maximum values of −2 to 2 m.The latter correction had a negligible impact (less than 0.02 m a −1 w.e.) on the region-wide mass balance because the Mont Blanc lies in the middle of the area covered by the two DEMs.After 3-D coregistration, 9 full years of elevation changes in the Mont-Blanc area are depicted (Fig. 8).
These satellite-derived elevation changes are compared to the mean elevation changes derived from GNSS measurements performed every year around 10 September by LGGE on nine transverse profiles (five on the Mer de Glace and four on the Argentière Glacier, Fig. 8).Due to the retreat of the front of the Mer de Glace, the lowest profile, Mottet, has been deglaciated since 2009 and could not be used in our comparison.The 2003-2012 elevation differences derived from satellite data are, on the average, only 0.3 m higher than those measured in the field (SD = 1.3 m, N = 8).To evaluate how the two satellite DEMs contribute to the dispersion of the elevation differences (±1.3 m), we directly extracted the DEM elevations at the location of the GNSS transverse profiles for each DEM separately.The mean elevation difference between the 2003 SPOT 5 DEM and the 2003 field data is 0.5 m (SD = 1.3 m, N = 8), and the mean elevation difference between the Pléiades DEM and the 2012 field measurements is 0.8 m (SD = 0.4 m, N = 8).As expected from its higher resolution, the Pléiades DEM is more precise than the SPOT 5 DEM with a SD three times lower.The fact that both satellite DEMs are higher than the GNSS profiles can be explained by glacier thinning between the DEM acquisition dates (around 20 August) and the dates of the field surveys (around 10 September) in late summer when strong ablation (and thus thinning) is still ongoing in the European Alps.
For each 50 m altitude interval, the histogram of the elevation changes is computed.The distribution is approximated by a Gaussian curve, which permits the calculation of the mean thickness change as the average of all the values less than 3 SDs from the mode of the Gaussian curve (Berthier et al., 2004;Gardner et al., 2012).Where no elevation change is available for a pixel, we assign to it the value of the mean elevation change of the 50 m altitude interval it belongs to, in order to assess the mass balance over the whole glacier area.Conservatively, the SD of the elevation differences at the eight transverse profiles (±1.3 m) is used as our error es-timate for the 2003-2012 satellite-derived elevation differences.For un-surveyed areas, this elevation change error is multiplied by a factor of 5, resulting in an error of ±8 m.The percentage of data voids equals 15 % for the whole Mont-Blanc area and range from 2 to 22 % for individual glaciers (Table 5).Elevation differences are converted to annual mass balances using a density of 850 ± 60 kg m −3 (Huss, 2013).
The resulting glacier-wide geodetic mass balance for Argentière Glacier (−1.12 ± 0.18 m a −1 w.e.) be-   5).Uncertainty for the glaciological mass balance is from Thibert et al. (2008).The glacier-wide mass balances for 10 selected glaciers are all negative (Table 5) as the region-wide mass balance of −1.04 ± 0.23 m a −1 w.e.A clear acceleration of the mass loss is observed for all glaciers that were measured both in 1979-2003and 2003-2012 (Table 5 (Table 5).These results reflect the strong mass loss that has occurred in the Mont-Blanc area over the last decade, in agreement with recent studies elsewhere in the European Alps (Abermann et al., 2009;Carturan et al., 2013;Gardent et al., 2014;Huss, 2012;Kropáček et al., 2014;Rabatel et al., 2013).

Summary and conclusions
So far, little work has been carried out based on Pléiades images over glaciers.Our evaluation over five different glacial environments demonstrates that Pléiades stereo images are a promising tool for the monitoring of glacier topography and elevation changes.Overall the precision of these DEMs (at the 1σ confidence level) is ca.±1 m, sometimes better (±0.5 m) for the flat glacier tongues, a result in agreement with a study on a small glacier in the Pyrénées (Marti et al., 2014).The coverage and precision of the accumulation areas is also promising.The higher precision on glaciers compared to the surrounding ice-free terrain implies that an error estimate performed on the ice-free terrain will be conservative.Vertical biases are greater (as much as 7 m) if no GCPs are available but can be greatly reduced through proper 3-D coregistration of the Pléiades DEMs with a reference altimetry data set on ice-free terrain.One or two accurate GCPs seem sufficient to reduce the vertical biases to a few decimetres.
There is a slight improvement of the DEM coverage when they are derived from a tri-stereo.We have shown 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 precision of the merged DEM.However, because glacier topography is often relatively smooth, a standard stereo coverage with a limited difference in incidence angles (typically a base-to-height ratio of about 0.35-0.45)provides a relatively comprehensive and cost-effective coverage of the glacier surfaces.
One strong advantage of DEMs derived from Pléiades (and from other optical stereo sensors) compared to DEMs derived from radar images (such as the SRTM and Tandem-X DEMs) is the absence of penetration into snow and ice.Thus, all measured elevation differences correspond to real ice and snow elevation changes.Given their accuracy, DEMs derived from Pléiades (or other similar optical sensor) could be used in the future to check the magnitude and spatial pattern of the penetration depth of the Tandem-X radar signal into snow and ice, if temporally concomitant acquisitions can be found in the image archives.As for all optical sensors, the main drawback of the Pléiades constellation is the need for clearsky conditions to obtain suitable cloud-free images.
Our results open some promising perspectives.In the future, the differencing of Pléiades DEMs acquired ∼ 5 years apart could make it possible to determine glacier-wide mass balances with an uncertainty of ±0.1 to ±0.2 m a −1 w.e.Such an error level is sufficiently low to check the cumulative glaciological mass balances measured in the field (e.g.Zemp et al., 2013) and explore the spatial variability of glacier-wide mass balances at the scale of a glaciated massif (Abermann et al., 2009;Soruco et al., 2009a).It is already possible to differentiate recent Pléiades and older DEMs to provide accurate glacier-wide and region-wide mass balance.In our study, Pléiades and SPOT 5 DEM differencing was used to measure the strongly negative mass balance of the entire Mont-Blanc area.Pléiades DEMs acquired at the beginning and end of the accumulation seasons could probably be used to map snow thickness if the ice dynamics component can be estimated.If this is confirmed, Pléiades will represent a good alternative to recently developed techniques based on lidar (Deems et al., 2013;Helfricht et al., 2014) and stereo-photogrammetry (Bühler et al., 2014), especially for remote areas where acquiring airborne data can be challenging.Still, the conversion of glacier elevation changes measured over short time periods (1 season or 1 year) to glacier-wide mass balances will remain a complicated task due to the lack of knowledge of the actual density of the material (ice-firn-snow) gained or lost.
Apart from their cost and their sensitivity to cloud coverage, the main limitation of Pléiades images is their relatively limited footprint, typically 20 km × 20 km for a single scene.No large-scale stereo mapping has yet been planned using these two satellites and the cost of covering all glaciers on Earth (> 700 000 km 2 ) would be very high.For mapping vast glaciated areas, the recently launched SPOT6 and SPOT7 satellites may prove to be a good compromise given their resolution (1.5 m) and wider swath (60 km).Like Pléiades, they benefit from a very broad radiometric range (12 bits), avoiding saturation in most cases and improving contrast on snow-covered areas.However, the accuracy of the DEMs that can be derived from these stereo images has yet to be demonstrated over glaciers, ice caps and ice sheets.

Figure 1 .
Figure 1.Study sites where Pléiades stereo pairs and tri-stereos were acquired.The background image is a MODIS mosaic from the Blue Marble Next Generation project.

Figure 2 .
Figure 2. Pléiades orthoimages for the five different study sites.Yellow dots locate the GNSS measurements used to evaluate the DEMs.For Tungnafellsjökull, the blue polygon marks the limits of the lidar DEM.For the Mont-Blanc area, the red rectangles separate two stereo pairs acquired the same day which are combined to cover the entire glacier complex.For Astrolabe Glacier, an enlargement of a 750 m by 600 m area in the upper part of the glacier (elevation of > 700 m above ellipsoid) is shown.

Figure 3 .
Figure 3. Map of the elevation differences between the lidar DEM (7-8 August 2011) and the Pléiades DEM (9 October 2013) of the Tungnafellsjökull Ice Cap.The margin of the ice cap is shown by a thick black line and snow patches are outlined with a thinner black line.On the ice cap, the elevation contour lines are drawn as thin dashed lines every 100 m (from 1000 to 1500 m a.s.l.).The study area has been divided into 3 × 3 tiles in which the median of the elevation differences on the ice-free terrain only is reported (in metres).Background: Pléiades image ( © CNES 2013, Distribution Airbus D & S).

Figure 4 .
Figure 4. Comparison of the 21 August 2003 SPOT 5 and 19 August 2012 Pléiades DEMs of the accumulation basin of the Mer de Glace.The upper panels show the SPOT 5 data and the lower panels the Pléiades data.From left to right, the panels show successively the satellite images, the DEMs with the 50 m elevation contour lines and the shaded relief images derived from the DEMs.Note the higher percentage of data voids and artefacts in the SPOT 5 DEM.

Figure 5 .
Figure 5.Comparison of the 14 December 2007 SPOT 5 and 6 February 2013 Pléiades DEMs of Astrolabe Glacier.(a) shows the Pléiades image.The two right panels show shaded relief images derived from the Pléiades DEM (b) and the SPOT 5-HRS DEM (c) with the 100 m elevation contour lines overlaid.

Figure 6 .
Figure 6.Elevation differences between the kinematic GNSS data (2 May 2013) and the Pléiades DEM (9 October 2013) on the Tungnafellsjökull ice cap.Black circles indicate the locations where elevation differences have been measured using repeated GNSS surveys (2 May 2013 and 18 September 2013).Yellow crosses show the location of most of the 19 GCPs used to generate the DEM (two are masked by the inset and the colour scale).Inset: distribution of these elevation differences with altitude, with repeat GNSS surveys shown as larger dots.Background: shaded relief image derived from the Pléiades DEM.

Figure 7 .
Figure 7. Elevation differences between the Pléiades DEMs of 19 August 2012 and 20 September 2013 over the Mont-Blanc area.The outlines of glaciers in August 2003 are shown as black lines.Yellow crosses show the location of GCPs.The southernmost triangle locates the summit of Mont Blanc (4810 m a.s.l.) and the northernmost triangle the summit of Aiguille Verte (4122 m a.s.l.).The grey background SPOT 5 image ( © CNES 2003, Distribution Spot Image) appears where no elevation change value is available.

Figure 8 .
Figure 8. Elevation differences between the SPOT 5 DEM from 19 to 23 August 2003 and Pléiades DEM from 19 August 2012 over the Mont-Blanc area.In yellow, the location of the transverse profiles where elevations are measured every year using differential GNSS.The field (noted GNSS) and satellite (SAT) 2003-2012 elevation differences averaged along these profiles are indicated.Inset: satellite-derived (SAT, small symbols) and field (GNSS, large symbols) elevation changes as a function of altitude for the Mer de Glace (blue) and the Argentière (grey) glaciers.Large symbols correspond to the field measurements.

Table 2 .
Influence of different processing parameter settings on the coverage and accuracy of the Pléiades DEMs.Statistics for the elevation differences between the Pléiades and lidar DEMs are computed for ∼ 3 000 000 data points on the ice-free terrain around the Tungnafellsjökull ice cap (Iceland).The parameter settings tested are as follows: type of terrain is mountainous (Mtn) or flat, DEM detail is low or high, data gaps are filled or not filled.All Pléiades DEMs are computed using 5 GCPs and with a final pixel size of 4 m.The table provides the horizontal shifts of the Pléiades DEMs and, after horizontal co-registration (i.e.correction of the mean horizontal shift between the Pléiades and the lidar DEMs on the ice-free terrain), statistics (mean, median, SD and NMAD) of the elevation differences (dh, Z Pléiades − Z lidar ) outside the ice cap (OFF ice).after horizontal co-registration.The parameter that was changed compared to the first row is in bold. *

Table 4 .
Statistics on the elevation differences between the Pléiades DEMs and the GNSS measurements for all study sites.When a tri-stereo was available, the statistics for the three different image combinations and for a merged DEM are given.For Agua Negra and Tungnafellsjökull, the statistics are also given separately on and off glaciers.The last column, N , indicates the number of points for which elevation differences are computed.

Table 5 .
(Berthier, 2005)lances for the entire Mont-Blanc area and its 10 largest glaciers, sorted by size.Where available, the 2003-2012 geodetic mass balances calculated in this study are compared to the 1979-2003 geodetic mass balances(Berthier, 2005).The glaciological mass balance is provided for Argentière Glacier.