Glacier-wide seasonal and annual geodetic mass balances from Pléiades stereo images: application to the Glacier d'Argentière, French Alps

Abstract

The increased availability of high-resolution and high-quality digital elevation models (DEMs) allows for the investigation of small-scale glaciological changes and improved precision in geodetic mass-balance estimates. However, high precision and careful methodological choices are required to retrieve glacier-wide mass changes at annual to seasonal time scales. Here, we used a 7-year time series of 12 DEMs of the Glacier d'Argentière, in the French Alps, derived from the Pléiades optical satellites to assess the ability of sub-metre stereoscopic satellite images to retrieve annual-to-seasonal mass balances. We found good agreement between the five annual and the four winter mass-balance values estimated using a geodetic method and those of in situ glaciological measurements: mean values via the geodetic method are −0.66 m w.e. and 1.47 m w.e. for annual and winter balances, respectively; mean absolute discrepancies are 0.25 m w.e. (annual) and 0.36 m w.e. (winter). Our study identified three main limitations of this methodology: (i) the intrinsic DEM precision; (ii) the lack of control over the satellite acquisition dates; and (iii) the density assumption. The consistency between the methods demonstrates the potential of short time-scale glacier mass-balance monitoring using very high-resolution satellite images.

1. Introduction

Changes in glaciers have important social and environmental consequences at local to global scales (Rasul and others, 2019; Hock and others, 2019). Glaciers have been globally losing mass for several decades as shown by in situ and remote-sensing measurements (e.g. Zemp and others, 2019; Hugonnet and others, 2021). The glacier mass balance is a relevant variable linked to climate (e.g. Vincent and others, 2009). The glacier-wide surface mass balance has historically been quantified by extrapolating in situ point measurements based on the glaciological method. It implies that large areas and/or inaccessible locations cannot be surveyed using in situ methods as they require important effort and means (Zemp and others, 2019). The geodetic method is a complementary approach which relies on digital elevation models (DEMs) usually acquired from airborne or spaceborne platforms both over single glaciers or large areas (e.g. Basantes-Serrano and others, 2018; Dussaillant and others, 2018; Pelto and others, 2019). While the glaciological method only measures surface mass balances, the geodetic mass balance includes all volume-changing processes (surface, internal and basal mass changes but also densification, e.g.; Thibert and others, 2008). The comparison between glaciological and geodetic mass balances is thus not completely straightforward. Geodetic mass balances are also used to validate and calibrate time series of glaciological measurements, and thus reduce the systematic errors in glaciological measurements (Zemp and others, 2013).

Geodetic mass balances with a time interval of several years have been widely quantified with satellite data at different resolutions and different spatial scales (Berthier and others, 2014; Paul and others, 2015). However, geodetic mass balances from spaceborne sensors are seldom calculated at annual and seasonal scales. Belart and others (2017) quantified the winter mass balance from optical sub-metre stereo images of Pléiades satellites over an Icelandic ice cap. Recently, other authors calculated geodetic mass balances of mountain glaciers at short time scales using high-resolution DEMs acquired from airborne laser scanning (Klug and others, 2018; Pelto and others, 2019). Two main limitations arise when calculating annual-to-seasonal geodetic mass balances. First, the volume-to-mass conversion factor (i.e. the density assumption) becomes uncertain for periods shorter than 5 years (Huss, 2013). Hence, the use of in situ density data and the delimitation of the different types of glacier surfaces are recommended (Pelto and others, 2019). The work of Belart and others (2017), focusing on winter mass balances, relies on a firn densification model and ice flow motion correction, which are not always accessible. Second, the precision of the DEMs must be sufficient to capture the signal of small elevation changes over a seasonal-to-annual period with limited uncertainties (Belart and others, 2017).

In this study, we tested whether high-resolution DEMs derived from commercial satellites provide the capability to measure the glacier-wide mass balance at annual and seasonal scales. We used a time series of 12 DEMs derived from Pléiades acquisitions from 2012 to 2020 covering the Glacier d'Argentière. This time series enabled the calculation of several annual and winter geodetic glacier-wide mass balances. By comparing cumulative and individual mass balances from both glaciological and geodetic methods, we assess to what extent Pléiades stereoscopic images can be used to quantify annual and seasonal mass balances.

2. Study area

Our study area is the Glacier d'Argentière (45°57′N, 6°59′E, 12 km² in 2015) located on the French side of the Mont-Blanc massif (Fig. 1). The glacier flows north-westwards from ~3600 to 1600 m above sea level (a.s.l.) and is ~9 km long. Part of the accumulation area is located on south-facing tributaries. The horizontal ice flow velocity at 2400 m a.s.l. above the serac fall is ~55 m a⁻¹ (Vincent and others, 2021). The lowermost part of the glacier tongue, at the north-western extremity, detached from the main body in 2011.

Fig. 1. Map of the study area, showing the Glacier d'Argentière (green outline delineation of 2015). AWS, Automatic Weather Station of GLACIOCLIM service. Background ortho-image: Pléiades 8 November 2018, © CNES.

