Slight glacier mass loss in the Karakoram region during the 1970s to 2000 revealed by KH-9 images and SRTM DEM

Abstract

An anomalously slight glacier mass gain during 2000 to the 2010s has recently been reported in the Karakoram region. However, to date, no investigations of the region-wide glacier mass balance in the Karakoram prior to 2000 have been reported, leaving a knowledge gap for assessing glacier responses to climate change. We calculated elevation and mass change using DEMs generated from KH-9 images acquired during 1973–1980 and the 1 arc-second SRTM DEM. We find a slight mass loss of −0.09 ± 0.03 m w.e. a⁻¹ (12 366 km²) for 1973–2000, which is less negative than the global average rate for 1971–2009 (−0.31 ± 0.19 m w.e. a⁻¹). Within the Karakoram, the glacier change patterns are spatially and temporally heterogeneous. In particular, a nearly stable state in the central Karakoram (−0.04 ± 0.05 m w.e. a⁻¹ during the period 1974–2000) implies that the Karakoram anomaly dates back to the 1970s. Combined with the previous studies, we conclude that the Karakoram glaciers as a whole were in a nearly balanced state during the 1970s to the 2010s.

1. INTRODUCTION

Mountain glaciers, which are a significant fresh water resource, play a pivotal role in regional hydrological cycles and the ecological environment (Kaser and others, 2010). Most glaciers have been retreating since the end of the Little Ice Age (Marzeion and others, 2014). In the recent decades, in particular, nearly worldwide glacier shrinkage and mass loss have been observed (Gardner and others, 2013). Certainly, the glaciers in High Mountain Asia are no exception and exhibit a high rate of area shrinkage (Li and others, 2008; Cogley, 2016). Glacier decline can cause flooding and geological disasters such as landslides and debris flows, which can seriously affect the security of the downstream areas (Benn and others, 2012; Rankl and others, 2014). A comprehensive view of the region-wide glacier changes is therefore essential to support hydrological management, disaster warning and glacier evolution projection.

The Karakoram is known as a home to giant glaciers, including many of surge type (Hewitt, 2005; Rankl and others, 2014). Usually, the surge-type glaciers experience more frequent terminus changes than the generic glaciers. Accordingly, most of the previous studies have paid attention to the surge-related length and area changes on a regional scale (Hewitt, 2005, 2007, 2011; Copland and others, 2011; Scherler and others, 2011; Bolch and others, 2012; Bhambri and others, 2013; Rankl and others, 2014; Herreid and others, 2015; Paul, 2015). However, of particular note is that neither area nor length changes can reflect unambiguously glacier response to climate changes, especially for heavily debris-covered glaciers (Scherler and others, 2011; Rankl and others, 2014). Hence, glacier mass balance, which is a direct and reliable indicator of climate change, must be investigated. Hewitt (2005) first reported a glacier thickening of ~5–15 m in the ablation zones of some large glaciers during 1997–2002 in the central Karakoram, based on field measurements. This finding was in sharp contrast to the situation in other High Asia regions and was therefore considered as an anomaly. Subsequently, using SPOT DEMs and SRTM DEM, Gardelle and others (2012a, 2013) confirmed the glacier anomaly from 2000 to the 2010s in the eastern and central Karakoram. However, several studies based on ICESat/GLAS trends have reported a slightly negative or nearly stable glacier mass balance in the entire Karakoram for 2003–09 (Kääb and others, 2012, 2015; Gardner and others, 2013). More recently, by subtracting SRTM DEM from TanDEM-X DEM and Cartosat-1 DEM, respectively, Rankl and Braun (2016) and Agarwal and others (2016) obtained similar results in the Shaksgam valley (central Karakoram, for 2000–12) and the Siachen glacier (eastern Karakoram, for 2000–07). At present, the mass-balance records of the Karakoram glaciers are mainly concentrated in the early 21st century, and the observation times cover ~10 years, which leaves scope for extending the timescale of the mass-balance monitoring. Herreid and others (2015) found that the debris-covered glacier area in the central Karakoram hardly varied between 1977 and 2014, and deduced that the Karakoram anomaly probably prevailed at times earlier than those previously identified. However, this conclusion is qualitative and remains to be verified. In order to determine the Karakoram glacier changes during the 1970s–2000 in a quantitative way, we adopted KH-9 stereo images from 1973–80 and the 1 arc-second SRTM DEM to measure the region-wide glacier mass balance.

2. STUDY AREA

