Geodetic mass balance of Abramov Glacier from 1975 to 2015

Abstract

Multi-decadal mass loss estimates are available for few glaciers of Central Asia. On Abramov Glacier (Pamir-Alay, Kyrgyzstan), comprehensive long-term glaciological measurements have been carried out from 1968 to 1999 and re-initiated in 2011. A climatological interpretation of this benchmark glacier in Central Asia requires bridging the gap between historical and renewed measurements. This is achieved here by computing the geodetic mass balance from 1975 to 2015 using previously unreleased Soviet aerial imagery and Pléaides stereo-imagery. During 1975–2015, Abramov Glacier lost 2.2 km² (8.2%) of its area. The mean annual thickness change was − 0.43 ± 0.14 m a⁻¹ for the period 1975–2015, corresponding to a volume change of − 0.45 ± 0.15 km³. The average specific geodetic mass balance amounts to − 0.38 ± 0.12 m w.e. a⁻¹. The 1975–2015 glacier mass loss lies within the range of glaciological and geodetic mass-balance estimates that were previously published for disparate and shorter time intervals since 1968. This study covers a much longer time period than earlier geodetic estimates and demonstrates the capacity to geodetically constrain glacier change at high spatial resolution in Central Asia using historic aerial imagery and Structure from Motion techniques. Therefore, it could serve as a benchmark for future studies of regional mass change.

Introduction

Glaciers are important indicators of climate change (Oerlemans, 2001; WGMS, 2008; Roe and others, 2017). In response to global atmospheric warming, glaciers experienced a world-wide retreat and a strong mass loss (Pörtner and others, 2019; Zemp and others, 2019). The current observed mass loss of glaciers in the two large Central Asian mountain ranges, the Tien Shan and the Pamir, is furthermore strongly impacting the runoff regimes of the two large river systems, Amu Darya and Syr Darya (Krasznai, 2019; Pritchard, 2019). Model projections for these rivers show that peak runoff is expected in this region in the middle of the century (Huss and Hock, 2018). Afterwards, a strong decline in water availability is expected until the end of the century. The expected runoff changes would affect people's livelihood downstream from the glacierised mountain ranges as meltwater from the glaciers is a major freshwater contributor to the total river runoff (Armstrong and others, 2019). A reduction of this resource might furthermore contribute to political instabilities (Bernauer and Siegfried, 2012; Krasznai, 2019).

In Central Asia, time series of in situ glaciological measurements are sparse and often discontinuous. Nevertheless, these direct measurements provide an essential source of data for assessing local and regional glacier changes (Hoelzle and others, 2017). Direct glaciological measurements are complemented by remotely sensed geodetic glacier changes, i.e. repeated measurements of glacier surface elevation (e.g. Kronenberg and others, 2016).

Only geodetic glacier mass balances can be derived for entire mountain ranges. Glaciological mass balance generally benefit from cross validation against geodetic estimates. There exists a considerable amount of region-wide to catchment or glacier specific geodetic measurements of glacier changes in Central Asia for the recent two decades (e.g. Gardelle and others, 2013; Brun and others, 2017; Shean and others, 2020). Some of these studies cover periods of several decades (Aizen and others, 2006, 2007; Surazakov and Aizen, 2010; Lambrecht and others, 2014; Holzer and others, 2015; Pieczonka and Bolch, 2015; Goerlich and others, 2017). These geodetic measurements of glaciers rely on declassified spy imagery that is photogrammetrically processed to produce digital elevation models (DEMs). Accuracy in such DEMs, however, generally does not reach the resolution and precision of DEMs from aerial photographs. Studies based on high-resolution and high-quality imagery are thus needed as reference. So far, no long-term geodetic glacier change estimate is available for glaciers located in the Pamir-Alay.

