A method to estimate surface mass-balance in glacier accumulation areas based on digital elevation models and submergence velocities

Abstract

Measuring surface mass-balance in the accumulation areas of glaciers is challenging because of the high spatial variability of snow accumulation and the difficulty of conducting annual field glaciological measurements. Here, we propose a method that can solve both these problems for many locations. Ground-penetrating radar measurements and firn cores extracted from a site in the French Alps were first used to reconstruct the topography of a buried end-of-summer snow horizon from a past year. Using these data and surface elevation observations from LiDAR and Global Navigation Satellite System instruments, we calculated the submergence velocities over the period between the buried horizon and more recent surface elevation observations. The differences between the changes in surface elevation and the submergence velocities were then used to calculate the annual surface mass-balances with an accuracy of ±0.34 m w.e. Assuming that the submergence velocities remain stable over several years, the surface mass-balance can be reconstructed for subsequent years from the differences in surface elevation alone. As opposed to the glaciological method that requires substantial fieldwork year after year to provide only point observations, this method, once submergence velocities have been calculated, requires only remote-sensing data to provide spatially distributed annual mass-balances in accumulation areas.

1. Introduction

Glacier surface mass-balance observations are widely used to assess climate change in various climate regimes because of the sensitivity of the surface mass-balance to climate forcings (e.g. IPCC, 2019; Zemp and others, 2019). Point surface mass-balances are primarily driven by meteorological variables such as precipitation and temperature (as opposed to the glacier-wide mass-balance that also depends on the glacier dynamics). Topographic-related factors such as shading, wind redistribution or accumulation due to avalanches may also influence the spatial pattern of accumulation (Machguth and others, 2006; Thibert and others, 2013; Sold and others, 2016; Vincent and others, 2017, 2018; Pulwicki and others, 2018). However, point surface mass-balance measurements require time-consuming fieldwork to collect data from stakes and pits across glaciers. For this reason, such measurements are conducted on only a few of the 200 000 mountain glaciers worldwide (World Glacier Monitoring Service (WGMS), 2017). Zemp and others (2015) report that only 37 glaciers worldwide are labelled as ‘reference’ glaciers by the WGMS, with ongoing mass-balance series covering more than 30 observation years.

Glacier mass-balance modelling mainly relies on melt models using detailed ablation measurements as a solid basis for parametrisation (e.g. Radić and Hock, 2006). Because of the high spatial variability, accumulation processes are more complex to take into account in the models and have attracted less attention (Machguth and others, 2006; Sold and others, 2013). The spatial variability of snow accumulation on glaciers remains largely uncaptured by the traditional snow pit and probing surveys (Pulwicki and others, 2018). This is problematic because winter accumulation has a large impact on the annual mass-balance, including via albedo feedback (Réveillet and others, 2017; Marshall and Miller, 2020). Better and spatially distributed estimates of snow accumulation will also help to constrain high-altitude precipitation models (Immerzeel and others, 2015). Finding alternative or complementary methods to estimate glacier surface mass-balance, particularly in the upper reaches of glaciers, is therefore highly relevant.

Remote-sensing elevation measurements (e.g. from airborne LiDAR, optical and radar satellites or photogrammetric UAV) are promising and have the advantage of covering wide areas. However, to use them to determine surface mass-balances, they must be complemented by field measurements or modelling to account for ice dynamics (e.g. Sold and others, 2013; Belart and others, 2017; Réveillet and others, 2020; Pelto and Menounos, 2021). Indeed, glacier flow leads to the emergence of ice in the ablation zone with a corresponding increase in the surface elevation height (Vincent and others, 2021) and a downward displacement of ice in the accumulation area with a corresponding decrease in the surface elevation height. Note that the total decrease in the surface elevation height, referred to as submergence, also includes the effect of firn compaction (Cuffey and Paterson, 2010).

Ground-penetrating radar (GPR) has been widely used in both polar and mountain environments to measure spatially-distributed snow accumulation, through time consuming manual identification of the buried isochronal internal reflection horizons (IRHs) (e.g. Eisen and others, 2004; Machguth and others, 2006; Sold and others, 2013, 2015; McGrath and others, 2015; LeMeur and others, 2018). These IRHs identified by GPR can also be used in accumulation areas to calculate submergence velocities, thus providing the measurements needed to complement the remote sensing data for mass-balance estimation.

Here, we test a method to obtain spatial estimates of surface mass-balance using digital elevation models (DEMs) that can be acquired regularly by remote-sensing techniques, together with initial estimates of submergence velocities. To do this, the spatial distribution of submergence velocities is first determined from data collected during an extended field campaign involving GPR and firn-core measurements. The surface mass-balance can subsequently be calculated annually by acquiring new DEMs alone, without additional glaciological ground measurements. In a recent study on Wolverine Glacier, Alaska, Zeller and others (2023) used a similar approach, using submergence velocities and DEM differences, to estimate spatially distributed seasonal and annual surface mass balances. Among three methods to estimate submergence velocities, they also used GPR measurements, but in an indirect approach different that the one described here. In their study, the submergence velocities were calculated using winter mass balance deduced from GPR measurements and DEM differences. Here, submergence velocities are calculated from direct elevation measurement, using GPR, of a buried past horizon.

We apply our method to the accumulation area of Mer de Glace, at the Col du Midi pass (Mont Blanc area, France) to estimate 2015–18 and 2012–21 surface mass-balances. After presenting the study area (Section 2), we detail the protocol used to acquire the in situ snow measurements, surface elevation measurements and GPR surveys (Section 3). In Section 4, we detail the method of estimating surface mass-balance from submergence velocities. In Section 5, we apply our method to the 2015–18 and 2012–21 periods. We discuss our results in Section 6.

2. Study site

The study site is located in the upper part of the accumulation area of the Mer de Glace glacier (Massif du Mont-Blanc, French Alps), at the Col du Midi pass (45°52′8 N, 6°53′12 E) between 3450 and 3530 m a.s.l. (Fig. 1). The western part of the pass descends slightly above an icefall. The eastern part, on which our study was carried out, is a vast, relatively flat plateau with an average slope of 2° and a maximum slope of 11°. The average ice thickness is ~180 m (Vallon and others, 1976) and the surface horizontal velocities are of the order of 10 m a⁻¹ (Réveillet and others, 2020). The study area covers approximately 0.6 km².

Figure 1. Location and topography of the Col du Midi study site in the accumulation area of the Mer de Glace glacier (Pléiades © CNES 2021 (year of acquisition), Distribution AIRBUS DS). The contour line (grey lines) interval is 20 m. The red pentagons are the 2019 drilling sites, the grey circles the GPR tracks and positions of GNSS measurements, the black asterisk the position of the fixed GNSS receiver. ACCU1 and ACCU2 are the regular accumulation measurement sites of the GLACIOCLIM Observatory. Additional elements from a previous study conducted at the same site (Réveillet and others, 2020) are also reported: drilling sites are represented by black triangles. Sites S1 to S5 correspond to locations where stake measurements were made in order to calculate submergence velocities. Red arrows represent annual surface flow velocities. The dashed black line is the outline of the zone scanned by the LiDAR. The orange dashed line represents the virtual GPR transect between the two boreholes at sites F2 and S2, schematised in Figure 4.