The Karakoram mountains (33°45′–37°30′N, 74°20′–78°15′E) are located in the northwest of the Tibetan Plateau and lie in a southeast–northwest orientation (see Fig. 1), with an average altitude of 5000 m a.s.l. In total, 15 peaks in the Karakoram mountains are over 7000 m a.s.l. (Xie and Liu, 2010). Among them, K2 (8611 m a.s.l.) is the second highest mountain in the world. Glaciers developed by the Karakoram mountains have a total area of ~ 18 000 km² (Bolch and others, 2012), and most of them are spread within the altitude range of 2300–8600 m a.s.l. (Hewitt, 2011). As mentioned above, many glaciers in the area are of giant size. Specifically, the Siachen, Baltoro, Biafo, Hispar and Batura glaciers have a total of drainage basin areas as large as 1400, 1270, 855, 785 and 710 km², respectively (Hewitt, 2013). Due to the frequent snow avalanches, the tongues of the glaciers in this area are generally covered by debris. However, Mayer and others (2006) reported that the debris thickness varies greatly, even within small zones. At the terminus of a giant glacier, the debris can be over 1.5 m thick.

Fig. 1. Overview of the study area. Snow and ice are marked in white. The yellow and purple rectangles denote the frames of the KH-9 images and the SRTM X-band DEM, respectively. The green dots mark the largest glaciers mentioned in context. The background is a Moderate Resolution Imaging Spectroradiometer (MODIS) image with a spatial resolution of 500 m (acquired by The Blue Marble Next Generation Project).

In winter, the climate of the Karakoram region is mainly dominated by westerlies (Gardelle and others, 2013), which bring approximately two-thirds of the high-altitude snowfall (Bolch and others, 2012). The Indian summer monsoon also has a small impact on the eastern margins of the Karakoram region (Bhambri and others, 2013). Furthermore, the jet stream plays an important role in adjusting the precipitation in the vertical direction (Hewitt, 2005; Mayer and others, 2006). The amount of precipitation at high altitudes is 5–10 times that of the low altitudes (Hewitt, 2005; Immerzeel and others, 2015). For example, the precipitation can reach 1200–2000 mm a⁻¹ (mainly in solid form) in the regions near the snowline (~5000 m a.s.l.) (Xie and Liu, 2010), but is only 150–200 mm a⁻¹ in the altitude range of 2000–3200 m a.s.l. (Mayer and others, 2006).

Several studies have demonstrated that glacier surges are frequent in the Karakoram region (Copland and others, 2011; Rankl and others, 2014). The updated glacier inventory in Rankl and others (2014) lists 101 surge-type glaciers in this area. Generally, surge-type glaciers have a greater length than normal glaciers (Barrand and Murray, 2006; Rankl and others, 2014). In fact, the surge-type glaciers should number more than the result of Rankl and others (2014), as the previous observations are limited in both the temporal and spatial scales. A glacier surge can easily block a stream valley and form a glacial lake, which can severely threaten the lives of those living downstream (Hewitt and Liu, 2010; Xie and Liu, 2010; Hewitt, 2011).

3. DATA AND METHODS

3.1. Glacier delineation with landsat imagery

Landsat ETM + L1 T images acquired ~2000 were employed to extract the glacier outlines. In order to minimize the influence of seasonal snow, we chose seven images collected during the ablation seasons (May–September) of 2001 and 2002. All Landsat L1 T images have been geometrically corrected (under UTM projection), and have a spatial resolution of 30 m (multispectral band). In this paper, we used thresholding of the band ratio, a semi-automated method, for the glacier mapping (Bolch and others, 2010). Specifically, a threshold of 2 was set for the TM3/TM5 band DN ratio to distinguish the clean ice and snow from other objects. The initial glacier mask was filtered using a 5 × 5 median filter to remove noise. Although the TM3/TM5 band DN ratio works well in shadows (Paul and Kääb, 2005; Andreassen and others, 2008), it tends to misclassify water bodies as glaciers (Bolch and others, 2010). Furthermore, debris-covered glaciers cannot be recognized. Hence, the post-editing including the elimination of water bodies and cloud, and the delineation of the debris-covered glacier areas was performed with the assistance of pan-sharpened band images (15 m).

3.2. DEM generation with KH-9 imagery

The KH-9 Hexagon mission, which operated from 1973 to 1980, recorded a large amount of surface information with a high spatial resolution (6–9 m) (Surazakov and Aizen, 2010). The US Geological Survey (USGS) declassified the KH-9 images in 2002 and distributed them to the scientific community at a low cost. Because the KH-9 images have an excellent stereo imaging ability (base-to-height ratio: 0.4), they are thus suitable to extract the valuable early DEMs. In this paper, a total of 14 KH-9 images acquired between 1973 and 1980 (Table 1) were employed to generate the new DEMs.

