Monitoring the annual geodetic mass balance of Bordu and Sary-Tor glaciers using UAV data

Abstract

The geodetic mass balance of a glacier corresponds to glacier-wide volume changes, converted to mass changes using density assumptions. It is typically calculated by differencing multi-temporal digital elevation models. In this study, we show how the annual geodetic mass balance of a glacier can be derived from uncrewed aerial vehicle (UAV) data. The presented workflow is applied to two small- to medium-sized glaciers in the Kyrgyz Tien Shan (Central Asia): Bordu glacier and Sary-Tor glacier. The obtained geodetic mass balance is compared with the glaciological mass balance derived from a network of ablation stakes and snow pits. A previously calibrated mass-balance model is used to correct for the difference in acquisition dates. The results show that the determined geodetic mass balance matches closely with the glaciological mass balance. Besides, for both glaciers the geodetic mass balance does not seem to be particularly sensitive to the assumptions regarding volume-to-mass conversion. Therefore, our results demonstrate that UAVs can serve as a valuable instrument to quantify the annual geodetic mass balance and to validate the glaciological mass balance. The conventional glaciological mass-balance estimation often relies on interpolation and extrapolation methods, whereas UAVs offer the potential for direct data acquisition over the entire glacier surface.

1. Introduction

Glaciers are regarded as important indicators of climate change since their variations directly reflect changes in temperature and precipitation. They are therefore intensively monitored on the global scale (Zemp and others, 2013). One of the most important components of glacier monitoring is the surface mass balance (SMB), which is the direct response of a glacier to atmospheric conditions. It is controlled by the processes adding mass to the surface, such as snowfall, and those removing mass from the surface, such as melt as well as accumulated change in glacier extent and geometry. Traditionally, the SMB of glaciers has been evaluated with an in situ network of ablation stakes and snow pits, resulting in point measurements, followed by an interpolation and extrapolation to obtain the glacier-wide mean-specific mass balance (Zemp and others, 2013). This method, which is called the glaciological method, requires intensive fieldwork, and inaccessible areas can often not be covered. If these areas are numerous and large, and if the SMB pattern is heterogeneous, uncertainties in the obtained mean-specific mass balance can be substantial.

An alternative approach for glacier mass-balance monitoring is the geodetic mass balance, which corresponds to glacier-wide volume changes (Hugonnet and others, 2021; Beraud and others, 2023). Geodetic mass balances are typically determined by differencing multi-temporal digital elevation models (DEMs). Hence, to derive the geodetic mass balance, repeated and accurate mapping of a glacier's surface topography is required, which is usually done over multiyear to decadal periods. The assumption is that bedrock elevations are (virtually) constant which implies that surface-elevation changes correspond to ice thickness changes. By integrating these surface-elevation changes over the entire glacier area, ice volume changes are determined, which are then generally converted to mass changes using density assumptions (Huss, 2013). To date, geodetic mass balances are mainly used to retrieve glacier-wide, region-wide and worldwide glacier mass changes (Brun and others, 2017; Goerlich and others, 2017; Pelto and others, 2019; Shean and others, 2020; Bhattacharya and others, 2021; Hugonnet and others, 2021). However, geodetic mass balances are also used to calibrate and correct cumulative mass balances from the glaciological method (Fischer, 2011; Zemp and others, 2013), and to calibrate mass-balance models (Huss and Hock, 2015). Besides, it is important to mention that with the glaciological method only the SMB is measured, whereas the geodetic mass balance includes all volume changing processes (therefore also basal and internal mass balance as well as compaction at depth) (Belart and others, 2017; Beraud and others, 2023).

Lately, several remote-sensing platforms have been used to infer the geodetic mass balance, such as satellites and airplanes. In recent years, uncrewed aerial vehicles (UAVs) are increasingly used to monitor glaciers as well. They are well suited for surveying small to medium-sized glaciers. In cryospheric sciences, UAVs are mainly used to create high-resolution DEMs (Groos and others, 2019; Yang and others, 2020), to delineate glacier areas (Fugazza and others, 2015; Rossini and others, 2018), to derive surface velocities (Immerzeel and others, 2014; Ryan and others, 2015; Kraaijenbrink and others, 2016; Rossini and others, 2018; Benoit and others, 2019) and to monitor the mass balance (Cao and others, 2021; Vincent and others, 2021; Van Tricht and others, 2021b).

Geodetic mass-balance studies can be particularly useful in regions with sparse direct glacier observations, like Central Asia, where accessibility is poor and climate conditions are harsh (Brun and others, 2017; Shean and others, 2020). Yet the region is particularly vulnerable to the impacts of climate change, and it is one of the mountain regions where glaciers are an important fresh water source (Dikikh, 1993; Aizen and others, 1995; Dikikh and Hagg, 2004; Kaser and others, 2010). Especially during the dry season, the glacial meltwater contributes significantly to the water supply for hydropower production, irrigation and (daily) consumption (Sorg and others, 2012; Farinotti and others, 2015; Chen and others, 2019). It is thus of the utmost importance to study these glaciers.