Geodetic surveys of glacier volume and mass change depend on the availability of DEMs. Today, spaceborne data provide a large pool of data for DEM acquisition. Before the onset of the satellite era, DEMs have been derived by means of terrestrial or aerial stereo photogrammetry (Finsterwalder, 1954). In recent years, with advances in 3-D computer vision algorithms and hardware improvements, Structure from Motion (SfM) photogrammetry has gained popularity and found numerous applications in Geosciences (Smith and others, 2015; Carrivick and others, 2016). Although traditional stereo photogrammetry requires metadata information about the calibration of cameras and ground control information, SfM only requires numerous overlapping 2-D images for input which can later be georeferenced using ground control information (Lowe, 2001; Westoby and others, 2012; Micheletti and others, 2015). Various authors have applied ground-based SfM to derive high-resolution DEMs of glaciers and to calculate geodetic volume changes and mass balances (e.g. Piermattei and others, 2015, 2016; Marcer and others, 2017), yet SfM has not extensively been applied to process historical aerial imagery acquired over steep, mountainous terrain and on glaciers. Girod and others (2018) used oblique imagery acquired in 1936 to derive high-resolution and high-quality elevation data over glaciers in Svalbard. Mölg and Bolch (2017) evaluated surface elevation changes over the tongue of Zmuttgletscher in Switzerland on DEMs created from aerial imagery between 1946 and 2005. Mertes and others (2017) derived high-resolution glacier DEMs from archival imagery but did not study glacier volume or mass changes. SfM bears a large potential to calculate DEMs for the early or pre-satellite era. Where the quality and availability of the aerial imagery is good, SfM allows generating DEMs that might serve as reference data sets for other studies.

Here, we used Soviet-era aerial photographs in combination with SfM photogrammetry to produce a 4 m resolution DEM of Abramov Glacier for the year 1975. In combination with recent high-resolution elevation data from Pléiades (Gleyzes and others, 2012; Berthier and others, 2014), we performed a detailed assessment of volume and mass change of Abramov Glacier located in the Pamir-Alay range. The high detail of the aerial photographs allowed for computation of an accurate DEM for a time period where previously no such data existed. We put our results in context with local and regional geodetic mass-balance estimates.

Study area and data

Abramov Glacier

Abramov Glacier is a valley glacier in the Pamir-Alay range in Kyrgyzstan, Central Asia. The glacier drains into the Koksu river and lies in the head watershed of Vakhsh river flowing south-westwards into the Amu Darya River Basin (Fig. 1).