Table 1. Details of the data adopted in this study

To generate the new DEMs, the first step was to correct the image geometric distortion based on the reseau crosses within the images (Surazakov and Aizen, 2010). The second step was to perform histogram equalization and locally adaptive filtering in order to enhance the local contrast (Pieczonka and others, 2013). A frame camera model in the ERDAS 2013 LPS (Leica Photogrammetry Suite) module was then adopted to extract the new DEMs. In order to achieve fine triangulation results, at least 30 GCPs (Table 2) were measured on regions with distinct features (such as angle of a road/river, road/river junctions and mountain ridges) and it was necessary to ensure that they were distributed as evenly as possible. The unit-weight triangulation RMSEs of all the stereo pairs were better than 1 pixel, and the RMSE residuals of the GCPs were ~10 m in both the horizontal and vertical directions (Table 2). Finally, the raw KH-9 DEMs were resampled to a 30 m resolution and were filtered by a 5 × 5 median filter. The new KH-9 DEM was employed to produce ortho-rectified KH-9 images with the same resolution (30 m), which were treated as base maps for revising the glacier inventory in 2001/02 in order to obtain the most extensive glacier outlines during our observation period.

Table 2. The triangulation information in the ERDAS LPS 2013 and the statistics of the nonglacial elevation difference

3.3. KH-9 DEM and SRTM DEM differencing

The SRTM C-band DEM acquired in February 2000 is considered a significant dataset for glacier research (Nuth and Kääb, 2011), and its relative accuracies in the horizontal and vertical directions are better than 45 and 10 m, respectively (Marschalk and others, 2004). In 2014, the USGS started to freely distribute the 1 arc-second SRTM C-band DEM covering 80% of the Earth's land surface to the public. We adopted the 1 arc-second SRTM C-band DEM as the benchmark DEM and subtracted the new KH-9 DEMs from it. In order to obtain accurate elevation changes, the horizontal offsets between the two DEMs should be as small as possible. A robust co-registration method with a sub-pixel precision was therefore applied to the two DEMs (Nuth and Kääb, 2011). After co-registration, we discerned a systematic spatial trend (related to the horizontal position) in the elevation difference map, and removed it by a second-order surface fitting (Bolch and others, 2011; Pieczonka and others, 2013). In addition, corrections of the potential errors related to altitude, slope, aspect and maximum curvature were carried out by establishing the relationship between the elevation differences and the related terrain parameters using polynomial fitting (Rolstad and others, 2009; Fischer and others, 2015). The statistics of the raw and corrected elevation differences over stable regions are listed in Table 2 (details can be seen in Table S1, Fig. S1–S2, Supplement). The series of corrections mentioned above resulted in a standard deviation (SD) improvement of ~15% for each elevation difference map (Table 2). We also used the normalized median absolute deviation (NMAD) as an additional index to describe the precision of the elevation difference (Pieczonka and Bolch, 2015). Höhle and Höhle (2009) showed that the NMAD is a more robust statistical description than the SD because the NMAD is less sensitive to outliers. Finally, following the method proposed by Pieczonka and Bolch (2015), with the assumption that glacier thinning reaches the maximum at the glacier front and gradually decreases with increasing altitude, we removed the outliers characterized by the peaks of surface thinning in the accumulation areas. The reason is that the outliers may represent the results of erroneous matching since the optical images have poor textures due to the brightness saturation in the snow-covered regions.

3.4. C-band penetration depth calculation

