Geodetic glacier mass balance (1975–1999) in the central Pamir using the SRTM DEM and KH-9 imagery

Abstract

Multiple studies on regional glacier mass balance in the Pamir Mountains have been conducted using the geodetic method, but they were rarely extended to the period before 2000. In this study, we used KH-9 imagery acquired in 1975 to generate the historical DEM for the central Pamir, and then obtained the glacier elevation change by comparing this with the SRTM C-band DEM. The penetration depth of the C-band radar was corrected for different glacier surfaces, i.e. 2.96, 1.68 and 0 m for firn/snow cover, bare ice and debris-covered areas, respectively. The final results suggest that the central Pamir glaciers, overall, experienced a near-zero mass balance of −0.03 ± 0.24 m w.e. a⁻¹ for 1975–99. Due to glacier surge activity, the elevation change patterns of individual glaciers were highly variable, and their mass balances varied from −0.12 ± 0.26 to 0.63 ± 0.20 m w.e. a⁻¹. The mean mass budgets of the surge-type glaciers and non-surge-type glaciers were 0.03 ± 0.14 and −0.05 ± 0.28 m w.e. a⁻¹, respectively. Concurring with previous studies, we conclude that the central Pamir glaciers may have been in a state of approximately balanced mass budget or slight mass deficit from the mid-1970s to the mid-2010s.

1. INTRODUCTION

Glacier mass balance, as one of the most direct and reliable indicators of climate change, can quantitatively reveal whether there has been mass gain or loss in a glacier system (Bolch and others, 2011; Gardelle and others, 2013). Most of the previous studies have suggested that glaciers around the world have experienced obvious mass loss since the 1960s or 1970s (Berthier and others, 2010; James and others, 2012; Fischer and others, 2015; Zhou and others, 2018b), and that this condition was aggravated after the 1990s or 2000s, especially for high mountain regions in Asia (Bolch and others, 2012; Yao and others, 2012; Gardner and others, 2013). Meltwater from glaciers in High Asia has not only contributed to a sea-level rise of 0.041 ± 0.009 mm a⁻¹ for the period of 2000–16 (Brun and others, 2017), but, more importantly, it plays a crucial role in adjusting river runoff, while a proportion of the meltwater is temporarily stored in proglacial lakes (Gardelle and others, 2013; Kääb and others, 2015; Brun and others, 2017). In addition, drastic glacier changes can induce geological hazards (such as landslides, debris flows and glacier lake outburst floods), bringing potential security risks for downstream residents (Mergili and others, 2013; Komatsu and Watanabe, 2014; Lamsal and others, 2016). For the Pamir Mountains, in particular, which hold a considerable number of alpine glaciers, the question of how glaciers have changed in this region is significant and deserves to be of concern.

Over the entire Pamir, field-based measurements of mass balance (i.e. the glaciological method) have mainly concentrated on several small glaciers that are easy to reach, revealing a highly variable pattern of mass change in both space and time. For example, it has been reported that the Muztagh Ata Glacier in the eastern Pamir experienced a positive mass balance of 0.25 m w.e. a⁻¹ for 2005–10 (Yao and others, 2012), whereas the Abramov Glacier in the Pamir Alay showed significant mass loss at a rate of −0.46 ± 0.06 m w.e. a⁻¹ for 1971–94 (Barandun and others, 2015). Since the field-based measurement technique is time-consuming and labor-intensive, it is not applicable to glacier change detection at a relatively large spatial scale, despite its advantage of high accuracy. As a complementary approach to investigate glacier mass balance, the geodetic method, which is primarily based on the differencing of multi-epoch topography products, has been widely and successfully used (Bolch and others, 2008; Kääb and others, 2012, 2015; Gardelle and others, 2013; Li and others, 2017). Geodetic-based studies have suggested that the Pamir glaciers, as a whole, have generally been in a state of mass loss from 2000 to the mid-2010s (Table 1). For example, the use of ICESat laser altimetry products revealed a surface lowering of −0.13 ± 0.22 m a⁻¹ for 2003–09 (Gardner and others, 2013) and −0.48 ± 0.14 m a⁻¹ for 2003–08 (Kääb and others, 2015). The discrepancy between these results is mainly due to the fact that the surface elevation of the Pamir glaciers is highly variable at the inter-annual scale (Brun and others, 2017). Based on the ASTER DEMs, updated results have suggested that there was only slight mass loss (−0.04 ± 0.07 m w.e. a⁻¹) for 2000–16 (Brun and others, 2017). On a local scale, although a negative mass balance of −0.15 ± 0.12 m w.e. a⁻¹ for 1971/76–2013/14 has been reported in the eastern Pamir (Zhang and others, 2016), most studies have suggested that the glaciers within this region have been in a balanced or positive state of mass budget from the 1970s to the mid-2010s (Table 1), e.g. −0.01 ± 0.30 m w.e. a⁻¹ for 1973–2013 in Holzer and others (2015) and 0.12 ± 0.07 m w.e. a⁻¹ for 1999–2014 in Lin and others (2017). The difference in the penetration depth correction of SRTM C-band radar largely accounts for the inconsistency of the existing results for the eastern Pamir (Table 1). In addition, over the central and western Pamir, previous studies mainly focusing on the changes since 2000 have reported seemingly divergent results for similar observation periods, e.g. 0.14 ± 0.14 m w.e. a⁻¹ for 2000–11 (Gardelle and others, 2013) versus −0.12 ± 0.06 m w.e. a⁻¹ for 2000–14 (Lin and others, 2017). Nevertheless, according to the latest results of Brun and others (2017) on a 1° × 1° geographic grid, it appears that the glaciers have likely been in a nearly stable state during 2000–16 in the central Pamir. Meanwhile, the western Pamir may have experienced a certain degree of mass wastage (see Fig. 2 in Brun and others (2017) for details). In addition, especially for the Fedchenko Glacier situated in the central Pamir, long-term geodetic-based results have suggested that the glacier-wide mass balance has increased from −0.28 m w.e. a⁻¹ for 1928–2000 (Lambrecht and others, 2014) to −0.34 m w.e. a⁻¹ for 2000–16 (Lambrecht and others, 2018).