In this context, our study shows how the annual geodetic mass balance can be derived for individual glaciers, using UAV surveys. The 2021/22 geodetic mass balance is derived for two Central Asian glaciers: Bordu glacier and Sary-Tor glacier. The derived geodetic mass balance is compared with the results from a previously calibrated mass-balance model (Van Tricht and others, 2021c) and with glaciological mass-balance values derived from ablation stake and snow pit data. The outputs are corrected for differences in acquisition dates to ensure an appropriate comparison. Drawing on our experience in on-site UAV glacier surveys, we elaborate on the advantages, disadvantages and pitfalls of this technique or methodology.

2. Bordu and Sary-Tor glaciers

In the Tien Shan, a Central Asian mountain range stretching over 3000 km from the east of Xinjiang, China, to the west of Tashkent, Uzbekistan, several glacier monitoring records date back to the 1950s and 1960s (Ushnurtsev, 1991; Dyurgerov and Mikhalenko, 1995; Bolch, 2015; Goerlich and others, 2017; Hoelzle and others, 2017; Barandun and others, 2018; WGMS, 2022). However, many monitoring programs were abandoned in the nineties to be only reinitiated in the last decade (Hoelzle and others, 2017; Satylkanov, 2018). In the context of this relaunch of measurement programs, the Bordu and Sary-Tor glaciers became part of a (re-)continued glacial measurement network, managed by the Tien Shan High Mountain Research Centre (TSHMRC).

In the summers of 2021 and 2022, we performed UAV surveys on Bordu glacier (41.81°N, 78.17°E) and Sary-Tor glacier (41.83°N, 78.18°E), two neighbouring glaciers located in the Inner Tien Shan on the north-western slopes of the Ak-Shyirak massif. The glaciers are situated within the Kumtor Gold Company concession, occupying the valleys adjacent to the Davydov glacier (to the east) and Glacier No. 354 (to the west). Both glaciers have an elevation range between 3900 and 4700/4800 m above sea level (a.s.l.) and consist of one large glacier tongue with smaller tributaries (Fig. 1). Ice thickness measurements were conducted on both glaciers revealing a maximum thickness of 159 m and an ice volume of 0.135 ± 0.001 km³ in 2013 for the Sary-Tor glacier (Petrakov and others, 2014). The Bordu glacier exhibited a maximum thickness of 148 m, accompanied by an estimated volume of 0.291 ± 0.058 km³ in 2017 (Van Tricht and others, 2021a). Between 1984 and 1989, SMB measurements were conducted on the Sary-Tor glacier and used to estimate the total mass balance in the entire Ak-Shyirak (Ushnurtsev, 1991; Dyurgerov and Mikhalenko, 1995). In 2014, these SMB measurements were reinitiated by the TSHMRC as part of the CHARIS Project (Contribution to High Asian Runoff from Ice and Snow) (Satylkanov, 2018). In 2015, SMB measurements were initiated on the Bordu glacier as well, and both measurement series continue today (Table 1). The mass-balance data are provided annually in altitudinal belts of 100 m (see supplementary Table 2). Due to numerous SMB and ice thickness measurements, as well as thermal regime (Van Tricht and Huybrechts, 2022) and internal accumulation studies (Golubev, 1976; Dyurgerov and Mikhalenko, 1995; Kronenberg and others, 2016), the Bordu and Sary-Tor glaciers are among the better monitored glaciers in the Tien Shan mountains. The availability of a continuous measurement series spanning seven recent years makes it possible to assess and validate novel techniques for glacier monitoring.

Fig. 1. Bordu and Sary-Tor glaciers located in Kyrgyzstan (Central Asia). The black outlines on the left map show the glacier extents in 2022. The former glacier tributary on the north-eastern margin of the Sary-Tor glacier is not connected to the main glacier tongue and is thus not considered in this study. The GCPs and GVPs used for the 2022 DEM are indicated with blue and red dots. The background image is a Sentinel-2 satellite image from 26 July 2021.

Table 1. Characteristics of both glaciers

3. Geodetic mass balance from UAV

3.1. Data acquisition

This study encompasses two fieldwork campaigns primarily dedicated to conducting UAV surveys on the two glaciers. During the first campaign in 2021, the glacier surfaces were almost snow-free at the time of the surveys, thereby providing optimal conditions for surveying. However, during the campaign in 2022, a significant portion of the glaciers was snow covered due to recent snowfall, making it more difficult to conduct the UAV surveys. The snow cover led to a diminished contrast, thereby reducing the number of discernible key points in the acquired images (Gindraux and others, 2017). Additionally, the snow cover posed challenges in terms of glacier accessibility, particularly in the higher areas prone to crevasses. In 2021, the processed area encompassed the entire extent of both glaciers. In 2022, because of the aforementioned reasons, the processed area was slightly smaller, covering 97% of the glacier area for the Sary-Tor glacier and 96% of the glacier area for the Bordu glacier.