The 1 arc-second SRTM X-band DEM, acquired simultaneously with the SRTM C-band DEM, has been used to quantify the C-band microwave penetration depth. At present, the penetration depths are primarily calculated by comparing the 1 arc-second SRTM X-band DEM with the 3 arc-second SRTM C-band data (Gardelle and others, 2012a, b, 2013; Junfeng and others, 2015). Paul (2008) and Gardelle and others (2012b) showed that the distinction between the two DEMs’ initial pixel spacing sizes can cause systematic bias in the elevation difference map. We therefore employed the 1 arc-second SRTM C-band DEM to calculate the penetration depth. In Gardelle and others (2012a) and Gardelle and others (2013), the penetration depths of the Karakoram region were listed as 3 and 3.4 m, respectively. To check the accuracy of the penetration depth based on the 1 arc-second SRTM C-band DEM, we also performed the calculation based on the 3 arc-second SRTM DEM (Gardelle and others, 2012b). Note that, in this paper, the study areas are classified into two categories (glacier and non-glacier areas), unlike the approach used in Gardelle and others (2012b) (bare ice-free, snow-covered ice-free, debris-covered glacier and snow-covered glacier areas). As shown in Figure 2a, for the glacier areas, the penetration depths increase as the altitude rises. However, the penetration depth based on the 1 arc-second SRTM C-band DEM is systematically lower than that of the 3 arc-second DEM. The average difference is ~0.5 m. For the non-glacier areas, below 5000 m a.s.l., the penetration depths are close to zero, while above 5000 m a.s.l., the penetration depths begin to increase with altitude and reach 5 m at 6000 m a.s.l.. Gardelle and others (2012b) obtained a similar but slightly higher trend and attributed the anomaly to the penetration in snow-covered ice-free areas. The average difference of the non-glacier areas is ~0.37 m. In addition, the relationship between maximum curvature and elevation difference indicates that there is systematic bias caused by the distinction between the pixel sizes of the two DEMs employed for differencing (Fig. 2b), which could account for the relatively smaller results when using the 1 arc-second SRTM data. In this study, the average differences in the glacier areas based on the 1 arc-second and 3 arc-second SRTM C-band DEM are 2.7 and 3.2 m, respectively. Although the latter is equal to the average of the penetration depths in Gardelle and others (2012a) and Gardelle and others (2013), we can still utilize the average difference between glacier areas and non-glacier areas to approximately estimate the region-wide mean penetration depth. This is because the residual average differences in the non-glacier areas likely include uncorrected vertical offsets at high altitudes. In addition, considering the fact that the average differences between glacier areas and non-glacier areas are approximate when using the 1 arc-second and 3 arc-second SRTM C-band DEM, respectively (2.66 and 2.53 m) (Figs 2c and d), our strategy thus can be considered to be reasonable. The final average penetration depth is 2.6 m (the average of the two results mentioned above). Within the error bars, this also agrees with the result of 2.4 ± 0.3 m in the Karakoram region calculated by Kääb and others (2012) based on the ICESat/GLAS data. This confirms that our result is reliable. The correction was then realized by adding the mean penetration depth to the region-wide average glacial elevation change.

Fig. 2. (a) The relationship between elevation difference and altitude; (b) the relationship between terrain maximum curvature and the elevation difference between the 1 arc-second SRTM-X data and the 3 arc-second SRTM-C data; (c) distribution and statistics of the SRTM-X data minus the 3 arc-second SRTM-C data over glacier and stable areas; (d) distribution and statistics of the SRTM-X data minus the 1 arc-second SRTM-C data over glacier and stable areas.

3.5. Glacier mass-balance calculation and accuracy assessment

When calculating the region-wide glacier mass balance, an assumption that glaciers experience similar changes within a certain altitude range is often made (Berthier and Toutin, 2008; Berthier and others, 2010; Gardelle and others, 2013; Fischer and others, 2015). After filling all of the data gaps (including data voids mostly resulting from elimination of outliers and initial data voids in the SRTM DEM) by using the average elevation change at the corresponding altitude range, the glacier volume change can be calculated by summing the products of the area and average thickness change of each altitude interval (100 m). For surge-type glaciers, the volume changes are usually calculated separately (Gardelle and others, 2012a, 2013). A glacier average density of 850 ± 60 kg m⁻³ has been used to convert the volume change into mass budget (Huss, 2013). In addition, in order to present the mass balance of integral years, we performed seasonal corrections. The average accumulation rate (0.13 m w.e. per winter month) measured on the Biafo glacier (35°53′N, 75°41′E) was employed (Wake, 1989; Gardelle and others, 2013). Note that for the images acquired in mid-June (the central Karakoram in Fig. 1), we corrected the result with the accumulation of 2 months (from mid-December to mid-February), assuming that the glacier mass loss between mid-June and mid-September was equal to the mass gain between mid-September and mid-December.

We used the SD of the mean elevation differences in each altitude interval to represent the error of the elevation change in an individual grid point (Gardelle and others, 2013). Considering the spatial correlation among the DEM pixels, we assumed a DEM autocorrelation distance of 600 m, which corresponds to 20 pixels spacing (Rolstad and others, 2009; Nuth and Kääb, 2011). Subsequently, according to the equations in Gardelle and others (2013), we can obtain the error of the mean elevation change for each altitude band and calculate the uncertainty of volume change (ΔV) based on the standard principles of error propagation. In order to estimate the uncertainty of mass balance (Δm), the glacier density uncertainty (Δρ _i ) was considered (±60 kg m⁻³). In addition, we also took the error of glacier inventory (ΔS) into account. Given that Guo and others (2015) presented a detailed and reliable assessment on the positioning accuracy of the glacier boundaries depicted from Landsat images, we applied their glacier area error of 5.5% in the Indus river basin to approximately represent the area error in the Karakoram. According to the error propagation, we can estimate the initial uncertainty of the mass balance (Δm ₀) using Eqn (1), where V and ρ _w denote the change of glacier volume and the water density (1000 kg m⁻³), respectively. Moreover, by considering the C-band penetration depth uncertainty (Δm _p) of 1.5 m (equivalent to mass balance of 1.28 w.e.) and assuming a seasonal correction uncertainty (Δm _s) of 100% (0.13 m w.e. per winter month), the ultimate mass-balance uncertainty (Δm) can be estimated by using Eqn (2).