Fig. 1. Overview map of glaciation (turquoise) in Central Asia. The study site, Abramov Glacier, is located in the red box. Source data from Natural Earth (http://www.naturalearthdata.com). Glacier extent source data: RGI version 6.0 (RGI, 2017).

Abramov Glacier features an abundance of glaciological in situ data because a research station was built in 1967 by the Central Asian Research Hydrometeorological Institute. The research station delivered continuous mass-balance observations from 1968 to 1999. The break-up of the former Soviet Union led to a lack of funding and the monitoring programme was interrupted in 1999 when the research station was attacked and destroyed (WGMS, 2001). In 2010/11, the monitoring of several glaciers in Central Asia, including Abramov Glacier, was re-initiated (Hoelzle and others, 2017, 2019).

From 1850 to 2015, a retreat of the glacier tongue of ~2 km has been observed (WGMS, 2017). Barandun and others (2015) report a retreat of the glacier tongue of ~1 km from 1975 to 2015 and an area reduction of 5%. Abramov Glacier advanced by ~400 m between 1972 and 1973. The event was monitored in a very detailed manner and is considered a surge (Emeljyanov and others, 1974; Suslov and others, 1980; Glazirin and others, 1987, 1993; Glazirin and Braun, 2005).

Kuzmichenok and others (1992) report an area of 23.3 km² for the year 1986. They report a mean ice thickness of 110 m and a maximum ice thickness of 247 m derived by radio echo soundings. They estimated a total glacier volume of 2.57 km³.

Barandun and others (2015) re-analysed the mass-balance data for the years 1968–98 based on field measurement data and modelled the mass balance for the period of missing field data between 1995 and 2011. Combining modelled and re-analysed data, they reported a mean specific mass balance of − 0.44 ± 0.10 m w.e. a ⁻¹ for 1968–2014 (Barandun and others, 2015). Barandun and others (2018) report a mean specific geodetic mass balance of $-0.39 \pm 0.16{\, \rm \, m\, w.e.\, a^{-1}}$ for 2003–2015.

Aerial imagery

Historical aerial photographs acquired during the time of the Soviet Union were available for this study (Mandychev and others, 2017). The aerial imagery was acquired with a Soviet camera AFA-TES system developed in the late 1940s. The camera system featured a Russar 63V lens with aperture values ranging from f/6.8 to f/16 and shutter speeds ranging from 1/70 s to 1/850 s (Iljin, 1969). The focal length was 101.39 mm. The flying altitude was 3800 m above ground and results in a scale of 1:32 000. The images are panchromatic and overall of good quality. We did not acquire the original aerial images but instead obtained contact prints of size 18 cm × 18 cm. For the year 1975, the whole glacier area is covered by the flight tracks. The along-track overlap is ~70% and cross-track overlap amounts to ~45%. For 1975, a total of 23 aerial photographs were used for analysis. The target imagery was acquired in mid July 1975 and not at the end of the ablation season; some snow is present on the glacier surface. The upper part of the accumulation area features low texture details and low contrast which can challenge SfM processing (Micheletti and others, 2015). Due to their age the prints show moderate signs of wear and tear such as artificial scratches but the glacier surface is not affected.

Satellite stereo images and Pléiades DEM

A total of three pairs of Pléiades stereo images (Berthier and others, 2014) were acquired over Abramov Glacier on 20th August 2015 and 1st September 2015 (Fig. 2). From this, two DEMs with a spatial resolution of 4 m were created using the open-source software Ames Stereo Pipeline (Shean and others, 2016). The panchromatic and the four multi-spectral channels of the Pléiades satellite were used to compute the DEM. Although they have a coarser resolution, the multi-spectral channels were added in the processing as they were less saturated than the panchromatic channel. This approach reduced the area of voids in the accumulation area. This DEM, DEM ₂₀₁₅, was used as a reference. Data voids of the DEM ₂₀₁₅ amount to a total of 0.64 km ² (2.6% of the total glacier area) and are mainly found in the accumulation area. The DEM ₂₀₁₅ was generated in WGS84 and set to UTM zone 42 North. Altitudes are defined above the EGM96 geoid. The standard deviation over stable terrain can be used as a measure of accuracy if a reference DEM is available (Fisher and Tate, 2006).

Fig. 2. Overview map of Abramov Glacier with footprints of the 2015 Pléiades stereo-pairs (red and orange dotted lines) and aerial imagery coverage of 1975 (blue dotted line). Glacier area derived from RGI 6.0 are illustrated in turquoise. The SRTM DEM (Farr and others, 2007) serves as shaded relief.

We used the DEM and orthoimagery of 1st September 2015 in the subsequent analyses, despite a snow-covered tongue. Fresh snow was not corrected for, even though present on the reference DEM ₂₀₁₅. The snow cover on DEM ₂₀₁₅ was not considered in the uncertainty calculation as it was small (i.e. an estimated 10–30 cm of snow and too small to be quantified by comparing the two Pléiades DEMs from directly before and after snowfall) and its influence on geodetic mass balance is negligible (estimated to be <$1{\, \rm \, cm\, w.e.\, a^{-1}}$). This DEM of 1st September 2015 was already used in Barandun and others (2018) to calculate geodetic mass balance of Abramov Glacier from 2003 to 2015. The Pléiades stereo pair of 20th August 2015 does not cover the entire glacier area of Abramov Glacier, but was used for glacier outline delineation.

Ground control measurements

Ground control measurements on ten reference points located on crests and summits around the glacier were performed during the field seasons of 2014 and 2015. A Topcon dGPS (Topcon GB-1000 receiver with PG-A1 antenna) was used in 2014 and measurements in 2015 were performed with a Trimble device (Trimble 5700 GPS and Trimble R8 GNSS receiver). Furthermore, during the acquisition phase of the Pléiades stereo pairs in late August 2015, 7649 GPS track points were recorded on the glacier. Location measurements from a permanently installed GPS receiver (~1.5 km from the glacier tongue) were used for differential correction of the measured GPS points of both years. The precision of the measured points is 2–3 cm.

Methods

The Abramov archival imagery was processed using the open-source photogrammetric pipeline MicMac (Rupnik and others, 2017). Our SfM workflow is visualised in Figure 3 and follows the recommendations by Girod and others (2018). Using MicMac has the main advantage of having control over intermediary steps and quality parameters in the processing chain which commercial SfM software often does not permit (Jaud and others, 2016).

Fig. 3. Workflow of SfM processing pipeline applied in this study, RM stands for reliability mask. Outputs are illustrated in turquoise.

Scanning imagery and internal orientation

The contact prints of size 18 cm × 18 cm were scanned at varying resolutions of 11 μm (2400 dpi) to 21 μm (1200 dpi) as 16 bit TIFF images for the purpose of SfM processing. The scanner model was a commercial scanner model (Canon MG7700 series) which resulted in a ground sampling distance of 89 cm per pixel. All aerial imagery was visually checked for apparent scanning errors such as blurring induced by not equally pressing the image on the scanning plate or different levels of contrast due to faulty scanning properties. Images were pre-processed in Adobe Photoshop. All four corner grid markers were highlighted with a small identification mark. In case the corner grid markers were located in very dark or very bright parts of the imagery and not clearly identifiable, we drew lines through the reseau crosses and reconstructed the corner grid mark position (Fig. 4).

Fig. 4. Example of aerial image pre-processing: reseau crosses on the aerial imagery are highlighted in red. The corner grid markers are marked manually (white dot), however the corner grid marker on the lower right was difficult to place because of shadows in the imagery. Inset b shows a zoom of the lower right corner of a. Inset c shows the applied reconstruction method by drawing two lines (orange) through the reseau crosses which allowed us to find the position of the corner grid mark.

The scanning process can result in slightly different image dimensions and introduce distortions, and thus the interior orientation needs to be standardised (Girod and others, 2018). As no calibration data of the corner grid marker position was available, their relative positions were measured manually and identified on each image. Iljin (1969) specified the distance in between the reseau crosses on the film as 10 mm. The images were then transformed and resampled to a virtual resolution in a standardised geometry using the MicMac command ReSampFid.

Relative orientation

Tie points on the overlapping images with different acquisition angle were identified in a next step. A total of 253 638 tie points were found for all images using the MicMac command Tapioca. The bundle adjustment estimates the camera location, intrinsic and extrinsic parameters of the camera and outputs a sparse point cloud. This step involves the actual process referred to as SfM: non-linear cost functions are solved to optimise the bundle of projection rays that connect the cameras with the 3-D points and estimate the orientation and internal parameter of the camera (Lowe, 2001).

The known focal length of 101.39 mm and the image size of 18 cm × 18 cm served as an additional input in the calibration estimation. The calibration using the radial lens model Fraser Basic (cf. Pierrot-Deseilligny, 2015) was initiated by using the MicMac command Tapas and the focal length was set to a fixed value in order to facilitate a stable calibration. Tapas obtained a residual error of 1.5 pixels.

Absolute orientation and creation of DEM and orthoimage

The orientation is transformed from an arbitrary one into an absolute reference system by georeferencing. The use of ground control points (GCPs) ensures an optimal absolute orientation. GPS field measurements from the end of ablation season 2014 and 2015 were mostly located in very steep terrain or in low contrast areas and not easily identified due to lower resolution of the 1975 imagery. Therefore, it was decided to use the Pléiades orthoimage and DEM as a reference to collect GCPs (Papasodoro and others, 2015; Belart and others, 2019). The criteria for the GCP selection was stability over time and an homogeneous distribution over the area of interest (James and Robson, 2012). Furthermore, clear identification needs to be possible on both the reference and target imagery. To limit the uncertainty of GCP placement, large round and flat boulders and bedrock features of sizes ranging from 5 to 10 m in diameter were favoured. A set of 18 GCPs was manually identified and used to perform the absolute orientation (Fig. 5). The uncertainty of GCP placement was estimated to be 3 pixels.

Fig. 5. Locations of GCPs for georeferencing. The small insets on the right show close ups of two GCPs and the corresponding orthoimagery. The GCP location on the aerial imagery 1975 is shown in a and c and the GCP placement of 2015 is shown in b and d. A shaded relief image of the Pléiades DEM from 2015 and the SRTM DEM (Farr and others, 2007) serve as background.

Using these GCPs, the model of the 3-D geometry is refined by rerunning the bundle adjustment using the MicMac tool Campari. A dense 3-D point cloud is created using semi-global image matching algorithms with the MicMac command Malt Ortho. Malt Ortho is optimised for matching the images into a regular grid in the ground geometry of the orthoimages (Rupnik and others, 2017). The image matching was performed using an 11 × 11 pixel moving window. From the image matching a DEM with a resolution of 4 m was produced. Using the MicMac module Tawny, an orthophoto mosaic with a resolution of 0.5 m was generated. The DEM and corresponding orthoimages were projected into the same reference system as the reference dataset from 2015.

Reliability map

The SfM chain in MicMac produces a correlation mask which is a measure for the quality of the terrain extraction process. The values of the correlation mask were normalised to values between 0 and 100 and by setting an empirical threshold of 40 (based on visual interpretation of the DEM and orthoimagery), values over 40 were classified with the label high confidence. Correlation mask values below 40 were classified as low confidence.

Determination of stable terrain and co-registration

Before performing a pixel-by-pixel subtraction between the two DEMs, potential systematic errors introduced by the acquisition or processing technique need to be corrected accordingly (Rolstad and others, 2009; Zemp and others, 2013; Paul and others, 2015). The DEMs were checked for horizontal and vertical offsets. This analysis was performed over stable ground.

Abramov Glacier is surrounded mostly by steep and potentially moving terrain and the demarcation of stable terrain proved to be difficult. Terrain with slope angle lower than 20^° was extracted in the 2015 reference DEM. The glacier areas derived from the Randolph Glacier Inventory (RGI) version 6.0 (RGI, 2017) were excluded. The remaining area was classified as non-moving terrain. The final stable terrain used for co-registration covers an area of 14.3 km ².

Here, we applied the well-established co-registration approach of Nuth and Kääb (2011). It provides a robust solution for co-registration also in cases when stable terrain is limited to <10% of the scene (Nuth and Kääb, 2011). The 2015 DEM served as a reference, therefore it was important to assess its absolute accuracy. As expected, due to the use of GCPs, no large shift was observed by applying this co-registration method. We shifted the DEM and corresponding orthoimagery by the vectors of 0.1 m in the x direction, 6.0 m in the y direction and 0.0 m in the z direction. The standard deviation of vertical differences across stable terrain between the 1975 and the 2015 DEMs after co-registration was 43.2 m.

Delineation of glacier outlines

Glacier outlines were digitised based on the orthophotos of 1975 and 2015. Glacier outlines based on Landsat imagery from Barandun and others (2018) served as an additional cross-reference in the digitising process. The Pléiades stereo-pairs of 20th August 2015 exhibited less fresh snow and were used to delineate glacier outlines. The area derived from Landsat outlines (Barandun and others, 2018) differs from our digitised outlines by 4.7% (1975) to 1.4% (2015) respectively, which is likely a conservative estimate of the uncertainty of our glacier delineation as we worked with high-resolution imagery. Our estimate agrees well with previously published area uncertainty estimates for high-resolution datasets (Paul and others, 2013).

Accuracy assessment of DEM ₂₀₁₅

The horizontal and vertical accuracy of the reference DEM ₂₀₁₅ was assessed by comparing dGPS field measurements in 2014 and 2015 against the horizontal and vertical position of the DEM ₂₀₁₅ off-glacier and on-glacier.

For a total of ten dGPS points off-glacier, the mean of the differences in the x direction (0.51 ± 0.99 m (1σ)), y direction (− 1.24 ± 0.72 m (1σ)) and z direction (0.03 ± 4.03 m (1σ)) were calculated. Additionally, a total of 7649 GPS track points acquired during the field campaign from 24th August to 27th August 2015 were analysed on-glacier in order to get a comprehensive assessment of the vertical DEM accuracy. A mean of difference in elevation of 3.53 ± 1.79 m (1σ)) was found for all points lying on the glacier surface of Abramov Glacier. The data points in the upper part of the accumulation area and near DEM data gaps indicate higher differences between GPS and DEM elevations.