In situ measurements of winter snow accumulation conducted by the French Service National d'Observation GLACIOCLIM (https://glacioclim.osug.fr) indicate a mean value of 2.3 m w.e. per winter (period 1995–2019) on the ACCU2 site at the centre of the study area (Fig. 1). Strong interannual variability (Std dev. (SD) of ± 0.5 m w.e.) is observed, but with no trend over the 25 years period. The mean annual point surface mass-balance at that site is also 2.3 m w.e. (SD of ± 0.7 m w.e.), but was quantified for only 12 of the 24 years because it was often impossible to find the accumulation stake in the late summer. The identical values of the reported mean winter snow accumulation and mean annual point surface mass-balance is simply a coincidence, due to the fact that the annual point balances are reported for only half the years over the period. It does not mean that there is no ablation at this site. The mean winter snow accumulation for the corresponding 12 years is 2.5 ± 0.5 m w.e.

Réveillet and others (2020) discussed the spatio-temporal evolution of accumulation at the Col du Midi pass based on repeated acquisitions of DEMs from terrestrial LiDAR measurements over the 2014–15 period. They used a sparse network of point measurements of submergence velocities, determined from stakes or firn core measurements, to quantify the accumulation based on the differences in the DEM data. They found a large spatial variability in submergence velocity in an area smaller than 2 km² (mean of 4.5 m a⁻¹, SD of 1.5 m a⁻¹ and a maximum difference of up to 5.9 m a⁻¹).

3. Data

The main field measurement campaign was carried out on 6 February 2019, including GPR, GNSS and firn core measurements. A LiDAR DEM was acquired previously in October 2015.

3.1. Radar measurements

GPR measurements were collected using a MALÅ ProEx control unit equipped with a 250 MHz shielded antenna using a trigger interval of 1 s, a 128-fold stacking, and a time window of 250 ns. This time window corresponds to a sounding of up to 26.9 m at a wave propagation speed of 0.215 m ns⁻¹ (Section 5.2.). The GPR system was mounted on a sledge and pulled along a 5.1 km route located within the study site, at a speed <1 m s⁻¹ (Fig. 1). The GPR acquisition was continuously interfaced with a GB1000 Top Con standalone GNSS with a horizontal and vertical accuracy of ~2–3 m depending on the quality of the satellite constellation.

Radargrams were processed using the ReflexW-2D software package (Sandmeier Scientific Software). All radargrams were corrected to time zero, taken as the first break in the first wavelet (Yelf and Yelf, 2006). Continuous and contrasting IRHs were manually picked along the entire 5.1 km profile. The background was removed when necessary, especially where the IRH was close to the surface and over-shadowed by the direct wave. Energy decay, bandpass filter and deconvolution techniques were also applied to improve the definition of deep IRHs.

In regions where the signal is noisy due to numerous reflectors, either annual horizons or ice layers (refrozen melt or rain water), IRHs are difficult to interpret from echograms. We solve this by using intersecting radar tracks at sites S2, S3, S5 and F2 which ensured consistency in identifying the annual horizons. The speed of wave propagation through the firn layer, and subsequent time-to-depth constant conversion factor, was determined using the drilled firn cores.

The accuracy of the GPR is contingent upon the resolution of the GPR signal (Hubbard and Glasser, 2005). The resolution is theoretically estimated to be about a quarter of the wavelength, λ/4; we instead assume λ/2 as it is considered to be more realistic and includes manual picking uncertainties (Le Meur and others, 2018). In this study, given a frequency of 250 MHz and a velocity of 0.215 m ns⁻¹, the uncertainty of λ/2 corresponds to ±0.43 m. This value represents 2 ns of the two-way travel time temporal scale (see left axis in Fig. 6). We discuss GPR uncertainty (including the uncertainty of the wave velocity) in Section 5.2.

3.2. GNSS measurements

For accurate positioning of GPR measurements, particularly of the vertical (z) component, discrete differential GNSS measurements were collected using a fixed and a mobile receiver. The fixed receiver, a Leica 1200 GNSS unit, was installed on a nearby reference point secured on bedrock <1 km from GPR's locations (black asterisk in Fig. 1). Every ~40 m, the GPR's antenna's location was acquired using a Leica 1200 GNSS unit (the mobile receiver) with dual frequencies comprising a total of 146 points measurements along the radar profiles. Occupation times were typically 1 min with 1 s sampling and the number of visible satellites (GPS and GLONAS) was always >7. Leica Infinity software was used for data post-processing. The differential GNSS positions have intrinsic vertical and horizontal accuracies of ±0.01 m. However, due to the vertical displacement caused by the weight of the radar equipment on fresh snow, an additional ±0.1 m is included in the vertical uncertainty estimates. In addition to radar measurements, these surface point XYZ values were used to generate the 2019 surface DEM. The position of ground control points needed to LiDAR measurements have been also obtained from GNSS measurements.

3.3. LiDAR measurements

We used an Optech ILRIS-LR long-range Terrestrial Laser Scanner with a near-infrared wavelength of 1064 nm suited to a snow-covered surface over scanned distances of 3 km. At the study site, the scanned distances range from 400 to 1800 m. The beam size of the scanner is 27 ± 7 mm at 100 m (factory data) and increases with distance (from 108 ± 28 mm at 400 m to 486 ± 126 mm at 1800 m). Measurements were made from the southernmost Aiguille du Midi tourist terrace (Fig. 1). Data was acquired with an averaged sampling step of 0.20 m at 900 m (ranging from 0.10 to 0.42 m depending on the distance). The LiDAR point cloud data were subsequently processed with Polyworks software to generate a DEM. Six ground control points were used for the georeferencing of the point cloud. The 1 m resolution DEM was obtained with a horizontal and vertical uncertainty related to the alignment and georeferencing of 3 and 10 cm, respectively. More details on the LiDAR acquisition process can be found in Réveillet and others (2020).

3.4. Satellite measurements

A map of the surface elevation changes between 19 August 2012 and 15 August 2021 was calculated by subtracting two DEMs derived from Pléiades stereo images. The 4-m DEMs were generated using the Ames Stereo Pipeline (Beyer and others, 2018). The DEMs were coregistered on stable terrain following Berthier and others (2007), masking out glacierised areas using a glacier inventory from Paul and others (2020). Spatially-coherent elevation differences were next corrected using a polynomial fit across-track (e.g. Gardelle and others, 2013) and a spline fit in the along-track direction (Falaschi and others, 2023). Uncertainty in the rate of elevation change at a single pixel is assumed to equal the standard deviation of the rate of elevation difference on the stable terrain for the same range of slope (i.e. 0 to 5°), so 0.12 m a⁻¹.

3.5. Firn cores

Six firn cores were drilled within the area of GPR acquisition. Four of these cores were drilled to a depth of 3–4.5 m using a hand drill. They were used to determine locally snow accumulation, identified on the basis of grain size and dust observations, between the GPR campaign (6 February 2019) and the previous end-of-summer 2018 surface. These cores, S1, S3, S4, and S5 (Fig. 1), match the locations of previously installed stakes (Réveillet and others, 2020). The accuracy of such direct measurements of snow accumulation obtained from the cores is estimated to be ± 0.1 m (Thibert and others, 2008).