The observations consist of images acquired by repeated UAV surveys with a DJI Phantom 4 RTK (P4RTK) quadcopter, equipped with a 20-megapixel camera. For flight planning and UAV piloting, DJI GS Pro software was employed. The flight plans were designed to minimise variations in ground sampling distance by flying parallel to the main surface slope. In general, the actual operational height (flight height above the operational area) ranged from ~180 to 200 m. To optimise the coverage at this height, the study area was divided into smaller sections (as shown in Fig. 2). Across all flight plans, sufficient overlap was maintained to facilitate the attachment and subsequent processing into one contiguous DEM.

Fig. 2. Take-off locations (red dots) for the survey on the Bordu glacier in 2022 and flight areas (coloured polygons). The right image shows the take-off location of number 5 and the flight area for the highest elevations of the Bordu glacier.

However, the flight height of 180–200 m above the surface could not be achieved for the higher parts. Factors like snow coverage and crevasses rendered these areas inaccessible, necessitating the execution of the surveys of the upper sections from a take-off location farther down and along the glacier margin (e.g. take-off spot number 5 in Fig. 2). In this case, the aircraft's height above the starting point was set at maximum 500 m (restricted by the software). This restriction of flying at a maximum of 500 m above the starting altitude poses limitations on the applicability of our workflow, as it prevents the surveying of areas at higher elevations, where achieving sufficient overlap becomes challenging.

Throughout the surveys, several orange plastic squares of 30 × 30 cm were strategically positioned as ground control points (GCPs) across the surface (Fig. 1). The precise positions of these GCPs were determined using a GNSS Trimble Geo7x, offering horizontal and vertical accuracies of ~0.1–0.2 m. Furthermore, several well-defined boulders were measured to augment the GCP density and ensure sufficient coverage. In this manner, we ensured a uniform distribution of GCPs within the ablation area and the lowest part of the accumulation area, except where crevasses or fresh snow impeded this (Fig. 1).

However, an advantage of the P4RTK is that the number of required GCPs is lower since it stores satellite observation data (GNSS data) that can be used for Post Processing Kinematics (PPK). This allows to correct the position of the UAV after the survey to achieve centimetre accuracy (in theory). As a result, some of our GCPs served as ground validation points (GVPs), meaning they were not used for georeferencing, but to verify the accuracy of the constructed DEM instead.

The GNSS data of the images were processed using the RTKLIB, an open-source program package designed for GNSS Positioning. To enhance the accuracy of the observation data, the geolocation of images captured by the P4RTK drone, before PPK, were corrected using data from the nearest operating UNAVCO permanent base station (Kumtor – KMTR). This base station is situated at a distance of 10 km from the glaciers.

Owing to the larger distances between the take-off point and the survey areas for the Bordu glacier, a smaller portion of each battery's flight time could be spent on the flight plan. Additionally, more time was required to reach the Bordu glacier from the basecamp, and to access the take-off locations. Consequently, in both 2021 and 2022, two survey-days were necessary for the Bordu glacier (specifically, on 10 and 11 August 2021, and 6 and 7 August 2022), while a single day was sufficient for the Sary-Tor glacier (12 August 2021 and 4 August 2022). The elevation difference between the areas surveyed on different days (1–2 days difference) is estimated to be negligible, as there was minimal or no surface melting due to the cold weather conditions and absence of snowfall during these days. Consequently, we did not apply any correction to account for the difference in acquisition dates.

3.2. Data processing

All acquired images from the UAV surveys were processed to generate DEMs using the photogrammetric workflow implemented in Pix4D. The resulting DEMs were established on a common local coordinate system (UTM44N) and referenced to the Earth's geoid (EGM96). To evaluate the accuracy of the DEM reconstructions, the GVPs, along with stable terrains outside the glacier boundaries, were employed. This is a typical approach in geodetic mass-balance studies (Dussaillant and others, 2018; Geissler and others, 2021; Hugonnet and others, 2021). The ice-free areas are assumed to have remained unchanged between two consecutive years, which is a valid assumption for exposed bedrock areas. Furthermore, the orthorectified images were combined into an orthophoto (geometrically correct combined image), which is a procedure to remove the distortions from the images using the obtained DEM.

To determine the geodetic mass balance, elevation differences between both years are calculated by subtracting the acquired DEMs. However, prior to this subtraction, some pre-processing steps are required. First, the pairs of DEMs were resampled to a common spatial resolution using bi-linear resampling. For both glaciers, the DEM with the higher resolution (2022) was resampled to match the resolution of the DEM with the lower resolution (2021). Subsequently, because small misalignments between DEMs can have a large impact on the derivation of elevation changes, the DEMs were co-registered by matching distinguishable and stable areas in Pix4D across the two consecutive years. For instance, in the (stable) pro-glacial areas, additional control points or boulders were selected that were easily recognisable on both DEMs. The orthophotos were used to delineate the glaciers.

3.3. Elevation changes and geodetic mass balance

After acquiring the data (section 3.1) and constructing, resampling and co-registering the DEMs (section 3.2), the DEMs are subtracted from each other to obtain elevation differences (Δh _xy). To ensure negative elevation changes where the surface elevation decreased over time, the oldest DEM is subtracted from the more recent one in each case (Eqn (1)). Next, the total volume change is calculated (Eqn (2)). Subsequently, the geodetic mass balance mb_geo (m w.e. a⁻¹) is determined (Eqn (3)).