Calculation of thickness change

The SfM-derived DEM ₁₉₇₅ was subtracted from the reference DEM ₂₀₁₅ on a pixel basis in ArcGIS to obtain the thickness change dh map. The values of thickness change were calculated over the glacier surface and divided by the number of years to obtain a mean annual thickness change.

Post-processing: filtering and void filling

Post-processing the dh map was necessary as the image matching process introduced noise and data voids, mostly due to low texture quality in the accumulation area (Fig. 6d). We filtered outliers by binning the thickness change values into 50 m elevation bins and calculating statistical parameters such as mean thickness change, standard deviation and median for each bin. In areas classified as low confidence we excluded values deviating more than one standard deviation from the median in each 50 m bin (Fig. 6e). We then filled these excluded areas with the median calculated for each 50 m elevation bin (Fig. 6f). Median filling was preferred over filling data voids with a mean value, because the median is less susceptible to outliers than the mean (McNabb and others, 2019).

Fig. 6. Data products of DEM ₁₉₇₅: (a) Orthophoto DEM ₁₉₇₅, (b) Hillshade DEM ₁₉₇₅, (c) reliability map DEM ₁₉₇₅, (d) raw elevation change, (e) filtered elevation change and (f) void-filled elevation change. RGI 6.0 outlines are illustrated in blue. The SRTM DEM (Farr and others, 2007) serves as shaded relief outside the area of interest.