(1) $$\Delta m_{\rm 0} = \sqrt {\left( {\Delta V\displaystyle{{\rho _i} \over {\rho _{\rm w} S}}} \right)^2 + \left( {\Delta \rho _i \displaystyle{V \over {\rho _{\rm w} S}}} \right)^2 + \left( {\Delta S\displaystyle{{\rho _i V} \over {\rho _{\rm w} S^2}}} \right)^2} $$

(2) $$\Delta m = \sqrt {\Delta m_0 ^2 + \Delta m_{\rm p} ^2 + \Delta m_{\rm s} ^2} $$

4. RESULTS

4.1. Geodetic mass balance

As shown in Figure 3, within the western Karakoram, the surge glaciers experienced significant thickness changes. Most of the surge-type glaciers experienced obvious thinning and thickening in their downstream and middle-stream parts, respectively. More precisely, for the small surge-type glaciers (in a quiescent or post-surge stage) in the north (Fig. 3a), the mass gain in the middle-stream part could not cover the mass loss in the downstream part, while for the large surge-type glaciers (in a surge stage) in the south (Fig. 3b), for example Toltore and Passu, the mass gain in the middle-stream parts surpassed the mass loss in the downstream parts. Similar thickness change characteristics were also found in the central and eastern Karakoram (see Figs S3–S7 in the supplementary material).

Fig. 3. Glacier thickness change maps in the western Karakoram (the location is shown in Fig. 1) for 1973–2000: (a) north part; (b) south part. The background is Landsat ETM+ images. The black curves denote glacier outlines. Note that the magenta triangles and circles represent the glaciers in a surge and quiescent or post-surge stage, respectively.

Figure 4 shows the regional mass balance of the sub-regions, and the corresponding statistics are shown in Table 3. For the period of 1973/74–2000, overall, the most negative region-wide mass balance of −0.11 ± 0.04 m w.e. a⁻¹ was observed in the western Karakoram, followed by the eastern Karakoram (denoted by E-K-E in Fig. 4) with a value of −0.07 ± 0.05 m w.e. a⁻¹. However, a less negative mass balance of −0.04 ± 0.05 m w.e. a⁻¹ occurred in the central Karakoram (C-K-W), being in a slight mass wastage or nearly stable status. Furthermore, on a more local scale, the north part of the western Karakoram (W-K-N) had a mass balance of −0.17 ± 0.06 m w.e. a⁻¹, which is more negative than that of the south part (W-K-S) with a rate of −0.08 ± 0.05 m w.e. a⁻¹ . A similar case was also found in the central Karakoram, for example −0.14 ± 0.05 m w.e. a⁻¹ in the north side (C-K-W-N) versus −0.03 ± 0.05 m w.e. a⁻¹ in the south side (C-K-W-S). For the period of 1979/80–2000, the central and eastern Karakoram areas (C-K-E and E-K-W) showed a mass balance of −0.16 ± 0.07 m w.e. a⁻¹ and −0.13 ± 0.08 m w.e. a⁻¹, respectively, which was found to be more negative than that of the other sub-regions (W-K, C-K-W and E-K-E) for 1973/74–2000 (Fig. 4). The overall mass balance was −0.09 ± 0.03 m w.e. a⁻¹ for the glacier areas (12 366 km²) covered by the new KH-9 DEMs, which is equivalent to a mass budget of −1.15 ± 0.31 Gt a⁻¹. If we extrapolate the measured regional mass balance to all the glaciers in the Karakoram region (a total area of 22 000 km²), the total glacier mass loss for 1973–2000 would be −51.48 ± 17.16 Gt.

Fig. 4. Regional mass balance of the sub-regions in the Karakoram for 1973–2000. Also note that considering the coverage and acquisition times of the KH-9 images (see Fig. 1), the central and eastern Karakoram regions are divided into two parts (east and west), respectively. For example, the west part of the central Karakoram is represented by C-K-W. In addition, for the W-K and C-K-W regions, we also separately divided them into north and south parts. Regarding the overlapping areas between images, they are exploited in Section 5.1 for consistency checks.