(1)$$\Delta h_{xy} = h_{2022, xy}-h_{2021, xy}$$

(2)$$\Delta V = {\rm \;}r^2\mathop \sum \limits_{xy}^{A_{\rm t}} \Delta h_{xy}$$

(3)$${\rm m}{\rm b}_{{\rm geo}} = \displaystyle{{\Delta V\;f_{\rm V}} \over {\rho_{\rm w} A_{\rm t}( {t_2-t_1} ) }}\;$$

x and y represent the x ^th and y ^th spatial component. h_xy is the local surface elevation (m a.s.l.). A _t is the average glacier area of 2021 and 2022, ΔV is the change in volume (m³), which is obtained by summing the individual observed elevation change values Δh _xy across the glacier, multiplied by the square of the pixel size r ² (Eqn (2)). f _V denotes the volume-to-mass conversion value (kg m^–3). ρ _w is the density of water (1000 kg m⁻³). t ₂ is equal to the survey date in 2022 while t ₁ is the survey date in 2021.

While the elevation changes can directly be calculated by subtracting the DEMs from each other (Eqn (1)), the volume-to-mass conversion value varies in time and space and is more difficult to determine (Huss, 2013). In large-scale geodetic mass-balance studies, a volume-to-mass conversion value of 850 kg m⁻³ is commonly employed to convert volume changes to mass changes (Huss, 2013; Brun and others, 2017; Shean and others, 2020). This value is primarily applicable for longer periods, stable mass-balance gradients and large volume changes (Huss, 2013). At smaller timescales, it is recommended to account for different surface types (ice, snow, firn) (Pelto and others, 2019; Beraud and others, 2023).

As our objective is to develop a UAV-based method that can easily be applied to other glaciers with limited data availability, our primary aim is to ensure simplicity for determining the geodetic mass balance. To achieve this, we initially rely on the volume-to-mass conversion factor proposed by Huss (2013) in the standard configuration. Nonetheless, in order to assess the influence of a more intricate approach, we also adopt an alternative method to determine the volume-to-mass conversion value as part of the uncertainty and sensitivity analysis in section 5. This method involves incorporating surface densities and ELA estimates, integrating elements from Geissler and others (2021), Pelto and others (2019) and Beraud and others (2023).

3.4. Uncertainty and error estimate

The overall uncertainty of the geodetic mass balance $\sigma _{{\rm m}{\rm b}_{{\rm geo}}}$ is calculated using the method applied in Pelto and others (2019) and Beraud and others (2023) (Eqn (4)). The assumption is that the uncertainty has two independent sources: (i) the volume change uncertainty σ _ΔV and (ii) the uncertainty related to the volume-to-mass conversion factor $\sigma _{f_{\Delta {\rm V}}}$. For the latter uncertainty, we assume a range of ±60 kg m⁻³ following Huss (2013).

(4)$$\Sigma _{{\rm m}{\rm b}_{{\rm geo}}} = \sqrt {{\left({\displaystyle{{\sigma_{\Delta {\rm V}}f_{\Delta {\rm V}}} \over {\rho_{\rm w}A_{\rm t}( {t_2-t_1} ) }}} \right)}^2 + {\left({\displaystyle{{\sigma_{\,f_{\Delta {\rm V}}}\Delta V} \over {\rho_{\rm w}A_{\rm t}( {t_2-t_1} ) }}} \right)}^2}$$

The uncertainty associated with the volume change resulting from the observed elevation changes is calculated as the square root of the quadratic sum of a term related to the std dev. of elevation changes over the stable terrain (random error $\sigma _{\Delta h_{xy, {\rm random}}})$ and a term related to the mean elevation changes over the stable terrain (systematic error $\sigma _{\Delta h_{{\rm systematic}}})$ (Eqn (5)). The random error is assumed to be spatially variable based on the slope of the terrain (see section 4.2 for details).

(5)$$\sigma _{\Delta {\rm V}} = \sqrt {\left({\sqrt {\displaystyle{{\pi L^2} \over {5A_{\rm t}}}} \times } \right.{\left({\left. {r^2\mathop \sum \limits_{xy}^{A_{\rm t}} \sigma_{\Delta h_{xy, {\rm random}}}} \right)} \right)}^2 + {( {A_{\rm t}\;\sigma_{\Delta h_{{\rm systematic}}}} ) }^2} \;$$

The cumulative random error (first term under the square root of Eqn (5)) is scaled by a correction factor $\left({\sqrt {\pi L^2/5A_{\rm t}} } \right)$, to account for the fact that the mean error will tend to approach zero for large areas characterised by uncorrelated random error (Rolstad and others, 2009; Brun and others, 2017; Shean and others, 2020). Following Pelto and others (2019) and Menounos and others (2019), the determination of L, the decorrelation length, facilitated through semi-variance analysis, yields a value of 250 m for both glaciers. This value is lower than the value of 500 m which is typically used for lower resolution satellite data (Brun and others, 2017; Shean and others, 2020), and slightly lower than 300 m which was found in the study of Menounos and others (2019) for higher resolution data. Using a decorrelation length of 250 m results in correction factors of 0.19 for the Bordu glacier and 0.27 for the Sary-Tor glacier. The delineation uncertainty, a factor commonly accounted for as well in lower resolution data (e.g. Pelto and others, 2019; Beraud and others, 2023), is assumed to be insignificant here due to the very high resolution of the data at centimetre scale.

