GNSS reflectometry from low-cost sensors for continuous in situ contemporaneous glacier mass balance and flux divergence

Abstract

Recent advances in remote sensing have produced global glacier surface elevation change data. Parsing these elevation change signals into contributions from the climate (i.e. climatic mass balance) and glacier dynamics (i.e. flux divergence) is critical to enhance our process-based understanding of glacier change. In this study, we evaluate three approaches for direct, continuous measurements of the climatic mass balance, flux divergence and elevation change at a site on Gulkana Glacier in Alaska using low-cost global navigation satellite system (GNSS) sensors, GNSS interferometric reflectometry (GNSS-IR), banded ablation stakes with time-lapse cameras and combinations thereof. Cumulative climatic mass balance over the season was 4.85 m and the three approaches were within 0.08 m through early July before the snowpack melted, and within 0.28 m through mid-August. The flux divergence increased from 0.52 ± 0.03 cm d⁻¹ before June 3 to about 0.73 cm d⁻¹ after June 27. We demonstrate a single GNSS system fixed atop an ablation stake can measure contemporaneous climatic mass balance, flux divergence and elevation change based on the antenna's position and GNSS-IR techniques. The ability of these systems to measure glacier mass balance and flux divergence offers unique opportunities for year-round observations on mountain glaciers in the future.

1. Introduction

Monitoring glacier change globally from remote sensing (Berthier and others, 2023) informs present-day and future estimates of glacier contributions to sea-level rise, water resources, natural hazards, climate and culture (Huss and Hock, 2018; Zemp and others, 2019; Rounce and others, 2023). Existing remotely sensed, geodetic observations measure spatially distributed total elevation change (e.g. Hugonnet and others, 2021; Jakob and Gourmelen, 2023), but this signal compounds the elevation change caused by climatic mass balance, firn densification and ice dynamics. Disaggregating total elevation change observations into climate-driven and ice flow-driven signals is critical to validate models and understand the discrete drivers of glacier change (Barandun and Pohl, 2023), but remains a challenge at large-scales (e.g. Miles and others, 2021).

The climate-driven signal is referred to as the climatic mass balance, defined as the sum of surface and internal accumulation and ablation (Cogley and others, 2011). Estimating the climatic mass balance at any point or area on a glacier requires the separation of total elevation change into its process-based components (i.e. the flux divergence and climatic mass balance) through mass continuity. In the ablation zone, mass continuity is captured by:

(1)$$\displaystyle{{{\rm d}h} \over {{\rm d}t}} = \displaystyle{{\;{\dot{b}}_{{\rm clim}}} \over \rho }-\nabla \cdot q$$

where (dh/dt) is the total elevation change in m a⁻¹, $\dot{b}_{{\rm clim}}$ is the climatic mass balance in mass per unit area per year, $\nabla \cdot q$ is the flux divergence in m a⁻¹, and ρ is the bulk density of snow, firn or ice. This formulation of $\dot{b}_{{\rm clim}}$ assumes that basal mass balance is negligible and recognizes that the density of seasonal snow on the surface is difficult to monitor and estimate, especially as the snow compacts over time. The flux divergence acts vertically, resulting from longitudinal strain caused by gradients in velocity and the advection of ice thickness gradients, defined as:

(2)$$\nabla \cdot q = h\cdot \left({\displaystyle{{{\rm d}u_x} \over {{\rm d}x}} + \displaystyle{{{\rm d}u_y} \over {{\rm d}y}}} \right) + u\cdot \left({\displaystyle{{{\rm d}h} \over {{\rm d}x}} + \displaystyle{{{\rm d}h} \over {{\rm d}y}}} \right)$$

where h is the ice thickness, u is the vertically averaged ice velocity, (du _x/dx) and (du _y/dy) are the along- and cross-flow velocity gradients, and (dh/dx) and (dh/dy) are the along- and cross-flow ice thickness gradients, respectively.

Remote-sensing techniques often derive the climatic mass balance by removing estimates of the flux divergence from the total elevation change (Miles and others, 2021; Pelto and Menounos, 2021; Van Tricht and others, 2021); however, the flux divergence estimates suffer from large errors and uncertainties in modeled ice thickness products and velocity data. The exact nature and origin of these errors is difficult to quantify without reliable ice thickness and velocity data (e.g. Pelto and Menounos, 2021) or reliable measurements of flux divergence (e.g. Vincent and others, 2021; Jourdain and others, 2023; Zeller and others, 2023). Regardless, these methods estimate the flux divergence at seasonal or annual timescales and consequently do not capture subseasonal changes. Understanding glacier changes as they continuously evolve thus requires contemporaneous measurements of the climatic mass balance, total elevation change, and flux divergence to adequately parse remotely sensed elevation change signals at the appropriate timescales.

Recent advances have enabled the continuous measurement of the climatic mass balance in situ over the ablation season with time-lapse cameras (Landmann and others, 2021; Cremona and others, 2023) and laser rangers (Wang and others, 2023) and/or the accumulation season with a cosmic ray sensor (Gugerli and others, 2019); however, none of these also measure contemporaneous flux divergence. Flux divergence has conventionally been measured using global navigation satellite system (GNSS) receivers mounted to ablation stakes (Anderson and others, 2004; Howat and others, 2008; Bartholomaus and others, 2011; Flowers and others, 2016; Armstrong and Anderson, 2020). While GNSS systems have been deployed for year-round measurements (e.g. Garbo and Mueller, 2024), systems mounted to stakes for direct estimates of flux divergence have only been deployed for ~30–80 days. These studies have found corresponding strain rates to vary substantially throughout the melt season (Anderson and others, 2004), show little variability (Howat and others, 2008), or display extension, compression or no trend depending on the distance being considered (Flowers and others, 2016; Armstrong and Anderson, 2020). The stark variability of glacier dynamics over portions of the ablation season in previous studies highlights the need for measurements over the entire ablation season (i.e. >80 days) and emphasizes the attention required with regards to parsing total elevation change signals into climatic mass balance and flux divergence. Large errors could arise from scaling subseasonal measurements (which may be anomalous) to longer time spans, which reinforces the importance of continuous data as validation for remote-sensing products with discrete acquisition times.