Two deeper cores were drilled to identify the depth of the end-of-summer 2015 horizon, using the Fast Electromechanical Lightweight Ice Coring System (FELICS) with a 58 mm core diameter (Ginot and others, 2002). The first core was drilled down to a depth of 21 m at site S2 (Fig. 1), corresponding to the site of stake 2 in Réveillet and others (2020). The second core was drilled upstream at site F2, down to a depth of 8 m. Snow accumulation is much lower at site F2 than at site S2 (Réveillet and others, 2020). Because of the low-density surface snow, sampling with the electromechanical drill began at 0.4 m depth. The depth accuracy of the S2 and F2 cores is therefore estimated to be ± 0.2 m.

Both the S2 and F2 cores were stored frozen in core boxes and brought back to the Institut des Géosciences de l'Environnement (IGE, Grenoble, France) cold rooms. The density along the S2 and F2 cores was measured in a cold room (−15 °C) and stratigraphic parameters (icy and dirty layers) were logged. The density was calculated on the raw cores, without subsampling, by measuring the length, circumference and weight of fragments (lengths from 5 to 83.5 cm). Next, the cores were analysed for major ionic composition using Ion Chromatography. Unfortunately, the ionic chemical analyses do not show clear seasonal peaks and thus do not help to separate annual layer and identify summer horizons. Therefore, the chemical profiles shown in supplementary information B are not used in our analysis (Figs S4 and Fig. S5).

3.6. Point surface mass-balance measurements

Point surface mass-balance measurements have been conducted each year since 1995 at site S2 by the French Service National d'Observation GLACIOCLIM using the classic glaciological method (Cuffey and Paterson, 2010, page 128; Six and Vincent, 2014). At the end of winter accumulation season, snow accumulation is measured using a hand drill, as described in Section 3.5 for the cores at points S1, S3, S4 and S5. The previous summer's snow horizon is determined based on grain size and/or snow colour observations. In order to determine winter mass-balance in m w.e., density is also quantified in the field by measuring the length, circumference and weight of each corresponding section of core extracted. A 3 cm diameter wooden stake, made of 2-m long poles linked together by a small chain, is then left in the drill hole. The section of the stake remaining above-surface is measured after the melt season in September. Assuming a firn density of 550 ± 30 kg m⁻³ at the end of the summer period, based on previous campaigns, this calculation produces each year the annual point surface mass-balance with an uncertainty of ± 0.21 m w.e. (Thibert and others, 2008; Six and Vincent, 2014).

3.7. Data interpolation and uncertainties

Although the manually picked IRH is continuous along the entire profile of the radargrams, in the calculations using GPR data we only take into account the locations where surface differential GNSS measurements were recorded (146 sites). This allows elevation variations at discrete locations to be analysed with a high level of accuracy. To spatially align all of the datasets to one identical grid, interpolation was implemented using a universal kriging method (linear variogram; slope = 1, anisotropy = 1.0; 10 m × 10 m resolution, xy). We use a 10 m gridcell size as this generates the most spatially complete raster with the highest resolution. All the uncertainties are calculated using the independent uncertainty propagation method (i.e. root mean square of the quadratic sum) and assuming independence of each term.

4. Method

4.1. General principles

The temporal change in the glacier surface elevation is obtained from Cuffey and Patterson (2010, p. 332):

(1)$$$\displaystyle{{\partial S} \over {\partial t}} = \displaystyle{{\dot{b_s}\;} \over {\rho _s}} + W_s-U_s\displaystyle{{\partial S} \over {\partial x}}-V_s\;\displaystyle{{\partial S} \over {\partial y}}\;$ $$

where S is the glacier surface elevation, U_s, V_s, W_s are the components of ice flow velocity at the surface in directions x, y and z, ∂S/∂x and ∂S/∂y are the surface elevation gradients in directions x and y, $\dot{b}$_s is the annual surface mass-balance (in m w.e. a⁻¹) and ρ_s is the density of snow and firn in the accumulation area.

The term $W_s-U_s{{\partial S} \over {\partial x}}-V_s\;{{\partial S} \over {\partial y}}\;$ is called the submergence velocity (V_sub). It integrates both the downward displacement of the ice due to its dynamics and the effects of the compaction of the firn. By rearranging Eqn (1), the surface mass-balance can therefore be calculated using the submergence velocities and changes in glacier surface elevation:

(2)$$\dot{b}_s = \left({\displaystyle{{\partial S} \over {\partial t}}-\;V_{sub}} \right)\rho _s\;$$

Figure 2 shows schematically the principle of the method used to calculate the spatial distribution of the surface mass-balance from Eqn (2). First, the elevations of the glacier surface at dates N and N + i are used to calculate the surface elevation change ∂S/∂t from the difference ΔS. Second, the submergence velocity V_sub is obtained from the position of the buried surface of the year N snow layer after i years (time interval Δt). This method is an Eulerian approach that calculates the submergence velocity for fixed points and produces a mean value for the considered layer.

Figure 2. Schematic representation of the key parameters for estimating the surface mass-balance by the method proposed in this paper. From top to bottom, the three lines correspond to the glacier surface in a reference year N (brown line), the surface in year N + i (orange line), and the horizon of the year N surface layer, moved downward under the effect of submergence, at the date N + i. The variables ∂S/∂t, Vsub and ${{\dot{b_s}\;} \over {\rho _s\;}}\Delta t$ are defined below (Table 1).

Table 1. Specific variables used in this study to estimate annual surface mass-balance between 2015 and 2018 (illustrated in Fig. 2 and Fig. 4)

The principles of our methodology are the following:

1. a) Acquire a high-resolution DEM (e.g. terrestrial or aerial LiDAR, drone survey, high-resolution satellite stereo images) at the end-of-summer season of year N. This is the reference surface. When covered by snow, it will become the end-of-summer horizon of year N.

2. b) Perform tight grid GPR surveys after i years to estimate submergence velocities. In order to determine the corresponding IRHs, drill firn cores to identify the buried end-of-summer horizon of year N. In addition, the firn cores provide density measurements.

Once steps a and b have been completed, this approach makes it possible to determine the spatial distribution of the surface mass-balance in a glacier accumulation area based on DEM differences only, assuming that the submergence velocities are stable over the period considered (this assumption will be discussed later).

4.2. Detailed methodology for this study (2015–18 period)

Our objective is to quantify the surface mass-balance over the three hydrological years 2015/16, 2016/17 and 2017/18 using the surface elevation changes over the 2015–18 period (Fig. 3). The problem is somewhat complicated due to the variability of the temporal scale of our datasets. Fieldwork took place in February 2019 and not in late summer 2018. Consequently, the surface elevation changes obtained are between October 2015 and 6 February 2019, and the submergence velocities were calculated from the elevation of the end-of-summer 2015 snow horizon measured in February 2019. For these reasons, some corrections were made to take into account: (i) the snow accumulation on the surface from the end-of-summer 2018 to February 2019 and (ii) the downward displacement of the end-of-summer 2018 and 2015 horizons over the same period due to submergence (Fig. 4). The variables used in Figure 4 are summarised in Table 1.

Figure 3. Schematic of the study's workflow. The yellow blocks correspond to field/remote sensing measurements for each year. The red and green blocks correspond to the intermediate calculated variables.

Figure 4. Schematic representation of a vertical section of the firn along a virtual GPR transect between the two boreholes at sites F2 and S2. Radargrams close to the drilling sites are represented with picked end-of-summer 2018 (orange) and 2015 (red) horizons, as well as the observed stratifications in the two cores (ice layer in grey, dirty layers in red). Arbitrary 6 February 2019 surface and buried past summer horizons are represented, as well as a schematic position of the October 2015 surface. The different variables calculated in this work, interpolated and plotted in Figure 5a to 5f, are represented by arrows (see Table 1 for a description of each variable).