3.5. Model data and comparison

3.5.1. Glaciological mass balance

The glaciological mass-balance data are based on stake readings from ablation stakes at the end of the mass-balance season, as well as measurements of snow and firn depth from snow and firn pits or probings in the accumulation area (supplementary Table 1). Subsequently, these measurements are interpolated to provide mean-specific mass balance values in elevation bands of typically 100 m covering the entire glaciers (Van Tricht and others, 2021c) (Fig. 3) (supplementary Table 2). By considering the area of the 100 m bands, the glacier-integrated mean-specific mass balance is calculated. The interpolation of sparse measurements in the accumulation area relies on (local) understanding of long-term accumulation patterns. Nevertheless, given that the maximum mass balance in the accumulation area is relatively small (typically ranging from 0.2 to 0.6 m w.e. a⁻¹ due to the aridity of the climate), errors introduced by uncertainties in the accumulation area are expected to have a minor impact on the mean-specific mass balance. The SMB is usually determined using a stratigraphic system, i.e. mass difference between formation dates of two summer surfaces, which is usually mid-September. For the Bordu and Sary-Tor glaciers, detailed glaciological mass-balance data are available for recent years (Table 1 in the paper and supplementary Tables 1 and 2). In this study, we compare the derived geodetic mass balance with the 2021/22 glaciological mass balance, which was −1.27 ± 0.15 m w.e. a⁻¹ for both glaciers (WGMS, 2022).

Fig. 3. Mean-specific glaciological surface mass balance (MB) in elevation bins covering the entire glaciers. The (areas of the) elevation bins are indicated with grey rectangles.

3.5.2. Mass-balance model

To account for the temporal difference between the geodetic mass balance and the glaciological mass balance, we introduce correction balances, denoted as mb_{cor_2021} and mb_{cor_2022}, which are modelled mass balances mb_mod between the dates of the UAV surveys and the end of the glaciological year. These correction balances are calculated using equations (Eqns (6) and (7))

(6)$${\rm m}{\rm b}_{{\rm cor}\_2021} = {\rm m}{\rm b}_{{\rm mod}, 30-09-2021}-{\rm \;m}{\rm b}_{{\rm mod}, {\rm survey}\_{\rm date}\_2021}$$

(7)$${\rm m}{\rm b}_{{\rm cor}\_2022} = {\rm m}{\rm b}_{{\rm mod}, 30-09-2022}-{\rm \;m}{\rm b}_{{\rm mod}, {\rm survey}\_{\rm date}\_2022}$$

The total correction balances mb_cor for both glaciers are then calculated based on the two annual corrections (Eqn (8))

(8)$${\rm m}{\rm b}_{{\rm cor}} = {\rm m}{\rm b}_{{\rm cor}\_2022}-{\rm m}{\rm b}_{{\rm cor}\_2021}$$

The mass balance is modelled using the previously calibrated simplified energy-balance model presented in Van Tricht and others (2021c) and Van Tricht and Huybrechts (2022). The input data for this model are limited to hourly precipitation and surface air temperature (from the nearby Tien-Shan-Kumtor weather station), as well as theoretical insolation values at the top of the atmosphere and a DEM. The model is run on a 25 × 25 m grid. The total energy flux in the model is determined by the flux of incoming shortwave solar radiation, the albedo of the surface and the atmospheric transmissivity (Nemec and others, 2009). Next to that, two parameters account for the other components of the radiation balance, i.e. the sum of the longwave radiation and the turbulent sensible heat exchange, linearised around the melting point. These two parameters were calibrated based on mass-balance observations (2015–2021) and they are kept fixed in space and time. The observations from the 2021/22 balance season were not used for the model calibration such that the modelled and observed glaciological mass balance are independent for this year. The uncertainty of the modelled annual mass balance was estimated to be ±0.27 m w.e. a⁻¹ for the Sary-Tor glacier and ±0.28 m w.e. a⁻¹ for the Bordu glacier (Van Tricht and others, 2021c).

For the workflow in this study, the geodetic mass balance is first compared to the independently modelled mass balance for precisely overlapping time periods. Subsequently, the geodetic mass balance of both glaciers is temporally corrected by adding the correction balances. This adjustment enables a meaningful comparison between the geodetic mass balance and the glaciological mass balance (Geissler and others, 2021).

4. Results

4.1. Elevation changes