Table 3. Glacier regional average elevation, volume, mass balance and mass budget change rates in the sub-regions

5. DISCUSSION

5.1. Internal consistency check

A triangulation method based on at least three datasets is often used to check the internal consistency in change monitoring (Nuth and Kääb, 2011). In this study, two of the new KH-9 DEMs partly overlapped in the eastern Karakoram (Fig. 1), which provided us with an opportunity to undertake an internal consistency check of the glacier mass change measurements based on the geodetic method. Three DEMs (two KH-9 DEMs and the SRTM DEM) were grouped into three pairs. Each pair corresponded to an elevation difference map. The same correction methods described above were applied to each elevation difference map. Note that, because of the same spatial resolution, correction of the curvature-related bias was not performed for the 1974–80 KH-9 DEM pair. Furthermore, no consideration of penetration depth was needed in the process of calculating the mass balance. Figure 5 shows the glacier thickness change patterns of the three periods. During 1974–80, significant thickening and thinning took place in the higher and lower reaches of the Rimo glacier (Fig. 5a), respectively, while during 1980–2000, the locations of significant thickening and thinning switched (Fig. 5b). Apparently, the Rimo glacier surged during the later period. From the thickness change map, the change magnitude of 1974–2000 is between those of 1974–80 and 1980–2000, which seems reasonable. Similar patterns can be observed on the non-surge glaciers (e.g. the Siachen glacier), whose thickness change rates varied slowly.

Fig. 5. Glacier thickness changes in the eastern Karakoram during the different time periods. Left: 1974–1980; center: 1980–2000; right: 1974–2000.

Regarding the evaluation of the internal consistency of glacier mass balance, we employed two kinds of strategies. The first is that we only used the pixels that are valid in all DEMs to calculate the average elevation change (it should be noted that the elevation change here is the result of the initial difference between two DEMs without penetration correction and thus does not represent the real glacier thickness change) and further obtain the uncertainty of the mass balance corresponding to the residual elevation. The final elevation differences (samples: 708 699) are 6.54 ± 14.24 m, −3.27 ± 17.24 m and 1.07 ± 14.25 m for the periods 1974–80, 1980–2000 and 1974–2000 (Table 4), respectively, which results in a residual elevation difference of 2.2 m, corresponding to an uncertainty of mass balance of 0.07 m w.e. a⁻¹. The second kind is that, using the method described in Section 3.5, we directly calculated the mass balance for each period. The results are 3.28 ± 0.36 m w.e. (samples: 998 163), −2.16 ± 1.63 m w.e. (samples: 977 346) and 0.70 ± 1.41 m w.e. (samples: 1 474 187) for 1974–80, 1980–2000 and 1974–2000 (Table 4). After triangulation, the residual mass balance is 0.42 m w.e., corresponding to an uncertainty of 0.02 m w.e. a⁻¹. To summarize, the uncertainty of the glacier mass balance derived from the DEM differencing method is acceptable for glacier change monitoring over periods of decades.

Table 4. Geodetic mass balance and elevation change in the eastern Karakoram (see Fig. 5) for 1974–80, 1980–2000 and 1974–2000

5.2. Comparison with previous studies

Bhambri and others (2013) reported an area change rate of −0.01% a⁻¹ (1973–2002) for the upper Shyok valley (eastern Karakoram), covering an area of 2977.9 km². This slight area change seems consistent with our result in the eastern Karakoram (E-K-E) with a slight mass loss (−0.07 ± 0.05 m w.e. a⁻¹ for 1974–2000). Moreover, Rankl and others (2014) reported that, during 1976–2012, many glaciers in the northwestern margin of the Karakoram obviously receded, while most of the glaciers in the southeastern margin did not. Only a few small glaciers retreated significantly. Furthermore, the advancing glaciers were concentrated in the central Karakoram. These length change patterns are in good agreement with our mass-balance patterns. In terms of mass-balance monitoring, only two studies in the Karakoram region have an observation period partly overlapping with ours. Based on a hydrological method, Zaman and Liu (2015) reported lower and upper bounds for mass balance of −0.23 and 0.22 m w.e. a⁻¹ for the Siachen glacier (Fig. 5 and S1) during 1986–1991. Considering a significant discrepancy in the observation periods, we did not make a direct comparison but just show our result at the Siachen glacier (−0.16 ± 0.08 m w.e. a⁻¹ for 1980–2000). Additionally, Gao and others (2010) derived a negative mass budget of −0.16 m w.e. a⁻¹ in the upstream of the Yarkant River for 1961–2006 via a degree-day glacier ablation model, which is numerically consistent with our result of −0.16 ± 0.07 m w.e. a⁻¹ (1980–2000) in the central Karakoram. What is more, compared with other regions for similar study periods, our Karakoram rate of −0.09 ± 0.03 m w.e. a⁻¹ is much less negative than that of the Central Tien Shan, with a value of −0.35 ± 0.34 m w.e. a⁻¹ for the period 1975–99 based on the same data sources (Pieczonka and Bolch, 2015). According to the study of Vaughan and others (2013), a global average mass balance of −0.31 ± 0.19 m w.e. a⁻¹ for 1971–2009 was also obtained (see Table 4.2 and Table 4.5 in Vaughan and others (2013)), which is more negative than that in the Karakoram. From the above, we can conclude that the Karakoram glaciers have experienced only slight mass loss, compared with the serious mass loss around the world.

