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, Reference Kaser, Großhauser and Marzeion2010). Most glaciers have been retreating since the end of the Little Ice Age (Marzeion and others, Reference Marzeion, Cogley, Richter and Parkes2014). In the recent decades, in particular, nearly worldwide glacier shrinkage and mass loss have been observed (Gardner and others, Reference Gardner2013). Certainly, the glaciers in High Mountain Asia are no exception and exhibit a high rate of area shrinkage (Li and others, Reference Li2008; Cogley, Reference Cogley2016). 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, Reference Benn2012; Rankl and others, Reference Rankl, Kienholz and Braun2014). 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, Reference Hewitt2005; Rankl and others, Reference Rankl, Kienholz and Braun2014). 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, Reference Hewitt2005, Reference Hewitt2007, Reference Hewitt2011; Copland and others, Reference Copland2011; Scherler and others, Reference Scherler, Bookhagen and Strecker2011; Bolch and others, Reference Bolch2012; Bhambri and others, Reference Bhambri2013; Rankl and others, Reference Rankl, Kienholz and Braun2014; Herreid and others, Reference Herreid2015; Paul, Reference Paul2015). 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, Reference Scherler, Bookhagen and Strecker2011; Rankl and others, Reference Rankl, Kienholz and Braun2014). Hence, glacier mass balance, which is a direct and reliable indicator of climate change, must be investigated. Hewitt (Reference Hewitt2005) 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 (Reference Gardelle, Berthier and Arnaud2012a, Reference Gardelle, Berthier, Arnaud and Kääb2013) 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, Reference Kääb, Berthier, Nuth, Gardelle and Arnaud2012, Reference Kääb, Treichler, Nuth and Berthier2015; Gardner and others, Reference Gardner2013). More recently, by subtracting SRTM DEM from TanDEM-X DEM and Cartosat-1 DEM, respectively, Rankl and Braun (Reference Rankl and Braun2016) and Agarwal and others (Reference Agarwal2016) 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 (Reference Herreid2015) 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, Reference Xie and Liu2010). 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 km2 (Bolch and others, Reference Bolch2012), and most of them are spread within the altitude range of 2300–8600 m a.s.l. (Hewitt, Reference Hewitt2011). 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 km2, respectively (Hewitt, Reference Hewitt2013). Due to the frequent snow avalanches, the tongues of the glaciers in this area are generally covered by debris. However, Mayer and others (Reference Mayer, Lambrecht, Belo, Smiraglia and Diolaiuti2006) 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.
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, Reference Gardelle, Berthier, Arnaud and Kääb2013), which bring approximately two-thirds of the high-altitude snowfall (Bolch and others, Reference Bolch2012). The Indian summer monsoon also has a small impact on the eastern margins of the Karakoram region (Bhambri and others, Reference Bhambri2013). Furthermore, the jet stream plays an important role in adjusting the precipitation in the vertical direction (Hewitt, Reference Hewitt2005; Mayer and others, Reference Mayer, Lambrecht, Belo, Smiraglia and Diolaiuti2006). The amount of precipitation at high altitudes is 5–10 times that of the low altitudes (Hewitt, Reference Hewitt2005; Immerzeel and others, Reference Immerzeel, Wanders, Lutz, Shea and Bierkens2015). For example, the precipitation can reach 1200–2000 mm a−1 (mainly in solid form) in the regions near the snowline (~5000 m a.s.l.) (Xie and Liu, Reference Xie and Liu2010), but is only 150–200 mm a−1 in the altitude range of 2000–3200 m a.s.l. (Mayer and others, Reference Mayer, Lambrecht, Belo, Smiraglia and Diolaiuti2006).
Several studies have demonstrated that glacier surges are frequent in the Karakoram region (Copland and others, Reference Copland2011; Rankl and others, Reference Rankl, Kienholz and Braun2014). The updated glacier inventory in Rankl and others (Reference Rankl, Kienholz and Braun2014) lists 101 surge-type glaciers in this area. Generally, surge-type glaciers have a greater length than normal glaciers (Barrand and Murray, Reference Barrand and Murray2006; Rankl and others, Reference Rankl, Kienholz and Braun2014). In fact, the surge-type glaciers should number more than the result of Rankl and others (Reference Rankl, Kienholz and Braun2014), 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, Reference Hewitt and Liu2010; Xie and Liu, Reference Xie and Liu2010; Hewitt, Reference Hewitt2011).
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, Reference Bolch, Menounos and Wheate2010). 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, Reference Paul and Kääb2005; Andreassen and others, Reference Andreassen, Paul, Kääb and Hausberg2008), it tends to misclassify water bodies as glaciers (Bolch and others, Reference Bolch, Menounos and Wheate2010). 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, Reference Surazakov and Aizen2010). 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.
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, Reference Surazakov and Aizen2010). The second step was to perform histogram equalization and locally adaptive filtering in order to enhance the local contrast (Pieczonka and others, Reference Pieczonka, Bolch, Junfeng and Shiyin2013). 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.
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, Reference Nuth and Kääb2011), and its relative accuracies in the horizontal and vertical directions are better than 45 and 10 m, respectively (Marschalk and others, Reference Marschalk, Roth, Eineder and Suchandt2004). 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, Reference Nuth and Kääb2011). 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, Reference Bolch, Pieczonka and Benn2011; Pieczonka and others, Reference Pieczonka, Bolch, Junfeng and Shiyin2013). 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, Reference Rolstad, Haug and Denby2009; Fischer and others, Reference Fischer, Huss and Hoelzle2015). 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, Reference Pieczonka and Bolch2015). Höhle and Höhle (Reference Höhle and Höhle2009) 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 (Reference Pieczonka and Bolch2015), 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, Reference Gardelle, Berthier and Arnaud2012a, Reference Gardelle, Berthier and Arnaudb, Reference Gardelle, Berthier, Arnaud and Kääb2013; Junfeng and others, Reference Junfeng2015). Paul (Reference Paul2008) and Gardelle and others (Reference Gardelle, Berthier and Arnaud2012b) 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 (Reference Gardelle, Berthier and Arnaud2012a) and Gardelle and others (Reference Gardelle, Berthier, Arnaud and Kääb2013), 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, Reference Gardelle, Berthier and Arnaud2012b). 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 (Reference Gardelle, Berthier and Arnaud2012b) (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 (Reference Gardelle, Berthier and Arnaud2012b) 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 (Reference Gardelle, Berthier and Arnaud2012a) and Gardelle and others (Reference Gardelle, Berthier, Arnaud and Kääb2013), 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 (Reference Kääb, Berthier, Nuth, Gardelle and Arnaud2012) 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.
(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, Reference Berthier and Toutin2008; Berthier and others, Reference Berthier, Schiefer, Clarke, Menounos and Rémy2010; Gardelle and others, Reference Gardelle, Berthier, Arnaud and Kääb2013; Fischer and others, Reference Fischer, Huss and Hoelzle2015). 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, Reference Gardelle, Berthier and Arnaud2012a, Reference Gardelle, Berthier, Arnaud and Kääb2013). A glacier average density of 850 ± 60 kg m−3 has been used to convert the volume change into mass budget (Huss, Reference Huss2013). 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, Reference Wake1989; Gardelle and others, Reference Gardelle, Berthier, Arnaud and Kääb2013). 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, Reference Gardelle, Berthier, Arnaud and Kääb2013). 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, Reference Rolstad, Haug and Denby2009; Nuth and Kääb, Reference Nuth and Kääb2011). Subsequently, according to the equations in Gardelle and others (Reference Gardelle, Berthier, Arnaud and Kääb2013), 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−3). In addition, we also took the error of glacier inventory (ΔS) into account. Given that Guo and others (Reference Guo2015) 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 0) using Eqn (1), where V and ρ w denote the change of glacier volume and the water density (1000 kg m−3), 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).
$$\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} $$
$$\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).
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−1 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−1. However, a less negative mass balance of −0.04 ± 0.05 m w.e. a−1 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−1, 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−1 . A similar case was also found in the central Karakoram, for example −0.14 ± 0.05 m w.e. a−1 in the north side (C-K-W-N) versus −0.03 ± 0.05 m w.e. a−1 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−1 and −0.13 ± 0.08 m w.e. a−1, 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−1 for the glacier areas (12 366 km2) covered by the new KH-9 DEMs, which is equivalent to a mass budget of −1.15 ± 0.31 Gt a−1. If we extrapolate the measured regional mass balance to all the glaciers in the Karakoram region (a total area of 22 000 km2), the total glacier mass loss for 1973–2000 would be −51.48 ± 17.16 Gt.
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.
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, Reference Nuth and Kääb2011). 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.
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−1. 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−1. 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.
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 (Reference Bhambri2013) reported an area change rate of −0.01% a−1 (1973–2002) for the upper Shyok valley (eastern Karakoram), covering an area of 2977.9 km2. 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−1 for 1974–2000). Moreover, Rankl and others (Reference Rankl, Kienholz and Braun2014) 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 (Reference Zaman and Liu2015) reported lower and upper bounds for mass balance of −0.23 and 0.22 m w.e. a−1 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−1 for 1980–2000). Additionally, Gao and others (Reference Gao, Zhang, Ye and Qiao2010) derived a negative mass budget of −0.16 m w.e. a−1 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−1 (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−1 is much less negative than that of the Central Tien Shan, with a value of −0.35 ± 0.34 m w.e. a−1 for the period 1975–99 based on the same data sources (Pieczonka and Bolch, Reference Pieczonka and Bolch2015). According to the study of Vaughan and others (Reference Vaughan, Stocker, Qin, Plattner, Tignor, Allen, Boschung, Nauels, Xia, Bex and Midgley2013), a global average mass balance of −0.31 ± 0.19 m w.e. a−1 for 1971–2009 was also obtained (see Table 4.2 and Table 4.5 in Vaughan and others (Reference Vaughan, Stocker, Qin, Plattner, Tignor, Allen, Boschung, Nauels, Xia, Bex and Midgley2013)), 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 (Reference Kääb, Berthier, Nuth, Gardelle and Arnaud2012) obtained a slightly negative or near-zero glacier mass balance (−0.03 ± 0.04 m w.e. a−1) over the entire Karakoram (covering a glacier area of 21 750 km2) for 2003–09 based on autumn ICESat/GLAS observations. Subsequently, by applying various data processing strategies to ICESat data, Gardner and others (Reference Gardner2013) and Kääb and others (Reference Kääb, Treichler, Nuth and Berthier2015) reported an updated mass balance of −0.11 ± 0.18 m w.e. a−1 (including Hindu Kush regions with a total of area of 24 700 km2) and −0.09 ± 0.05 m w.e. a−1 (21 000 km2) 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 (Reference Rankl and Braun2016), who found a glacier mass balance of −0.08 ± 0.10 m w.e. a−1 for 2000–12 in the Shaksgam River valley of the central Karakoram (Fig. 4) with an area of 1107 km2. In addition, Agarwal and others (Reference Agarwal2016) showed a mass balance of −0.03 ± 0.21 m w.e. a−1 during 1999–2007 for the Siachen glacier (937 km2). However, by differencing SPOT5 DEMs and the SRTM DEM, Gardelle and others (Reference Gardelle, Berthier and Arnaud2012a, Reference Gardelle, Berthier, Arnaud and Kääb2013) reported a nearly stable or even slightly positive mass balance of 0.11 ± 0.22 m w.e. a−1 (2000–08) and 0.11 ± 0.14 m w.e. a−1 (2000–10) for the central and eastern Karakoram regions, respectively, covering a total area of ~9700 km2. 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−1) 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−1 for 1973–2000), which also confirms the conclusion of Herreid and others (Reference Herreid2015), in that the glaciers were almost stable between 1977 and 2014, when considering a nearly balanced condition from 2000 to 2010. In addition, Cogley (Reference Cogley2016) 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.
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 (Reference Archer and Fowler2004) 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 (Reference Fowler and Archer2005, Reference Fowler and Archer2006) 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, Reference Fowler and Archer2006). 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, Reference Tahir, Chevallier, Arnaud and Ahmad2011; Tahir and others, Reference Tahir, Chevallier, Arnaud, Ashraf and Bhatti2015, Reference Tahir, Adamowski, Chevallier, Haq and Terzago2016). 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, Reference Chen, Xu, Chen, Li and Liu2010). 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 (Reference Yao2012) 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 (Reference Mukhopadhyay and Khan2014) and Mukhopadhyay and others (Reference Mukhopadhyay, Khan and Gautam2015), 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, Reference Tahir, Chevallier, Arnaud, Ashraf and Bhatti2015, Reference Tahir, Adamowski, Chevallier, Haq and Terzago2016). 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, Reference Gautam, Hsu, Lau, Tsay and Kafatos2009; Prasad and others, Reference Prasad, Yang, El-Askary and Kafatos2009), which agrees with the observations obtained by Fowler and Archer (Reference Fowler and Archer2005, Reference Fowler and Archer2006) and the partial results obtained by Chen and others (Reference Chen, Xu, Chen, Li and Liu2010). Moreover, according to the climate model simulations, Kapnick and others (Reference Kapnick, Delworth, Ashfaq, Malyshev and Milly2014) 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, Reference Winiger, Gumpert and Yamout2005). Furthermore, Mukhopadhyay and others (Reference Mukhopadhyay, Khan and Gautam2015) 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−1 in the western, central and eastern Karakoram for 1973/74–2000). Moreover, Chen and others (Reference Chen, Xu, Chen, Li and Liu2010) 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−1) than that in the south part (−0.08 ± 0.05 and −0.03 ± 0.05 m w.e. a−1).
6. CONCLUSIONS
Our results suggest that the Karakoram glaciers had a slight negative mass balance of −0.09 ± 0.03 m w.e. a−1 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−1) than the eastern Karakoram (−0.07 ± 0.05 m w.e. a−1), and the weakest negative mass balance appeared in the central Karakoram (−0.04 ± 0.05 m w.e. a−1). 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−1) experienced a more negative mass balance than the south parts (−0.08 ± 0.05 and −0.03 ± 0.05 m w.e. a−1), 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.
SUPPLEMENTARY MATERIAL
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2016.142
ACKNOWLEDGEMENTS
We thank the USGS for supplying the C-band SRTM data, the Landsat images and the KH-9 images, some of which were purchased. We also thank the German Aerospace Center (DLR) for providing the SRTM X-band data. The work in the paper was supported by the Nature Science Foundation of China (41474007, 41404013, 41404011, 41574005, 41222027), the Hunan Provincial Foundation of China (No. 13JJ1006), the Project of Innovation-driven Plan of Central South University (Grant No. 2016CX004) and the Fundamental Research Funds for the Central Universities of Central South University (Grant No. 2016zzts083).
Open access