The development of open-source, low-cost GNSS systems designed for cryosphere applications provide unique opportunities for GNSS applications at unprecedented low costs and accessibility that have already been optimized for low-power sensing (Still and others, 2023; Garbo and Mueller, 2024; Pickell and Hawley, 2024). These GNSS systems further enable unique opportunities to measure surface melt and accumulation using GNSS interferometric reflectometry (GNSS-IR) (Roesler and Larson, 2018). Specifically, reflected signals from the surfaces around a GNSS antenna affect the interference pattern in the GNSS signal-to-noise ratio data, which can be used to estimate the distance between the surface and antenna for each satellite pass. The technique requires a relatively flat surface surrounding a GNSS receiver with an unobstructed view of the sky. GNSS-IR has thus been applied to ice sheets and seasonal snow depth (e.g. Larson and others, 2009; Larson and others, 2020; Steiner and others, 2023), lakes (Fagundes and others, 2021), tides (e.g. Larson and others, 2013; Peng and others, 2021) and soil moisture (Chew and others, 2016), but its suitability for mountain glaciers has scarcely been evaluated. While GNSS-IR has been employed to measure snow melt on a small glacier in Greenland (Dahl-Jensen and others, 2022), no studies have demonstrated the potential of GNSS-IR to measure ice ablation that is critical to understand mountain glacier mass balance.

In this study, we utilize low-cost, open-source GNSS systems with GNSS-IR to measure contemporaneous climatic mass balance and flux divergence over 4 months (129 days; 16 April 2023 to 23 August 2023) on Gulkana Glacier in Alaska. We deploy a collection of on- and off-ice GNSS systems with a monitored banded ablation stake to assess the accuracy and uncertainties of different daily climatic mass-balance estimates and GNSS-related measurements. We demonstrate three independent methods for measuring daily climatic mass balance with high accuracy. By leveraging GNSS-IR, we show that coincident, daily climatic mass balance and subseasonal flux divergence measurements can be obtained from a single GNSS system fixed to an ablation stake and discuss the potential for future applications including the accumulation season.

2. Study site