We chose the Glacier d'Argentière as our study region for two main reasons. First, it is one of the most studied glaciers in the French Alps, including a dense measurement network of ablation stakes and snow pits that produces reference data for our study. It is part of the network of the French Service National d'Observation (SNO) GLACIOCLIM (https://glacioclim.osug.fr/) which maintains a long-term glaciological in situ monitoring network including seasonal surface mass balances, surface flow velocities, ice thickness measurements and snout coordinates (Vincent and others, 2009; Rabatel and others, 2018). Second, sub-metric satellite data are for now only accessible through orders to commercial services, and the time series of Pléiades sub-metre stereo images acquired since 2012 through the CNES-KALIDEOS Alpes project (https://alpes.kalideos.fr/) over this area is quite unique and well suited to our study goal. Although we do not aim to study specific glaciological processes, data collected here will enrich and complement ongoing data collection on the Glacier d'Argentière.

3. Data and methods

3.1. Glaciological mass balance

The glacier-wide glaciological mass balances and the point mass-balance data are provided by the GLACIOCLIM service, and are based on stake measurements or on snow pits/probes, depending on the season. For each point, we converted the height variation into mass by a density value, which is measured in situ for snow and a constant value for ice. More details on the measurement method can be found in Vincent (2002). While the location and number of point measurements can change from one year to another, measurements are normally taken at approximately nine points in the accumulation area and at 30 points in the ablation area (Fig. 1).

The GLACIOCLIM service calculated annual glacier-wide mass balances using a non-linear analysis of variance applied to each accumulation/ablation site measurement and then to each glacier location (Vincent and others, 2018). The resulting time series of the glacier-wide mass balance is available on the GLACIOCLIM website (https://glacioclim.osug.fr/). This time series is regularly calibrated with long-term geodetic mass balances (~10 years) using DEMs acquired from aerial optical photogrammetry or LiDAR (Vincent and others, 2009). However, this calibration does not match our study period from 2012 to 2020. We consequently recalibrated all glacier-wide annual surface mass balances quantified with the glaciological method using the geodetic mass balance quantified here with Pléaides images between 19 August 2012 and 25 August 2019. We adjusted this geodetic mass balance to the glaciological dates (13 October 2012 to 4 October 2019) with a degree-day melt estimation.

We calculated the glaciological glacier-wide winter mass balances, which we extrapolated from mass balances measured on the stakes to the entire glacier using a hypsometric method (area-altitude distribution; Cuffey and Paterson, 2010). We computed the mean mass balance from the stakes for each 200 m hypsometric band. We then computed the glacier-wide winter mass balance by weighting each band mean mass balance by its area. Despite the representative hypsometric distribution of the stakes, ~5% of the study area did not contain measurement data. More than half of this area is the lowermost part of the glacier, which is disconnected from the main body. For this part of the glacier, we assigned the winter point mass-balance value at its median altitude from a linear extrapolation of all stake mass balances according to their altitude. For these interpolations, the coefficients of determination range from 0.58 to 0.78 depending on the years, and the slope of the linear regression is ~2 m w.e. 1000 m⁻¹. Furthermore, for some winter mass-balance extrapolation, a few accumulation spots of the glacier were not measured on the same dates as the other in situ measurements. We adjusted the corresponding winter stake mass balances over the 7–42 days difference using a degree-day melt calculation. The melt adjustments ranged from 0.00 to 0.72 m w.e. Precipitation was not included in this adjustment and was estimated at between 130 and 660 mm at Chamonix weather station for the considered periods.

3.2. Calculating the geodetic mass balance

We calculated changes in glacier volume (in m³) from the sum of changes in pixel elevation (in m) on the maps of interpolated elevation differences. We then converted the volume change into a mass change (in kg). Following Pelto and others (2019), we calculated a volume-to-mass conversion factor (i.e. a density assumption) for the annual mass balance based on the proportion of the different types of surface (ice, and snow/firn, as the latter ones are not distinguished in the images), as well as in situ measurements of winter snow density. We set a value of ρ _ice = 910 ± 60 kg m⁻³ for the density of ice covering the surface S _ice and an estimation of ρ _snow/firn = 550 ± 100 kg m⁻³ for annual snow/firn density covering S _snow/firn. In contrast to Pelto and others (2019), we calculated a mean density according to the following formula:

(1)$$\rho _{{\rm mean}} = \rho _{{\rm ice}} \times \displaystyle{{S_{{\rm ice}}} \over {S_{{\rm total}}}} + \rho _{{\rm snow/firn}} \times \displaystyle{{S_{{\rm snow/firn}}} \over {S_{{\rm total}}}}.$$

We assessed the distribution of ice and snow/firn from a manual delineation on Landsat satellite images of the surfaces at the end of the ablation season. We could not use a specific firn density value to quantify the annual mass balance due to the lack of density data measured in situ as well as the small proportion of the surface of the glacier covered by firn compared to that covered by ice and annual snow. For the five annual mass balances and the 2013–15 mass balance, the calculated mean density over the whole glacier ranged from 721 to 783 kg m⁻³ (mean 756 kg m⁻³, median 0.774 kg m⁻³). For multi-annual mass balances, we used an estimated constant density of 850 ± 60 kg m⁻³ (Huss, 2013). The density we used for snow for the quantification of the seasonal winter mass balances was the mean of the densities measured in situ. We obtained average densities ranging from 430 to 520 kg m⁻³ depending on the year. For all the density values except the multi-year density, we assumed an uncertainty of ±100 kg m⁻³, approximately three times the std dev. of the snow density measurements.

For the whole time series, we used a glacier outline based on 2015 Sentinel-2 satellite images (Paul and others, 2020; Fig. 1) that we manually adjusted. As the glacier lost ~1% of its surface area each year during the study period, some margins of the glaciers are not captured by the delineation (in the years before the inventory was performed), or areas within the outline include parts of the terrain that have been de-glacierized (in the years after the inventory was performed). To account for the bias induced by false null signals or the lack of observations of glacier change, we divided the total change in volume by the estimated mean surface area between the two dates (Fischer and others, 2015). To estimate the surface area at each date, we derived a surface change rate of ~−0.10 km² a⁻¹ from a linear regression (coefficient of determination of 0.99) using delineations taken from 2008, 2015 and 2018 (Gardent and others, 2014; Paul and others, 2020; Rabatel, unpublished data).

To assess our annual and seasonal geodetic mass balances (internal and basal mass balance included), we compared them to glaciological surface mass balances. However, there is a gap of up to 2 months between the dates of the Pléiades acquisitions and glaciological measurements, and considerable melt can occur during this period. To compare the mass-balance values captured at different dates, we made two different estimates of the geodetic mass balances, with one of them adjusted over the glaciological mass-balance period. The first estimate considers no correction. The second estimate of the melt was based on the degree-day calculation detailed in section 3.4 below.

3.3. Time series of Pléiades stereo-images and digital elevation model processing

We used multi-annual stereoscopic sub-metre optical images from the Pléiades constellation that have been acquired over the Mont-Blanc area. The panchromatic band products have a spatial resolution of 0.5 m (resampled), and 2 m for spectral bands in the blue, green, red and near-infrared wavelengths (ADS Pléiades user guide, 2021). The 12-bit encoding of the images provides a better contrast of texture-less areas such as snow-covered accumulation areas. Together with the specific tuning of the gain parameter of the acquisition, the risk of saturation is thus reduced and DEM gaps are limited (Berthier and others, 2014). The mean base-to-height ratio (B/H) of stereo acquisitions is 0.34 (std dev. of 0.09). We processed the stereo images with the Ames Stereo Pipeline (ASP) to compute DEMs using a block matching method (Moratto and others, 2010; Shean and others, 2016). We produced non-interpolated DEMs with a spatial resolution of 4 m.

We used a time series of 12 DEMs from 2012 to early 2020 (out of the 18 acquired, six were outside the seasonal boundaries). Except for the years 2014 and 2020 for which no images were acquired in autumn, we were able to estimate five annual mass balances (Table 1). We quantified four winter mass balances with an image near the end of the accumulation season, from 2017 to 2020. Due to a 50% gap in the DEM dated 13 May 2019, we also used the DEM dated 22 March 2019 to reduce the potential impact of the gap. The 2019 winter geodetic mass balance was thus the sum of the mass balances from 8 September 2018 to 22 March 2019 and 22 March 2019 to 13 May 2019.

Table 1. Details of DEMs

The DEM dated 19 August 2012 was the reference for all co-registrations. The date format is YYYYMMDD.

^a Symbols next to the dates identify the DEMs used for annual mass balances.

^b Symbols identify the DEMs used for winter mass balances.

We then 3-D co-registered the DEMs over stable areas with a co-registration algorithm (Nuth and others, 2011), using an open-source python-based workflow (Shean and others, 2016; https://github.com/dshean/demcoreg). For the whole study, the reference DEM was the Pléiades DEM dated 19 August 2012 (chosen because of its low snow cover and the small proportion of gaps on the glacier). We used a three-step co-registration procedure, where we filtered out snow-covered areas (i.e. considered as non-stable terrain) at the third step only. This choice permits keeping enough surface for horizontal registration and for tilt and undulation correction, for which snow bias was considered negligible. We defined the stable area as the terrain surrounding the glacier after removal of the following zones: slopes of more than 40 degrees, main forested areas, and, for the third step, seasonal snow cover determined by an automated classification (adjusted manually when needed, notably for winter acquisitions). First, we performed horizontal co-registration, where off-glacier snow-covered terrain was considered stable terrain (Nuth and others, 2011). Second, after this co-registration, we identified tilt and undulation errors in some DEMs, which we corrected by estimating the errors using the same reference DEM on stable areas. We corrected the tilt by fitting a plane in both x and y directions to the median of pixel rows/columns: the final modelled tilt is <0.2 mm km⁻¹ amplitude along both axes. The undulations are the effect of high-frequency satellite jitter and affect both along and across-track directions with superimposed sinusoidal altitude errors (Girod and others, 2017). This effect has been reported for several optical stereo satellites, including Pléiades (Deschamps-Berger and others, 2020). In our case, we observed only along-track undulations; their amplitude sometimes reaches several tens of centimetres. We used a Fourier transform to filter out undulations of more than 1 500 m of wavelength (Deschamps-Berger and others, 2020). The frequencies and amplitudes may vary between DEMs. Third, after these corrections, we re-adjusted vertically the DEMs (second registration) over more restricting stable terrain (i.e. excluding seasonal snow cover; Table 1 and Fig. 2 on the right).

Fig. 2. Left: Map of the Glacier d'Argentière showing the number of interpolations between consecutive DEMs (total 11 pairs). Right: Snow-free stable areas used by three co-registrations or more. Spring areas nearly fully cover summer areas.

The mean std dev. of DEM difference on stable area after co-registration for annual mass balances is 0.57 m, which is similar to the results obtained by Berthier and others (2014). The DEMs we used for winter mass balances at the end of the winter season have larger std dev. (all >1 m except one), likely due to seasonal changes in, for example, vegetation, and to the smaller stable area used to calculate the std dev.. The mean std dev. for all pairs is 1.30 m (median 0.65 m).

Due to the local topography and the brightness of the accumulation area, at medium or high altitudes, DEMs were not always correctly generated over these areas (Fig. 2, left panel). The main sources of gaps were first the shadows projected on the glacier due to the topography, and second, the season and the angle of the sun at image acquisition time. The lack of data over the accumulation areas can create a bias by under-representing areas with a given elevation change (McNabb and others, 2019). Data gaps may also occasionally be created by the two-step filtering procedure we apply to each DEM difference map: we used ±15 m a⁻¹ as a threshold to remove outliers (±15 m for winter mass balances), and we used a minimum of pixel group threshold (perimeter threshold of seven pixels threshold) to avoid error-induced isolated pixels. To fill these data gaps, we used interpolation methods assessed by McNabb and others (2019). We implemented these methods using McNabb's own publicly available Python code (https://github.com/iamdonovan/dem_voids) and pybob library.

3.4. Degree-day adjustment

As the geodetic and glaciological measurements are not simultaneous, we used meteorological data for degree-day adjustments to account for the melt during in-between days. We did not correct the mass balances with precipitation, which are presented along temperatures in Supplementary Figure S1. We applied the adjustments to the geodetic mass balances to fit the dates of the glaciological surveys. We calculated the melt M (m w.e.) as:

(2)$$M = {\rm DDF} \times \mathop \sum \limits_{i = 1}^n T_i^ + {\times} \Delta t, \;$$

with DDF the in situ calibrated degree-day factor (Six and others, 2014), $T_i^ +$ positive air temperature during a period of n time intervals Δt (Hock, 2003). Daily temperatures have been recorded since 2007 by the GLACIOCLIM AWS on the right-hand margin of the glacier at ~2400 m a.s.l. (Fig. 1). We adjusted the temperatures locally to the altitude of the glacier using a lapse rate according to the season (Six and others, 2014). We calculated the lapse rates from the mean of 2016–2019 seasonal lapse rates between four AWS: two of the French national meteorological institute (Météo-France) in the nearby Chamonix valley (c.a. 1040 and 1500 m a.s.l.), the GLACIOCLIM AWS and the Aiguille du Midi Météo-France AWS (3850 m a.s.l.). The lapse rates we obtained are −0.65°C 100 m⁻¹ in May–June and −0.53°C 100 m⁻¹ in September–October. We took from Six and others (2014) the degree-day factors for all mass-balance adjustment. For geodetic glacier-wide adjustments, this factor is 4.0 × 10⁻³ m w.e. °C⁻¹ d⁻¹ for snow and 5.3 × 10⁻³ m w.e. °C⁻¹ d⁻¹ for ice, we considered to remain constant even though it may change with altitude. The degree-day factor we used for each adjustment was based on the share of the different types of surface (ice, snow/firn; see section 3.2 for the surface classification).

3.5. Uncertainty assessment

3.5.1. Geodetic uncertainties

We calculated geodetic uncertainties mainly using the method presented by Pelto and others (2019), assuming random uncertainties. In our case, we did not assess systematic uncertainties, as we artificially set the median elevation difference to zero on the stable terrain. We assumed random uncertainties to have three independent sources: elevation change uncertainty (σh _ΔDEM), delineation uncertainty (σA) and volume-to-mass density conversion uncertainty (σρ).

σh _ΔDEM depends to a great extent on the area over which we calculated the uncertainty, with larger areas leading to smaller uncertainties (e.g. Rolstad and others, 2009; Menounos and others, 2019). In this study, we faced the challenge of small co-registration areas for some pairs (from 0.19 to 37.56 km²; Table 1), which could introduce systematic biases in the elevation difference. To account for this potential bias, we calculated both the uncertainty due to the co-registration (σdh _coreg) and to the spatial averaging of the difference in height on the glacier (σdh _gla). We calculated σh _ΔDEM as the quadratic sum of these terms:

(3)$$\sigma h_{\Delta {\rm DEM}} = \sqrt {\sigma {\rm d}h_{{\rm coreg}}^2 + \sigma {\rm d}h_{{\rm gla}}^2 } .$$

We calculated both terms following the simplified equation of Fischer and others (2015):

(4)$$\sigma {\rm d}h = \left\{{\matrix{ {\sigma h\sqrt {\displaystyle{{A_{{\rm corr}}} \over {5A}}} \;; \;\;{\rm if}\;A > A_{{\rm corr}}} \cr {\sigma h\;; \;\;{\rm if}\;A < A_{{\rm corr}}, \;} \cr } } \right.$$

where A is the stable area used for vertical co-registration (for σdh _coreg) or the glacierized area (for σdh _gla), σh =0.52 m (the mean value of the stable terrain std dev., for the three stable areas larger than 30 km²; Table 1), and A _corr = πL ², L is the correlation length. We calculated L with the skgstat library (Mälicke and others, 2019) as the mean of ten rounds of the correlation length using a spherical variogram with 2000 random samples (maximum lag distance of 12 km). For annual mass balances, the mean individual correlation length is 1335 m (median value of 866 m), ranging from 687 to 3747 m. We retained L = 1335 m for all the DEM differences.

σA, the delineation error, is the glacier perimeter multiplied by two times the resolution of the satellite images used for the delineation (10 m). σρ is the different density uncertainties, i.e. 100 kg m⁻³ for surface-weighted and seasonal densities, and 60 kg m⁻³ for multi-year mass balances (Huss, 2013).

From these variables, we calculated the following volume uncertainty σ _ΔV (in m³), considering p the share of non-voided glacier area and then interpolated and assigning an uncertainty factor of 5 to this area (Berthier and others, 2014), h _ΔDEM is the mean elevation change over the glacier, and A _mean is the mean estimated area of the glacier between the two mass-balance dates:

(5)$$\sigma _{\Delta {\rm V}} = \sqrt {{( {\sigma h_{\Delta {\rm DEM}} \times A_{{\rm mean}}( {\,p + 5 \times ( {1-p} ) } ) } ) }^2 + {( {\sigma A \times h_{\Delta {\rm DEM}}} ) }^2} .$$

We then expressed the mass-balance uncertainty, σ _ΔM (in m w.e. per unit area), as:

(6)$$\sigma _{\Delta {\rm M}} = \sqrt {{( {\sigma_{\Delta {\rm V}} \times \rho } ) }^2 + {( {\sigma \rho \times \Delta V} ) }^2} \times \displaystyle{1 \over {A_{{\rm mean}}}}.$$

To assess the uncertainty after the degree-day adjustment, we calculated the uncertainty of each date adjustment σ _{adjust. indiv.}. For this purpose, we used an uncertainty σ _gradient for the temperature gradient used to calculate the temperature over the glacier. σ _gradient are 0.01 and 0.02°C⁻¹ 100 m⁻¹ for spring and autumn, respectively, from the std dev. of the gradient over the two seasons (see section 3.1). We neglected the degree-day factor uncertainty (order of 10⁻³ m w.e. °C⁻¹ d⁻¹). We did not include precipitation in the adjustment, but we accounted for larger uncertainties when precipitation occurs. For that, we multiplied the precipitation P (in m w.e.) recorded in Chamonix occurring during the N days between the acquisitions, by a factor 2 to scale the orographic effect (Vincent, 2002; Six and others, 2014). We assumed a conservative snow precipitation density over-estimation at 400 kg m⁻³:

(7)$$\sigma _{{\rm adjust}{\rm .\;indiv} .} = \sqrt {{\left({\sigma_{{\rm gradient}} \times \sqrt N } \right)}^2 + {\left({P \times 2 \times \displaystyle{{400} \over {1000}}} \right)}^2} .$$

The combined uncertainty of the adjustments at the beginning and at the end of the mass balance, σ _{adjust. total}, is the quadratic sum of the two individual uncertainties σ _{adjust. indiv.}. The final uncertainty of adjusted geodetic mass balances is thus the quadratic sum:

(8)$$\sigma _{{\rm \Delta M, \;adjusted}} = \sqrt {\sigma _{\Delta {\rm M}}^2 + \sigma _{{\rm adjust}{\rm .\;total}}^2 } .$$

3.5.2. Glaciological uncertainties on the winter surface mass balances

We calculated the seasonal surface mass balances from stake measurements at different altitudes, we calculated therefore an uncertainty for each altitude bin (see section 3.1). This uncertainty σ _bin is the std dev. of stake mass-balance measurements in the altitude bin (2–23 stakes per bin). We neglected the uncertainty of empty bins filled with neighbouring data accounting for 2% of the overall glacier surface area. For the bin on the glacier tongue, the uncertainty is the std dev. of hypsometrically extrapolated surface mass balance applied to the whole disconnected area. The first estimated surface mass-balance uncertainty (σ _M.B.seas.raw in m w.e.) is then

(9)$$\sigma _{{\rm M}{\rm .B}{\rm .seas}{\rm .raw}} = \sqrt {\mathop \sum \limits_{{\rm bin}} {\left({\sigma_{{\rm bin}} \times \displaystyle{{A_{{\rm bin}}} \over {A_{{\rm total}}}}} \right)}^2} , \;$$

where A _bin is the area of each bin. For another mountain glacier, Thibert and others (2008) found an uncertainty due to the lack of representative spatial sampling of σ _sampling =  0.12 m w.e. The total surface mass-balance uncertainty (σ _{M.B.seas.glacio} in m w.e.) is then:

(10)$$\sigma _{{\rm M}{\rm .B}{\rm .seas}{\rm .glacio}} = \sqrt {\sigma _{{\rm M}{\rm .B}{\rm .raw}}^2 + \sigma _{{\rm sampling}}^2 } .$$

We adjusted some stakes by degree-day melt estimation: we considered the degree-day uncertainty for the corresponding bins. The stake mass-balance uncertainty σ _stake is expressed as (7). For each bin, with n _stakes,bin stakes, the adjustment uncertainty is then $\sigma _{{\rm bin, \;adjusted\;stakes}} = ( \sigma _{{\rm stake}}/\sqrt {n_{{\rm stakes, bin}}} )$. With p the proportion of adjusted stakes per bin, the total band uncertainty is:

(11)$$\sigma _{{\rm bin, adjusted}} = \sqrt {\left({\,p\sqrt {\sigma_{{\rm bin, adjusted\;stakes}}^2 + \sigma_{{\rm bin}}^2 } } \right)^2 + ( {( {1-p} ) \times \sigma_{{\rm bin}}^2 } ) ^2} .$$

3.5.3. Glaciological uncertainties on the annual mass balances

The calculation uncertainty of the annual mass balances follows the analysis of Wagnon and others (2020), detailed in their supplements for non-linear models. The uncertainty formula for a mass balance inside calibration dates is the following:

(12)$$\sigma _{{\rm M}{\rm .B}{\rm .an}{\rm .glacio}} = \sqrt {\displaystyle{{\sigma _{{\rm M}{\rm .B}{\rm .geod} .}^2 } \over N} + \mathop \sum \limits_i S_i^2 \sigma _\varepsilon ^2 } , \;$$

where σ _M.B.geod. is the error of the calibrating geodetic mass balance, N the number of years between the DEMs, S _i the relative area in proportion to the total glacier area of the altitude bins used in the mass balance calculation and σ _ɛ the std dev. of the residuals of the non-linear method. As we lacked calculations details of the calibration of GLACIOCLIM glaciological mass balances for the period including 2012–19, we used the results of a 2003–2019 geodetic mass balance, whose DEMs processed from aerial data have been used for the GLACIOCLIM glaciological mass-balance calibration. As all the mass balances are for the period 2003–2019, we applied the previous formula. We used an uncertainty of the calibrating geodetic mass balance of ±1.24 m w.e. From a previous assessment of the non-linear model on the Glacier d'Argentière (Vincent and others, 2018), σ _ɛ is ±0.37 m w.e. Instead of bins, the non-linear model of this glacier was calculated by the GLACIOCLIM service from 0.2 × 0.2 km cells around the stakes (Vincent and others, 2018). The uncertainty we obtained for all glaciological annual mass balances is ±0.31 m w.e. We considered the uncertainty σ _i for each year i to be independent and following a normal distribution, so the uncertainty σ of the cumulative mass balance follows the equation:

(13)$$\sigma = \sqrt {\mathop \sum \limits_i \sigma _i^2 } .$$

Using (13) and a constant annual uncertainty, the uncertainty of the cumulative 7-year glaciological mass balance is simplified to $0.31 \times \sqrt 7$ = ± 0.82 m w.e.

4. Results

All the annual glaciological mass balances have been calibrated with our 2012–19 DEM pair: the resulting adjusted geodetic mass balance is −5.41 ± 0.80 m w.e. The cumulative non-calibrated glaciological mass balance over this time span is −7.79 ± 0.82 m w.e. The mass-balance adjustment added to each annual glaciological glacier-wide mass balance is +0.34 m w.e. a⁻¹.

4.1. Changes in elevation, volume and mass

Due to unconstrained changes in the density of the type of surface during the ablation season, we are not able to apply our volume-to-mass conversion factor for the summer period. Consequently, we focus on the winter and annual estimates. Figure 3 shows the variability in hypsometric elevation change (dh) rates between years and Figure 4 shows the spatial distribution of dh for each year and season. On an annual scale, dh is mostly negative, both in terms of hypsometric distribution (Fig. 3) and averaged at the glacier scale (Fig. 4). A notable exception is the 2015–16 year (Figs 3, 4c) when thickening of the upper reaches contrasts with the slight thinning of the main tongue. For each year, the disconnected part of the glacier tongue has the highest thinning rates. Overall, at annual scales, different patterns of elevation changes are observed. For some years (2012–13, and 2015–16), we observe a strong thinning at the lowermost elevations, and a thickening at higher elevations and at the foot of the north-oriented rock faces surrounding the glacier (Figs 3, 4a, c). For other years, the thinning pattern is relatively homogeneous over most of the glacier surface (Figs 4b, d, f). The bins ~3000 m a.s.l., mainly composed of the largest accumulation area, show the most variable behaviour between years.

Fig. 3. Elevation change (m) for each hypsometric bin. The two different line styles represent annual and winter changes, respectively. The grey areas represent the uncertainty (mean of the std dev.). The blue bar plot represents the surface of each 50 m hypsometric bin (0.6 km² being 5% of the total 12 km² of the glacier). No point or bin surface mass balance can be estimated from these elevation changes due to glacier dynamics.

Fig. 4. Surface elevation changes used to calculate the time series (bi-linear interpolation of the height difference). Thumbnails are showing: annual variations (a–f), winter variations (g–j), 7-year variation (k). All thumbnails, except one on the bottom right, refer to the scale on the bottom left. The date format is YYYYMMDD.

In winter, the mean glacier elevation changes are positive, ranging from 2.33 to 3.52 m on average at the glacier scale (glacier-wide geodetic mass balances from 1.02 to 2.05 m w.e., Fig. 5). The winter elevation gain is either homogeneous with elevation (Figs 4g, i, 3) or shows larger gains in surface elevation at higher elevations (Figs 4h, j). The DEM resolution is fine enough to see the effect of the serac falls depositing ice on the upper part of the disconnected lower section of the glacier tongue (particularly in the winter season), as well as the snow avalanche deposits in the accumulation areas (Figs 4g, h).

Fig. 5. Time series of cumulative mass balances with void proportions. Both blue and orange plots are from Pléiades DEMs. The red vertical lines show the mass-balance uncertainties (cumulative except for multi-year geodetic mass balances). The vertical bars show the gap proportions over the glacier. The annotated numbers correspond to the annual and seasonal geodetic mass balances and to the difference with the multi-year geodetic mass balance, when available. The dotted line represents the summer mass balances derived from the winter and annual mass balances only, and they are not estimated directly. The date format is YYYYMMDD.

Figure 5 shows the time series of cumulative geodetic mass balances produced from the surface elevation changes. For each year, we also calculate the geodetic mass balance relative to the first date of the time series (19 August 2012), using a volume-to-mass conversion factor of 850 ± 60 kg m⁻³, even for periods shorter than 5 years (Huss, 2013). To facilitate the comparison, we superpose the plots, and both glaciological and geodetic time series start at 0 m w.e., despite the 55-day interval (Table 2). Winter geodetic mass balances are not used for the calculation of the total cumulative geodetic mass balance, as summer mass balances are required for this purpose. We compare the two geodetic methods (the sum of all the computed annual mass balances and the ‘long-term mass balance’ with respect to 19 August 2012 DEM) for the possible six end-of-season dates. For the six dates in 2013, 2015, 2016, 2017, 2018 and 2019, the difference in the mass balance calculated with the two methods is 0.04, 0.12, 0.08, 0.29, 0.36 and 0.40 m w.e., respectively (mean 0.21 m w.e.). The mean error of the first three dates (0.08 m w.e.) differs significantly from the mean of the last three dates (0.35 m w.e.) with a p-value of 0.04. The DEM dated 25 October 2017 has an on-glacier void percentage of 31.1% (Table 1); which may explain this significant break in the quantified differences, and its propagation afterwards.

Table 2. Details of the comparisons between glaciological and geodetic mass balances

Precipitation was recorded by Chamonix weather station.

^a Identifies dates for which stakes of accumulation areas of different dates are included with degree-day adjustment (see section 3.1). The date format is YYYYMMDD.

4.2. Comparison of geodetic and glaciological mass balances

Table 3 compares geodetic and glaciological mass balances, according to the dates and adjustments shown in Table 2.

Table 3. Resulting mass balances corresponding to Table 2

^a Only annual glaciological mass balances are calibrated.

^b Glaciological glacier-wide mass balances provided by the GLACIOCLIM service are geodetically calibrated, over a period different from our study one.

All the glaciological mass balances agree with non-interpolated geodetic ones (i.e. differences are within the uncertainties; Fig. 6). They also agree with the degree-day adjusted geodetic mass balances except for two annual values (2016 and 2018). The mean absolute difference between the annual glaciological and geodetic mass balances is 0.47 m w.e. (winter: 0.42 m w.e.) for the ‘degree-day adjusted’ and 0.25 m w.e. (winter: 0.36 m w.e.) for the ‘not adjusted’ methods. Degree-day adjustments to fit the dates of the glaciological surveys have a mean impact of 0.46 m w.e. on the mass balances, for a mean number of 21 days between field surveys and satellite acquisition dates.

Fig. 6. Geodetic versus glaciological annual and winter mass balances. The type of symbol distinguishes ‘degree-day adjusted’ from ‘not adjusted’ values. The shaded areas correspond to the uncertainties.

4.3. Influence of different assumptions and corrections

4.3.1. Interpolation method

We assess the impact of the mass-balances quantification depending on the method of interpolation used to fill the gaps in the dh maps. We compare quantification without interpolation and with interpolation methods, hereafter called ‘bilinear interpolation of the elevation difference’ and the ‘mean 300 m neighbours’ (i.e. using the mean on-glacier difference within a radius of 300 m to fill in the gaps).

For the mass-balance time series quantified using different spatial interpolation methods, Table 4 lists the mean difference between the two pairing methods (consecutive dates for the first one, and with 2012 as reference for the second one). The differences in the ‘mean 300 m neighbours’ are slightly smaller than those in the ‘bilinear interpolation of the elevation difference’. Note that the increase in the differences after 2017 already mentioned in section 4.1 is also visible here. Without spatial interpolation, the differences are similar or smaller, except for one year (2017) with a much larger difference than with the two spatial interpolation methods (Table 4).

Table 4. Absolute differences between the cumulative mass balances of the two geodetic pairing methods, i.e. multi-year and the cumulative annual mass balances, in m w.e

The mean of the absolute difference between the annual mass balances of this interpolation method and the bilinear interpolation of elevation difference is 0.02 m w.e., which is much smaller than the calculated uncertainties. For the four winter mass balances, this indicator is 0.06 m w.e., also smaller than the uncertainties. The mean of the absolute difference between the glaciological and geodetic annual mass balances quantified with or without the ‘degree-day adjustment’ are listed in Table 5. All the indicators point to a small improvement with the ‘mean 300 m neighbours’ interpolation method compared to the ‘bilinear interpolation of elevation difference’. Without interpolation, the results may be more sensitive to data gaps. Indeed, an unusually strong negative mass balance can be noted without interpolation for the long-term 13 October 2012 to the 25 October 2017 mass balance (the 2017 DEM contains 31% of data gaps), mitigated by the other interpolation methods.

Table 5. Mean of the absolute difference between glaciological and geodetic annual mass balances according to the different interpolation methods and choices of date adjustment, in m w.e.

4.3.2. Density assumptions

To assess the influence of the uncertainty related to density on the mass balances, we calculate the main results with the upper and lower bounds of the uncertainty ranges on the density values. We use the mass balance quantified from the ‘bilinear interpolation of the elevation difference’ and compare raw values (with no adjustment for date). The differences in the mass balances quantified considering the original and extreme density values range from 0.02 to 0.39 m w.e., with a mean difference of 0.07 and 0.34 m w.e. for annual and winter mass balances, respectively. Compared to the original mass-balance values, this represents an average 10 and 22% change in the annual and winter mass balances, respectively. While the mean of the absolute difference between the non-temporally adjusted geodetic mass balances and glaciological results is 0.25 m w.e. for annual mass balances and 0.36 m w.e. for the three first winter mass balances, they are respectively 0.26 and 0.61 m w.e. for the upper density bound, and 0.78 and 0.39 m w.e. for the lower density bound. The density values we use produce the best results, and the lowest density bounds produce the worst results.

5. Discussion

The goal of this study is to evaluate whether time series of Pléiades DEMs can be used to accurately measure the glacier-wide mass balance at annual and seasonal scales. We discuss three main challenges: DEM precision, the volume to mass conversion (i.e. the density assumption) and the time gap between the acquisition of the satellite images and the field measurements as we have no control over the dates the Pléiades data are acquired.

5.1. DEM precision

Among all sources of uncertainty of the methodology used, DEM quality has the most important impact on the overall uncertainty of mass-balance quantification. As the seasonal and annual elevation changes are typically around the order of magnitude of 1–5 m, DEMs with high precision are required. This is particularly challenging for calculating winter mass balances, as the DEMs are acquired in late spring when a lot of snow remains on stable areas outside the glacier used for 3-D co-registration. The snowpack on the off-glacier terrain complicates the vertical co-registration and may cause a vertical shift of up to several tens of centimetres (+0.86 m estimated from the 13 May 2019 DEM). In addition, the automatic classification of the spring snow-covered areas is difficult in areas under shadow, and the remaining snow-free areas outside the glacier might be too small to allow accurate vertical adjustment between the DEMs. The stable surface areas used in this work are listed in Table 1, with one critical date (22 March 2019) when snow covered most of the valley floor. Consequently, spring DEMs are more difficult to co-register than end-of-summer DEMs, with higher std dev. over stable area: over 1 m for the spring DEMs versus a mean of 0.57 m for the end-of-summer DEMs. The std dev. are, outliers excepted, sub-metric to metric, and on-glacier data gap proportions are ~10%, which is of the same order of magnitude as the results of Berthier and others (2014) and Belart and others (2017) for Pléiades and WorldView2 DEMs.

In addition, some parts of the glacier are more prone to data gaps (mainly areas under shadow; Fig. 2 left panel), underlining the need for interpolation methods that are not much sensitive to this source of uncertainty. However, all interpolation methods are not equal and may result in marked surface elevation differences variations. The two interpolation methods compared here produce similar results, in agreement with McNabb and others (2019), while other methods (tested but not shown here) are not so consistent. Errors of interpolation could be the reason for a constant discrepancy of ~0.2 m w.e. between the long-term and annual time series of geodetic mass balances, which appear on 25 October 2017. The DEM at this date contains major gaps (31% of the surface area of the glacier) due to the low-angled sun at this period of the year, leading to a high proportion of shadow. The efficiency of the two interpolation methods compared here is proven by the mitigation of this discrepancy: the difference is of 0.26 and 0.29 versus 0.72 m w.e. without interpolation (Table 4). This difference is mainly caused by the multi-year mass balance, and it does not impact the rest of the time series. This highlights the sensitivity of the mass-balance time series to large data gaps: a single large gap error can influence the rest of the time series. While not assessed here, the image processing parameters and the software that produces DEMs can change the quality of the results, including gaps. DEMs derived from satellite images may show tilt and jitter distortions that we correct following the approach used by Deschamps-Berger and others (2020). The tilt and jitter combined corrections on DEMs can here reach 1 m. The selection of stable terrain close to the glacier and all around it can limit the impact of these phenomena. However, our stable terrain is mostly at the west of the glacier (Fig. 2, on the right). The absolute change of geodetic mass balance ranges from 0.02 to 0.42 m w.e. (mean 0.15 m w.e.). For Pléiades images and with our registration method, these characteristics are often not negligible and require special attention.

Apart from the above-mentioned limitations, the precision of Pléiades derived DEMs permits the assessment of annual-to-seasonal elevation changes. Note that here, the claim is made for a ‘large’ glacier (12 km²) but might not be true for small glaciers (typically <0.1–0.5 km²) for which the approach should be also tested and validated.

5.2. Volume to mass conversion

The choice of the volume to mass conversion factor has large impacts on the quantification of the geodetic mass balances. It is also the least constrained source of uncertainty, as it is a sort of average of multiple processes. For the annual mass balances, the share of the different surface types (ice, snow and firn) needs to be taken into consideration due to the marked difference in density and the year-to-year variation in the proportion of these different surfaces. In our case, the end-of-summer images used to generate the DEMs had a mean snow-covered area on the glacier of 40% (range from 26 to 52%). Hereafter, we compare two methods to apply the volume-to-mass conversion with the share principle. For the first method, we use a bulk mean density value that we multiply with the spatially averaged dh, corresponding to the method used in this study. For the second method, we multiply the corresponding density value (called ‘local density’) to each surface-type area by the corresponding dh. Note that for this method the mass conservation principle is the least respected due to ice flow and vertical ice velocity. The mean difference in the cumulative mass balances with the long-term pairing method and the ‘local density’ is −0.19 versus 0.22 m w.e. with the method used here. The use of the ‘local density’ results in more negative mass balances than with our method. The mean difference obtained with the glaciological annual mass balances is −0.08 m w.e. (0.21 m w.e. for absolute difference) with the ‘local density’, versus 0.02 m w.e. (0.25 m w.e.) with the ‘default’ density value. Dividing the annual mass balances calculated with the ‘local density’ by the mean elevation change, we find annual volume to mass conversion factors ranging from 100 to 1350 kg m⁻³ (median 833 kg m⁻³, mean 792 kg m⁻³). These values are closer to the multi-year density (850 ± 60 kg m⁻³) though much more variable, compared to the method we use (mean densities from 721 to 783 kg m⁻³, mean 756 and median 774 kg m⁻³).

In addition, a suitable date to estimate the proportion of different types of surfaces needs to be chosen. Due to light snowfalls that frequently occur at high elevations in late summer, snow covered more than half of the glacier on 20 September 2013 and 28 September 2016, and it covered the entire surface of the glacier on 25 October 2017. Hence, the snow-free glacier surface visible on these dates is not representative of the end-of-summer surface state. The maximum discrepancy in the distribution of the different types of surfaces appears for the year 2017 between the Landsat 8 and the Pléiades image. While the glacier was entirely snow-covered on the Pléiades images, only 37% of the glacier surface area was covered with snow on the Landsat image acquired on 19 August 2017. With a glacier that is entirely snow-covered (with related density value), the annual geodetic mass balance of 2016–17 is −1.13 m w.e., versus −1.60 m w.e. when considering the different types of surfaces from the Landsat delineation (a difference of 0.47 m w.e). For the above-mentioned dates, satellite images show almost no snow on the ground around the lower part of the glacier (excluding steep faces at higher altitude). Therefore, there is strong evidence that the snowfalls that occurred a few days before the image acquisitions were insignificant compared to the annual glacier surface elevation variation. Thus, Landsat surface share distribution is preferred, and the height and density change are considered negligible. No snowfall adjustment is implemented due to the hardly estimable amount of precipitation at the glacier surface and unknown spatial distribution of snow.

Density values and associated uncertainties may be difficult to retrieve. Density can differ between climatic regions and previous studies recommended using in situ measured densities when available (Pelto and others, 2019). Density also varies with the season; in this study, we only use in situ density values to quantify the winter mass balance. Snow compaction in spring and summer leads to strong spatial (mainly altitudinal) variations in density, so no geodetic summer mass balance is calculated. The mass-balance change resulting from the choice of these density assumptions within our rough uncertainty range can reach 10 and 22% for the annual and seasonal mass balances, respectively.

5.3. Temporal discrepancy

Due to differences in the acquisition dates, geodetic and glaciological mass balances are calculated for different periods. As Pléiades acquisitions are ordered within a window of several weeks and with a cloud-cover threshold, discrepancies between image acquisition dates and field measurements can represent several tens of days and several tens of centimetres of melt and/or solid precipitation. This is a major limitation to validation at short time scales. Degree-day adjustment to fit the dates of glaciological surveys is one possible solution which does require only limited measurements and that can be widely implemented, but it produced poor results in the present study with similar or higher discrepancies than without date adjustment (Table 5, Fig. 6). Degree-day adjustments, with melt estimation, can reach more than 1 m (Table 2). Unlike the adjustments of annual mass balances, degree-days adjustments for winter mass balances are mostly below the uncertainties (Supplementary Fig. S2). Other strategies not applicable to our data have been used in the literature: for example, Klug and others (2018) used different cases such as stake trend extrapolation and snow depth on non-glacierized terrain with better resolved data (laser scanning) and for shorter date differences hence resulting in limited adjustments, and Pelto and others (2019) used GNSS surveys.

Supplementary Figure S3 shows the discrepancies between glaciological and geodetic mass balances against (i) the number of days not considered in common by the two mass balances, and (ii) the cumulative precipitation fallen in the valley for this same period. Regarding the non-adjusted geodetic mass balances, discrepancies are weakly linked to the number of days, both for annual and winter mass balances (dots in Fig. S3-a), but there is no obvious link with cumulative precipitation, especially for winter mass balances (Fig. S3-b). Note that we do not consider the phase of the precipitation. Figure S1 shows that valley precipitation may occur in-between the acquisition, but more often at temperatures above 0°C at the glacier AWS (2400 m a.s.l.), which is still too uncertain to determine the precipitation phase. For long periods, such as for the 2012–2013 mass balance, some of the precipitation events at low temperature are followed by several days with temperatures above 5°C. Hence, a part of these precipitation may have no significant impact on the surface mass balance, either because the phase is liquid or by melting shortly afterwards. Still, for a hypothetical constant factor of two between precipitation at Chamonix and at the Glacier d'Argentière AWS (Vincent, 2002; Six and others, 2014), accounting for all precipitation would help reconcile six out of nine mass balances. For mass balances where in-between precipitation adjustment is >0.10 m w.e., three out of four annual mass balances are positively affected: 2013, 2016 and 2018. For these three annual mass balances, which also are those with the largest difference from the glaciological mass balances, the improvement covers 60% (i.e. 0.41 m w.e.), 33% (0.19 m w.e.) and 25% (0.20 m w.e.) of the difference, respectively. However, due to the strong spatial variability of snow precipitation in mountainous areas, point precipitation data recorded at an automatic weather station in the valley can hardly be spatialized at the glacier-wide scale, and the estimated precipitation factor has only been used for uncertainties and this sensitivity analysis.

5.4. Geodetic and glaciological measurements

Beyond glaciological calibration and local versus regional coverage, several complementarities between in situ and satellite measurements should be emphasized. As detailed previously, the DEM differences revealed both large and fine-scale spatial patterns. The difference in surface elevation can vary from year to year over the entire glacier, with a changing relationship with respect to hypsometry, precipitation, avalanches, etc. (Figs 3, 4). In addition, fine-scale impacts of serac falls or debris cover may also be explored. Unravelling the link between ice dynamics and surface mass-balance variability is well beyond the scope of this study, but the recent progress both on the observations and theory opens the way towards monitoring surface mass balance from elevation changes (Van Trichet and others, 2021; Vincent and others, 2021). Snow depth mapping with Pléiades images was positively assessed along with features of avalanche deposits by Deschamps-Berger and others (2020). Another possible application is the study of the seasonal variations in glacier thickness, which is still poorly known.

6. Conclusions

We have assessed the ability of high-resolution DEMs to quantify annual and seasonal mass balance. We compared annual and seasonal geodetic mass balances estimated from a time series of 12 DEMs produced from stereoscopic sub-metre Pléiades satellite images covering the Glacier d'Argentière over a period of seven hydrological years (2012–19) with mass balances based on two methods (i) geodetic multi-year mass balances using a reference density, and (ii) calibrated glaciological measurements. DEM quality is the main source of uncertainty; we found that the mean vertical precision of the DEMs is ±0.6 m (median ± 0.6 m) for annual mass balances and ±2.2 m (median ± 1.6 m) for spring DEMs. The mean of absolute discrepancy with glaciological mass balances geodetically calibrated is 0.25 m w.e. (maximum 0.48 m w.e.) and 0.36 m w.e. (maximum 0.63 m w.e.) for annual and winter mass balance, respectively. We note that the two methods are not fully independent due to the geodetic calibration of the glaciological mass balances with our data over the 7 years, and absolute discrepancies should be considered carefully and in a relative way.

Retrieved annual and winter geodetic mass balances were in good agreement with the glaciological and geodetic multi-year mass-balances estimates. The trends captured by the geodetic method were consistent, and individual mass balances could be retrieved with limited differences from the glaciological ones (after geodetic calibration). There are two limitations of this method: (i) the satellite and ground data need to be acquired as close in time as possible, the comparison becomes less reliable as the time difference lengthens. Degree-day adjustments did not reduce these discrepancies, as they depend not only on temperature, but also on the number of days between the two acquisitions and the precipitation amount that might have fallen in between; (ii) it is important to delineate snow/ice-covered surface areas and the choice of the density values has a high impact, especially between a bulk volume-to-mass conversion factor (i.e. the density assumption) or a ‘local’ application of the surface density.

Due to the different methods resolutions, geodetic and glaciological mass balances do not represent the same thing. Geodetic mass balances include process-scale phenomena that are not captured by glaciological mass balances. Analysis of geodetic time series shows that it is possible to describe several processes that are difficult to observe using other methods at lower spatial resolution (e.g. topographically influenced snow accumulation patterns), which can be observed at a wide scale with Pléiades-derived DEMs. One primary advance is the first estimation of winter geodetic mass balances of this glacier, which can be used for hydrological modelling and other applications. Geodetic mass balances are of great value to constrain models and calibrate in situ time series of glaciological mass balances.

Our results confirm the ability of high-resolution optical stereo-images for the quantification of both annual and seasonal geodetic mass balances. The DEMs derived from Pléiades images are generally of good quality after careful 3-D co-registration. Indeed, for the end-of-ablation season DEMs, a median std dev. of 0.55 m is found on stable terrain compared with a reference DEM and an on-glacier data gap median of 6%. Our approach can largely be automatized to produce similar results on a regional scale.