Figure 5. 10 m × 10 m interpolation maps of: (a) Surface elevation change of the Col du Midi area between October 2015 and 6 February 2019 (dZ1). (b) Snow accumulation between the end-of-summer 2018 and 6 February 2019 (dH1). (c) Net snow accumulation between the surface and the end-of-summer 2015 horizon on 6 February 2019 (dH2). (d) Submergence velocities (dZ4/Δt). (e) Mean annual surface elevation change between the end-of-summer 2015 and the end-of-summer 2018 (∂S/∂t). (f) Mean annual surface mass-balance calculated between 2015 and 2018 (ḃ_s). The red pentagons show the snow or firn core locations and the grey dots the GNSS measurements on the GPR tracks. The contour line (grey lines) interval is 20 m.

From Eqn (2), the quantity ${{\dot{b_s}\;} \over {\rho _s\;}}\Delta t$ over the period Δt (2015–18) can be calculated from the difference in surface elevation ${{\partial S} \over {\partial t}}\Delta t$ and the submergence velocity V_sub, multiplied by the time interval Δt, assuming that the density ρ_s of the 2015–18 snow layer is known (Fig. 4). Given that the measurements were made on 6 February 2019, the measured surface elevation change dZ1 is affected by (i) the snow accumulation dH1 between the end-of-summer 2018 and 6 February 2019 and (ii) the downward motion dH3 of the end-of-summer 2018 horizon between the end-of-summer 2018 and 6 February 2019. Similarly, the estimated submergence velocities dZ4/Δt are affected by (i) the measured surface elevation change dZ1 between October 2015 and 6 February 2019 and (ii) the downward motion dH4 of the end-of-summer 2015 horizon between the end-of-summer 2018 and 6 February 2019.

Thus, we can write (see Fig. 4):

(3)$$\displaystyle{{\partial S} \over {\partial t}} = \displaystyle{{dZ1-dH1-dH3} \over {\Delta t}}$$

(4)$$V_{sub} = \displaystyle{{dZ1-dH2-dH4} \over {\Delta t}}$$

Then Eqn (2) becomes:

(5)$$\dot{b}_s = \left({\displaystyle{{dZ1-dH1-dH3} \over {\Delta t}}-\displaystyle{{dZ1-dH2-dH4} \over {\Delta t}}} \right)\;\rho _s$$

Thus, the annual surface mass-balance between 2015 and 2018 is calculated from:

1. (i) the surface elevation changes measured between the end-of-summer 2015 and 6 February 2019 (dZ1);

2. (ii) the depths of the end-of-summer 2018 and end-of-summer 2015 horizons (from dH1 and dH2 respectively) obtained from GPR measurements performed on 6 February 2019. This is used to obtain the submergence velocity.

3. (iii) dH3 and dH4, calculated from the calculated submergence velocity along with the time between end-of-summer 2018 and 6 February 2019.

4.3 Application to subsequent determination of SMB using DEMs differences only (period 2012–2021)

Finally, our method has been validated over the period 2012–2021. For this purpose, the average annual SMB have been reconstructed, based on Eqn (2), from the previous submergence velocities and the surface elevation changes between 19 August 2012 and 15 August 2021 calculated by subtracting two DEMs derived from Pléiades stereo images. These results have been compared directly to the in situ observations of SMB over the period 2012–2021.

5. Results

5.1. Surface elevation changes between 2015 and 2019 (dZ1)

We interpolated a DEM of the surface topography for the zone delimited by the radar tracks (Fig. 1) from the differential GNSS point measurements that were acquired on 6 February 2019. We compared the 2019 radar-DEM with the 2015 LiDAR-DEM (Fig. 5a). Over October 2015 to February 2019, there was a surface gain in the central part of the study area of 1–2 m. Other regions experienced an elevation loss by as much as 4 m. The accuracy of the calculated surface elevation change, determined from the ± 0.10 m accuracy of both the 2015 and 2019 surface DEMs is ± 0.14 m.

5.2. Snow accumulation between the end-of-summer 2018 horizon and 6 February 2019 (dH1)

The fresh snow accumulated between the end-of-summer 2018 and 6 February 2019, dH1, is estimated from GPR measurements. The end-of-summer 2018 horizon is well marked on the radargrams (orange line in Fig. 6 and Figs S1–S3) and cannot be confused with an IRH linked to ice layers given that none are observed above this end-of-summer horizon (see Fig. 7). The corresponding IRH was manually picked along the radar transects to determine the spatial variability of snow accumulation at the Col du Midi pass during the first part of the 2018/19 winter. We used the cores drilled in the field to determine the wave propagation speed in the upper part of the snowpack and obtain a time-to-depth conversion factor for the radar signal. For each of the six cores (S1–S5 and F2), we identified the depth of the previous end-of-summer 2018 horizon. These depths were compared to the end-of-summer 2018 IRH picked locally at the same sites, assuming different wave propagation speeds. A speed of 0.215 m ns⁻¹ led to the best match between these two measurements (Table 2) and was applied to all the GPR measurements. For the six coring sites, differences in horizon depths were always <0.21 m (with an average of 0.01 m), well within the uncertainty range of the measurements.

Figure 6. (a) Example of radargram obtained along the S2–S5 transect at 250 MHz. (b) Same radargram with manually picked end-of-summer 2018 (orange) and end-of-summer 2015 (red) horizons. The calculated uncertainties (10% for the red curve, 14% for the orange one) are also shown (thin red lines). Three other well marked IRHs are shown by light white lines. The blue squares on the top bar and on the insert map correspond to the points precisely measured by differential GNSS. (c) Enlarged details of the beginning of the radargram (blue rectangle) without and with picked IRHs and calculated uncertainties. Various radargram processing approaches are also displayed for the same part of the radargram: (d) Additional energy decay, background removal and deconvolution processes. (e) Additional energy decay, background removal, deconvolution and dewow processes. Examples of other radargrams are reported as supplementary information (Fig. S1, S2, S3).

Figure 7. Profiles of firn cores' densities from sites S2 and F2. Density measurements start at 0.40 m below the surface because the fresh snow at the surface is highly inconsistent. Icy layers are represented in grey and dirty layers in red.

Table 2. Comparison of the snow accumulated between the end-of-summer 2018 and 6 February 2019 measured directly in the field by drilling snow cores and determined from radar measurements for a wave propagation speed of 0.215 m ns⁻¹

The last column indicates the difference in snow accumulation, in metres of snow, between the radar and the core methods.

We estimated an uncertainty of ± 0.01 m ns⁻¹ for the wave propagation speed, taking into account the spatial variability of the firn density within the studied area. Combining this uncertainty with that of 2 ns defined for the radargram picking, the uncertainty is ± 0.57 m for the 2018 horizon depth (over an average 2018 horizon depth of 4 m). Given that the firn cores were used to identify the horizon depth, their ± 0.2 m uncertainty was incorporated with picked radar horizon uncertainty which resulted in a total uncertainty of ± 0.60 m at site S2.

Figure 5b presents the snow accumulation pattern obtained for the period between the end-of-summer 2018 and the date of the field campaign (6 February 2019). The snow accumulation pattern for the small (<1 km ²), flat study site has large spatial variability. Snow accumulation ranges from 2 m in the upstream area near the pass (site F2) and in the east steeper region where surface crevasses are present, to 6 m in the centre of the study site, a factor of three within a distance of ~400 m.