Table 1. Summary of the existing region-wide glacier mass balances reported in the Pamir Mountains

From the above, clearly, there is a knowledge gap about how the central and western Pamir glaciers reacted to climate change before 2000. Hence, the key purpose of this study was to extend the timescale of investigating glacier mass balance in the central Pamir using the early KH-9 (Keyhole-9, also known as Hexagon) images, and to further promote the understanding of the glacier changes of this region for a longer time period.

2. STUDY AREA

Due to being affected by the prevailing westerlies, the Pamir Mountains in central Asia receive a large amount of precipitation at high elevations, and are one of the most important glacier centers in the mid- to low-latitudes (Chevallier and others, 2014), covering a total area of ~12 000 km² (excluding the Pamir Alay) (RGI Consortium, 2017). Most of the glaciers concentrate in the central Pamir, which has an area of ~5400 km². The percentage of debris-covered glacier areas is ~8% in this region (Mölg and others, 2018). Of these glaciers, the Fedchenko Glacier, with a length of ~72 km and a total ice volume of ~124 km³, is one of the largest alpine glaciers in the world (Aizen and others, 2009; Lambrecht and others, 2014). In addition, meteorological records from two observation stations (Gorbunov Station and Altynmazar Station, at altitudes of 4169 m and 2782 m a.s.l., respectively; Fig. 1a) reveal a considerable difference in the vertical distribution of the precipitation. At Altynmazar Station, the mean annual precipitation is 159 mm, which is only one-seventh of that at Gorbunov Station (~1140 mm). Most of the precipitation occurs in winter and spring (from December to May), with an average of ~890 mm (Gorbunov Station) and 113 mm (Altynmazar Station), accounting for 78 and 71% of the annual precipitation, respectively. This indicates that the central Pamir glaciers are of the winter-accumulation type (Maussion and others, 2014; Sakai and others, 2015). Moreover, the mean annual air temperature is −6.9 and 3.5 °C at Gorbunov Station and Altynmazar Station, respectively, indicating relatively cold climate conditions in this region. Furthermore, a typical characteristic for the central Pamir is that glacier surge behavior frequently occurs (Kotlyakov and others, 2008; Sevestre and Benn, 2015; Wendt and others, 2017), reflecting the instability of the glacier systems.

Fig. 1. (a) Mosaic true-color Landsat image for the central Pamir (24 August 2000 and 16 September 2000). (b) The topography of the study area (SRTM DEM) shown as a hillshade map. The red lines with arrows denote major glaciers and their flow directions. The capital letters from A to Q represent the glacier names. A − Fedchenko Glacier, B − Grum Grjimailo Glacier, C − Tanimas-2 Glacier, D − Tanimas-3 Glacier, E − N. Kyskurgan Glacier, F − Kosinenko Glacier, G − Little Tanimas Glacier, H − Ayujilga Glacier, I − Mushketov Glacier, J − Bivachny Glacier, K − Vavilov Glacier, L − Shokalski Glacier, M − RGS (Russian Geographical Society) Glacier, N − Bears Glacier, O − Abdukagor Glacier, P − Yazgulem Glacier, Q − Ulugbeka Glacier, R − N. Tanimas Glacier. Note that there are various groups of Tanimas Glaciers, five draining into the Tanimas River and one (N) on the eastern slope.

3 DATA PROCESSING

3.1 Glacier outlines