DEM ₁₉₇₅

The final resolution of the SfM-derived DEM ₁₉₇₅ is 4 m. The on-glacier elevation values are distributed between 3606 and 4969 m a.s.l. Figure 6a–c show the different data products that resulted from the SfM processing within MicMac. The data coverage of the DEM ₁₉₇₅ is 97.4% over the glacier surface. The data voids (2.6%) in the reconstructed DEM ₁₉₇₅ are mainly found in the southeastern tributary of the accumulation area of Abramov Glacier above 4250 m a.s.l. (Fig. 6). Based on the correlation score image in MicMac and the resulting reliability mask pictured in Figure 6c, 44.9% of the glacier pixels were classified as high confidence data, 52.5% as low confidence data and 2.6% are no data values. Low confidence data are predominantly found in the accumulation area (Fig. 6c).

Calculation of volume change and mass balance

Total volume change ΔV in m³ was calculated by using Eqn (1):

(1)$$\Delta V = r^2 \sum\limits_{k = 1}^K \Delta h_{k}\comma\; $$

where K is the number of pixels covering the glacier at the maximum extent, Δh _k is the difference in elevation of the two DEMs at pixel k and r is the pixel size (4 m). The maximum extent in glacier area (1975 for observation period 1975 to 2015) was chosen in order to capture all pixels that were affected by a change in elevation (Koblet and others, 2010).