5.3. Net snow accumulation between the end-of-summer 2015 horizon and 6 February 2019 (dH2)

5.3.1. End-of-summer 2015 horizon depth at site S2

As the radar wave penetrates deeper through the firn, numerous reflectors from annual horizons or ice layers make an unambiguous identification of individual summer annual horizons almost impossible. We used physical measurements and observations (density, icy and dirty layers) on the 21 m firn core drilled at site S2 to identify the depth of the end-of-summer 2015 horizon (Fig. 7).

Unfortunately, these data alone proved insufficient. Consequently, we used the order of magnitude of the net annual accumulation reported by the GLACIOCLIM observatory for the hydrological years 2015/16, 2016/17 and 2017/18 to aid in identifying the end-of-summer 2015 horizon from several visible layers in the drilled core.

A few dirty layers were observed in the S2 core (Fig. 7): one at a depth of 9.34 m (1 cm thick), one between 15.28 and 15.34 m (6 cm thick), and two others at 18.63 and 21.46 m (both 1 cm thick). These dirty layers correspond to end-of-summer horizons in relation with impurities accumulated at the surface during summer months. However, as illustrated by the absence of a dirty layer for the end of-summer 2018 horizon (Fig. 7), end-of-summer horizons are not necessarily associated with a visible dirty layer (Sold and others, 2015). Hence, there was some doubt concerning the position of the end-of-summer 2015 horizon, making it necessary to conduct a sensitivity test to determine the correct depth of the 2015 horizon. We assumed that the 9.34 m dirty layer corresponds to the end-of-summer 2017 horizon. Then we successively assigned two depths (corresponding to 2 visible dirty horizons) to the 2015 end-of-summer horizon and compared their values with the order of magnitude of the 2015 horizon depth estimated by GLACIOCLIM mass-balance measurements (Table 3).

Table 3. In situ measurements at site S2 conducted since summer 2015: annual surface mass-balance measurements conducted by the GLACIOCLIM program

^a This study.

^b Estimation (see text).