Between the UAV survey dates, the Bordu glacier and Sary-Tor glacier lost a significant amount of ice, firn and snow, as indicated by the negative elevation changes (Fig. 4). These elevation changes reach values of up to −4 m in the frontal area of the Bordu glacier, and −3.5 m at the front of the Sary-Tor glacier (Fig. 4). With increasing altitude, the elevation changes approach zero, but also the heterogeneity increases. The latter is mainly a result of moving crevasses and snowdrift. The lower areas of the glaciers are mainly characterised by flatter and uncrevassed terrain which results in more uniform elevation changes. The large area with an elevation gain to the east of the Sary-Tor glacier tongue is due to the construction of a road by the Kumtor Gold Company. In the other non-glacierised areas, the elevation changes vary around 0 m, indicating a close correspondence between the pairs of DEMs for both glaciers, as expected for stable terrain. Part of the terrain considered as stable (purple rectangles) is used for the uncertainty analysis. In 2022, small parts of the upper reaches could not be surveyed (green areas on Fig. 4). Consequently, no elevation changes are available for these specific areas.

Fig. 4. Elevation differences between 2022 and 2021 for both glaciers. The background image is a Sentinel-2 satellite image from 26 July 2021. The black outline is the glacier outline in 2022. The green areas were not surveyed in 2022. The purple rectangles correspond to the part of the terrain outside the glacier outlines, considered as stable terrain for the uncertainty analysis.

4.2 Accuracy assessment

The accuracy of the elevation changes is assessed by examining the elevation changes of stable areas over the two consecutive years (Fig. 5). The std dev. of the elevation change within these stable areas (0.25 m for Bordu glacier and 0.22 m for Sary-Tor glacier) is used as the random uncertainty $( \sigma _{\Delta h_{xy, {\rm random}}})$. Additionally, the mean difference in elevation (0.05 m for Bordu glacier and 0.08 m for Sary-Tor glacier) is considered to reflect the systematic uncertainty $( \sigma _{\Delta h_{{\rm systematic}}})$.

Fig. 5. Histogram of the elevation differences between 2022 and 2021 for the areas considered as stable terrain.

By comparing the elevation of the GVPs with the obtained elevation of the DEMs (Fig. 1), we find an RMSE for 2021 of 0.04 m horizontally and 0.11 m vertically for the Sary-Tor glacier and 0.08 m horizontally and 0.09 m vertically for the Bordu glacier. For 2022, we find an RMSE for 2022 of 0.13 m horizontally and 0.16 m vertically for the Sary-Tor glacier and 0.10 m horizontally and 0.15 m vertically for the Bordu glacier. The lower values for 2021 are likely related to the better survey conditions (more contrast, more key points) (Gindraux and others, 2017). In an additional analysis, we find that even by using only three GCPs, the RMSE of the remaining GVPs is of the order of 0.1–0.2 m, close to the RMSE values when using a larger number of GCPs. In addition, no relation was found between the residuals of the GVPs as a function of distance to the closest GCP.

It should be emphasised that we could only place GCPs and GVPs in the ablation area and the lower portion of the accumulation area due to accessibility and time constraints. Consequently, the mentioned RMSE values based on the GVPs primarily reflect the conditions in the lower and flatter regions of the glaciers, while the higher and steeper areas may exhibit greater uncertainties in terms of elevation changes. To elaborate further on this, we conducted an additional analysis to determine the elevation change values in the stable areas as a function of slope and aspect. Through our analysis, we observed a correlation between the std dev. and slope, while no discernible relationship was found with respect to aspect. Notably, including only stable areas with slopes exceeding 30° demonstrated a significant increase in std dev., reaching a maximum value of 0.45 m for the areas with slopes within the range 40–50°. Consequently, we adopted an elevation change uncertainty (random uncertainty, $\sigma _{\Delta h_{xy, {\rm random}}}$) of 0.45 m for locations (x,y) characterised by slopes exceeding 30°.

4.3 Geodetic mass balance

Part of the upper areas of both glaciers could not be reached in 2022 with the UAV due to fresh snowfall (Fig. 4). However, these areas need to be considered as well to calculate the geodetic mass balance. This is often done by extrapolation. Here, Δh values for the unmeasured areas in 2022 (green areas in Fig. 4) are obtained using an extrapolation based on a fit between Δh and surface elevation for the measured areas of both glaciers (Fig. 6). This fit is then used to obtain Δh values for the unsurveyed areas depending on their elevation. To consider the ensuing larger uncertainty for these parts of the glaciers, we doubled the random uncertainty for these parts.

Fig. 6. Surface-elevation change (Δh) as a function of elevation (blue dots for every 10 × 10 m cells) and the best fit (red line) used to derive the elevation change in the unsurveyed areas.

Next, Eqn (2) is used to calculate the overall volume change of both glaciers, resulting in a value of −5.2 × 10⁶ m³ for the Bordu glacier and −2.2 × 10⁶ m³ for the Sary-Tor glacier. Then, using the volume-to-mass conversion value of 850 kg m⁻³, the geodetic mass balance for 2021/22 is calculated, resulting in −0.92 ± 0.09 m w.e. for the Bordu glacier (spanning from 10 August 2021 to 7 August 2022) and −0.87 ± 0.09 m w.e. for the Sary-Tor glacier (covering the period 12 August 2021 to 4 August 2022). The uncertainty range for both values is calculated using the standard values inserted into Eqns (4) and (5).

4.4 Comparison with modelled and observed mass balances