The volume change for the period 1975–2015 was converted into a specific mass balance using Eqn (2):

(2)$$B_{\rm geod} = {\Delta V\over \bar{A}} {\,f \Delta V\over \rho_{\rm water}}$$

where $\bar {A}$ is the averaged glacier area between the acquisition dates, fΔV is a density conversion factor of 850 kg m ⁻³ adopted from Huss (2013) and ρ_water is the density of water (1000 kg m ⁻³). An annual rate of geodetic mass balance B _{a geod} specified in m w.e. a⁻¹ was further calculated by dividing B _geod by the number of years n = 40.

Uncertainty analysis

Uncertainties in volume change and mass balance arise from DEM generation such as data voids in low contrast or snow-covered areas (Kääb and Funk, 1999; Bolch and others, 2019a; McNabb and others, 2019).

Random and systematic errors can be introduced during the geodetic calculation of glacier mass changes (Zemp and others, 2013). Here, we estimated only the random error. To assess the uncertainty in the DEMs and errors introduced in the subsequent volume and mass-balance calculation, the uncertainty approach of Fischer and others (2015) was adapted.

A _cor in m² describes a circular area in which errors introduced by the DEM differencing are spatially correlated and is calculated with Eqn (3) adapted from Rolstad and others (2009):

(3)$$A_{\rm cor} = \pi L^2\comma\; $$

where L describes the spatial correlation distance in m derived by semivariogram analysis on stable terrain. Within this range, a spatial autocorrelation between two points is established (Rolstad and others, 2009). A spatial correlation length of L = 827 ± 68 m (1σ) was obtained for the observation period 1975–2015. The calculated value for A _cor is 2.15 km ² for 1975–2015.

The uncertainty in surface elevation change on the scale of individual pixels is calculated with Eqn (4) as the correlated area is smaller than the maximum glacier area (Rolstad and others, 2009; Fischer and others, 2015):

(4)$$\sigma _{\Delta z} = \sqrt{\sigma^2 _{\Delta h} {A_{\rm cor}\over 5 A_{\rm max}}}\comma\; $$

where σ_Δh in m is defined as the standard deviation over stable terrain (43.2 m). A _cor is the correlated area described in Eqn (3) and A _max the glacier area at its maximum extent of the observation period (in our case 1975).

The uncertainty in volume change σ_ΔV in m³ is assessed with Eqn (5):

(5)$$\sigma _{\Delta V} = \sigma _{\Delta z} A_{\rm max}\comma\; $$

The uncertainty of the geodetic mass balance σ_geod in m w.e. is calculated with Eqn (6):

(6)$$\sigma _{\rm geod} = \sqrt{\lpar \Delta z\, \sigma _{\,f \Delta V}\rpar ^2 + \lpar f _{\Delta V}\, \sigma _{\Delta _z}\rpar ^2}\comma\; $$

where Δz specifies the mean elevation change in meter, fΔV is a density conversion factor of 850 kg m⁻³ and σ_fΔV the corresponding uncertainty in material density of 60 kg m⁻³ (Huss, 2013). σ_Δz is the uncertainty in surface elevation change derived with Eqn (4).

Results

Changes in glacier area and length

The total glacier area in 1975 was 26.53 km² of which 1.23 km² (2.6%) are data voids in the DEM ₁₉₇₅. The total glacier area in 2015 was 24.36 km². From 1975 to 2015, Abramov Glacier lost 2.17 km² ($-8.2\percnt$) of its surface area. The loss in area was most pronounced at the glacier tongue. The front experienced a retreat of ~1040 m between 1975 and 2015. The relatively small northeastern tributary disconnected from the glacier between 1975 and 2015 and receded ~100 m from 1975 to 2015.

Change in elevation, volume and mass balance from 1975 to 2015