^c data from the GLACIOCLIM observatory (available at https://glacioclim.osug.fr).

Bold values correspond to the estimated annual surface mass-balance for the hydrological years 2015/16, 2016/17 and 2017/18.

Table 3 reports the significant variability of the annual winter snow accumulation measured by the glaciological method at site S2, ranging from 4.78 to 8.29 m of snow between 2016 and 2018. The mean value for this period, expressed in water equivalent, is 2.69 m w.e., similar to the mean observed value over 1995–2019 of 2.3 ± 0.5 m w.e. The net annual accumulation determined in September 2016, 2017 and 2018 was 3.03, 1.05 m and 3.16 m w.e. respectively. On 6 February 2019, the end-of-summer 2015 horizon is consequently estimated at a depth equivalent to 8.79 m w.e. (1.55 + 3.16 + 1.05 + 3.03 m w.e., see Table 3) from the surface, or 7.24 m w.e. (3.16 + 1.05 + 3.03 m w.e.) below the end-of-summer 2018 horizon.

In the S2 firn core, we initially linked the dirty layer at 15.3 m to the end-of-summer 2016 horizon and the dirty layer at 18.6 m to the end-of-summer 2015 horizon (grey values in Table 4). However, the corresponding net accumulations calculated using the measured density profile differed from the GLACIOCLIM values. The surface mass-balance would be significantly over-estimated for the 2016/17 hydrological year (3.58 m w.e. compared to the 1.05 m w.e. of the GLACIOCLIM data) and underestimated for 2015/16. We therefore had to consider the possibility that the end-of-summer 2016 horizon was not linked to a dirty layer in the S2 firn core, as was the case for the 2018 horizon. In this hypothesis, the dirty layer at 15.3 m corresponds to the end-of-summer 2015 horizon (black values in Table 4) and the end-of-summer 2016 horizon is situated between 9.34 and 15.30 m. This depth was estimated to be 10.9 m (ice layer). As shown in Tables 3 and 4 (bold values), this hypothesis led to good agreement with the GLACIOCLIM values, even though the pluri-annual firn core method estimates an accumulation at site S2 (6.53 m w.e. between the end-of-summer 2018 and 2015 horizons, or 2.18 m w.e. a⁻¹) on the average 10% lower than the annual glaciological method.

Table 4. In situ measurements at site S2 conducted since summer 2015: accumulation deduced from the S2 core study

^a This study.

^b Estimation (see text).

Grey values in brackets refer to a sensitivity test on the position of the end-of-summer 2016 and 2015 horizons (see text). Bold values correspond to the estimated annual surface mass-balance for the hydrological years 2015/16, 2016/17 and 2017/18.

5.3.2. End-of-summer 2015 horizon depth over the whole studied area

The end-of-summer 2015 horizon for the study-site was identified from the GPR survey (Fig. 5c). Unlike the end-of-summer 2018 horizon, the end-of-summer 2015 horizon is visually less distinctive on the radargrams because of the attenuation of the radar signal with increasing snow depth (red horizon in Fig. 6 and supplementary material Fig. S1, S2, S3). In the firn core, this horizon was found at a depth of 15.3 ± 0.2 m at site S2. Using the mean wave propagation speed of 0.215 ± 0.010 m ns⁻¹ in the firn, a corresponding IRH was visible on the recorded radargrams at a mean depth of 15.75 ± 0.11 m (mean observation and SD of the four radargram tracks available at site S2). At this depth, the error propagation calculation leads to a total uncertainty of ± 1.47 m at site S2, or 10% of the depth. This IRH was assigned to the end-of-summer 2015 horizon. Manual picking of the IRH, associated with 10% uncertainty, was possible along almost all the radar tracks, except for the lower eastern zone (around points M9 and M15, see Fig. 1 and Fig. S3) where the radargrams were too noisy. The crossing of four radar tracks at site S2, where the end-of-summer 2015 horizon was initially identified, as well as at sites F2, S3 and S5, made it possible to validate the manual picking and to correct some drift caused by zones in a few radargrams where several IRHs could have been selected for the end-of-summer 2015. Figure 6 and Figs S1, S2 and S3 show examples of radargrams with the picked end-of-summer 2015 horizon, its uncertainty and different radar data processing approaches. Analyses conducted in the 8-m firn core of site F2 did not allow clear identification of annual horizons, but the depth of the picked 2015 IRH matched a summer dust layer observed at a depth of 4.1 m (Fig. 7).

5.4. Submergence velocities (dZ4/Δt)

The submergence velocities are calculated by differencing the elevations of the October 2015 surface DEM and the February 2019 end-of-summer 2015 horizon DEM and dividing by the time between the two timestamps (3.28 years; Fig. 5d).

The submergence velocities exhibit large spatial variability, with values ranging from −2 m a⁻¹ in the upper section of the Col du Midi pass to ~−5 m a⁻¹. Given an uncertainty of ± 1.47 m at site S2 for the radargram manual picking, ± 0.1 m for the vertical GNSS position and ± 0.1 m for the 2015 surface DEM, we calculated an uncertainty of ± 1.48 m for the elevation difference (−14.5 m) of the end-of-summer 2015 horizon between 2015 and 2019. This produced an uncertainty of ± 0.46 m a⁻¹ for the average submergence velocities or, ± 0.27 m w.e. a⁻¹ when using a mean firn density of 540 ± 20 kg m⁻³ (measured in the S2 core from the surface down to the 2015 horizon). This uncertainty quantified at site S2, where in situ density measurements were available, could be larger within the interpolated study area because of the variability of the snow density. Assuming an uncertainty of ± 1.2 m for the end-of-summer 2015 horizon depth (i.e. 10% of the 12 m mean value of end-of-summer 2015 horizon depth) and an uncertainty of ± 80 kg m⁻³ for density, we estimated a mean uncertainty of ± 0.38 m a⁻¹ or ± 0.34 m w.e. a⁻¹ for the spatially averaged submergence velocity. This value of ± 80 kg m⁻³ is the standard deviation from the mean density measured on the S2 core between the two 2015 and 2018 horizons; it is also the density difference observed between the S2 and F2 cores for the same 2015–2018 firn layer.

5.5. Surface elevation changes between 2015 and 2018 (∂S/∂t)

The 2018/19 winter snow accumulation was subtracted from the surface elevation measured in February 2019 at each GNSS measurement point to obtain the elevation of the buried end-of-summer 2018 horizon on 6 February 2019. A DEM of the end-of-summer 2018 horizon at that date was then interpolated. However, as illustrated in Figure 4 (dH3 arrow), the end-of-summer 2018 horizon measured on 6 February 2019 was deeper than it was at the end-of-summer 2018. To correct this change in elevation due to submergence, the submergence velocities determined previously were used to calculate the downward elevation change of the end-of-summer 2018 horizon over the 100 days leading up to 6 February 2019 (dH3 = V_sub Δt). It was then possible to calculate the mean annual surface elevation change between 2015 and 2018 (Fig. 5e).

This surface elevation change ranges between −0.5 and −2.1 m a⁻¹ across the study area, with spatial variations more pronounced where the observed 2019 accumulation is the lowest, i.e. in the upper area near the geographic pass and the lower area (East) where the slope increases.

5.6. Average annual surface mass-balance for 2015–18

To calculate the average annual surface mass-balance corresponding to the 2015–18 period over the area delimited by the radar survey, we subtracted the submergence velocities from the mean annual elevation change between 2015 and 2018.

The end-of-summer 2015 horizon DEM, at the end-of-summer 2018, was calculated from the 6 February 2019 end-of-summer 2015 horizon DEM by correcting for the change in elevation due to submergence between the end-of-summer 2018 and February 2019 (dH4 arrow in Fig. 4). To correct for this change in elevation, we used the previously calculated submergence velocities to calculate the downward displacement of the end-of-summer 2015 horizon over the 100 days leading up to 6 February 2019, just as we did for the end-of-summer 2018 horizon.

Figure 5f presents the interpolated 2015–18 average surface mass-balance with values ranging from 0.17 to 2.26 m w.e. a⁻¹ (mean 1.64 m w.e. a⁻¹). We assumed a mean firn density of 600 kg m⁻³ in this calculation. This value corresponds to the mean density measured in the S2 core for the firn between the two 2018 and 2015 horizons. Note that the density measured in the F2 core for the same period was a 13% lower (520 kg m⁻³). The uncertainty of the calculated average annual surface mass-balance was estimated to be ± 0.33 m w.e. a⁻¹ at site S2 where in situ density measurements were available, but could be larger elsewhere. Assuming a mean uncertainty of ± 1.2 m for the identification of the end-of-summer 2015 horizon on the radargrams and ± 80 kg m⁻³ for density, we estimated a mean uncertainty of ± 0.34 m w.e. a⁻¹ for the average annual surface mass-balance shown in Figure 5f.

These results are validated from in situ GLACIOCLIM measurements. These in situ measurements show an annual surface mass-balance of 3.03 ± 0.22, 1.05 ± 0.22 m w.e. and 3.16 ± 0.22 m w.e. for the hydrological years 2015/16, 2016/17 and 2017/18, respectively at the site S2 (Table 3). This gives a mean value of 2.41 ± 0.22 m w.e. a⁻¹ over the period 2015–2018. Using our method, we reconstructed a 3-year averaged annual surface mass-balance of 2.16 ± 0.33 m w.e. a⁻¹ at the same site, which is slightly lower than the in situ value but within the uncertainty range of the measurements.

5.7. Average annual surface mass-balance for 2012–21

The previous comparison of our reconstructed surface mass-balance at site S2 with GLACIOCLIM measurements over the 2015–2018 period are not completely independent given that GLACIOCLIM observations have been used to aid in identifying the buried end-of-summer 2015 horizon.

To strengthen the validation, we used an independent dataset over the period 2012–2021. For this purpose, the average annual surface mass-balance has been reconstructed from the previous submergence velocities and the surface elevation changes between 19 August 2012 and 15 August 2021 using two DEMs derived from Pléiades stereo images (Fig 8a), and compared directly to in situ observations of cumulative surface mass-balance over the period 2012–2021. Figure 8b shows the reconstructed average calculated surface mass-balance over the same period using a firn density of 550 ± 30 kg m⁻³ (based on GLACIOCLIM measurements). The average surface mass-balance varies from 0.67 to 3.15 m w.e. a⁻¹ (mean 2.17 m w.e. a⁻¹) over the studied area. At site S2, the measured surface elevation change is 0.08 ± 0.12 m a⁻¹. Using the previous value of −4.79 ± 0.46 m a⁻¹ for the submergence velocity at S2 site, this leads to a mean annual surface mass balance of 2.68 ± 0.30 m w.e. a⁻¹.

Figure 8. 10 m × 10 m interpolation maps of: (a) Mean annual surface elevation change between 19/08/2012 and 15/08/2021. (b) Mean annual surface mass-balance calculated between 19/08/2012 and 15/08/2021. The red pentagons show the snow or firn core locations and the grey dots the GNSS measurements on the GPR tracks.

The mean annual surface mass balance measured by the GLACIOCLIM observatory is 2.35 ± 0.21 m w.e. a⁻¹. The difference is within the uncertainty of measurements. Note also that the period 2012–2021 is not exactly the same between in situ observations (September) and the reconstructed values (August). However, GLACIOCLIM observations show that the ablation rate between August and September were very similar in 2012 and 2021. From these results, one can conclude that our method enables reconstruction of the surface mass balance over large areas with similar accuracy as in situ methods.

6. Discussion

6.1. Approximation of the estimated downward displacement of snow layers between end-of-summer 2018 and 6 February 2019

Because our objective is to quantify the surface mass-balance over the period from the end-of-summer 2015 to the end-of-summer 2018, the study would have been simpler if we had conducted the measurement campaign in late summer 2018. In our case, it was necessary to take into account the fact that from the end-of-summer 2018 to 6 February 2019, snow accumulated on the surface and at the same time the end-of-summer 2018 and 2015 horizons moved downward because of submergence (dH3 and dH4 variables in Fig. 4).

To account for this additional change in elevation when calculating surface elevation changes between 2015 and 2018 (Section 5.5), we used submergence velocities determined in Section 5.4 to calculate the downward displacement of the end-of-summer 2018 horizon over the 100 days leading up to 6 February 2019: dH3 = V_sub Δt. This assumes no temporal or seasonal variability for the submergence velocities, which is a reasonable assumption, at least at that site where stable submergence velocities throughout the course of a year have been previously observed (Réveillet and others, 2020).

The same calculation was applied in Section 5.6 for the end-of-summer 2015 horizon. This use of the submergence velocities calculated at the surface for a deeper layer implies that the two 2018 and 2015 horizons moved downward at the same speed during those 100 days (that is, the distances dH3 and dH4 in Fig. 4 are equal). This approximation neglects the compaction of the firn between the 2018 and 2015 horizons during this 100-day period. To estimate this compaction, we calculated the average core density between the two horizons with the given density profile (mean density of 600 kg m⁻³, SD of 80 kg m⁻³) and the average density between the horizons after shifting the density profile 1 m downward (mean density of 610 ± 80 kg m⁻³), corresponding to a compaction of approximately 2%. Applied to the approximately 11 m thick layer, this corresponds to a compaction of approximately 0.2 m, well below the ± 1.47 m uncertainty in the 2015 depth horizon calculation. We ran a firn compaction model adapted to mountain environments (Herron and Langway, 1980; Huss, 2013) and compared it to direct estimates of compaction. Output differences were less than uncertainty.

6.2. Comparison with previous submergence velocities measurements at the Col du Midi pass

Réveillet and others (2020) estimated an uncertainty of ± 0.5 m a⁻¹ (in metres of snow) for the submergence velocity data obtained from measurements performed in winter 2014/15. Table 5 and Figure 9 present a comparison with our submergence velocity data on eight sites (S1 to S5, M8, M9 and M12). Very close results (RMSE of 0.22 m w.e. a⁻¹) were obtained by both methods, although the observation periods are different (31 October 2014 to 27 May 2015 for the core method and 23 October 2015 to 6 February 2019 for the radar method).

Figure 9. Correlation between the submergence velocities (m w.e. a⁻¹) calculated using the radar method (this study) and the core method (Réveillet and others, 2020).

Table 5. Comparison of the submergence velocities, expressed in m w.e. a⁻¹, calculated at the same sites by the core method (see Réveillet and others, 2020) and the radar method proposed in this study

^a This study

^b Reveillet and others, 2020

The values of the radar method are the values interpolated by the kriging method. The mean density observed in the S2 firn core (540 ± 80 kg m⁻³) between the end-of-summer 2015 horizon and the February 2019 surface was used to convert into w.e. the values calculated with the radar method. For the core method, a density of 420 kg m⁻³, generally observed for winter snow at that site, was used.

In reality, the submergence velocities determined in our study are integrated values and not purely local values. They take into account the displacement of the glacier over the 3-year 2015–18 period and correspond to a mean value between the considered point and a few metres upstream. Réveillet and others (2020) showed that horizontal velocities are low in that area (between 5 and 14 m a⁻¹). The spatial variability of the submergence velocities observed in our study, within a few dozen metres, is of the same order as the estimated uncertainties (<0.30 m w.e. a⁻¹), and consequently does not question the validity of the comparison. For a detailed comparison of the principles of the two methods, see supplementary information C (Figs S6 and S7).

6.3. Comments on the proposed method to measure the spatial distribution of surface mass-balance using DEM differences

Our study shows that GPR measurements are an effective way to obtain a spatially distributed estimate of submergence velocities. Once these are known, and assuming that their temporal stability is confirmed by future studies, it is possible to calculate subsequent surface mass-balance using DEM differences alone.

6.3.1. Advantages and disadvantages of the method

The main difficulty with the method presented here is that it relies on an extended initial field campaign to determine submergence velocities. This fieldwork requires significant logistical and technical resources including a helicopter for transporting equipment and snow cores and analytical resources for processing the cores and radargrams. This is far more complex than the conventional glaciological accumulation measurements that can be conducted faster and with minimal logistics.

However, once this first step has been achieved, the advantage of this method is to obtain, for the following years, the surface mass-balances from differences in DEMs without conducting in situ measurements each year. In addition, spatially distributed surface mass-balances are obtained, as opposed to the glaciological method that provides only point measurements.

Given that Réveillet and others (2020) have shown that the submergence velocities are constant throughout most of the year (February – October), it is also possible to use our approach to measure seasonal surface mass-balances. For instance, the elevation difference between a DEM performed at the end of a summer season and another at the end of the following winter accumulation period could be used, after correction for the submergence velocities, to estimate the winter surface mass-balance. However, in all cases it will be necessary, at the time of acquisition of the new DEM, to know or estimate the value of the surface layer density involved in the calculation of Eqn (2). If no measurement is made at the site, this estimate should be based on previous measurements or on modelling.

This method seems to be well suited to relatively homogeneous accumulation areas, with low slopes, regardless of their size. A few tests conducted on the accumulation area of the Argentière glacier (Mont Blanc area) show that multi-year IRHs are also visible at this site (Jourdain and others, unpublished data). Sold and others (2015) also report such observations on the Findelengletscher glacier in Switzerland. The use of this method on steep or highly variable terrain is more uncertain, mainly because of the difficulty of using radar measurements on this type of terrain. First, the presence of crevasses often makes the radar signal very noisy, and second, radar measurements on steep terrain can be very challenging, except when using an airborne system (drone or helicopter).

Although not tested in our study, use of the method in polar regions may be possible, provided that the annual accumulation is sufficient to distinguish the buried annual horizons with the GPR. The method would likely be operational in Greenland (see for example the GPR measurements of Dunse and others (2008)). In Antarctica, the situation is probably more complex. In coastal regions, although a clear seasonal signal is recorded in firn core chemistry (Goursaud and others, 2017), this annual signal is rarely linked to a clear IRH (Le Meur, personal communication, 2021). Finally, the method is unlikely to be readily applicable on the Antarctic plateau where the very low accumulation and post deposition processes prevent the clear identification of seasonal snow horizons (Legrand and others, 1988; Gautier and others, 2016).

In any case, the density of radar measurements must be adapted to the area to be covered and the spatial variability that can be expected a priori or from previous studies. The same applies to the interpolation grid, which must be adapted to the density of measurements and the topography of the area. The time interval should also be adjusted to the annual accumulation rate. Finally, the frequency of the radar wave may have to be adjusted to take into account the annual accumulation rate of the site and the depth of the target horizon.

6.3.2. Uncertainties

The mean uncertainty of 0.34 m w.e. a⁻¹ in the surface mass-balance estimated by our approach comes from two main sources: (1) the uncertainty in the identification and detection of the deep N horizon for the year N + i (quantification of the submergence velocities) and (2) the uncertainty in the density used to convert into water equivalent.

1. 1. The uncertainty associated with the identification of the deep N horizon for the year N + i depends on the GPR resolution, the accuracy of the wave propagation velocity used to convert the two-way travel time into a depth and the accuracy of the depth of the firn core. This uncertainty is not uniform over the study area and depends on the depth considered. In our case, it is of the order of 10% of the depth, except at low depths where it becomes higher (12% at 5 m, 22% at 2 m). For the mean observed depth of the 2015 horizon within the studied area (12 m), the calculated uncertainty is ± 1.2 m.

   • • First, the identification of the depth of the target horizon relies on the dating of the snow layers based on the firn core studies. Assuming a well-marked seasonal variability in the firn cores (seasonality of chemical or physical parameters), the horizon can be dated from this seasonality and the accuracy of the depth of the horizon is of the order of a few tens of cm.

   • • Next, the identification of the corresponding IRH on the GPR profiles relies on the determination of the wave propagation velocity applied for time-to-depth conversion. It can be measured in the field by GPR common-midpoint measurements (Hubbard and Glasser, 2005) or, if the firn density is known, determined based on the relationships between firn density and wave propagation speed found in the literature (e.g. Kovacs and others, 1995). An alternative approach, used in our study, is to use the firn cores to constrain the radar survey. Direct observations in the firn cores make it possible to identify the successive annual past summer horizons. The wave propagation velocity is then adjusted in order to match well-marked IRHs with these annual summer horizons. However, firn core observations are not easy to acquire and depend on the quality of the snow profiles at each site. Although time-consuming, we recommend obtaining several firn cores in the field to better constrain the radar wave velocity. This approach with multiple firn cores has the advantage of taking into account the specificities of the site (temperate firn, impurity content, etc.) and the variability of the firn density within the studied area. For example, at site F2 where the accumulation is lower than at site S2, the density is higher at a given depth. However, it is lower if we consider an equivalent accumulation period. Indeed, on the period from the end-of-summer 2015 to end-of-summer 2018, density is 600 kg m⁻³ at site S2 but only 520 kg m⁻³ at site F2.

   • • Finally, the accurate picking of the selected IRH along the set of radar tracks relies on the quality of the radargrams. Although the choice of a deeper horizon improves the accuracy of calculated submergence velocities by incorporating a longer time period, the difficulty in accurately identifying the radar IRH increases with depth because of the attenuation of the radar signal in the snowpack. In our study, we might have increased our accuracy by considering the end-of-summer 2017 horizon instead of the 2015 horizon (but no surface DEM was available for end-of-summer 2017). The best compromise is likely different for each accumulation site, depending on the accumulation rate, the physico-chemical characteristics of the firn, as well as the radar equipment used.

2. 2. The use of the appropriate density to convert into water equivalent is the second limitation of the method. This source of uncertainty is not specific to the proposed method but is associated with any surface mass-balance measurement. Given that firn density is not homogeneous even within a single firn layer in the studied area, a mean value must be estimated. In our study, we used the value measured in the S2 firn core, and estimated an uncertainty of ± 80 kg m⁻³ based on the variability observed with respect to the F2 firn core.

   • • Firn cores are the best way to acquire in situ density measurements. In the case of spatially distributed measurements, the density must be estimated from one or more local measurements. The GLACIOCLIM measurements routinely conducted at the end of the accumulation season show that the density is quite homogeneous within a winter snowpack in the accumulation area of a glacier (variability of 5%). However, based on only two firn cores drilled at sites S2 and F2, our study shows that on a multi-year snowpack, density variability of more than 10% is observed in the firn even on a small spatial scale (<1 km²). This variability is probably linked to wind exposure and variability of snow accumulation and is the main source of uncertainty in the calculation of spatially distributed surface mass-balances. Measurements of the density at more numerous sites should reduce this uncertainty.

   • • However, this method aims to estimate operationally the surface mass-balance using new surface DEMs alone, once the submergence velocities have been determined. At this stage, an assumption on the density must be made for Eqn (2). The density measured in the firn cores during the initial extended field campaign can be used, assuming a temporal stability. Another approach to estimate this density is to use a firn densification model, as in Sold and others (2015), provided that it has been previously calibrated with in situ measurements.

6.3.3. Temporal stability of the submergence velocities

Finally, the feasibility of this approach relies on the stability of submergence velocities over time. From the definition of the submergence velocities (see section 4.1), note that the submergence velocity depends on the ice-flow velocities of glacier (U_s, V_s, W_s), the slopes ∂S/∂x and ∂S/∂y and the density of firn ρ_s when it is converted to m w.e. From previous studies, one can assume that these terms do not show significant changes at the annual scale. First, the components of ice flow velocity and surface slopes show few changes in accumulation areas from one year to the next (e.g. Stocker-Waldhuber and others, 2019). Second, the annual changes of the density of firn in accumulation areas are small at the decadal scale (Huss, 2013; Huss and others, 2015; Vincent and others, 2020). Consequently, we can assume that the vertical profile of density does not change significantly from one year to another.

Here, we showed that the mean submergence velocity of the 2015–19 period is similar to the one measured in 2014/15. We can therefore assume the stability of the submergence velocity at the Col du Midi pass over at least a five-year period (2015–19). Previous studies (Vincent and others, 1997, 2007, 2020; Zeller and others, 2023) have shown that changes in submergence velocities are small in accumulation areas on a decadal scale and appear to offer a good way of assessing the long-term average surface mass-balance.

Over several decades however, Vincent and others (2020) observed a gradual decrease in the submergence velocities between the 1997–2004 and 2016–17 periods because the studied glacier was not in a steady state. They noted that although the average difference is close to the uncertainty of submergence velocities, differences are systematic and the decrease reached 21% over the whole 20-year period, i.e. an average decrease of around 1% per year. It would therefore seem useful, in our approach, to plan decadal field campaigns to measure submergence velocities.

7. Conclusion

The classic way to determine the point surface mass-balance in the accumulation zone of a glacier is to dig pits or drill firn cores. However, this method requires major efforts to conduct annual field campaigns and only provides point measurements. This is often insufficient to provide a good representation of the spatial mass-balance variability as required to constrain surface mass-balance models or estimate precipitation at high elevations.

Here, we have shown that, in the accumulation zone, distributed surface mass-balances can be reconstructed over a large area from elevation changes alone provided that the submergence velocities have first been measured. The data required by our method are: (i) GPR measurements carried out to identify buried past horizons for the calculation of the submergence velocities and (ii) annual DEM differences that can be obtained from remote sensing if sufficiently accurate. With this data, the annual surface mass-balances can be estimated by subtracting the submergence velocity from the annual elevation changes.

In our experimental studies carried out at the Col du Midi pass, this method was used to reconstruct the spatially distributed surface mass-balances over a surface area of 0.6 km². For this, we used first DEMs derived from LiDAR and GNSS measurements performed in 2015 and 2019 respectively, and second DEMs derived from 2012 and 2021 Pléiades stereo images. In addition, we performed GPR measurements in 2019 to calculate the submergence velocities by comparing the elevation change over time of the end-of-summer 2015 horizon. Our method assumes that the submergence velocities are almost constant over a number of years. This is supported by our detailed observations performed on stakes. Comparison between the reconstructed point surface mass-balances and the observed values shows close agreement with a mean uncertainty of ± 0.34 m w.e. a⁻¹. Our results show that the surface mass-balance can then be calculated for various time scales (seasonal, annual, multi-year) as soon as a new surface DEM is available, whatever the technique used (LiDAR, UAV or satellite data), along with an estimate of the snow density. The main difficulty in this approach is the identification of the buried reference horizon and the estimate of the snow density.

Finally, the advantage of our method is two-fold. First, it can provide surface mass-balances distributed over large areas in accumulation zones where the classic glaciological method would require substantial fieldwork each year and would not provide the spatial variability of surface mass-balances. The method is well suited to relatively homogeneous accumulation areas, with low slopes, regardless of their size. Although not tested in our study, it is possible that the method could be used in polar areas, provided that annual accumulation is sufficient. Its application to steep or highly variable terrain is more uncertain, mainly because of the difficulty of using radar measurements in this type of terrain. Secondly, this method could be very useful in remote areas where in situ surface mass-balance measurements are not possible each year.