To allow a comparison of the obtained geodetic and the glaciological mass balances, a temporal correction is required. Therefore, the cumulative mean-specific mass balance is modelled between 1 October 2020 and 30 September 2022 with the same parameters from Van Tricht and others (2021c). The mass-balance model is thus not re-calibrated by incorporating the 2021/22 mass-balance data.

The results show a clearly negative evolution over time with a cumulative mass loss of almost −2.5 m w.e. over two balance years (Fig. 7). The modelled mean-specific mass balance for the overlapping period of the geodetic mass balance corresponds to −0.99 m w.e. for the Bordu glacier and −1.02 m w.e. for the Sary-Tor glacier (Fig. 7). This matches rather well with the derived geodetic mass balance (−0.92 ± 0.09 and −0.87 ± 0.09 m w.e., respectively). Next, the temporal correction of the geodetic mass balance is determined using Eqns (6–8) by taking into account the difference in acquisition dates between the measurements of the glaciological mass balance (30 September) and the surveys of the geodetic mass balance (see section 3.5.2), using the model output (Fig. 7). As such, we find mb_cor to correspond to −0.20 m w.e. for the Bordu glacier and −0.22 m w.e. for the Sary-Tor glacier.

Fig. 7. Modelled evolution of the cumulative mean specific mass balance of both glaciers and procedure to derive the temporal correction for the geodetic mass balance.

By incorporating the correction balances mb_cor into the calculated geodetic mass balance, we find a geodetic mass balance of −1.12 ± 0.09 m w.e. a⁻¹ for the Bordu glacier and −1.09 ± 0.09 m w.e. a⁻¹ for the Sary-Tor glacier. These values align relatively well with the mean-specific mass balance of −1.27 ± 0.15 m w.e. a⁻¹ for 2021/22. They are slightly less negative, yet fall within the uncertainty range for both glaciers. There is thus a good correspondence between the glaciological mass balance, derived through interpolation of stake and snow pit measurements, and the geodetic mass balance derived from the UAV observations, which provided coverage of the (almost) entire glacier surface.

5. Uncertainty and sensitivity analysis of the geodetic mass balance

To investigate the influence of employing the conventional volume-to-mass conversion value (Huss, 2013) versus a more detailed approach, we also implement an alternative method for determining the volume-to-mass value (Pelto and others, 2019; Geissler and others, 2021; Beraud and others, 2023). In this approach, a surface density value ρ _xy is assigned to each pixel encompassing the glacier (Eqn (9)). Below the ELA – 50 m, we adopt the density of ice (910 kg m⁻³) (Clarke and others, 2013). Conversely, above the ELA + 50 m, we assume an average density of firn (550 kg m⁻³). In line with the methodology of Geissler and others (2021), an altitude-dependent transition between ice and firn density is employed for intermediate elevations. The ELA, which is assumed to be close to the highest observed snowline, is derived from Sentinel-2 satellite imagery captured at the end of the 2021/22 mass-balance year (4420 m a.s.l. for the Bordu glacier and 4460 m a.s.l. for the Sary-Tor glacier).