From 1975 to 2015, a mean surface lowering of − 17.03 ± 5.50 m occurred on Abramov Glacier (Fig. 6d). Owing to the ~1 km front position retreat, the lowering is strongest at 3624 m a.s.l. on the tongue, about 1.1 km upstream of the glacier front in 1975. The tongue of Abramov Glacier shows maximum thinning of up to 151.32 ± 5.50 m. A high level of noise (standard deviation outside of glacier area is 81.8 m) is visible in the terrain surrounding the glacier (remarkably in higher regions at the Southern borders of the DEMs, Fig. 6d).

The volume change for 1975–2015 amounts to − 0.45 ± 0.15 km ³. This corresponds to 18% of the glacier volume determined by Kuzmichenok and others (1992) in 1986. The calculated mean annual geodetic mass-balance rate for the years 1975 to 2015 is $-0.38 \pm 0.12{\, \rm \, m\, w.e.\, a^{-1}}$.

Figure 6e, f illustrates the result of the void filling approach. In total, 14.8% of the total glacier area was interpolated using the median fill method. Figure 6f shows that an elevation change of close to zero occurred in the accumulation area. Relatively few locations had positive elevation changes (9.8% of the glacier area in 1975; cf. Fig. 6f), generally situated at the uppermost, very steep sections of the glacier. Steep terrain in the accumulation area with low contrast is often associated to high uncertainty in DEMs derived from optical images (Barandun and others, 2018). Mean elevation change in the accumulation area is close to zero, frequent deviations from the mean, including positive values, are related to high levels of noise.

Discussion

Different studies calculated mass balance of Abramov Glacier, however calculations over different time periods challenge a direct comparison (Fig. 7 and Table 1). Only the glaciological measurements exceed the time frame of our estimate, albeit with a large data gap of 16 years. Geodetic mass-balance estimates are only published for the recent period (after the year 2000).

Fig. 7. Mean specific mass balance (geodetic and glaciological/modelled) results for different studies and different time periods for Abramov Glacier. The lines illustrate the mean specific annual mass-balance estimates and boxes display the uncertainty range for the respective publications.

Table 1. Summary of existing glaciological and geodetic mass-balance estimates for Abramov Glacier

Mass-balance estimates and corresponding uncertainty are specified in m w.e. a⁻¹.

Comparison to other geodetic mass-balance estimates of Abramov Glacier

Brun and others (2017) calculated geodetic mass balance of $-0.38 \pm 0.10{\, \rm \, m\, w.e.\, a^{-1}}$ for the period 2002–2014 using DEMs derived from multi-temporal Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) imagery. Their average annual mass balance is identical to our study ($-0.38 \pm 0.12{\, \rm \, m\, w.e.\, a^{-1}}$). Based on high-resolution Pléiades data (same data source as the reference DEM ₂₀₁₅ used in this study) and Satellite Pour l'Observation de la Terre (SPOT)5 DEMs from 2003 and 2011, Barandun and others (2018) calculated an average geodetic mass balance of $-0.39 \pm 0.16{\, \rm \, m\, w.e.\, a^{-1}}$ for the period 2003 to 2015 and $-0.36 \pm 0.26{\, \rm \, m\, w.e.\, a^{-1}}$ for the period 2011–2015. The same study reconstructed a slightly lower mass loss of $-0.26 \pm 0.20{\, \rm \, m\, w.e.\, a^{-1}}$ from 1998 to 2016 combining geodetic measurements, numerical modelling, snowline observations and glaciological observations. The dh map in Barandun and others (2018) shows a pronounced surface lowering mainly on the tongue and almost no change in the accumulation area from 2003 to 2015. Our results show that surface elevation change for the period 1975–2015 in the accumulation area are close to zero. Shean and others (2020) reprocessed thousands of ASTER and very high-resolution imagery (mostly Worldview) DEMs over High Mountain Asia (HMA) from 2000 to 2018. For Abramov Glacier they report a specific annual mass balance of − 0.26 ± 0.12 m w.e. a ⁻¹ from 2000 to 2018. The studies of Brun and others (2017) and Shean and others (2020) used stacks of multiple DEMs to derive geodetic mass balance, thus their time stamp is not clearly defined and their results are difficult to compare with a geodetic mass balance derived from only two DEMs (Barandun and others, 2018). Gardelle and others (2013) suggested a mean annual geodetic mass balance of − 0.03 ± 0.14 m w.e. a ⁻¹ for Abramov Glacier by differencing SPOT5 imagery from 2011 and SRTM data from 2000. The dh dataset of Gardelle and others (2013) shows significant positive changes of elevation in the accumulation area of Abramov Glacier and at the front of the tongue. Barandun and others (2018) question the findings of Gardelle and others (2013) and pointed out a possible positive bias in elevation change because of an underestimated penetration of the SRTM X-band radar signal into snow and ice (Dehecq and others, 2016; Lambrecht and others, 2018).

