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 (Reference Vincent2002). 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, Reference Vincent2018). 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, Reference Vincent, Soruco, Six and Le Meur2009). 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, Reference Cuffey and Paterson2010). 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 (Reference Pelto, Menounos and Marshall2019), 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 (Reference Pelto, Menounos and Marshall2019), 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, Reference Huss2013). 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, Reference Paul2020; 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, Reference Fischer, Huss and Hoelzle2015). 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, Reference Gardent, Rabatel, Dedieu and Deline2014; Paul and others, Reference Paul2020; 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, Reference Berthier2014). 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, Reference Moratto, Broxton, Beyer, Lundy and Husmann2010; Shean and others, Reference Shean2016). 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, Reference Nuth and Kääb2011), using an open-source python-based workflow (Shean and others, Reference Shean2016; 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, Reference Nuth and Kääb2011). 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, Reference Girod, Nuth, Kääb, McNabb and Galland2017). This effect has been reported for several optical stereo satellites, including Pléiades (Deschamps-Berger and others, Reference Deschamps-Berger2020). 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, Reference Deschamps-Berger2020). 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 (Reference Berthier2014). 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, Reference McNabb, Nuth, Kääb and Girod2019). 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 (Reference McNabb, Nuth, Kääb and Girod2019). 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, Reference Six and Vincent2014), $T_i^ +$ positive air temperature during a period of n time intervals Δt (Hock, Reference Hock2003). 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, Reference Six and Vincent2014). 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 (Reference Six and Vincent2014) 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 (Reference Pelto, Menounos and Marshall2019), 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, Reference Rolstad, Haug and Denby2009; Menounos and others, Reference Menounos2019). 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 (Reference Fischer, Huss and Hoelzle2015):

(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, Reference Mälicke and Schneider2019) 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, Reference Huss2013).

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, Reference Berthier2014), 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, Reference Vincent2002; Six and others, Reference Six and Vincent2014). 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 (Reference Thibert, Blanc, Vincent and Eckert2008) 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 (Reference Wagnon2020), 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, Reference Vincent2018), σ _ɛ 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, Reference Vincent2018). 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.