Gulkana Glacier is a mountain glacier located in the eastern Alaska Range (Fig. 1) that spans an elevation from 1253 to 2461 m and has an area of ~15.5 km² (U.S. Geological Survey, Benchmark Glacier Program, 2016). Gulkana Glacier experiences an interior climate with an average annual temperature of about −3℃ and average glacier-wide winter accumulation of ~0.9 m w.e. (O'Neel and others, 2019). The main trunk of the glacier flows south and spans nearly 7 km in total length from the top of the accumulation area. The western portion of the glacier features three separate accumulation basins of varying areas, which contribute relatively small fluxes into the main trunk (~20%). Gulkana Glacier is one of the ‘benchmark glaciers’ in Alaska with a long-term mass-balance program carried out by the U.S. Geological Survey (USGS) since 1966. Currently, an automatic weather station and a network of seven ablation stakes are installed and maintained biannually at the beginning and end of the melt season by the USGS.

Figure 1. Map of Gulkana Glacier showing deployed equipment (main panel) and Gulkana Glacier location in the eastern Alaska Range (inset). A 2021 outline is plotted atop an orthophoto from 30 August 2016 (McNeil and others, 2019). The GNSS base station, fixed GNSS, floating GNSS and monitored banded ablation stake sites are shown.

3. Methods

We collected daily mass balance and horizontal and vertical velocity data from 16 April 2023 to 23 August 2023, using two GNSS systems and a monitored banded ablation stake at long-term monitoring site AB on Gulkana Glacier (Fig. 1). Site AB was selected for this study as it is located in the middle of the ablation area. We also deployed GNSS systems in the middle of the accumulation area at site D and off-glacier on stable terrain at the Nunatak site as a base station. Poor weather conditions and the failure of our steam drill did not enable us to install a GNSS system on an ablation stake at site D, thus only horizontal velocity data are available at this site.

3.1 Field setup

Two GNSS systems, a time-lapse camera and a banded ablation stake (Fig. 2) were deployed to measure the climatic mass balance, total elevation change and flux divergence. The setup was designed to include redundancy to support an evaluation of how well the different methods measure each component.

Figure 2. The field setup upon deployment in the spring on 16 April 2023 (left) and retrieval in the fall on 23 August 2023 (right).

The monitored banded ablation stake used a trail camera that captured daily images from ~1 m. The ablation stake was spray painted in 10 cm bands alternating unpainted bands with rotating black, red and yellow colors to distinguish daily changes. Additional equipment and installation details are available in the supplement (see Supplementary Text S1).

Two ‘Cryologger Glacier Velocity Tracker (GVT)’ GNSS systems were installed at the same location as the banded ablation stake (Garbo, 2024; Garbo and Mueller, 2024). The Cryologger GVT is an open-source, low-cost instrument optimized for precise positioning in remote environments that is built entirely with off-the-shelf products and uses publicly available software. This system features an Arduino microcontroller that enables users to easily adjust desired sampling modes along with a multi-frequency, multi-constellation u-blox ZED-F9P GNSS receiver. The GNSS systems were operated daily for 6 h from 1:00 to 7:00 AKST at a sampling rate of 1 Hz and tracked all available satellite constellations (GPS, GLONASS, Galileo and BeiDou). An additional system was installed on off-glacier stable terrain as a base station to estimate our GNSS position accuracy. Detailed information about the design, assembly and software for the Cryologger GVT can be found on GitHub (Garbo, 2024). In the following, we refer to the Cryologger GVT generally as our ‘GNSS system(s)’.

Our setup included one GNSS system that was placed to ‘float’ on the glacier surface using a tripod and a second that was ‘fixed’ atop the banded ablation stake (Fig. 2). The floating system measured the glacier surface elevation each day, thus enabling the measurement of elevation change (dh/dt), similar to remotely sensed observations. The fixed system was unaffected by surface processes (i.e. accumulation and ablation) and instead measured the elevation change due to glacier dynamics, i.e. the flux divergence. Both systems also experience elevation changes caused by downhill glacier motion. We thus applied a slope correction to both systems using a 2 m resolution DEM from 19 September 2021 (McNeil and others, 2019) to account for downhill glacier motion (Fig. 3, see Supplementary Text S2), which effectively converted the measurements from a Eulerian to Lagrangian reference. Note the elevations measured by both the fixed and floating systems can be affected by ice-bed separation, caused by the evolving storage of liquid water in the subglacial hydrologic system at the ice-bed interface, on short timescales.

Figure 3. Schematic outlining all components of the field setup system for the three independent methods used to measure the climatic mass balance, total elevation change and flux divergence. Monitored banded ablation stake (green), fixed GNSS (teal) and floating GNSS (pink) are depicted. Initial positions are shown on the left (time 0) and final positions are shown on the right (time 1). The climatic mass balance $( \dot{b}_{{\rm clim}}/\rho )$, total elevation change (dh/dt) and flux divergence ($\nabla \cdot q$) can be calculated from GNSS elevation positions (z) and GNSS-IR reflector heights (H) (see Eqns (3–6)). Note the distance to the glacier bed is not drawn to scale.

3.2 Monitored banded ablation stake image detection

The position of the surface on the ablation stake was estimated daily from images using the tkinter module (Lundh, 1999) in Python for all days with suitable images (see Section 4.1). For each monitored banded ablation stake image, a GUI was used to identify the glacier surface on the stake relative to the nearest color band. We used the known ablation stake diameter as a reference length in the image to determine the glacier surface distance below or above a color band. Since each color band is a defined distance from the end of the ablation stake, we used these measurements to derive daily changes.

The uncertainty in the monitored banded ablation stake mass balance accounted for both detection error and a measurement uncertainty arising from errors in manual delineation of the glacier surface, stake band position and conversions from image pixel to actual distance. The uncertainty on any given day was estimated to be within 2 cm (see Supplementary Text S3), which is similar to other studies using monitored banded ablation stakes (Landmann and others, 2021; Cremona and others, 2023).

3.3 GNSS processing

Cryologger GVT GNSS data were stored locally on a microSD card in proprietary u-blox (.ubx) files. u-blox files were converted to standard RINEX 3.03 (Receiver Independent Exchange) files using the Emlid Studio desktop application (Emlid Tech Kft.). The antenna model and type were manually added to the RINEX file header following the NOAA Antenna Calibration codes. RINEX files were then converted to position using the Canadian Spatial Reference System Precise Point Positioning (CSRS-PPP) tool in kinematic mode to resolve small changes in position over time. The mean value of the CSRS-PPP solution was used to estimate one precise position each day.

GNSS position accuracy was estimated from the base station system. The base station system was processed with the CSRS-PPP tool in static mode. The output resulted in up to ~22 000 position solutions during a 6 h logging period. The exact number of observations each day varied based on satellite coverage and due to partial GNSS system failures associated with a bug in the u-blox ZED-F9P firmware and SparkFun Artemis Processor hardware that caused random premature data logging stoppages (Garbo, 2024). Thus, GNSS accuracy was analyzed with respect to the number of position solutions (Section 4.2). While position accuracy can be further improved by incorporating the base station in GNSS processing (see Supplementary Text S4), the improved positioning accuracy reduces data availability to only those dates with sufficient data from both the base and rover GNSS systems. Furthermore, measurement uncertainties are dominated by stake tilt, tripod migration and/or bed separation, which are all sources of error that will not be eliminated by improved accuracy in GNSS positioning.

3.4 GNSS-IR processing

GNSS-IR was performed using the Python gnssrefl software (Larson, 2024) to measure the vertical distance between the GNSS antenna phase center and the reflecting surface (reflector height). To that end, signal-to-noise ratio data were obtained from RINEX files using the rinex2snr function. The azimuth limit was set from 170° to 260°, which coincided with the uninterrupted sky view given that surrounding steep valley walls would block signals from the east and northwest. Limits were also set for reflector heights (0.4–8.0 m) and elevation angles (5°–25°), and the peak-to-noise ratio was set to 3.0. We executed the method with the gnssir function and consolidated data daily by taking the mean of accepted satellite tracks using the daily_avg function, removing solutions with fewer than three satellite tracks. GNSS-IR thus provided a daily estimate of the distance between the surface and antenna for an approximate footprint of 100–1000 m² around the system, depending on the antenna height and satellite elevation angle. Individual solutions (from which the daily average was taken) can be used to estimate uncertainty in reflector height for a given day. Further details on these processing methods are available in Roesler and Larson (2018).

The fixed system reflector height provides an independent measurement of the climatic mass balance at daily resolution (Fig. 3). The floating system reflector height allows us to correct the total elevation change measurements to account for the tripod sinking into the snow in the early summer and for tripod tilting on uneven local surface topography in the late summer (Fig. 3).

3.5 Continuity equation outputs

The combined monitored banded ablation stake and GNSS systems measure the climatic mass balance, total elevation change and flux divergence with built-in redundancy (Fig. 4). Specifically, we demonstrate three independent methods to isolate the climatic mass balance with daily resolution: (1) a monitored banded ablation stake, (2) the change in reflector height of a fixed GNSS system using GNSS-IR and (3) the difference between the GNSS elevations of the fixed and floating systems. We do not parse density (ρ) in our definition of height change attributable to climatic mass balance (e.g., $( \dot{b}_{{\rm clim}}/\rho )$ in Eqn (1)) to avoid introducing uncertainties in density assumptions that would be required to convert observed height changes to units of mass (Huss, 2013). Our results for climatic mass balance are thus reported as changes in height.

Figure 4. Schematic outlining input data, processing steps and final products for the continuity equation components. We show three independent field methods for daily climatic mass-balance measurements. The fixed system alone with a reference DEM can be used to obtain climatic mass balance, flux divergence and total elevation change at centimeter-level accuracy. Note that total elevation change (dh/dt) and flux divergence ($\nabla \cdot q$) measurements require a DEM to apply a slope correction and can be influenced by bed separation over short time scales.

The fixed and floating system elevations also isolate flux divergence and total elevation change, respectively. Alternative methods for flux divergence and total elevation change also exist; flux divergence can be derived as the residual of the monitored banded ablation stake climatic mass balance and the floating GNSS system total elevation change and total elevation change can be calculated from the fixed GNSS system position and GNSS-IR reflector height. Based on the system setup (Fig. 3), the climatic mass balance $( \dot{b}_{{\rm clim}}/\rho )$ is derived as:

(3)$$\displaystyle{{{\dot{b}}_{{\rm clim}}} \over \rho } = ( {( {z_{{\rm flo}, 1}-z_{{\rm flo}, 0}} ) -( {z_{{\rm fix}, 1}-z_{{\rm fix}, 0}} ) } ) -\Delta z_0-( {H_{{\rm flo}, 1}-H_{{\rm flo}, 0}} ) $$

(4)$$\displaystyle{{{\dot{b}}_{{\rm clim}}} \over \rho } = H_{{\rm fix}, 1}-H_{{\rm fix}, 0}$$

where H represents GNSS-IR reflector heights and z represents GNSS elevation positions for the fixed (fix) and floating (flo) systems at an initial (t = 0) and final (t = 1) position. The total elevation change (dh/dt) and flux divergence $( {\nabla \cdot q} )$ can similarly be determined:

(5)$$\displaystyle{{{\rm d}h} \over {{\rm d}t}} = ( {z_{{\rm flo}, 1}-z_{{\rm flo}, 0}} ) -\Delta z_{{\rm slope}} + z_{{\rm ice}/{\rm bed}}-( {H_{{\rm flo}, 1}-H_{{\rm flo}, 0}} ) $$

(6)$$\nabla \cdot q = \Delta z_{{\rm slope}}-z_{{\rm ice}/{\rm bed}}-( {z_{{\rm fix}, 1}-z_{{\rm fix}, 0}} ) $$

where Δz _slope is the elevation change resulting from motion down the glacier slope and z _ice/bed is the ice-bed separation. While bed separation can be ignored for seasonal measurements, it must be accounted for over short timescales in the ablation season. All GNSS system measurements for flux divergence and total elevation change require a reference DEM to correct for elevation changes due to their advection down glacier.

3.6 Velocity

Daily GNSS positions were used to estimate cumulative and weekly velocities. Cumulative velocity was calculated as the horizontal distance traveled since equipment deployment divided by the elapsed time (Eqn (7)). To account for data gaps, weekly velocity was calculated as the distance traveled between two dates that are at least 6 days apart divided by the elapsed time (Eqn (8)) and reported for the average of the two dates,

(7)$$v_{{\rm cumulative}} = \displaystyle{{d_{{\rm cumulative}}} \over {t_2-t_ 0}}$$

(8)$$v_{{\rm weekly}} = \displaystyle{{d_{{\rm weekly}}} \over {t_2-t_1}}$$

where d _cumulative is the distance traveled over the entire instrument deployment period, d _weekly is the distance traveled during the multi-day (>6 day) period, t ₀ is the study period start date, t ₁ is the weekly velocity start date, and t ₂ is the velocity end date.

Three primary sources of error affect the velocities: (1) GNSS positioning accuracy (see Section 4.2), (2) ablation stake tilt and (3) tripod migration while on ice. GNSS positioning accuracy affects velocities derived from all GNSS systems, while the ablation stake tilt only affects the fixed GNSS systems. Tripod migration refers to the slight downhill sliding of the tripod resulting from the tripod not being secured to the surface and was identified from the change in horizontal distance between the fixed and floating systems. This only occurred when the tripod was on ice, as the tripod sank or melted into the snow during the initial half of the study period thereby preventing sliding. Thus, this bias only applies to the floating system at site AB in the latter half of the study period. For ablation stake tilt, we assumed 5° of tilt uncertainty and used the fixed GNSS reflector height to calculate uncertainties in horizontal antenna positions. For tripod migration, we applied an uncertainty of 2.5 cm d⁻¹ (see Section 4.4) only in the downhill direction. Uncertainty in the distance was calculated from error in the cardinal directions summed in quadrature. Uncertainties in velocity (σ _v) are derived from uncertainties in the distance (σ _D) (Eqn (9)),

(9)$$\sigma _{\rm v} = \displaystyle{{\sigma _{\rm D}} \over {t_2-t_1}}$$

4. Results

4.1 Monitored banded ablation stake

The monitored banded ablation stake measured 4.85 m of height change between April 17 and August 23 (Fig. 5). Daily height change measurements were available for 104 of the 129 days (81% of days), although the camera successfully recorded images each day. Data gaps resulted from substantial accumulation that buried the camera, camera sinkage into snow or minor accumulation that prevented visibility of a color band edge, or significant camera tilt (Fig. S3).

Figure 5. Time series of monitored banded ablation stake results for site AB at Gulkana Glacier showing glacier surface type (top), cumulative climatic mass balance (middle) and daily climatic mass balance and air temperature (bottom) in 2023. Winter accumulation has completely melted on July 8 (blue line). Daily mass balance uncertainty (see Section 3.2) is shown during periods with data. Air temperature is from the USGS Nunatak automatic weather station (Baker and others, 2019: ver 2.0, 2024) and adjusted to site AB using a lapse rate of 6.5 ℃ km⁻¹.

The daily measurements revealed a significant increase (p<0.001) in melt rate after mid-June, corresponding to periods of higher air temperatures and after the snowpack likely ripened. Specifically, <20% of height change occurred during the first half of the ablation season, i.e. 0.85 m of change between April 17 and June 20 (64 days) compared to 4.00 m of change between June 20 and August 23 (65 days). Early spring surface lowering may be affected by snow compaction, but we expect this impact to be small considering the minor lowering during periods of negative air temperature (Fig. 5). Daily melt rates in the mid and late summer consistently exceeded 5 cm per day, with a maximum daily melt of 19.5 cm occurring on July 31. Our cumulative measurement agrees within 0.5 cm of the glaciological measurement that we measured at the end of the study period.

4.2 GNSS positioning accuracy

The accuracy of the GNSS system was estimated from our Cryologger GVT base system as a function of the number of observations (Fig. 6). With at least 1000 observations (~16 min), position accuracy was within 2 cm in the vertical and horizontal directions (Fig. 6). For days with at least 5000 observations (~80 min), the maximum error was 1.8 cm vertically and 0.8 cm in the XY direction. The mean absolute error was 0.6 cm vertically and 0.4 cm in the XY direction (0.3 cm in either the X or Y directions individually). We took a conservative approach and estimated the error of the GNSS systems in both the vertical and horizontal directions to be 1.2 cm (2–3 times the mean) for our analyses moving forward and note that errors are independent in time.

Figure 6. GNSS base station error in XY plane and Z direction error as a function of the number of observations. Error is measured as the difference between the known fixed base station position and the measured position for each day.

An analysis of absolute positions between our base station and the known USGS Nunatak base station revealed no bias in the vertical direction and ~4 cm of bias in the easting and northing direction (see Supplementary Text S5). However, because absolute GNSS position is only used to apply slope corrections, this bias does not impact our results.

Lastly, the fixed GNSS system continually melted out over the ablation season such that the antenna was ~5 m above the glacier surface at the end of the study period. Ablation stake tilt could cause a small amount of downward bias. We thus assumed a 5° ablation stake tilt uncertainty to incorporate this bias into the fixed GNSS system climatic mass-balance estimates.

4.3 GNSS-IR

GNSS-IR reflector heights from April 17 to August 23 showed that the height of the fixed system relative to the glacier surface continuously increased throughout the melt season (Fig. 7). The initial decrease in floating system height through June 1 was from the floating system tripod sinking into the snow. When the snow completely melted on July 8, as shown in camera imagery, the floating system reflector height was close to its initial height, as expected. After this time, the height continuously increased to a maximum of 25 cm higher than the initial setup, which indicates the tripod legs were on local topographic maxima. An increase in floating system height of 15 cm was measured in the field (Figs S1 and S2).

Figure 7. GNSS-IR results for the fixed and floating systems at Gulkana Glacier site AB from 17 April 2023 to 23 August 2023 show the relative height of the fixed system increases over time as the surface melts, while the floating system remains relatively constant. Shading shows the 95% confidence interval, which increases in late-summer due to surface roughness. Dashed lines represent the initial reflector height above the glacier surface.

The fixed system reflector height increased over time as the ablation stake melted out, consistent with observations from the monitored banded ablation stake of limited melt in the first half of the summer followed by a rapid increase through August (Fig. 5). Both the fixed and floating systems measured a snow event around June 3 that increased snow depth by 36 ± 2 cm, buried the monitored banded ablation stake camera, and temporarily decreased the GNSS system reflector heights (Fig. 7).

Mean uncertainty was calculated as the 2-sigma standard deviation of reflector height plus 2 cm to account for any systematic errors (Kristine Larson, University of Colorado Boulder, personal written communication, 19 February 2024). We further incorporated a surface heterogeneity uncertainty after ice was exposed to account for the decrease in GNSS-IR signal-to-noise ratio amplitude on rough surfaces (Larson and others, 2009); we empirically estimated this to be 10% of the melt signal after seasonal snow melt-out on July 8. Before July 8, the mean daily uncertainty in GNSS-IR measurements was 4.68 and 4.78 cm for the fixed and floating systems, respectively. After bare ice was exposed starting July 8, the mean daily uncertainty was 10.5 and 13.2 cm for the fixed and floating systems, respectively. We note that GNSS-IR uncertainty is not correlated with reflector height but is instead driven by the reflecting surface quality (i.e. if the surface is smooth or not). We observe no evidence for bias based on the surface slope or GNSS-IR footprint, as filtering the considered azimuth and elevation angle ranges did not impact results (see Supplementary Text S6).

4.4 Velocity

The weekly and cumulative velocities from the on-ice floating GNSS systems showed clear seasonal variations of a speedup and slowdown in June, the timing and magnitude of which varies for ablation zone site AB and accumulation zone site D (Fig. 8). At site AB, the velocity was ~19 m a⁻¹ from April 23 through June 3. The velocity doubled to 45 m a⁻¹ from June 3 to June 11, when it temporarily slowed, before increasing further to reach its peak velocity of 83 m a⁻¹ on June 19. The velocity then rapidly decreased through June 26 ultimately steadying at 30 m a⁻¹ (~1.5 times faster than the spring velocity) through August.

Figure 8. Time series of weekly (a) and cumulative (b) surface velocities at site AB (ablation area) and site D (accumulation area) from the GNSS systems in 2023. Shaded regions represent uncertainties from GNSS positioning, ablation stake tilt and tripod migration when applicable.

At site D in the accumulation area, located nearly 3.5 km up-glacier from site AB, the velocity was steady at ~53 m a⁻¹ until early June. From June 15 to June 26, the velocity almost doubled from 59 m a⁻¹ to its peak of 106 m a⁻¹. The velocity remained relatively constant for almost three weeks before slowing again. Ultimately, the velocity at site D dropped below its initial velocity and steadied at ~45 m a⁻¹. The average velocity for the entire 2023 summer, between April 19 and August 23, was 63.6 m a⁻¹ at site D and 32.1 m a⁻¹ at site AB. Further discussion about the timing and magnitude of velocity changes is in the Supplementary material (Supplementary Text S7).

Site AB provided an opportunity to evaluate uncertainties in velocities between the fixed and floating systems. The cumulative velocity uncertainties were <0.8 m a⁻¹ after 2 weeks, but weekly velocity uncertainties increased at site AB in the late-summer (Fig. 8). The velocities of both systems were similar through mid-July (cumulative velocities of 28.9 and 30.3 m a⁻¹). However, cumulative fixed and floating system velocities differed by 12.9% (28.2 and 32.1 m a⁻¹, respectively) by the end of the study period. The weekly velocity derived from the fixed system became noisy during this time period, as small changes in position (only ~0.5 m per week) were encumbered by motion associated with the potential tilt of the antenna mounted on the stake ~5 m above the ground (see Fig. 2, Fig. S8). Additionally, the floating system tripod appeared to slowly migrate downhill as this unanchored system slid on the ice. From mid-June until the end of the study period, we observed the floating system travel ~1.3 m farther than the fixed system in the downhill direction (Fig. S8), which amounts to a plausible estimate of 2.2–2.9 cm d⁻¹ (i.e. up to ~30% error in the late-summer floating GNSS weekly velocity signal).

4.5 Flux divergence and total elevation change

The elevation of the fixed and floating GNSS systems diverged over the study period as the fixed system remained at an almost constant height, while the floating system decreased in elevation (Fig. 9). The near constant elevation of the floating system suggests the rate of uplift from flux divergence (and bed separation on short timescales, see Section 5.2) was almost identical to the decrease in elevation from moving down the glacier. The decrease in elevation of the floating system was due both to the surface melting and the down glacier motion. The rest of this section will refer to the slope-corrected changes.

Figure 9. 2023 Gulkana Glacier site AB fixed and floating system absolute elevation (before applying the slope correction). The fixed system (teal) maintains its elevation throughout the summer despite the stake flowing down glacier (dashed line) due to the opposite vertical motion (upward) caused by flux divergence.

On seasonal time scales, the slope-corrected elevation change of the fixed GNSS system directly measures the flux divergence ($\nabla \cdot q$) and the slope-corrected elevation change of the floating GNSS system directly measures the total elevation change (dh/dt). The cumulative total elevation change signal measured by the floating system was −3.83 m by the end of the study period. This is more than 1 m smaller in magnitude than the −4.85 m for climatic mass balance measured by the ablation stake (Fig. 5) due to the impact of flux divergence (Eqn (1)). The cumulative total elevation change measured from the fixed GNSS system shows agreement within 11 cm of the total elevation change derived from the floating GNSS system, with only minor discrepancies occurring in late summer.

A mean seasonal flux divergence of −2.7 m a⁻¹ (−0.91 m) from the fixed GNSS system aligns with the −2.6 m a⁻¹ (−0.92 m) as calculated from standard glaciological measurements and GPS position, corresponding to a 4% difference between values. Note a slight mismatch in measurement timescales: the glaciological measurements spanned April 17 to August 23 while the first and last dates with fixed system GNSS-IR results were April 20 and August 19, respectively. In general, a reliable flux divergence signal appeared by May 3, after only 2 weeks in the early ablation season (Fig. S6). Deriving the flux divergence as the residual of the floating GNSS system total elevation change and monitored banded ablation stake climatic mass balance yields comparable (−2.5 m a⁻¹) results to the fixed GNSS system.

Determining subseasonal flux divergence is complicated, as emergence (i.e. flux divergence in the upward direction) is difficult to parse from the temporally dynamic signal of bed separation as the subglacial hydrologic system evolves over the course of the melt season (see Section 5.2). While bed separation can be difficult to quantify, we estimate an emergence of 0.52 cm d⁻¹ (1.9 m a⁻¹) before June 3 and an emergence of 0.69–0.77 cm d⁻¹ (2.5–2.8 m a⁻¹) after June 27, which corresponds to an ~40% increase in emergence velocity.

4.6 Climatic mass balance

The daily climatic mass balance was estimated from three independent methods: the monitored banded ablation stake, GNSS-IR and the difference between the fixed and floating systems (Fig. 10). The monitored banded ablation stake accurately estimates the climatic mass balance within 2 cm except for a few days with camera failures (Section 4.1). The climatic mass-balance estimates from GNSS-IR and the fixed-floating difference agree well with the stake measurements, with the added benefit of being able to capture accumulation (e.g. changes in snow depth around June 3). However, the fixed-floating difference estimates are only accurate when the floating system reflector height is used to correct for the height that the tripod sinks into the snow (Fig. S7).

Figure 10. Gulkana Glacier site AB climatic mass-balance method comparison. These include the monitored banded ablation stake (black), the reflector height of the fixed GNSS system (teal) and the difference between the fixed and floating GNSS systems adjusted with floating system reflector height (red). The glacier surface transitions from snow to ice at site AB around July 8 (blue line).

On the last day with a GNSS-IR measurement for the fixed system (August 13), there was a 0.28 m (6.1%) difference compared to the monitored banded ablation stake. On the final day with fixed and floating GNSS positions (August 18), we observed a −0.19 m (−4.0%) difference between the fixed-floating estimate and the monitored banded ablation stake. Before the surface snow melted out (July 8), the mean difference between the GNSS-IR and monitored banded ablation stake estimates was 0.08 and 0.05 m between the fixed-floating difference estimate and the monitored banded ablation stake. The less negative climatic mass balances from GNSS-IR, compared to the banded ablation stake (Fig. 10), is likely due to the spatial footprint of GNSS-IR that included both topographic highs and lows, while the ablation stake appeared to be in a topographic low (Fig. 2). For the fixed-floating system climatic mass balance, the floating system appeared to be on a lower relative surface (Fig. S2), which increased the perceived climatic mass balance in late summer compared to the other two methods. Still, these discrepancies are small (<0.3 m) and only appear in late summer over the final 2 weeks of the monitoring period. We note that even without the apparent topographic highs and lows, we would expect some discrepancy between methods given the spatial footprint of GNSS-IR versus the banded ablation stake that measures a single point. The observed discrepancy between point ablation stake and area-averaged GNSS-IR measurements is within the bounds of measurement uncertainty.

5. Discussion

5.1 GNSS-IR measures summer ablation and accumulation

Our results demonstrate that GNSS-IR can be used to obtain daily, high-accuracy measurements on mountain glaciers. The applicability of GNSS-IR may be limited by the surrounding topography, but even for a moderately sized (10¹ km² magnitude) glacier like Gulkana, GNSS-IR yields successful results in the south-flowing ablation area and the sheltered northeast-flowing accumulation area. GNSS-IR thus appears well suited to be applied at other mountain glacier sites, as our findings indicate the method only requires a narrow azimuth range with an unobstructed view of the near-horizon.

GNSS-IR has a few distinct advantages over daily ablation stake measurements. The footprint of reflected signals used for GNSS-IR yields a spatially distributed (100–1000 m²), rather than strictly point scale, mass-balance estimate. GNSS-IR is thus more robust to variability in local topography, which could introduce bias in ablation stake measurements; although this assumes the topography is not so rough that it leads to errors in GNSS-IR. We do not have the information needed to determine if the distributed GNSS-IR or point stake measurements are more accurate in our data (i.e. if the stake is in a topographic low). Still, GNSS-IR methods are shown to accurately measure ablation throughout the melt season. We expect this to be especially true in the accumulation area, where perennial snow cover results in a relatively smooth, consistent surface even at the end of the ablation season. Additionally, GNSS-IR is able to measure accumulation events, whereas these events are logged as data gaps by the monitored ablation stake. Specifically, two independent GNSS-IR methods for the climatic mass balance measured 0.36 ± 0.02 m of accumulation on June 3. Such quantification of accumulation is essential for energy-balance modeling and understanding the nuances of glacier mass balance, especially at high-temporal resolutions. Furthermore, banded ablation stakes must be installed at the end of the accumulation season, while GNSS-IR could theoretically be installed at any time as long as the GNSS receiver remains sufficiently above the surface (~0.4 m). In summary, GNSS-IR derived from a low-power GNSS system can be used to measure the climatic mass balance over the ablation season and has potential to measure accumulation over the winter in areas, although accumulation rates must be considered to ensure the system does not get buried.

5.2 Bed separation and emergence

Subseasonal total elevation change and flux divergence signals are affected by bed separation caused by changes in basal water storage (Fig. 11). We refer to the combination of emergence (i.e. flux divergence in the upward direction) and bed separation as uplift. During the spring, the velocity remains constant as little melt occurs (Fig. 10), and any melt that does occur likely refreezes in the cold snowpack such that there is no meltwater runoff from the glacier (Fig. 11a). We thus infer there is no water reaching the glacier bed and all the slope-corrected uplift can be attributed to flux divergence, which we calculate to be 0.52 cm d⁻¹ (1.9 m a⁻¹) from April 20 to June 3.

Figure 11. Time series of streamflow discharge (a), weekly velocity (b) and slope-corrected cumulative uplift from the fixed GNSS system (c) at site AB. Uplift is parsed into bed separation (d) and emergence (e) for four scenarios based on assumptions about bed separation during the drainage and efficiency phases (see Section 5.2). Streamflow discharge is from a stream gauge located ~2 km downstream from the terminus of Gulkana Glacier (U.S. Geological Survey, National Water Information System, 2016).

In June there is a sudden onset of increased velocity and uplift (Fig. 11). These observations are a hallmark of bed separation, which is caused by an increase in basal water pressure due to surface meltwater reaching the bed and an inefficient subglacial drainage network (e.g. Iken and others, 1983; Flowers and others, 2016). The maximum stream discharge occurs directly after peak velocity, on June 21 and 22 (Fig. 11a), which we surmise is the release of stored basal water that continues to drain until June 27. The glacier also slows, and uplift stagnates during this time, which we attribute to the bed separation decreasing as the basal water pressure is reduced. The stagnated uplift on June 23 reflects a slight lag following the onset of decreased velocity, consistent with previous literature that shows the onset of velocity decrease before bed separation (Harper and others, 2007; Howat and others, 2008). The observed rates of change in fixed GNSS system uplift (Fig. 11c) have similar patterns compared to previous observations (Anderson and others, 2004; Howat and others, 2008; Flowers and others, 2016).

To provide a first-order estimate of the flux divergence after the efficient drainage system is achieved, we make assumptions about the rate of bed separation and remove this from the uplift term. Four clear phases of basal hydrology are identified based on the uplift rate: spring, build-up, drainage and efficiency (Fig. 11). During each phase, flux divergence and bed separation are assumed to be constant. For the spring phase, we assume uplift is solely from emergence (0.52 cm d⁻¹), as stated above. During the subsequent build-up phase (June 3 to June 23), we assume emergence remains constant at 0.52 cm d⁻¹ (1.9 m a⁻¹) and all additional uplift is from bed separation. After the drainage phase, we provide a range of estimates for emergence based on four assumed bed separation scenarios. These scenarios include: (i) a full decrease during the drainage phase (‘full decrease’), (ii) a constant decrease over both the drainage and efficiency phases (‘constant decrease’), (iii) a partial decrease that accounts for half the bed separation during the drainage phase and the remainder during the efficiency phase (‘partial decrease’), and (iv) the same partial decrease scenario but extending the efficiency phase by 30 days (‘partial decrease + 30 day’).

Our full decrease and constant scenarios provide minimum (0.0 cm d⁻¹) and maximum (0.40 cm d⁻¹) bounds for bed separation during the efficiency phase, respectively. These correspond to a minimum emergence of 0.55 cm d⁻¹ and maximum emergence of 0.95 cm d⁻¹ during the efficiency phase. The minimum scenario reveals that emergence during the efficiency phase is at least 6% greater than the emergence during the spring phase. The partial decrease scenarios, which are more physically realistic (see Supplementary Text S8), estimate emergence during the efficiency phase in the range of 0.69–0.77 cm d⁻¹ (2.5–2.8 m a⁻¹). These emergence estimates correspond to an increase of ~40% from early to late summer, during which velocity increases by 70% (i.e. mean velocity of 20.8 m a⁻¹ before June 3 and 35.4 m a⁻¹ after June 26).

We note that quantifying the flux divergence is encumbered by a lack of information about the glacier bed, basal hydrology and strain rates. We also rely on key assumptions about constant emergence rates throughout the spring and build-up phases (an alternative set of assumptions yields comparable results as detailed in Supplementary Text S8). While the systems were deployed well before the onset of the velocity speed-up, they were removed before velocities returned to their initial magnitudes. To fully constrain the seasonal evolution of bed separation, GNSS systems would need to be deployed and maintained until bed separation has undoubtedly returned to zero. Additionally, differences between surface and bed slopes could be a source of error in distinguishing flux divergence from bed separation, particularly during sliding events.

5.3 A fixed GNSS system for all continuity equation components

The fixed GNSS system measures the climatic mass balance, flux divergence and total elevation change at subseasonal timescales via traditional GNSS processing and GNSS-IR. Thus, a fixed GNSS system enables in situ glacier data acquisition that is relevant for standalone studies and remote-sensing validation. In particular, the direct measurement of flux divergence in early and late summer offers information about the dynamic contribution to glacier thinning that is notoriously difficult, if not impossible, to obtain from remote sensing alone. Our results show a relatively marginal (~40%) increase in flux divergence despite substantial velocity changes (~70%), suggesting that increased velocity due to glacier sliding has little effect on changes in longitudinal strain. The large subseasonal variability in velocity changes at sites AB and D further suggests a complex relationship between changes in velocity and strain, reaffirming the importance of direct flux divergence measurements and urging caution when deriving flux divergence from velocity data. A denser network of GNSS systems across the glacier is required to concretely investigate the relationship between velocity and flux divergence.

GNSS-IR also facilitates contemporaneous measurements of climatic mass balance and flux divergence, offering insight into the ice dynamic and climate-related surface height change at a particular site on the glacier. GNSS-IR yields an area-averaged surface reflector height around the stake, thus providing a representative reflector height around the stake that is robust to local melt hotspots or extrusions. GNSS-IR is thus a robust method to measure all components of the continuity equation (Eqn (1)) using a single GNSS receiver, which vastly reduces the time, cost and potential points of failure from these field measurements in harsh environments.

6. Conclusions

In this study, we measured contemporaneous daily climatic mass balance and subseasonal flux divergence over 4 months on Gulkana Glacier in Alaska using low-cost, open-source GNSS systems with GNSS-IR and banded ablation stakes with time-lapse cameras. We developed and deployed three independent methods for estimating the climatic mass balance which were consistently within 5% of one another despite measuring different spatial footprints. GNSS-IR was also able to capture an accumulation event at the start of the melt season that buried the monitored ablation stake equipment. Furthermore, results show a subseasonal increase in flux divergence temporally aligned with changes in velocity, which is also affected by subglacial hydrology and bed separation.

Findings show that contemporaneous climatic mass balance, flux divergence and elevation change can be obtained from a single GNSS system fixed atop an ablation stake. Thus, the data provide critical high-temporal resolution calibration and validation data for modeled and remote sensing-derived products alike, especially for efforts that seek to separate flux divergence signals from total elevation change. As one of the first studies leveraging GNSS-IR on a mountain glacier, we see potential in GNSS-IR for future in situ investigations of mountain glaciers, as it reduces the barrier for otherwise difficult simultaneous climatic mass balance and flux divergence measurements. While our study deployed systems during the ablation season, we see the potential of GNSS-IR to measure winter snow depths as well.