Overall, our mass balance 1975–2015 is in the same range as mass-balance estimates published for a ~15 years subset (~2000 to 2015: Brun and others, 2017; Barandun and others, 2018; Shean and others, 2020) of our 40-year time frame (Fig. 7). Thus, we cautiously conclude that mass balance for the two time periods 1975 to ~2000 and ~2000 to 2015 are similar. This would indicate that the rate of mass loss of Abramov Glacier was relatively stable. This assumption is partly confirmed by re-analysis of glaciological measurements (see the following subsection and Barandun and others, 2015) which indicates no general long-term trend in mass balance over the time frame 1968–2014. At this point we can only speculate to why Abramov's rate of mass loss appears to be relatively stable. Kronenberg and others (in press) have repeated Soviet firn studies in the accumulation area of Abramov Glacier. The results suggest that accumulation has increased since the 1970s, potentially offsetting increased melt.

Comparison to glaciological and modelled mass-balance estimates

Glaciological mass-balance measurements generally benefit from cross validation against geodetic estimates. Barandun and others (2015) (re)analysed the glaciological mass balance for the period 1968–1995, 2012–2014 and reconstructed the mass balance from 1995 to 2011. Their mass balance of − 0.44 ± 0.10 m w.e. a ⁻¹ from 1968 to 2014 lies in the same range and is valid for a similar time frame as in our study. Re-analysed glaciological measurements (Barandun and others, 2015) in combination with improved modelling (Barandun and others, 2018) provide a mass balance of − 0.37 ± 0.17 m w.e. a ⁻¹ from 1975 to 2015. The 1975–2015 area change of 2.17 km ² and glacier retreat of ~1040 m are in line with the findings by Barandun and others (2015).

Comparison to regional glacier changes

Considering longer time frames, similar to the one of the current study, Holzer and others (2015) found a nearly balanced mass-balance budget for Muztagh Ata Glacier basin in the eastern Pamir using declassified spy imagery for the years 1975 to 2013. Using declassified satellite imagery, Zhou and others (2019) reported an almost balanced mass budget (− 0.03 ± 0.24 m w.e. a ⁻¹) between 1975 and 1999 for Central Pamir. Bolch and others (2019b) reported an averaged mass loss of − 0.26 m w.e. a ⁻¹ for the Western Pamir. Barandun (2018) reported an averaged mass loss of − 0.37 ± 0.42 m w.e. a ⁻¹ for the Western Pamir. The results of Pieczonka and others (2013) suggest that the overall mass budget of − 0.33 ± 0.15 m w.e. a ⁻¹ for 1976–2009 for the glaciers in the Aksu-Tarim catchment further north in the Central Tien Shan lies in the same range as here calculated for Abramov. Our results, together with those from the literature, draw a picture of moderate multi-decadal mass loss (mass balance of typically − 0.2 m to − 0.4 m w.e. a ⁻¹) in the Western Pamir, Pamir Alay and Central Tien Shan. The literature furthermore indicates that over a similar time period glacier mass budgets in Central and Eastern Pamir were almost balanced, despite heterogeneous response for individual glaciers within each region (Scherler and others, 2011; Farinotti and others, 2015; Kääb and others, 2015; Brun and others, 2017; Wang and others, 2017; Shean and others, 2020).

Conclusions

In this study, we derived a DEM and outlines of Abramov Glacier for the year 1975 based on Soviet aerial images. A high-resolution Pléiades DEM acquired in September 2015 served as a reference for the GCPs and for the calculation of area, volume and mass changes between 1975 and 2015. For 1975–2015, there was a reduction of $-8.2\percnt$ in glacier area and a mean annual thickness change of − 0.43 ± 0.14 m a ⁻¹. We calculated a mean specific mass balance of − 0.38 ± 0.12 m w.e. a ⁻¹. This mass balance is in the same range as glaciological and geodetic mass-balance values published in recent studies for the time period ~2000 to 2015. Calculated mass balance is very similar to the re-analysed glaciological mass balance 1968–2014 of the glacier. In the future, the historic 1975 DEM that we generated from declassified spy imagery could serve as an input for a time series of Abramov Glacier combining all available DEMs, modelled and glaciological mass-balance values. Our results illustrate that SfM is a useful approach for processing historical imagery from Soviet times. Our high-resolution map of elevation change could also help to validate similar estimates using spy satellite imagery.