In this study, we did not generate our own glacier boundaries, but directly used the Randolph Glacier Inventory (RGI) Version 6.0 dataset, which was released on 28 July 2017 by the Global Land Ice Measurements from Space (GLIMS) initiative (http://www.glims.org/RGI/rgi60_dl.html). Considering that the Glacier Area Mapping for Discharge from the Asian Mountains (GAMDAM) inventory (with a mapping uncertainty of 15%, see Nuimura and others (2015)) is an important part of the RGI dataset, we assume that the uncertainty of the glacier outlines used in this study are at the order of 15%. In order to obtain the glacier extent that corresponds to the observation period of this study, the RGI v6.0 dataset was further modified by referring to two Landsat/ETM+ images (path/row: 151/033 and 152/033; date: 24 August 2000 and 16 September 2000) and ortho-rectified KH-9 images (acquired on 13 July 1975). The adjustments that we made included changes in the terminal positions of the glaciers due to glacier surge behavior and in lateral outlines because of surface thickening and thinning (King and others, 2017). In addition, given the relatively low contrast of the KH-9 images, the final elevation difference map was used to help determine the glacier boundaries, especially for the debris-covered areas (Maurer and others, 2016).

3.2 KH-9 DEM generation

The KH-9 images declassified by the US Geological Survey in 2002 record land surface information between 1971 and 1986 at a relatively fine scale (a spatial resolution of 6–9 m), and have become a very valuable data source for reconstructing historical topography (Surazakov and Aizen, 2010; Pieczonka and others, 2013). In this study, we used the KH-9 triplet images to extract the KH-9 DEM via the Automated Terrain Extraction (ATE) routine contained in the ERDAS Leica Photogrammetry Suite (LPS) module. In the process of collecting the ground control points (GCPs), we used the Landsat/ETM+ images and the 1 arc-second SRTM DEM as the horizontal reference and the vertical reference, respectively (Pieczonka and others, 2013; Zhou and others, 2017, 2018b). A total of 43 GCPs were chosen at easily distinguished areas such as mountain ridges, river bends and road intersections. To independently examine the model error of the aerial triangulation, seven check points (CPs) were employed. Finally, the root-mean-square error (RMSE) of the stereo model was ~0.39 pixels, satisfying the constraint condition of being less than one pixel. The RMSEs of the GCP and CP residuals are ~10 m in the three directions (Table 2). For the extracted KH-9 DEM, statistics of the vertical accuracy with the GCPs and CPs as reference are given in Table 2. In addition, the DEM extraction report shows that most of the mass points (75.2%) used to interpolate the elevation in a regular grid are of good quality, with a correlation score above 0.7. The rest of the points (24.8%) have scores ranging from 0.5 to 0.7. Despite this, at high altitudes (e.g. in the accumulation zones), there are some clearly erroneous pixels characterized by irregular bulges in the hillshade map.

Table 2. The aerial triangulation information and vertical accuracy check for the generation of the KH-9 DEM

* The SD, MAE and RMSE, respectively, denote the std dev., mean absolute error and root-mean-square error.

3.3 DEM differencing and outlier elimination for the accumulation zones

As one of the most important topography data sources, the 1 arc-second SRTM C-band DEM taken in February 2000 was used as a reference in this study. Once the KH-9 DEM was produced, we followed the data processing strategy of Zhou and others (2017, 2018b) to check and correct some potential offsets and biases between the two DEMs. To be specific, we first co-registered the KH-9 DEM to the reference SRTM DEM using the classic method proposed by Nuth and Kääb (2011). Subsequently, we examined and corrected the tilt by fitting a second-order trend surface (Fig. S1), which is something that is often required for the DEMs derived from the early aerospace imagery (Bolch and others, 2011, 2017; Pieczonka and others, 2013; Li and others, 2017; Zhou and others, 2017, 2018b). In addition, with regard to some potential biases associated with terrain parameters (e.g. elevation, slope, aspect and maximum curvature), we first checked the maximum curvature-dependent bias, given that it is likely to be the cause of the elevation-dependent bias (Paul, 2008; Gardelle and others, 2012). We then examined the elevation-, slope- and aspect-dependent biases in sequence (Fig. S2). The residual elevation-dependent bias may result from an uneven spatial distribution of GCPs in the x-y-z planes (Nuth and Kääb, 2011). The final elevation difference in the off-glacier areas is shown in Figure S3. The overall processing flow resulted in an improvement of the std dev. of the elevation difference in ice-free areas of 15.95% in total (Table 3).

Table 3. Statistics of the elevation differences in the ice-free areas for the raw and corrected difference maps

After performing the above corrections, similar to the findings of previous studies utilizing the KH-9 DEMs (Pieczonka and Bolch, 2015; Maurer and others, 2016; Bolch and others, 2017; Zhou and others, 2017, 2018b), a remarkable surface lowering (more than 50 m) was found in the upper parts of the glacier accumulation zones (Fig. S4). Given that: (1) the KH-9 images have a relatively poor contrast in these snow-covered areas; and (2) the steep topography at high altitudes can easily induce the geometric distortion of stereo images, this phenomenon can be ascribed to erroneous matching of optical images in the process of producing the KH-9 DEM, and can thus be treated as outliers that need to be eliminated. In this study, with a prerequisite hypothesis that glacier thinning gradually decreases with the increase in altitude following a non-linear trend (Schwitter and Raymond, 1993), we applied the simple and effective method proposed by Zhou and others (2018b), who adopted a quadratic function to approximately reflect the non-linear variation. To be more specific, the first step was to normalize the elevation in the glacier areas (Eqn (1)), and the next step was to determine a proper range (Δh _range) of glacier elevation changes by adjusting an empirical coefficient (A) (Eqn (2)).

(1)$$\gamma = \displaystyle{{E_{\max} -E_{\rm g}} \over {E_{\max} -E_{\min}}}, $$

(2)$$\Delta h_{{\rm range}} = A\gamma ^2,$$

where E _max,E _min and E _g represent the maximum elevation, the minimum elevation and the elevation at a given point in the glacier areas. γ is the normalizing factor. It should be noted that the setting of the empirical coefficient primarily depends on the performance of the outlier removal in the accumulation zones. After many trials, the coefficient was specified as −190, which results in the elevation changes in the upper part of the accumulation zones being limited to the range of ±10 m after the outlier elimination. However, this also causes somewhat excessive elimination for effective elevation changes in local ablation zones, especially for the surge-type glaciers. We filled these gaps using the original values. The final elevation difference map after removing the outliers is shown in Figure S4, where the proportion of the data gaps reaches 42%.

3.4 Radar penetration correction

Previous studies have shown that the penetration depth of C-band radar has become one of the most significant uncertainty sources when comparing SRTM C-band DEMs with other topography data obtained by optical or laser sensors (Gardelle and others, 2013; Kääb and others, 2015). In this study region, based upon the differencing operation between the SRTM C-band DEM and the X-band DEM, Gardelle and others (2013) and Lin and others (2017) have reported an average penetration depth difference of 1.80 and 1.88 m, respectively. In this study, given that the penetration depth mainly depends on the physical characteristics of the glacier surface (Rignot and others, 2001), in order to better elaborate the penetration depth difference for different glacier surfaces, we used the Landsat-7 image acquired on 16 September 2000 (close to the end of the ablation season) to separate the debris-covered areas, bare ice and firn/snow-covered areas based on the band ratio method. To be specific, the ratio of band 3 (red channel) to band 5 (mid-infrared channel) with a threshold of 2 was used to distinguish the ice and snow from debris-covered areas. We further applied a single-band thresholding approach (band 3: digital number>165) to identify the perennial snow-covered areas. After producing the classification map of glacier surfaces, we generated the mean penetration depth difference of each 100 m altitude interval for every category and performed the correction for each class using an overall average. It should be noted that we eliminated some pixels where the terrain slopes were >30°, in order to avoid introducing elevation bias induced by the steep topography. Moreover, prior to computing the penetration depth difference for the glacier areas, a mean residual (0.67 m) in the ice-free areas, possibly reflecting a systematic vertical offset, was subtracted.

In addition, given that: (1) the X-band radar can penetrate the glacier surface (mainly referring to the dry snow layer) to as deep as 3–6 m (Dehecq and others, 2016; Zhao and Floricioiu, 2017); and (2) Wendt and others (2017) and Lambrecht and others (2018) found a maximum penetration depth of the X-band radar of ~6 m (with an average of 3–4 m) (TanDEM-X data) in the upper part of the Fedchenko Glacier's accumulation zones, there is no doubt that the C-band radar penetration depth we obtained is underestimated. Nevertheless, considering that this region generally receives a large amount of winter snowfall that is dry and cold, especially at high altitudes, the SRTM X-band DEM likely represents the glacier surface elevation at the end of the ablation season in 1999. This also means that the glacier surface elevation in summer of 1999 can be approximately recovered by performing the correction of the penetration depth difference for the SRTM C-band DEM. In others words, the impact on the final mass balance of the C-band radar absolute penetration depth being underestimated may not be significant, as the KH-9 images were only acquired in summer. Based on this, we assume that the seasonal change caused by the difference of the data acquisition times (mid-July versus mid-February) is negligible, given the long observation period.

3.5 Mass-balance calculation

To calculate the glacier mass balance, we first generated the histogram of the glacier elevation change with elevation bins of 100 m. Due to the large number of data gaps caused by outlier elimination at high altitudes, in order to obtain the volume change of the whole glacier area, we assigned the average elevation change to these voids in a given elevation band. The volume change was then converted to the mass change using a factor of 0.85 ± 0.06 (Huss, 2013). For each surge-type glacier, the mass balance was separately calculated. In addition, given that the difference in the strategies of filling data voids may lead to discrepant results, we calculated the mass balance for individual glaciers by replacing these holes with zeros based on the above method. To further test the influence of different settings of the empirical coefficient (see Section 3.2) on the mass-balance results, we used a simple approach of directly calculating the average mass balance for the entire region when the coefficient was specified as −130, −160, −190, −220 and −250, respectively.

3.6 Uncertainty evaluation

With regard to the accuracy assessment of the glacier mass balance, we used the std dev. (σ _z) of the elevation differences for each elevation band (i) in stable regions to represent the uncertainty of the glacier elevation change of each pixel (Gardelle and others, 2013). When calculating the average elevation change for a given elevation band, we needed to consider the influence of spatial autocorrelation between DEMs for the sake of evaluating the uncertainty (E _z) (Eqn (3)) (Rolstad and others, 2009). Based on the elevation difference map over the off-glacier areas, we fitted an experimental anisotropic semi-variogram using a spherical model (Fig. S5), and obtained a spatial autocorrelation distance (D) of 720 m (Eqn (4)).

(3)$$E_{{\rm z},i} = \displaystyle{{\sigma _{{\rm z},i}} \over {\sqrt {N_{{\rm eff},i}}}}, $$

(4)$$N_{{\rm eff},i} = \displaystyle{{N_{{\rm t},i}R} \over {2D}},$$

where N _t, N _eff and R represent the total number of observations, the number of independent measurements and the spatial resolution (30 m). Subsequently, the volume change uncertainty for a given glacier or region (E _v) was calculated by Eqn (5).

(5)$$E_{\rm v} = \sqrt {\sum {{(E_{{\rm z},i}S_i)}^2}}, $$

where S denotes the glacier area in a given altitude interval. To further calculate the uncertainty of the mass balance, with an assumption of a 15% glacier area error (E _S) (Nuimura and others, 2015) and an uncertainty for the conversion factor of 0.6 ($E_{{\rm f}_{\rm c}}$), the initial mass-balance uncertainty ($E_{{\rm b}\_0}$) can be obtained by Eqn (6) based on the standard error propagation law.

(6)$$E_{{\rm b}\_0} = \sqrt {{\left( {\displaystyle{V \over {S_{\rm t}}}\cdot E_{{\rm f}_{\rm c}}} \right)}^2 + {\left( {\displaystyle{{\,f_{\rm c}} \over {S_{\rm t}}}\cdot E_{\rm V}} \right)}^2 + {\left( {\displaystyle{{V\cdot f_{\rm c}} \over {S_{\rm t}^2}} \cdot E_{\rm s}} \right)}^2}, $$

where V and S _t denote the volume change and the total glacier area for a glacier or region. f _c is the conversion factor (0.85). In addition, since the penetration depth correction has a systematic effect on the final mass change, the uncertainty of the mass balance can thus be further expressed as $E_{{\rm b}\_1}$ (Eqn (7)).

(7)$$E_{{\rm b}\_1} = \sqrt {E_{{\rm b}\_0}^2 + E_{{\rm b}\_{\rm p}}^2}, $$

where $E_{{\rm b}\_{\rm p}}$ represents the uncertainty of the penetration correction, which was set to 1.28 m w.e. in this study (assuming a penetration depth error of 1.5 m) (Gardelle and others, 2013). Moreover, to further account for the uncertainty caused by the data gaps and the outlier elimination, we introduced a scale factor (1/(1 − R _void)) and an additional factor ($E_{{\rm b}\_{\rm c}}$) to calculate the final mass-balance uncertainty (E _b) by using Eqn (8).

(8)$$E_{\rm b} = \displaystyle{2 \over {1-R_{{\rm void}}}}\cdot E_{{\rm b}\_1} + E_{{\rm b}\_{\rm c}},$$

where R _void denotes the proportion of the data gaps (Table 4). It should be pointed out that the scale factor was conservatively multiplied by 2, given that almost half the pixels were removed, which likely led to larger uncertainty for the glacier mass change. In particular, for the uncertainty induced by the setting of the empirical coefficient, we chose the maximum range of the change of the mass-balance results in different scenarios as the additional factor (0.06 m w.e. a⁻¹ in the present study).

Table 4. The average elevation change and the mass balance for individual glaciers and the whole region between 1975 and 1999. Method 1: assuming zero change for missing data. Method 2: using average elevation changes at the same altitude bands to fill the gaps. Note that this study uses the results obtained by Method 2

* The SLE represents the snow line elevation, which was obtained by reference to the Landsat images (16 September 2000) and the SRTM DEM.

† denotes the surge-type glaciers. Note that the SG1–3 glaciers were used in this study to represent three unnamed glaciers.

4. RESULTS

4.1 Penetration depth difference

The final results show that the average penetration depth difference is 0.06 m for debris-covered areas (accounting for 18.5% of the selected areas), 1.68 m for bare ice areas (46.6%) and 2.96 m for firn/snow-covered areas (34.9%). In particular, based on a visual analysis of the available optical images, we speculated that the 0.06 m residual in the debris-covered areas is possibly caused by wet snow, as there was heavy snowfall in this area before 5 February 2000 (not shown), and most of the snow remained on 13 February 2000 (Fig. 2). Hence, we did not take this term into account for the correction. The penetration difference estimate we obtained is close to that of Holzer and others (2015) in the bare ice areas (i.e. 1.68 versus 1.50 ± 0.90 m), but is lower than that in the firn/snow areas (2.96 versus 4.30 ± 0.90 m). The latter represents the average of the penetration depth estimation over three regions (Karakoram, Hindu Kush and Jammu Kashmir), which were calculated by extrapolating the ICESat trend to the moment of SRTM acquisition, theoretically reflecting an absolute penetration depth. In addition, our penetration depth difference results are more or less higher than the estimates of Zhang and others (2016) for different glacier surfaces in the eastern Pamir (i.e. 0.79 and 2.41 m for the bare ice and the firn/snow covered areas, respectively) (Table 1). This can be attributed to the regional difference in glacier characteristics. Finally, the area-weighted average penetration depth difference is 1.82 m for the whole region, which is in line with the existing estimates in Gardelle and others (2013) and Lin and others (2017). In contrast, over the entire Pamir, Kääb and others (2015) found an average penetration depth of 5–6 m, which may be overestimated. This is because ICESat data samples are relatively sparse in the Pamir, and frequent glacier surges probably induce elevation bias, as stated by Kääb and others (2015). The updated glacier mass balance by Brun and others (2017) (−0.04 ± 0.19 m w.e. a⁻¹ for 2000–08) indicates that the results of Kääb and others (2015) (−0.41 ± 0.12 m w.e. a⁻¹ for 2003–08) may be biased; thus, obtaining the penetration depth from ICESat data is problematic.

Fig. 2. (a) The Landsat false-color image (RGB: bands 5/4/3) acquired on 16 September 2000. (b) The classification map of the glacier surface. (c) The Landsat false-color image (RGB: bands 5/4/3) on 13 February 2000. (d) The area-altitude distribution and the penetration depth difference for each category.

4.2 Glacier elevation change

From Figure 3, we can clearly see that some glaciers experienced significant surface thickening and thinning in their lower and upper parts, respectively (e.g. Bivachny Glacier, Kosinenko Glacier, Vavilov Glacier, Shokalski Glacier, Yazgulem Glacier and Bears Glacier), implying the emergence of surging in our observation period. In particular, at the forefront of the Bears Glacier, a striking thinning of ~−130 m can be observed (Fig. 4b). Such a pattern can be ascribed to two catastrophic surge events that occurred in 1973 and 1989, respectively. The magnitude of the earlier surge was greater than that of the later (Harrison and others, 2015), resulting in continuous thinning at the advanced fronts after 1973. A similar condition emerges at the terminus of the Kosinenko Glacier. In contrast, for the N. Kyskurgan Glacier, its accumulation and ablation zones, respectively, showed significant thickening and thinning (i.e. 0.50 ± 0.31 and −1.87 ± 0.31 m a⁻¹ on average) (Table 4; Fig. 4b). This pattern indicates that the glacier surged before 1975 and was in a post-surge or quiescent phase for this study period. Similar conditions can be found on the nearby glaciers and the N. Tanimas Glacier. On the RGS Glacier, the ablation zone as a whole exhibited a moderate thickening (~10–15 m). This may have been caused by a micro-surge at the beginning of the 1970s (Kotlyakov and others, 2008). Moreover, according to previous studies, there are four additional small surge-type glaciers in this region (i.e. Ravak Glacier and the SG1–3 glaciers in Fig. 3) (Kotlyakov and others, 2008; Gardelle and others, 2013). Although from the elevation change pattern only the Ravak Glacier clearly exhibits the surge features (Fig. 3), the statistical results suggested that the others may also have experienced micro-surges (Table 4). For the non-surge-type glaciers, the most representative Fedchenko Glacier showed an average thinning of −0.33 ± 0.33 m a⁻¹ in the ablation zone and a near-zero change (0.02 ± 0.33 m a⁻¹) in the accumulation zone (Table 4; Fig. 4b). Detailed information of the remaining glaciers is given in Table 4.

Fig. 3. Glacier elevation changes in the central Pamir for 1975–2000. The purple dotted lines represent the profiles in Figure 4b along the central flowline. The black lines denote the glacier boundaries. The cyan triangles represent the surge-type glaciers. The background is the Landsat true-color images, and the white areas (not covered by the elevation change map) represent the extent of the data gaps.

Fig. 4. (a) The area-elevation distribution and average elevation change with a 100 m elevation interval for all of the non-surge-type glaciers. (b) The elevation changes of some representative glaciers along the central flowline.

4.3 Glacier mass balance

With regard to the mass change, from 1975 to 1999, the central Pamir glaciers (covering a total glacier area of 2020.3 km²), as a whole, showed a nearly balanced or slightly negative mass balance of −0.03 ± 0.24 m w.e. a⁻¹, which is equivalent to a mass deficit of −0.06 ± 0.48 Gt a⁻¹. Of these glaciers, the Tanimas-3 Glacier and N. Kyskurgan Glacier showed the most negative mass balance, with a value of −0.12 ± 0.26 m w.e. a⁻¹, while the most significant mass gain, at a rate of 0.63 ± 0.20 m w.e. a⁻¹, occurred at the Ravak Glacier. For the Fedchenko Glacier, the mass balance was −0.09 ± 0.28 m w.e. a⁻¹, corresponding to a mass wastage of −0.05 ± 0.15 Gt a⁻¹. The mass change results for the other glaciers are listed in Table 4. In addition, for the surge-type glaciers (468.9 km²) and non-surge-type glaciers (1651.4 km²), their mass balances were 0.03 ± 0.14 and−0.05 ± 0.28 m w.e. a⁻¹, respectively.

4.4 Sensitivity to outlier elimination and data gap filling

With regard to the outlier elimination for the accumulation zones, when the empirical coefficient was set as −130, −160, −190, −220 and −250, the corresponding region-wide mean mass balance was −0.01, −0.02, −0.04, −0.05 and −0.07 m w.e. a⁻¹, respectively. The maximum difference between these results is 0.06 m w.e. a⁻¹, which is within the error bars. In addition, the estimate of −0.04 m w.e. a⁻¹ calculated based on the average of the whole area agrees with that (−0.03 m w.e. a⁻¹) obtained by the different calculation approach with respect to the surge-type glaciers. For the processing of the data gaps, the two results based on the different gap-filling strategies are basically consistent (Table 4). The region-wide mass balance (filled with zeros) was −0.01 ± 0.24 m w.e. a⁻¹, and the non-surge-type and surge-type glaciers showed a mass balance of −0.02 ± 0.28 and 0.03 ± 0.14 m w.e. a⁻¹, respectively. For individual glaciers, the differences between these two results are generally <0.03 m w.e. a⁻¹, but reach as high as 0.06 m w.e. a⁻¹ for the Kosinenko Glacier and N. Tanimas Glacier. Even so, all the differences are within the error bars, suggesting good compatibility.

5. DISCUSSION

5.1 Comparison with previous studies

For the Fedchenko Glacier, Lambrecht and others (2014) reported an ice loss of −4.1 km³ (equivalent to a mass budget of −0.39 m w.e. a⁻¹) for 1958–2000 in its main trunk. This is more negative than our estimate of −0.20 m w.e. a⁻¹ for 1975–99 in similar areas. Similarly, in its ablation zone, volume change rates for each 100 m elevation band in this study are smaller than theirs on the whole (Fig. 5), but the variation trends are basically consistent, particularly in areas below 3800 m elevation. This seems to suggest a relatively high rate of mass loss in the lower part of the Fedchenko Glacier between 1958 and 1975. Nevertheless, there is no other evidence to support this intuitive speculation due to the lack of both in situ measurements and remote-sensing data for that period. Above an elevation of 3800 m, the two rates are almost equal in areas between 3800 and 4300 m, and exhibit a greater fluctuation above 4300 m. This local variation may be caused by inaccurate matching in the course of generating the KH-9 DEM and/or the ‘terracing’ effect within the earlier DEM (for 1958) obtained by interpolating from the topographic maps (Racoviteanu and others, 2007). Other reason for the above discrepancies is the difference in the observation periods. Detailed comparison results are given in Table S1. Furthermore, in comparison with the glacier-wide results after 2000 (e.g. −0.27 ± 0.05 m w.e. a⁻¹ for 2000–11 and −0.51 ± 0.06 m w.e. a⁻¹ for 2011–16) reported by Lambrecht and others (2018), this study period experienced less mass loss (−0.20 m w.e. a⁻¹). The findings reveal an intensified mass wastage at the Fedchenko Glacier, especially after 2011, as reported in Lambrecht and others (2018).

Fig. 5. Comparison of the mean volume change rate (for each altitude interval of 100 m) in the Fedchenko Glacier's ablation zone (with an altitude range of 2800–4850 m) between Lambrecht and others (2014) and this study.

As for the region-wide glacier mass change, it is important to note that we did not make a comparison of our geodetic results (in local regions) before 2000 with those of the entire Pamir after 2000, given the spatio-temporal heterogeneity of the glacier change pattern within this area. For the central Pamir, our estimate of −0.03 ± 0.24 m w.e. a⁻¹ lies in between two previous results (i.e. 0.14 ± 0.14 in Gardelle and others (2013) and −0.12 ± 0.06 m w.e. a⁻¹ in Lin and others (2017)), but is closer to the result of Lin and others (2017), despite a substantial overlap of error bounds. However, considering the pattern of glacier elevation change from Brun and others (2017) (not affected by the radar penetration depth), we can deduce that the central Pamir glaciers have most likely been in an approximately balanced state over the past four decades (from the mid-1970s to the mid-2010s). Moreover, by comparing the mass changes of the surge-type and non-surge-type glaciers (i.e. 0.19 ± 0.22 and 0.12 ± 0.11 m w.e. a⁻¹, respectively), Gardelle and others (2013) concluded that the surge behavior did not significantly affect the region-wide mass balance over a short time period (~10 years). This agrees with our findings prior to 2000, with a change of 0.03 ± 0.14 (surge-type) and −0.05 ± 0.28 m w.e. a⁻¹ (non-surge-type). Coincidentally, in two observation periods (1975–99 and 2000–11), the differences between these two results are very close, i.e. 0.07 m w.e. a⁻¹ in Gardelle and others (2013) and 0.08 m w.e. a⁻¹ in this study. From the above, it appears that the conclusion of Gardelle and others (2013), i.e. the surge behavior did not significantly affect the region-wide mass balance, is valid for a longer period, at least from our observation results.

5.2 Potential reasons for the glacier changes

In order to understand the causes of the glacier mass changes in the central Pamir, we used the meteorological records for 1935–94 from Gorbunov Station and Altynmazar Station (Williams and Konovalov, 2008). The relevant data have been subjected to rigorous quality checks (e.g. correcting the precipitation bias for the gauge type) and homogeneity adjustment (only for the precipitation data). More details can be found in the user guide (http://nsidc.org/data/G02174). In this study, we analyzed the variation trends of annual precipitation and mean air temperature for two sub-periods (1935–75 and 1975–94), as shown in Figure 6. Overall, it seems that both precipitation and temperature at high altitudes (i.e. Gorbunov Station) show an opposite tendency compared with low altitudes (Altynmazar Station) within the two sub-periods. For the period 1975–94, the low altitudes show signs of increased precipitation and decreased temperature, whereas there is rising temperature and decreased precipitation at high altitudes, especially for 1990–94. Moreover, it should be stated that the mean annual air temperature at low altitudes showed an increasing trend for the whole period, and the magnitude of the warming (k = 0.0041, p > 0.05) appears to have been lower than that (k = 0.0085, p > 0.05) at high altitudes. From the above, we intuitively speculate that the rising air temperature and decreased total precipitation at high altitudes to a great extent account for the slight mass loss, despite the increased precipitation and decreased temperature at low altitudes.

Fig. 6. Mean annual air temperature and precipitation for 1935–94. The red dashed line represents the regressed linear trend, and ‘k’ and ‘p’ denote the slope of the regressed trend and the significance level.

In addition, on a larger spatial scale, Zhou and others (2018a) generated a unified air temperature and precipitation dataset (including a gap-filled station dataset and a derived grid dataset) for 1951–2010 by compiling and merging meteorological observation materials from 457 stations in central Asia. Analyzing the datasets suggests that central Asia, overall, experienced statistically significant warming around the 1970s, while the mean annual precipitation did not show prominent change over the whole period (Zhou and others, 2018a). This, to some extent, supports the above speculation.

6. CONCLUSIONS

Based on the early KH-9 stereo images and the 1 arc-second C-band SRTM DEM, this study evaluated the glacier mass change for 1975–99 in the central Pamir by employing the approach of DEM differencing. To alleviate the impact of the C-band radar penetrating into the glacier surface, we separately estimated the penetration depth differences for different glacier facies by subtracting the SRTM C-band DEM from the X-band DEM. The obtained penetration corrections were zero for debris-covered areas, 1.68 m for clean ice and 2.96 m for firn/snow areas. Finally, the overall mass balance was found to be −0.03 ± 0.24 m w.e. a⁻¹ for 1975–99, indicating a nearly balanced state or very slight mass wastage. This may be attributed to the decreased precipitation and increased air temperature at high altitudes. Individual glaciers have shown a highly heterogeneous pattern of mass change, due to the glacier surge activity, with the mass balance varying from −0.12 ± 0.26 to 0.63 ± 0.20 m w.e. a⁻¹. For the largest Fedchenko Glacier, it showed a negative mass balance of −0.09 ± 0.28 m w.e. a⁻¹, corresponding to a mass budget of −0.05 ± 0.15 Gt a⁻¹. In addition, for the different types of glacier (surge and non-surge), the mean mass balances were 0.03 ± 0.14 and −0.05 ± 0.28 m w.e. a⁻¹, respectively. Given the natural defects of the KH-9 images and the uncertainty sources caused by the data processing, we tested the effect of outlier elimination and data gap filling on the final mass balance. The results suggest that our findings are quite robust for the chosen geodetic method, given the large uncertainty range. Finally, integrating our results with those of previous studies indicates that, over the past four decades (between the mid-1970s and the mid-2010s): (1) the central Pamir glaciers have been in a nearly stable state or have shown only slight mass loss; (2) the collective geodetic mass balance of glaciers that surged during the study period was similar to that of glaciers that did not surge.