(9)$$\rho _{xy}( \rm kg\,m^{{-}3}) = \left\{{\matrix{ {910} \hfill & {{\rm for}\;{\it h} < {\rm ELA}-50} \hfill \cr {550} \hfill & {{\rm for}\;{\it h} > {\rm ELA} + 50} \hfill \cr {730-3.6( {\it h}_{\it xy}- ELA) } \hfill & {{\rm for}\;{\it h}_{\it xy} \,{\rm between\ ELA} \pm 50} \hfill \cr } } \right.$$

Subsequently, the volume-to-mass conversion value is computed by determining a weighted average of the surface densities, with the local elevation change serving as the weighting factor (Eqn (10)).

(10)$$f_{\Delta {\rm V}} = \displaystyle{{\mathop \sum \nolimits_{xy}^{A_{\rm t}} ( {{\rm abs}( {\Delta h_{xy}} ) \times \rho_{xy}} ) } \over {\mathop \sum \nolimits_{xy}^{A_{\rm t}} {\rm abs}( {\Delta h_{xy}} ) }}$$

Using this approach, a volume-to-mass conversion value of 859 kg m⁻³ is obtained for the Bordu glacier and 871 kg m⁻³ for the Sary-Tor glacier. These values are relatively close to the typical value of 850 kg m⁻³ (Huss, 2013). By employing these values to determine the geodetic mass balance, we obtain −0.93 ± 0.09 m w.e. for Bordu glacier and −0.89 ± 0.09 m w.e. for Sary-Tor glacier, which is only slightly more negative than with the standard setup (Table 2).

Table 2. Geodetic mass balance (m w.e.) for both glaciers using different setups

A critical component of the alternative approach lies in the estimation of the ELA. Therefore, we also determine the geodetic mass balance using different ELA values, i.e. once for ELA–100 m and once for ELA + 100 m. Due to the steepness of the glaciers near the ELA, the difference in geodetic mass balance adopting different ELA values is limited, ranging from −0.88 to −0.96 m w.e. for the Bordu glacier and −0.85 and −0.91 m w.e. for the Sary-Tor glacier (Table 2). However, it should be noted that when the ELA is situated in areas characterised by more gentle slopes, uncertainties in the ELA's position could have a more pronounced influence on the results.

Aside from the ELA estimate, uncertainties in the alternative approach are mainly related to the surface density estimate. While it is reasonable to assume that the density of ice remains relatively stable at 910 kg m⁻³ (Clarke and others, 2013), the density of firn density entails greater uncertainties. Therefore, we also compute the geodetic mass balance employing firn densities of 450 and 650 kg m⁻³ in Eqns (9) and (10). Notably, these alternative density assumptions yield only minimal adjustments to the geodetic mass balance (Table 2).

In addition to uncertainties linked to the volume-to-mass conversion value, further uncertainty arises from the elevation changes. As such, we also calculate the geodetic mass balance for both glaciers by excluding the unsurveyed areas. In this case, the geodetic mass balance only slightly decreased with −0.01 m w.e. for both glaciers (Table 2). This can be attributed to the fact that the unsurveyed areas are situated at the uppermost extents of the glaciers and encompass only a minor proportion of the overall area (Fig. 4).

Furthermore, it is important to acknowledge that within the accumulation zone, firn compaction occurs, which commonly leads to a gradual lowering of the base of the yearly snow layer (Belart and others, 2017). This process can potentially result in an underestimation of changes in snow/firn volume. Nevertheless, due to the lack of direct data and the limited accumulation observed on both glaciers (attributable to the arid climatic conditions), we opted not to integrate the consideration of firn compaction into our methodology. This decision introduces an additional uncertainty into our analysis. However, given the minimal accumulation occurring on both glaciers, it is plausible that the impact of firn compaction on the geodetic mass balance is minimal. For glaciers situated in more maritime environments, the role of firn compaction might be more important, necessitating its incorporation.

An additional uncertainty is related to the fact that the UAV surveys and the glaciological measurements were not performed simultaneously. To consider this, we temporally corrected the geodetic mass balance using model output. For a fully independent comparison, it is necessary for both measurements to have been obtained simultaneously. Furthermore, in 2022, part of the glacier surfaces was covered with a fresh layer of snow. Although we assume the effect hereof to be minimal due to the thinness of the snow layer, it did result in a slight increase of the surface and might therefore have led to an underestimation of the volume changes for the balance year of 2021/22. As such, this may constitute a contributing factor to the observed difference between the geodetic mass balance and the mean specific mass balance determined through the glaciological method.

In general, the conducted analysis reveals that for the two studied glaciers, the conducted analysis reveals that for the two studied glaciers, the difference of the calculated geodetic mass balance using the conventional volume-to-mass conversion value and a more complex approach using surface densities and ELA estimates to determine the volume-to-mass conversion value is minimal (Table 2). However, it is important to highlight that this may not always be the case. For example, in years or for glaciers experiencing substantial elevation changes in the accumulation area, the influence of firn density could be more significant. In such cases, differences in volume-to-mass conversion values may be more pronounced, leading to varying estimates of the geodetic mass balance.

When the alternative approach is used, the geodetic mass balance is predominantly influenced by the determined ELA, with only minor sensitivity to the employed firn density. This can mainly be attributed to the fact that the most significant changes in volume occur in the ablation area, while accumulation rates are low. For glaciers characterised by higher accumulation rates such as, for example, Abramov glacier in the Pamir Alay (Kronenberg and others, 2022), uncertainties related to firn density may be more pronounced. Minimising the uncertainty associated with the density estimation can be achieved through in situ measurements of surface densities, particularly in firn areas, as highlighted by Pelto and others (2019). However, in practice this approach is limited by the challenges and logistical constraints of conducting fieldwork in demanding and high-altitude terrains.

6. Conclusions

Geodetic mass balances are of great value to constrain models and validate or calibrate in situ time series of glaciological mass balances. They are increasingly derived using satellite observations over multi-annual timescales. The results of this study highlight the ability of UAV surveys to quantify the annual geodetic mass balance with uncertainties of the order of ±0.09 m w.e a⁻¹. A close concordance between the geodetic mass balance and the glaciological mean-specific mass balance was observed for 2021/22.

The study further emphasises the advantage of using a RTK or PPK UAV system, for which fewer GCPs are required to achieve an accurate high-resolution DEM. This finding is especially important for glaciers or glacier areas that are difficult to access, where placing GCPs is more difficult. In such high-altitude regions, accumulation areas or areas with many crevasses, UAV surveys now provide an alternative to obtain accurate results. However, it should be noted that the approach has certain limitations, including the limited battery capacity, as well as the restricted flight altitude, which necessitates sufficiently high take-off locations. Using a fixed-wing UAV could potentially offer advantages in addressing these challenges.

With our study, we underline that UAVs can clearly be an add-on to classical fieldwork programs to determine the geodetic mass balance of small- to medium-sized glaciers. Our approach can easily be applied to other glaciers. Given the high accuracy of UAV data, it should also be possible to determine the seasonal geodetic mass balance, if surface densities, crucial for estimating the volume-to-mass conversion value, can be prescribed accurately.