1. Introduction
Snow accumulation and snow/ice melt on glaciers are directly influenced by atmospheric conditions, such that long-term monitoring of these processes provides valuable information on climate change (Kaser and others, Reference Kaser, Cogley, Dyurgerov, Meier and Ohmura2006; Gardner and others, Reference Gardner2013; Marzeion and others, Reference Marzeion, Cogley, Richter and Parkes2014; Huss and Hock, Reference Huss and Hock2018; Zemp and others, Reference Zemp2019). These processes directly control the glacier surface mass balance (SMB), a key indicator of glacier behavior and its direct response to changing climate conditions (Vincent and Moreau, Reference Vincent and Moreau2016). Traditionally, SMB has been estimated from point-based in situ observations by repeatedly measuring melt with ablation stakes and snow accumulation from pits or cores (e.g. Thibert and others, Reference Thibert, Blanc, Vincent and Eckert2008). Although this technique provides the longest continuous measurements, spanning multiple decades, its application is limited to a small number of glaciers worldwide, known as “reference glaciers”, which are used by the World Glacier Monitoring Service (WGMS) for global assessments of impacts of climate change on glacier mass balance (Zemp and others, Reference Zemp2019). These in situ measurements are labour intensive, as they require systematic field campaigns. As a result, they generally provide a temporal resolution limited to the number of field visits, typically a few measurements per year (at most on a monthly basis), each measurement representing the cumulative melt over that time interval (e.g. Francou and others, Reference Francou, Ribstein, Saravia and Tiriau1995; Rabatel and others, Reference Rabatel2013; Mölg and others, Reference Mölg, Ceballos, Huggel, Micheletti, Rabatel and Zemp2017). These limitations make it challenging to accurately constrain models that describe surface melt on the basis of atmospheric conditions. Various SMB models have been developed, ranging from empirical and statistical approaches that rely solely on temperature (Hock and Jensen, Reference Hock and Jensen1999; Braithwaite, Reference Braithwaite2002) to physics-based models that include all energy fluxes on the glacier surface. The widely used degree-day model (Réveillet and others, Reference Réveillet2018; Thibert and others, Reference Thibert, Dkengne Sielenou, Vionnet, Eckert and Vincent2018) has been successful in a range of contexts, but is highly dependent on the ability to constrain the degree-day factor from observations, particularly its spatial and temporal evolution over the course of the melt season (Hock and Jensen, Reference Hock and Jensen1999; Réveillet and others, Reference Réveillet, Vincent, Six and Rabatel2017).
Automated methods for quantifying SMB with higher temporal resolution have been developed in recent years. Ice ablation can be measured using pressure transducers (Bø ggild and others, Reference Boggild, Olesen, Andreas and Jorgensen2004); draw-wire sensors (Hulth, Reference Hulth2010); thermistor strings (Carturan and others, Reference Carturan, Cazorzi, Fontana and Zanoner2019), and time-lapse cameras combined with graduated stakes (Landmann and others, Reference Landmann, Künsch, Huss, Ogier, Kalisch and Farinotti2021), while snow accumulation is monitored with sonic depth sensors (Kinar and Pomeroy, Reference Kinar and Pomeroy2015); snow pillows (Johnson and others, Reference Johnson, Gelvin, Duvoy, Schaefer, Poole and Horton2015) and cosmic ray sensors (Howat and others, Reference Howat, de la Peña, Desilets and Womack2018; Gugerli and others, Reference Gugerli, Salzmann, Huss and Desilets2019). Although these methods provide high temporal resolution, their small spatial footprints limit their ability to capture spatial variability (Gutmann and others, Reference Gutmann, Larson, Williams, Nievinski and Zavorotny2012). Furthermore, the installation and maintenance of these instruments and their telemetry systems in harsh glacial environments remain costly and logistically challenging.
A widely used technique for SMB estimation is the geodetic method, which uses remote sensing data to quantify changes in surface elevation using satellite, airborne, or unmanned aerial vehicle carriers. This approach allows retrieving the glacier-wide SMB by differencing digital elevation models (DEMs) (Beraud and others, Reference Beraud, Cusicanqui, Rabatel, Brun, Vincent and Six2023; Berthier and others, Reference Berthier2023), allowing large-scale monitoring, even in remote regions. However, the method has three main shortcomings: (i) its temporal resolution is low, as it depends on the frequency of satellite revolutions or the number of aerial surveys (Berthier and others, Reference Berthier2014; Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017; Van Tricht and others, Reference Van Tricht, Huybrechts, Van Breedam, Vanhulle, Van Oost and Zekollari2021); (ii) as these systems measure changes in surface elevation rather than SMB directly, field measurements are still required to correct for the effects of ice dynamics (flux divergence) (Belart and others, Reference Belart2017; Réveillet and others, Reference Réveillet, MacDonell, Gascoin, Kinnard, Lhermitte and Schaffer2020; Kneib and others, Reference Kneib2024) and snow/ice densification (Fischer, Reference Fischer2011); and (iii) its accuracy on a homogeneous snow-covered surface is reduced due to the lack of distinguishable features on the surface that are needed for the correlation algorithm used to build DEMs from stereo images (Rolstad and others, Reference Rolstad, Haug and Denby2009).
Continuous Global Navigation Satellite System (GNSS) measurements on glacier surfaces are commonly used to study glacier dynamics (e.g. Anderson and others, Reference Anderson2004; Harper and others, Reference Harper, Humphrey, Pfeffer, Fudge and O’Neel2005; Bartholomew and others, Reference Bartholomew, Nienow, Sole, Mair, Cowton and King2012; Le Bris and others, Reference Le Bris2025), but can also be used to monitor changes in surface elevation using GNSS Interferometric Reflectometry (GNSS-IR) (Larson and others, Reference Larson, Gutmann, Zavorotny, Braun, Williams and Nievinski2009). This technique has been applied in a variety of contexts to monitor the environment around the GNSS antenna, including water level, snow accumulation, and soil moisture (Larson, Reference Larson2016). In glaciology, GNSS-IR has been proven to be effective in measuring snow depth (Larson and others, Reference Larson, Wahr and Munneke2015; Reference Larson, MacFerrin and Nylen2020; Dahl-Jensen and others, Reference Dahl-Jensen, Citterio, Jakobsen, Ahlstrøm, Larson and Khan2022) and even firn densification when combined with GNSS elevation changes (Larson and others, Reference Larson, Wahr and Munneke2015; Shean and others, Reference Shean2017); however, its application for measuring ice-surface ablation has been explored only in a few studies (Shean and others, Reference Shean2017, Wells and others, Reference Wells2024). Since GNSS-IR can detect changes in both snow and ice surfaces, it has the potential to become a valuable tool for systematically measuring SMB across various temporal scales (from daily to multi-year), reducing the number of instruments deployed on glaciers, and possibly enabling extraction of SMB estimates from past GNSS acquisitions. However, the evaluation of GNSS-IR for glacier applications and its associated uncertainties remains limited (Pickell and others, Reference Pickell, Hawley and LeWinter2025). Most studies have been conducted on the Greenland and Antarctic ice sheets (Shean and others, Reference Shean2017; Larson and others, Reference Larson, MacFerrin and Nylen2020; Dahl-Jensen and others, Reference Dahl-Jensen, Citterio, Jakobsen, Ahlstrøm, Larson and Khan2022; Steiner and others, Reference Steiner, Schmithüsen, Wickert and Eisen2023) rather than on mountain glaciers (Wells and others, Reference Wells2024). Moreover, GNSS-IR measurements tend to exhibit higher uncertainties on ice surfaces due to greater roughness, whereas performance is generally better over snow, which provides a smoother and more planar reflective surface (Larson and others, Reference Larson, Wahr and Munneke2015, Dahl-Jensen and others, Reference Dahl-Jensen, Citterio, Jakobsen, Ahlstrøm, Larson and Khan2022). However, the specific conditions that amplify uncertainties related to scale of surface roughness (ranging from centimetre to metre scales), its temporal variability, and potential mitigation strategies remain largely unexplored (Nievinski and Larson, Reference Nievinski and Larson2014a; Pickell and others, Reference Pickell, Hawley and LeWinter2025). In addition, other influencing factors, including surrounding topography (e.g. slope, obstructions), antenna height, and signal penetration depth, have not been thoroughly investigated (Gutmann and others, Reference Gutmann, Larson, Williams, Nievinski and Zavorotny2012).
In this study, we assess the feasibility of using GNSS-IR to continuously measure snow accumulation and glacier surface ablation over a 3-year period on an alpine glacier. We use measurements from 13 GNSS stations on Glacier d’Argentière in the French Alps (Togaibekov and others, Reference Togaibekov, Gimbert, Gilbert and Walpersdorf2024) that we compare with independent measurements from a range of techniques with different temporal resolutions and spatial coverage, as well as with predictions using the degree-day model.
2. Data and methods
2.1. Field setup
Glacier d’Argentière is located in the Mont-Blanc massif in the French Alps (Fig. 1). It has an area of approximately 10 km2, spanning an elevation range of 3400 to 1600 m a.s.l. (Vincent and Moreau, Reference Vincent and Moreau2016). In 2008, the glacier had a total length of approximately 10 km, with the equilibrium-line altitude (ELA) lying at about 2900 m a.s.l. (Vincent and others, Reference Vincent2021). The annual SMB ranges roughly from 2 m of water equivalent per year (m w.e. a−1) in the accumulation area to about −10 m w.e. a−1 close to the snout (Vincent and others, Reference Vincent, Soruco, Six and Meur2009). The lower part of the glacier is largely covered by rock debris, particularly in the vicinity of two longitudinal superficial moraines (Fig. 1b). Our study site is located in the ablation zone at approximately 2,400 m a.s.l. in the vicinity of these moraines, where the glacier has a relatively shallow slope (surface slope of approximately 10%) and a maximum thickness of about 250 m (Rabatel and others, Reference Rabatel, Sanchez, Vincent and Six2018) (Fig. 1a). Glacier d’Argentière is one of the 61 WGMS “reference glaciers” (Zemp and others, Reference Zemp2019) on which multidecadal SMB measurements have been conducted (Vincent and others, Reference Vincent2021).
Observational setup at Glacier d’Argentière. (a) Aerial picture of the study area in the ablation zone of the glacier. The red rectangle indicates the zone shown in the upper-right inset map. The orange areas next to the GNSS stations indicate the area covered by GNSS-IR measurements for antenna height of 2 m. The blue isolines in this inset show the glacier thickness in metres. (b) Ablation zone of Glacier d’Argentière. The photograph was taken by a drone in July 2022 when the glacier was not snow covered. The superficial moraines can be seen at the centre of the glacier. The terminus of the glacier is highly crevassed. (c) Reflective environment of a glacier valley. Mountains reflect satellite signals and block part of the sky obstructing GNSS signals. Ground reflection varies depending on the snow depth, and once the snow has melted, signals can be reflected directly from the ice surface (dashed lines). The photographs of (d) survey GNSS stations, (e) an ablation stake, (f) the SmartStake, and (g) the automatic weather station (AWS). Profile 4 corresponds to an elevation of 2400 m a.s.l. in the ablation zone of the glacier.
Five GNSS stations were deployed in the ablation zone of the glacier in April 2019, followed by seven more in February 2020 (red dots in Fig. 1a), as part of a study on the processes governing basal sliding (Vincent and others; Reference Vincent2022; Togaibekov and others, Reference Togaibekov, Gimbert, Gilbert and Walpersdorf2024; Roldán-Blasco and others, Reference Roldán-Blasco2025). We also include GNSS data from site ARGG, located near the ELA, which belongs to a different observatory network. These GNSS stations are equipped with multifrequency Leica GR25 receivers and Leica AS10 antennas (Fig. 1d), continuously recording GPS, GLONASS, and Galileo signals at 1 Hz. GNSS antennas are mounted on aluminium masts, anchored up to 6 m deep in the ice. The distance between neighbouring stations ranges from 50 to 200 m. To ensure the antenna masts remain upright and powered without interruption, we maintain the GNSS network every few months, and sometimes as frequently as every other week in summer, when melt causes the masts to tilt. We log all changes in antenna height during fieldwork. Unlike the ablation stakes, which are installed at the same locations every year, the GNSS stations have not been relocated throughout the duration of the study.
During the study period of 2019–2021, four to nine ice ablation measurements were made annually from the ablation stakes (green dots in Fig. 1a). The ablation stakes are 10 m long and consist of five 2 m long wooden sticks tied together with metallic chains (Fig. 1e). Annual SMB measurements in the ablation zone have an estimated uncertainty of ±0.14 m w.e. a−1 (Thibert and others, Reference Thibert, Blanc, Vincent and Eckert2008).
Ice ablation has also been automatically and continuously monitored with a SmartStake equipment (Fig. 1f) since 2019, located 30 m northwest from the GNSS station ARG6 (cyan dot in Fig. 1a). The SmartStake measures the length of a string anchored several meters deep in the ice and wound around an encoded wheel on a tripod at the surface. As the glacier surface melts, the tripod moves downwards, causing the string to wind onto the wheel, providing a continuous record of the ice-surface melt. The technique is similar to the draw-wire sensor (Hulth, Reference Hulth2010), although it is more compact and autonomous. The SmartStake has a temporal sampling interval of 30 minutes with a millimetric precision and becomes operational after the seasonal snow cover has melted on the glacier, typically in June or July depending on the amount of snow that accumulated during the previous winter season, until the first snowfall events in late September to early October.
Snow depth and density are measured once a year at the end of winter (typically in May) at five points near the GNSS network along the so-called “Profile 4” (blue line in Fig. 1a) using snow pits and probes. An automatic weather station (AWS) is also installed in the glacier’s ablation zone (Fig. 1g). We also use air temperature data obtained at 30-minute intervals from the SAFRAN meteorological reanalysis in the French Alps (Vernay and others, Reference Vernay2022).
2.2. The GNSS-IR technique
Unlike GNSS positioning, which relies on carrier phase and pseudorange observables, GNSS-IR uses signal-to-noise ratio (SNR) data produced by the constructive and destructive interference between direct and reflected signals (multipath effect) from a satellite (Roesler and Larson, Reference Roesler and Larson2018). This interference manifests as periodic SNR oscillations in each satellite track, which are especially prominent at low elevation angles e (less than 25 degrees). While interference introduces noise into the position time series, it can be used to estimate the reflector height HR through the following model:
\begin{equation}
\text{SNR}(e) = A(e) \sin \left( \frac{4\pi H_R}{\lambda} \sin(e) + \phi(e) \right),
\end{equation} where λ is a carrier wavelength, ϕ a phase constant, and A the amplitude of the SNR data (Nievinski and Larson, Reference Nievinski and Larson2014a, Reference Nievinski and Larson2014b). The frequency term 
$\frac{4\pi H_R}{\lambda} \sin(e)$ in Eq. 1 is primarily influenced by the additional path length traveled by the reflected signal (Fig. 1c) and is estimated using the Lomb-Scargle method, which is designed for detecting periodic signals in unevenly spaced data (Lomb, Reference Lomb1976). For a planar horizontal reflector, the frequency of interference between direct and reflected signals, as observed in SNR data, is a function of the sine of the elevation angle 
$\sin(e)$. Thus, the only unknown variable in the frequency term is the reflector height HR (Nievinski and Larson, Reference Nievinski and Larson2014a; Roesler and Larson, Reference Roesler and Larson2018).
We estimate HR using the gnssrefl software (Larson, Reference Larson2024), applying the parameter settings that are most effective for our survey: an elevation angle range of 5∘ to 25∘, a quality control parameter (ediff) of 2 which allows tracking of GNSS signals at elevation angles between 7∘ up to 23∘, and a peak-to-noise ratio of 2.8. The 1 Hz GNSS data are decimated to a 10-second interval to reduce processing load (Tabibi and others, Reference Tabibi, Geremia-Nievinski and van Dam2017). The software estimates a daily averaged HR from all rising and setting multi-GNSS satellites vehicles (SV) (GPS, GLONASS, and Galileo) over a 24-hour period. Due to topographical obstructions, low elevation angle signals are available only in the northwest quadrant (between 271∘ and 330∘) as indicated with orange sectors in the inset maps of Fig. 1a. To maximize the number of SNR tracks N, we employ SNR data from all available multi-GNSS frequencies, including GPS (L1, L2, L2C), GLONASS (L1, L2), and Galileo (E1, E5a, E5b, E6), which have shown effectiveness in the French Alps (Tabibi and others, Reference Tabibi, Geremia-Nievinski and van Dam2017).
A possible optimal antenna height for GNSS-IR measurements on the glacier would be on the order of 2 m, which has a footprint of approximately 1800 m2 in the case of our setup (orange sector in Fig. 1a, which represents 16
$\%$ of the full circular area). A higher antenna covers a larger area. However, maintaining a 2 m antenna height requires frequent adjustments to compensate for surface ablation or snowfall, leading to discontinuities. Although we log changes in antenna height during fieldwork, we do not rely solely on these values, especially with tilting antenna masts (a 10∘ tilt at 2 m height introduces a 3 cm error). Instead, we extrapolated the trend of the closest measurements in time to adjust for these discontinuities. Independent SMB measurements (such as the SmartStake equipment in our study) can serve as a reference to correct for discontinuities in the HR estimates. We then invert the continuous reflector height time series and translate it by setting the first value to zero. To avoid introducing uncertainties associated with density assumptions required for converting height changes to mass, the seasonal-scale SMB is reported in terms of height change (in metres).
2.3. Degree-day model
The degree-day model 
$M_d^{DD}$ is based on an empirical relationship between positive temperature and surface melt (Braithwaite and Olesen, Reference Braithwaite, Olesen and Oerlemans1989):
\begin{equation}
M_d^{DD} =
\begin{cases}
DDF_{snow/ice} \times T, & \text{for } T \gt 0\,^{\circ}\text{C}, \\
0, & \text{for } T \leq 0\,^{\circ}\text{C},
\end{cases}
\end{equation} where DDFsnow and DDFice are the degree-day factors for snow and ice (m w.e.
${}^{\circ}$ C−1 d−1), respectively, and T is daily mean air temperature above 
$0^{\circ}\mathrm{C}$. The output of the degree-day model is expressed in metres water equivalent (m w.e.) which provides a standardized measure of mass changes in terms of water. Réveillet and others Reference Réveillet, Vincent, Six and Rabatel(2017) determined glacier-wide DDFsnow of 0.0042 m w.e.
${}^{\circ}\mathrm{C}^{-1}\,\mathrm{d}^{-1}$ and DDFice of 0.0052 m w.e.
${}^{\circ}\,\mathrm{C}^{-1}\,\mathrm{d}^{-1}$ with associated uncertainties of about 10% based on long-term in situ measurements at Glacier d’Argentière (Vincent and Moreau, Reference Vincent and Moreau2016). We adjust DDFsnow and DDFice in order to best fit the GNSS-IR results at site ARG6, the site with the closest location to the in situ measurements, and then discuss these values compared to those found in Réveillet and others Reference Réveillet, Vincent, Six and Rabatel(2017). We integrate the melt time series from the degree-day model, considering periods when the air temperature exceeded 
$0^{\circ}\,\mathrm{C}$, to calculate cumulative SMB over the melt season. To compare the outputs of the degree-day model with the SMB derived from GNSS-IR, we convert the latter to m w.e., as required by the definition of the DDFsnow and DDFice. This conversion is performed assuming a constant density of 0.45 g cm−3 for snow, as measured from pits on May 25, 2021, and of 0.91 g cm−3 for ice. Since we focus on the period after the maximum snow height (typically in early May) when snow height decreases primarily due to melt rather than snow compaction, it is reasonable to use end-of-winter snow density from snow pit measurements to convert snow height to m w.e.
3. Results
3.1. GNSS SNR oscillation quality
The SNR pattern varies throughout the measurement period, influencing the quality of the reflector height HR estimates (Fig. 2). The number of retained SNR tracks N is higher when the glacier is covered with snow (two top panels in Fig. 2), and decreases as the melting period progresses, increasing the measurement uncertainty (1σ in Fig. 3a). The sinusoidal pattern of the SNR signals becomes noisier and sometimes exhibits irregular frequencies leading to more scattered HR estimates during the late ablation period (Fig. 2e–l).
Examples of detrended SNR data from the GPS, GLONASS, and Galileo constellations at the ARG6 site for six different cases in 2021, along with histograms showing the distribution of all the estimated reflector heights. The top two panels (a–d) display examples of SNR patterns for a mid-winter snow surface on March 4 (a,b) and a late-spring snow surface on May 20 (c,d). The middle two panels show SNR patterns before (e,f) and after (g,h) manually lowering the antenna height by 4.4 m on July 17 (DOY 198). The bottom two panels show SNR patterns before (i,j) and after (k,l) manually lowering the antenna height by 2.2 m on August 10 (DOY 222).
(a) Number of SNR tracks used for daily HR estimation with associated 1 σ error at the GNSS site ARG6. The vertical purple dashed lines correspond to the days of the year (DOY) shown in Figure 2: 1—DOY 063, 2—DOY 140, 3—DOY 191, 4—DOY 199, 5—DOY 219, and 6—DOY 223. (b–e) Site ARG6 before and after manually lowering the antenna height by 4.4 m on July 17, 2021 (DOY 198) (b,c) and lowering it by 2.2 m on August 10, 2021 (DOY 222) (d,e).
Significant frequency changes in SNR signals are observed when the antenna height is manually adjusted during fieldwork (Fig. 3b–e), leading to changes in the number of fringes, which is then directly related to the antenna height above the reflecting surface (Fig. 2e–l). We observe that variations in antenna height impact the uncertainties in the HR estimates. Interestingly, reducing the antenna height from approximately 6 m to less than 2 m (Fig. 3b, c) decreased σ (Fig. 2e–h), while lowering it from around 3 to 1 m (Fig. 3d, e) more than doubled σ (Fig. 2i–l). Potential explanations for these variations are discussed in Section 4.1.2.
3.2. Daily to seasonal variations
3.2.1. Overview
The analysis reveals that the cumulative glacier surface melt over the three years is approximately 20 m at 2400 m a.s.l. in the ablation zone, with noticeable annual variations (Fig. 4). As expected, the temporal extent of the melt period (light blue boxes in Fig. 4b) correlates with the duration of the period with positive air temperature (Fig. 4a). In 2019, high air temperatures and an extended ablation period resulted in the highest recorded surface melt, reaching nearly 10 m; whereas, the slightly cooler temperatures and shorter ablation periods of 2020 and 2021 led to approximately 2 m less surface melt. Surface ablation at higher elevations near the ELA at site ARGG (300 m above the Profile 4 area) is predictably lower.
Time series of different observations at Glacier d’Argentière from 2019 to 2021: (a) Air temperature and liquid precipitation, including the duration and mean of positive air temperature for each year; (b) Cumulative SMB derived from 13 GNSS sites, ablation stakes, and SmartStake measurements. The relative positions of the GNSS sites and the SmartStake along the flow direction are schematically shown in the top-right corner. Light blue boxes indicate the total annual surface mass loss corresponding to the ablation period. Horizontal black bars represent the lowest boundary reference for snow depth Hs.
The end of the ablation period varies yearly, ranging from September to November, and is consistent with air temperature and the occurrence of the first snowfalls (Fig. 4a). For example, in 2020, the air temperature fell below 
$0^{\circ}\,\mathrm{C}$ in September, and the SMB remained close to zero until December with limited precipitation. In contrast, during 2019 and 2021, the ablation period extended until November, followed by a rapid increase in surface elevation, coinciding with precipitation at temperatures below 
$0^{\circ}\,\mathrm{C}$, indicating snowfall. At these times, snowfall events (dotted vertical lines in Fig. 4) are captured well by repetitive short-term increases in SMB, often followed by slow decreases likely due to snow compaction. We use the end of the ablation period as a reference to distinguish between snow and ice surfaces (black horizontal bars in Fig. 4b), considering all measurements above this level as snow depth Hs. This reference indicates that snow cover remains until late June or early July, depending on the amount of snow accumulation throughout winter. A snow-depth measurement of 
$H_s=2.71\,\mathrm{m}$, taken from a pit near the site ARG6 on 20 May 2021, agrees with the GNSS-IR results within 1 cm (cyan star in Fig. 4b).
3.2.2. Comparison against SmartStake measurements
Zooming in on the ice ablation period (July–September) and examining the reflector height HR together with SmartStake data (flipped and brought to the level of the HR estimates to facilitate comparison, Fig. 5a, c), the GNSS-IR estimates of reflector height appear to be smoother during the first half of the ice ablation period (July to mid-August) and more scattered in the second half (mid-August to September), when the roughness of the ice surface increases. We define these periods as the low-roughness ice period (red dots) and the high-roughness ice period (blue dots), with the transition occurring in mid-August. However, the boundary between these two periods is also influenced by the height of the GNSS antenna (discussed later). In our specific case, we mark the boundary based on a field work day in mid-August when the antenna height was adjusted (Fig. 5a, c). As observed in the SNR data (Fig. 2), both low antenna heights (around 1 m) and high antenna heights (above 3 m) affect GNSS-IR measurement uncertainties; however, the increased data scatter in Fig. 5a–d indicates that this impact is particularly significant during the high-roughness period.
Evaluation of uncertainties in daily melt from the GNSS site ARG6 (
$M_d^{ARG6}$), compared to SmartStake measurements (
$M_d^{SS}$). (a, c) Reflector height time series and inverted SMB values from SmartStake during the ice ablation periods in 2020 (a) and 2021 (c). The discontinuity in elevation are due to antenna height manual lowering during fieldwork (vertical dashed red lines). (b, d) Variations in σ of the difference between 
$M_d^{ARG6}$ and 
$M_d^{SS}$ as a function of different moving average time windows. In (a–d), low-roughness (early ice) and high-roughness (late ice) ice periods are represented in red and blue, respectively. (e) Scatter plot of daily ice melt Md comparing GNSS-IR at ARG6 with SmartStake measurements for the 2020 and 2021 ablation seasons. (f) Corresponding histogram of Md differences between GNSS-IR and SmartStake measurements.
To evaluate the uncertainties of the GNSS-IR technique, we analyze the standard deviation (1 σ) of the differences between daily melt at station ARG6 
$M_d^{ARG6}$ and at SmartStake 
$M_d^{SS}$, using SmartStake data as a reference (Fig. 5b, d). On the daily time scale, σ increases as the summer ablation period progresses, consistent with an increase in surface roughness that causes noisy SNR oscillations (Fig. 2). In 2020, σ ranges from 0.036 m for low-roughness ice to 0.131 m for high-roughness ice, while in 2021 it ranges from 0.095 m for low-roughness ice to 0.184 m for high-roughness ice. σ being higher in 2021 compared to in 2020 is likely attributed to increased surface roughness (see Section 4.1.1). As expected, for any given period of the year, we observe that averaging daily melt over a longer period results in a lower σ. A 5-day moving average significantly reduces σ by a factor of 5 for low-roughness ice and by a factor of 6 for high-roughness ice (Fig. 5b, d). However, for high-roughness ice, a 14-day time window is required to achieve an accuracy level of 
$\sigma=0.02\,\mathrm{m}$. Therefore, we use a 5-day time window for low-roughness ice and a 14-day time window for high-roughness ice. After applying this filtering, the mean difference between SMB quantified with GNSS-IR versus with the SmartStake for the ice ablation periods in 2020 and 2021 averages to 0, with a standard deviation σ of 0.022 m (Fig. 5e,f).
3.2.3. Comparison against degree-day model predictions
We further compare the SMB measured from the GNSS-IR measurements with that predicted using the degree-day model (Fig. 6), and we assess the ability of the GNSS-IR measurements to constrain the degree-day factors for ice DDFice and snow DDFsnow. The degree-day factors for ice and snow are tuned to best fit the cumulative SMB derived from GNSS-IR measurements at site ARG6 (Fig. 6a–f). For DDFsnow, we use GNSS-IR estimates obtained daily during the snow ablation period (May to June), when temperatures remain consistently positive to constrain the degree-day model, while for DDFice, we use 
$M_d^{ARG6}$, averaged over 5-day periods for low-roughness ice and 14-day periods for high-roughness ice, as described in the previous section. Visual inspection reveals that the predictions from the degree-day model generally agree with the measurements from GNSS-IR but tend to diverge as the ablation season progresses (Fig. 6a–f). In particular, the degree-day model significantly overestimates the SMB values quantified from GNSS-IR toward the end of the ablation season, typically in late September. Since the estimated DDF varies slightly from year to year, we calculate the mean value to determine this coefficient for snow ablation DDFsnow and ice ablation DDFice. We find that DDFsnow is 0.0042 m w.e.
${}^{\circ}\,\mathrm{C}^{-1}\,\mathrm{d}^{-1}$, which is equal to the glacier-wide value previously reported by Réveillet and others Reference Réveillet, Vincent, Six and Rabatel(2017). In contrast, we find that DDFice is 0.0062 m w.e.
${}^{\circ}\,\mathrm{C}^{-1}\,\mathrm{d}^{-1}$, which is 0.001 m w.e.
${}^{\circ}\,\mathrm{C}^{-1}\,\mathrm{d}^{-1}$ higher than the glacier-wide value obtained by Réveillet and others Reference Réveillet, Vincent, Six and Rabatel(2017).
Comparison of SMB derived from GNSS-IR measurements and degree-day model predictions. (a–f) degree-day model predictions using the best-fit 
$DDF_{snow/ice}$ values based on cumulative SMB at site ARG6 during the snow, low-roughness ice (early ice), and high-roughness ice (late ice) ablation periods in 2020 and 2021. (g) Time series of daily melt simulated by the degree-day model and derived from the GNSS-IR measurements at site ARG6 in 2020 and 2021. (h, j) Scatter plots comparing GNSS-IR-derived daily melt 
$M_d^{ARG6}$ with degree-day model estimates 
$M_d^{DD}$ for (h) the snow surface and (j) the ice surface. (i, k) Histograms showing the distribution of daily melt differences between the respective measurements.
Using the obtained values of DDFsnow and DDFice, we model the daily SMB 
$M_d^{DD}$ and compare it to the daily melt derived from GNSS-IR at the site ARG6 
$M_d^{ARG6}$ (Fig. 6g). By analysing daily melt rather than the cumulative time series, we minimize the influence of accumulated errors. The degree-day model accurately reproduces the daily melt derived from GNSS-IR with a mean difference of −0.0001 ± 0.013 m w.e. for snow ablation and 0.006 ±0.024 m w.e. for ice ablation (Fig. 6h–k).
3.3. Spatial patterns
To assess spatial variations in SMB throughout the season, we divide the season into three periods: snow accumulation, snow ablation, and ice ablation. These periods are defined based on the dates of in situ measurements at ablation stakes and snow pits, as they provide independent SMB measurements across the study area for comparison with GNSS-IR estimates (Fig. 7). Snow and ice ablation (Fig. 7a–c) are calculated by simply taking the difference between extrema in the cumulative time series over the periods of interest (Fig. 4b). To compare GNSS-IR measurements with independent snow-pit observations conducted at the end of spring, we estimate snow depth from the GNSS-IR time series as the difference between the end of the ablation period (November 11, 2020) and the day of manual pit measurements (2021 May 20) for each GNSS site (Fig. 7d).
Spatial pattern of (a,b) ice ablation in 2020 and 2021, respectively, (c) snow ablation in 2021, and (d) snow deposition in 2021, derived from the GNSS-IR and in situ measurements. In situ measurements consist of either ablation stake (a,b) or snow pit (d) measurements, represented by the large circles that follow the corresponding color scale. The orange rectangle in (a) indicates the zone shown in the bottom-left inset map. The orange circular sector next to the GNSS site ARG6 indicates the area covered by GNSS-IR measurements for antenna height of 2 m. Blue numbers represent in situ measurement values, while red numbers correspond to GNSS-IR values. The color scales vary to accommodate the differences in the corresponding quantity in each panel and are adjusted to the range of GNSS-IR–derived values; as a result, some ablation stake circles that fall outside this range appear in black or white. Orange contour-lines in (d) illustrate glacier surface topography. SS in the captions of maps (a) and (b) stands for SmartStake.
GNSS stations located along the central flow line of the glacier show higher ice ablation compared to those at the edges of the area delimited by the GNSS network (Fig. 7a, b). In addition, we find that ice ablation as measured by the network of ablation stakes (large blue-colored scaled circles in Fig. 7a, b) is significantly more heterogeneous, with variability reaching up to 0.5 m between stakes only 10 m apart (inset map in Fig. 7a). The GNSS-IR technique, on the other hand, averages out this heterogeneity, providing smoothed SMB values over a larger and more representative area (see Section 4.3).
Unlike ice ablation, GNSS-IR-derived snow ablation shows a relatively homogeneous spatial pattern, with a maximum difference of less than 0.4 m across the entire network (Fig. 7c). The accumulation of snow throughout the study area varies from 1.86 to 2.71 m (Fig. 7d). GNSS stations along the central transect between two longitudinal moraines tend to record slightly higher values than those near the edges of our GNSS network. Two out of five manual snow depth measurements were conducted near the GNSS stations AR6D and ARG6, showing a strong agreement within 0.01 m (Fig. 7d). Unfortunately, station AR6G was destroyed by an avalanche in February 2021. The avalanche significantly increased the snow depth at pit 5, nearly doubling it compared to other GNSS sites.
4. Discussion
4.1. Uncertainties in GNSS-IR SMB measurements
Most previous studies were conducted on ice sheets (Shean and others, Reference Shean2017; Siegfried and others, Reference Siegfried, Medley, Larson, Fricker and Tulaczyk2017; Larson and others, Reference Larson, MacFerrin and Nylen2020; Dahl-Jensen and others, Reference Dahl-Jensen, Citterio, Jakobsen, Ahlstrøm, Larson and Khan2022; Hoffman and others, Reference Hoffman, Maclennan, Lenaerts, Larson and Christianson2025; Pickell and others, Reference Pickell, Hawley and LeWinter2025) or in relatively flat areas using Plate Boundary Observatory (PBO) GNSS stations (Larson and others, Reference Larson, Gutmann, Zavorotny, Braun, Williams and Nievinski2009; Gutmann and others, Reference Gutmann, Larson, Williams, Nievinski and Zavorotny2012; Larson and Nievinski, Reference Larson and Nievinski2013), where antennas receive signals omnidirectionally. In contrast, our study area is characterized by complex mountain topography, which reduces signal reception at low elevation angles by a factor of six, limiting the reception of low elevation signals to the glacier flow direction (northwest) (Fig. 1a). To increase the number of SNR observables under conditions of limited sky view, it is advantageous to use all available frequencies from multi-GNSS constellations (Tabibi and others, Reference Tabibi, Geremia-Nievinski and van Dam2017). However, we acknowledge that the multi-GNSS mode can introduce additional sources of error, such as unaccounted phase center variations (PCV) and increased noise in certain frequencies. For example, legacy GPS signals tend to be noisier than modernized signals such as L2C and L5 (Tabibi and others, Reference Tabibi, Geremia-Nievinski and van Dam2017). In addition to the previously mentioned sources of error, tropospheric delay may also introduce uncertainties, particularly for high antenna setups (Williams and Nievinski, Reference Williams and Nievinski2017). However, these errors remain minor compared to the other sources of uncertainties, which we discuss in the following sections.
4.1.1. The effect of surface roughness
GNSS-IR requires a surface around the antenna that acts as a specular reflector, which is typically the case when the surface is smooth, such as a snow-covered surface (Larson and others, Reference Larson, Gutmann, Zavorotny, Braun, Williams and Nievinski2009). However, as the melt season progresses, the small-scale roughness of the surface in the ablation zone of the glacier increases due to the presence of supraglacial debris and crevasses causing scattering of signals and leading to increased measurement errors from ±5 cm in the early melt period to ±14 cm in the late melt period. These values align with previous studies, which reported an uncertainty of ±5 cm on ice surfaces in Greenland (Dahl-Jensen and others, Reference Dahl-Jensen, Citterio, Jakobsen, Ahlstrøm, Larson and Khan2022) and ±13 cm in Alaska (Wells and others, Reference Wells2024). We also observed greater uncertainties in 2021 than in 2020, presumably due to glacier melting, which progressively releases more rock debris stored within its body each year. We also acknowledge that large-scale roughness (on the order of tens of metres), such as nearby moraines, might introduce measurement errors; however, since it is difficult to disentangle the effects of large-scale and small-scale roughness, we refer to their combined influence as random roughness (Nievinski and Larson, Reference Nievinski and Larson2014a; Pickell and others, Reference Pickell, Hawley and LeWinter2025). Nevertheless, our results demonstrate that applying a simple moving average smoothing technique can reduce the uncertainties on the melt estimates by a factor of six, depending on the surface roughness, resulting in an uncertainty of around ±2 cm which is lower than values reported in previous studies (Gutmann and others, Reference Gutmann, Larson, Williams, Nievinski and Zavorotny2012; Dahl-Jensen and others, Reference Dahl-Jensen, Citterio, Jakobsen, Ahlstrøm, Larson and Khan2022; Steiner and others, Reference Steiner, Schmithüsen, Wickert and Eisen2023; Wells and others, Reference Wells2024). A more thorough spectral analysis of SNR observables could further enhance precision beyond the simple moving average approach (Tabibi and others, Reference Tabibi, Geremia-Nievinski and van Dam2017).
On the other hand, the increase in uncertainty due to higher surface roughness can serve as an indicator of when the snow cover has melted, similar to recent studies on the evolution of surface roughness in the interior of the Greenland ice sheet (Pickell and others, Reference Pickell, Hawley and LeWinter2025). As surface roughness increases, the multipath modulation power decreases, reducing the visibility and range of interference fringes (Fig. 2) resulting in higher uncertainty in reflector height estimates (Nievinski and Larson, Reference Nievinski and Larson2014a, Tabibi and others, Reference Tabibi, Nievinski, van Dam and Monico2015). We also observed a sudden increase in uncertainties in early July, which marks the end of the snow ablation period (Fig. 3a). At this point, the GNSS-IR-derived SMB crosses the horizontal bar that indicates the ice level at the end of the previous melt season, confirming the boundary between the snow-covered and snow-free glacier surface. Furthermore, GNSS-reflected signals may offer more detailed information on the surface conditions of glaciers, including the presence of crevasses, moulins, supraglacial debris, and fine-scale roughness at the snow-ice interface (Gutmann and others, Reference Gutmann, Larson, Williams, Nievinski and Zavorotny2012; Steiner and others, Reference Steiner, Schmithüsen, Wickert and Eisen2023; Pickell and others, Reference Pickell, Hawley and LeWinter2025).
4.1.2. The effect of antenna height
The height of the antenna above the reflecting surface determines the frequency of SNR oscillations, with a higher antenna resulting in a higher frequency of SNR data (Eq. 1), favourable for observing multiple cycles of the modulation pattern, which is essential for accurately estimating HR (Nievinski and Larson, Reference Nievinski and Larson2014a). The effect of changing antenna height was particularly evident after field work on DOY 223 in 2021 (Fig. 3d, e), when the antenna height at site ARG6 was manually lowered to less than 1 m, resulting in low-frequency SNR data (Fig. 2k). The reduction in observed fringes, combined with the increased noise due to the higher surface roughness at this time of year, led to the highest uncertainties recorded during the study period of approximately ±0.14 m (Fig. 2k,l).
Since low elevation GNSS signals are only available from the glacier flow direction where the slope between the antenna location and the reflection point is negative (inset maps in Fig. 1a), GNSS-IR measurements tend to overestimate HR as the reflected signals originate from the down-valley direction relative to the antenna. This effect is particularly impactful for higher antennas, which have a larger footprint, thus increasing the range of the footprint and amplifying the associated error. Larson and Nievinski Reference Larson and Nievinski(2013) estimated that for an antenna height of 2 m, the error due to a surface slope of up to 8∘ (6∘ in our study area) ranges between 2 cm and 5 cm and can increase by a factor of three for an antenna height of 6 m. This effect was observed before and after fieldwork on DOY 195 (Fig. 2e–h), where the error decreased from 0.113 to 0.057 m as the antenna height at the site ARG6 was lowered from approximately 6 to 2 m (Fig. 3b, c). Although the surface slope effect contributed to the measurement uncertainty (Nievinski and Larson, Reference Nievinski and Larson2014a; Pickell and others, Reference Pickell, Hawley and LeWinter2025), our results indicate that the primary issue with an excessively high antenna is random surface roughness, as the larger footprint captures a wider area of rough surface on the glacier. We observe this effect when an antenna height of nearly 4 m exhibits lower uncertainties on the snow surface (Fig. 5a, b) compared to a similar antenna height on a bare ice surface (Fig. 5c, d). It is important to note that part of the mismatch (5–30
$\%$) between the GNSS-IR and SmartStake measurements may not be due solely to uncertainty, but to physical differences arising from one measurement being spatially integrative while the other being a point-scale one.
4.2. Using GNSS-IR to further constrain the degree-day model predictions
The continuous and high temporal resolution GNSS-IR measurements enable better constraints on the degree-day factor 
$DDF_{snow/ice}$, in particular its variation over the course of the melting season. During the early melt season (May–June) when the glacier is covered with snow, we assume a constant snow density of 0.45 g cm−3, based on pit measurements from May 20, 2021, which represents the high snow density conditions of the early ablation period. The accuracy of the chosen snow density is further supported by the fact that the best-fit DDFsnow is equal to the value previously reported in Réveillet and others Reference Réveillet, Vincent, Six and Rabatel(2017). The degree-day model for snow melt in the early ablation period provides an accurate estimate of the SMB, with daily melt mean differences from GNSS-IR estimates below ±15 mm (Fig. 6i).
During the rest of the melt season from July to September, when the glacier is snow-free, we observe mixed results. The mean daily melt difference between the GNSS-IR and degree-day model is close to zero, with a standard deviation of 0.024 m w.e., indicating that the degree-day model overall accurately simulates the SMB (Fig. 6). However, as the ablation season progresses, simulations with the degree-day model begin to deviate from the SMB quantified with the GNSS-IR, especially after mid- to late September, when the degree-day model significantly overestimates the observed daily melt likely due to a varying value of DDFice during the late ablation period. This suggests that the proportion of the role of air temperature among the various variables taken into account in the surface energy balance changes over time, potentially decreasing toward the end of the melt season. In contrast, the role of surface albedo and/or turbulent fluxes, likely linked to the increase in surface roughness (Braithwaite, Reference Braithwaite1995), appears to become more pronounced during the late ablation period. This offers valuable insights for future research, as the GNSS-IR technique can provide both high-temporal resolution data and information on ice surface roughness.
4.3. Interpreting spatial variations in GNSS-IR-derived SMB
The spatial extension of the GNSS network provides access to the analysis of complementary mechanisms driving the SMB variations, beyond its temporal evolution discussed in the previous sections. Indeed, the 12 GNSS sites can capture the spatial variability of both snow accumulation and surface melt. Moreover, the GNSS-IR technique quantifies an average reflector height over a certain surface (Fig. 7) rather than in a single point like traditional in situ measurements. We suggest that ablation stake measurements are strongly influenced by local preferential melt and consequently may be biased by the operator’s judgment when determining the average surface level around the stake during field surveys, which likely led to highly heterogeneous measurements at neighboring stakes (inset map in Fig. 7a). In contrast, the GNSS-IR approach, which is free from these biases, provides more representative estimates of SMB and allows better retrieving spatial variations in glacier surface melt (Fig. 7a, b). As a result, we are able to infer a higher ice melt along the central flow line of the glacier compared to the glacier borders, which would not be observed using traditional stake measurements. We interpret this as a higher surface impurity content along the central line, such as rock, sand, and dust (Fig. 1b), reducing the albedo and thus affecting the shortwave radiation budget (Warren, Reference Warren1982). Another contributing factor may be the higher exposure to sunlight along the central line, while the glacier borders remain in the shadow of surrounding mountains for a longer period.
Our GNSS network offers an unprecedented spatial coverage of snow depth measurements. The observed increase in snow depth towards the uppermost GNSS stations (Fig. 7d) is likely due to wind-driven snow redistribution, a well-documented factor influencing snow distribution patterns (Machguth and others, Reference Machguth, Eisen, Paul and Hoelzle2006, Dadic and others, Reference Dadic, Mott, Lehning and Burlando2010). The saddle-shaped area located between two longitudinal superficial moraines, where the central transect of the GNSS sites is located, probably also experiences increased snow accumulation since it is less exposed to wind erosion. In contrast to snow accumulation, the GNSS-IR measurements show that snow melts uniformly throughout the study area, with melt rates varying between sites only within the uncertainty of the GNSS-IR measurements (Fig. 7c).
5. Conclusions
We measure SMB variations associated with snow accumulation and snow/ice ablation using the GNSS-IR technique over three consecutive years at an environmentally and topographically challenging site, the ablation zone of Glacier d’Argentière in the French Alps. We find that GNSS-IR is a reliable method for autonomously measuring daily SMB in alpine glaciers, with a precision of ±2 cm d−1. This precision is achieved by averaging signals over 3–5 days during periods when the surface of the glacier is smooth (e.g. snow-covered surface or when the bare ice surface is relatively smooth), and over 6–14 days when surface roughness increases due to strong melt and the presence of debris. Thus, we identify the roughness of the bare ice surface as the main source of measurement uncertainty. In contrast, the snow-covered surface provides a planar reflection of GNSS signals, increasing the measurement precision to ±1 cm d−1, and eliminating the need for smoothing filters. Both excessively low and high GNSS antenna heights reduce the precision of measurements, particularly toward the end of the ablation period, further amplifying the impact of surface roughness. Therefore, we recommend maintaining the antenna height between 1.5 and 2.5 m.
Compared to in situ measurements, the high temporal resolution and spatial averaging provided by the GNSS-IR technique are particularly advantageous (i) to constrain the degree-day factors DDFsnow and DDFice and its potential variation through time; and (ii) to retrieve spatial patterns in ice ablation that would be difficult to retrieve otherwise with the more traditional stake measurements, which are point measurements and likely biased by local preferential melt. We observe that ice ablation is highest along the central flowline, likely due to reduced topographical shading and lower albedo. In contrast, snow accumulation and snow ablation exhibit a relatively uniform distribution. Application of the GNSS-IR technique to retrieve SMB more systematically, through revisiting past records or acquiring new ones, has the potential to help us better understand the different factors that control snow accumulation and snow and ice melt, and provide more reliable long-term SMB timeseries.
Data
The GNSS data used in this study are archived in the Oreme repository (Walpersdorf and others, Reference Walpersdorf2023a, Reference Walpersdorf2023b, Reference Walpersdorf2023c). Other data are available on the Zenodo repository platform (Vincent, Reference Vincent2021, Togaibekov, Reference Togaibekov2025).
Acknowledgements
This work is a part of the SAUSSURE project (ANR-18-CE01-0015-01) supported by the French National Research Agency (ANR). The GPS equipment of the temporary network was provided by the French national mobile parc GPSmob and permanent GPS data are from the French national geodetic network Rénag, both integrated in the French national research infrastructure Epos-France. We also thank many people who conducted the GNSS fieldwork on the Glacier d’Argentière, particularly Laurent Ott, Agnès Helmstetter, Christian Sue and Martin Champon. We thank Christian Vincent and Delphine Six for providing in situ glaciological data obtained within the framework of the French Service National d’Observation GLACIOCLIM (UGA-OSUG, CNRS-INSU, IRD, INRAE, Météo-France, IPEV). We also thank Sajad Tabibi for sharing his insights on the GNSS-IR technique. Finally, we acknowledge Romain Biron for his involvement in the development of the SmartStake.
Open access