5.3. The long-term glacier change anomaly in the Karakoram

We also made a comparison with other geodetic measurements for 2000 to the 2010s (Table 5). Kääb and others (2012) obtained a slightly negative or near-zero glacier mass balance (−0.03 ± 0.04 m w.e. a⁻¹) over the entire Karakoram (covering a glacier area of 21 750 km²) for 2003–09 based on autumn ICESat/GLAS observations. Subsequently, by applying various data processing strategies to ICESat data, Gardner and others (2013) and Kääb and others (2015) reported an updated mass balance of −0.11 ± 0.18 m w.e. a⁻¹ (including Hindu Kush regions with a total of area of 24 700 km²) and −0.09 ± 0.05 m w.e. a⁻¹ (21 000 km²) in the Karakoram, respectively. These results suggest that during the period from 2003 to 2008/09, the glaciers underwent a slight mass loss in the whole Karakoram. At a local scale, a similar result was obtained by Rankl and Braun (2016), who found a glacier mass balance of −0.08 ± 0.10 m w.e. a⁻¹ for 2000–12 in the Shaksgam River valley of the central Karakoram (Fig. 4) with an area of 1107 km². In addition, Agarwal and others (2016) showed a mass balance of −0.03 ± 0.21 m w.e. a⁻¹ during 1999–2007 for the Siachen glacier (937 km²). However, by differencing SPOT5 DEMs and the SRTM DEM, Gardelle and others (2012a, 2013) reported a nearly stable or even slightly positive mass balance of 0.11 ± 0.22 m w.e. a⁻¹ (2000–08) and 0.11 ± 0.14 m w.e. a⁻¹ (2000–10) for the central and eastern Karakoram regions, respectively, covering a total area of ~9700 km². The above results demonstrate a consensus that the Karakoram glaciers were basically stable from 2000 to the 2010s. Combining with our result of weak mass loss (−0.09 ± 0.03 m w.e. a⁻¹) for the period 1973–2000, we can conclude that the glaciers in the Karakoram as a whole were approximately in balance during the past four decades (1970s–2010s). Particularly, this can be shown in the central Karakoram, in which the south part (C-K-W-S) exhibited a slightly negative but insignificant (or near-zero) mass balance (−0.03 ± 0.05 m w.e. a⁻¹ for 1973–2000), which also confirms the conclusion of Herreid and others (2015), in that the glaciers were almost stable between 1977 and 2014, when considering a nearly balanced condition from 2000 to 2010. In addition, Cogley (2016) also reported, if only in terms of the rate of area change, that the ‘Karakoram anomaly’ is a long-lived feature. Therefore, we arrive at the conclusion that the Karakoram anomaly (at least for the central Karakoram) actually dates back to at least as far as the 1970s.

Table 5. Comparisons with region-wide glacier mass-balance measurements in the Karakoram

5.4. Factors of glacier changes

Glacier changes are mainly determined by temperature and precipitation changes. In the Upper Indus Basin (UIB), by analyzing the observations of 17 climate stations, Archer and Fowler (2004) found that the winter precipitation increased with a rate of 7% per decade from 1961 to 1999, while the annual and summer precipitation also showed a significantly rising trend. Meanwhile, Fowler and Archer (2005, 2006) demonstrated that the mean and minimum temperature in summer had a decreasing trend during the same period. Specifically, the average summer temperature dropped by 1.2°C (Fowler and Archer, 2006). Subsequently, similar trends of the precipitation and temperature were observed for the period from 2000 to the 2010s, which accounts for a slight expansion of snow cover during 2000–12 in the UIB (Tahir and others, 2011; Tahir and others, 2015, 2016). Nevertheless, in contrast to the above observations on the south side of the Karakoram mountains, meteorological records (1961–2000) from Tashkurghan station in China (on the north side of the central Karakoram) demonstrated that annual average temperature exhibited a significantly warming trend (with a rate of 0.25°C per decade) and the annual precipitation trend showed no significant change (Chen and others, 2010). This may indicate that climate conditions to the north and south of the Karakoram divide are different. In addition, based on the Global Precipitation Climatology Project (GPCP) datasets, Yao and others (2012) reported an increasing winter precipitation in Pamir and Karakoram ranges for 1979–2010 and attributed it to the strengthened westerlies. This seems to be inconsistent with the findings of Mukhopadhyay and Khan (2014) and Mukhopadhyay and others (2015), who found, by analyzing the same GPCP datasets, a significantly increasing trend of summer precipitation instead of winter precipitation. However, considering the increasing summer and winter precipitation observed at lower elevations, it seems that both the two trends of precipitation derived from the GPCP datasets may be reasonable. The stable or slightly expanded summer snow cover also confirms the patterns of increasing summer precipitation at high altitudes (Tahir and others, 2015, 2016). Furthermore, regarding the atmospheric temperature at higher elevations, the Microwave Sounding Unit (MSU) datasets (1979–2008) indicated that there was an apparent warming trend in winter but a stable or slightly declining trend in summer for the Himalaya-Karakoram-Tibetan Plateau (HK-TP) region (Gautam and others, 2009; Prasad and others, 2009), which agrees with the observations obtained by Fowler and Archer (2005, 2006) and the partial results obtained by Chen and others (2010). Moreover, according to the climate model simulations, Kapnick and others (2014) also reported a stable or increased snowfall. From all of these observations, we can draw the conclusion that in the whole Karakoram, the annual and seasonal precipitation had increased during the 1960s to the 2010s while the summer temperature had slightly dropped, despite an exception in Tashkurghan climate station. Such climate change likely results in strengthened mass accumulation and weakened ablation, which illustrates that there is a nearly balanced state for the Karakoram glaciers during the past four decades (from the 1970s to the 2010s) in combination with the possible warming in the north part of the Karakoram.

Within the Karakoram, the heterogeneity of the observed glacier change patterns possibly can be ascribed to the discrepancy of climate settings and their changes. For example, in the north part of the Batura glacier (in the western Karakoram, see Fig. 3b) the annual precipitation is ~600 mm at ~5000 m a.s.l., while in the south part, the annual precipitation is as high as 1500–1800 mm (Winiger and others, 2005). Furthermore, Mukhopadhyay and others (2015) also reported that the precipitation patterns in different seasons were spatially heterogeneous in the interior of Karakoram ranges. In summer, the amount of precipitation in the central and eastern Karakoram is higher than that in the western Karakoram, whereas an opposite phenomenon occurred in winter. This may be related to the distinction of sub-regional mass changes in the east-west direction (e.g. −0.11 ± 0.04, −0.04 ± 0.05 and −0.07 ± 0.05 m w.e. a⁻¹ in the western, central and eastern Karakoram for 1973/74–2000). Moreover, Chen and others (2010) also revealed a significant increase in summer and autumn temperature, which illustrates that the north part has higher rate of glacier mass wastage (−0.17 ± 0.06 and −0.14 ± 0.05 m w.e. a⁻¹) than that in the south part (−0.08 ± 0.05 and −0.03 ± 0.05 m w.e. a⁻¹).

6. CONCLUSIONS

Our results suggest that the Karakoram glaciers had a slight negative mass balance of −0.09 ± 0.03 m w.e. a⁻¹ from 1973 to 2000, which is in contrast to the severe mass loss seen globally. Within the Karakoram mountains, the patterns of glacier mass balance are heterogeneous both in space and time. During 1973/74–2000, the western Karakoram had a more negative mass balance (−0.11 ± 0.04 m w.e. a⁻¹) than the eastern Karakoram (−0.07 ± 0.05 m w.e. a⁻¹), and the weakest negative mass balance appeared in the central Karakoram (−0.04 ± 0.05 m w.e. a⁻¹). Furthermore, on a more local scale, for central and western Karakoram, both the north parts (−0.17 ± 0.05 and −0.14 ± 0.06 m w.e. a⁻¹) experienced a more negative mass balance than the south parts (−0.08 ± 0.05 and −0.03 ± 0.05 m w.e. a⁻¹), which possibly can be attributed to the differences of local climate changes, especially the rising summer temperature in the north part of the Karakoram. Compared with previous studies, we found the Karakoram anomaly can be dated back to the 1970s, particularly for the central Karakoram. The increasing winter precipitation and the decreasing summer temperature may be the main reasons.