1. Introduction
Since the beginning of this century, most glaciers have retreated worldwide at an increasing rate (Dehecq and others, Reference Dehecq2019). However, there are important regional differences in glacier changes and behavior in High Mountain Asia (HMA). The Karakoram Range has been identified as a region where glaciers have remained stable or even advanced since the late 1990s (Hewitt, Reference Hewitt2005). Gardelle and others (Reference Gardelle, Berthier, Arnaud and Kääb2013) indicated that the ‘Karakoram Anomaly’ should include the Pamir regions and be renamed as the ‘Pamir-Karakoram Anomaly’ because of the similar mass gain during early 21st century in the central Karakoram (0.10 ± 0.16 m w.e. a−1) and western Pamir (0.14 ± 0.13 m w.e. a−1). Recent studies have shown that anomalous glacier behavior probably occurs in the nearby west Kunlun Mountains and Pamir (Kääb and others, Reference Kääb, Treichler, Nuth and Berthier2015; Brun and others, Reference Brun, Berthier, Wagnon, Kääb and Treichler2017; Berthier and Brun, Reference Berthier and Brun2019; Farinotti and others, Reference Farinotti, Immerzeel, de Kok, Quincey and Dehecq2020), the extent of mass gain could even reach to the central of the Tibet Plateau from 2003 to 2009 (Neckel and others, Reference Neckel, Kropáček, Bolch and Hochschild2014). Changes in ice flow velocity were shown to be unrelated to the surge-type glaciers in the Karakoram and the west Kunlun Mountains, but there was a signal of increased ice deformation and sliding due to glacier thickening and associated mass gaining (Dehecq and others, Reference Dehecq2019). This accelerating ice velocity contrasted with the other parts of HMA, where the glacier flow decreased. Farinotti and others (Reference Farinotti, Immerzeel, de Kok, Quincey and Dehecq2020) suggested more mass gain in the west Kunlun Mountains compared to the results in the Karakoram. Negi and others (Reference Negi, Kumar, Kanda, Thakur and Singh2021) also reported that there was no glacier anomaly in the eastern Karakoram (−0.10 ± 0.07 m w.e. a−1) in the early 21st century. However, few studies have focused on the exact extent of this anomaly.
The Kunlun Mountains in the northern sector of the Tibetan Plateau are located in the center of HMA. The remoteness and harsh environmental conditions of the Kunlun Mountains have prevented the collection of in situ glacier observations. This impedes our understanding of spatial and temporal changes of glaciers in this region. Satellite data obtained from a wide range of sensors provide data for studying changes in these remote glaciers (Paul and others, Reference Paul2015). So far, most studies about glacier changes in the Kunlun Mountains concentrate on the west Kunlun Mountains. Using Landsat 7 TM/ETM+ and the Chinese glacier inventory data, Shangguan and others (Reference Shangguan2007) revealed that 278 glaciers (with a total area of 2711.57 km2) in the west Kunlun Mountains reduced in area by only ~0.40% from the 1970s to 2001. Comparing two Chinese glacier inventory datasets, Bao and others (Reference Bao, Liu, Wei and Guo2015) observed the glacier area (283 with a total area of 2822.04 km2) in the west Kunlun Mountains decreased by 3.4 ± 3.1% from 1970 to 2010. Based on ICESat/GLAS data and Shuttle Radar Topography Mission (SRTM) elevation data, glacier surface elevation decreased for glaciers located below 5800 m a.s.l. but increased for glaciers situated above 6000 m a.s.l. (Wu and others, Reference Wu, Wang, Guo and Wu2014), and these glaciers gained mass during 2003–2009 (Bao and others, Reference Bao, Liu, Wei and Guo2015) in the west Kunlun Mountains. Recently, Wang and others (Reference Wang2018), using multiple DEM sources, estimated average glacier mass change of −0.06 ± 0.13 m w.e. a−1 in the west Kunlun Mountains between 1970s and 1999. However, Zhou and others (Reference Zhou, Li, Li, Zhao and Ding2018) inferred a nearly stable or slight mass gain state in the west Kunlun Mountains from 1970s to 2000. Differences between these studies are related partly to the variations in the extent of the regions considered in the analysis. The presence of surge-type glaciers in west Kunlun, confirmed by observations of ice velocity and terminal changes (Yasuda and Furuya, Reference Yasuda and Furuya2015; Chudley and Willis, Reference Chudley and Willis2018), further complicates studies of mass change in the region. To our knowledge, studies of glacier change across the entire Kunlun Mountain region have not been carried out to date. Hence, there is a need to systematically analyze the spatial distribution of glacier change in the entire Kunlun Mountains to investigate the extent of the mass-balance anomaly.
In order to find the eastern limit of the Kunlun-Pamir-Karakoram Anomaly, this study quantifies changes in glacier area and surface elevation across the Kunlun Mountains. We examine spatial and temporal changes in glacier area between three glacier inventories (1970s to 2016) and calculate surface elevation and mass changes for over 70% of the glaciers in the region between 2000 and 2013.
2. Study area
The Kunlun Mountains (33°–38°N, 74°–100°E) are situated in the northern part of the Tibetan Plateau (Fig. 1) and span ~2000 km from west to east. The main ridge is above 5500 m a.s.l. The climate in the west and the central Kunlun Mountains is controlled mostly by westerlies nearly all year, while the climate in the east Kunlun Mountains is affected by both westerlies in the wintertime and, occasionally, by Asian monsoon systems in the summertime (Shi, Reference Shi2002; Tian and others, Reference Tian2007; Yao and others, Reference Yao2013; Ma and others, Reference Ma, Guan, Guo, Gan and Xie2017; Thompson and others, Reference Thompson2018). Although the whole mountain range is situated in a cold and dry region, annual precipitation can reach ~200 mm at 5000 m a.s.l., and 400 mm at 5900 m a.s.l. (Zhang and Li, Reference Zhang and Li1999). The high terrain of the Kunlun Mountains results in a large number of glaciers. According to the second Chinese Glacier Inventory (CGI-2) (Liu and others, Reference Liu, Guo and Xu2012), there are 8452 glaciers with a total area of 11 099 ± 400 km2 in the Kunlun Mountains, with numerous large glaciers in the west and central Kunlun and small ones in the east (Fig. 1). Most of the glaciers are characterized as the continental type.
Fig. 1. Distribution of studied glaciers in the Kunlun Mountains and meteorological stations (pink triangles) around the study area. Pink rectangles are the outlines of the TanDEM-X DEMs used to calculate the elevation change, covering nearly 70% of the glacier area in the Kunlun Mountains, and black dashed lines are the dividing lines for the west, the central and the east Kunlun Mountains for the estimation of glacier area change.
3. Data and methods
3.1 Data sources
3.1.1 Landsat OLI images and previous glacier inventories
Landsat 8 Operational Land Imager (OLI) images (N = 38) collected between 2014 and 2017 are used to delineate recent glacier extents, with the majority of the imagery (79%) collected between 2015 and 2016 (Supplementary Table S1). Two Chinese Glacier Inventories from the 1970s (CGI-1) and 2006–2010 (CGI-2) are compared with the inventory created in this work to assess glacier area change in the Kunlun Mountains from the 1970s to 2016. Different dates are assigned to glacier boundaries delineated for the eastern, central and western parts of the Kunlun Mountains in CGI-1 and CGI-2 (Supplementary Table S2). CGI-1 was compiled with topographical maps at 1:50 000 and 1:100 000 scales, aerial photos and field validation, while CGI-2 was based on Landsat TM images with a band-ratio threshold method together with visual interpretation (Guo and others, Reference Guo2015).
3.1.2 SRTM and TanDEM-X DEM
In this study, the non-void filled NASA SRTM DEM (Version 1; https://dx.doi.org//10.5066/F7K072R7) and the TanDEM-X DEM 90 m from TanDEM-X mission are used to calculate glacier surface elevation change between 2000 and 2013. The SRTM DEMs are processed from InSAR data acquired from 11 to 22 February 2000, with WGS84, coordinate system set to EGM96 geoid (Farr and Kobrick, Reference Farr and Kobrick2000). The nominal absolute vertical accuracy of SRTM DEM is ~15 m at the 90% confidence level, equal to an RMSE ~10 m (Wan and others, Reference Wan, Wang, Hou, Huai and Liu2020). The TanDEM-X DEMs used in this study are processed from InSAR data acquired in winter 2013 (Supplementary Table S3). The height of TanDEM-X DEM is ellipsoidal height, which is the sum of the ellipsoidal height and orthometric height.
Absolute vertical accuracy of TanDEM-X DEM is below 10 m at the 90% confidence level (https://geoservice.dlr.de/web/dataguide/tdm90/). There are slightly different observation periods within and between regions due to TanDEM-X coverage from different time periods and regions. This may cause some uncertainty that is difficult to directly quantify, although it can be considered small for an observational period of ~13 years (Braun and others, Reference Braun2019).
3.2 Methods
3.2.1 Glacier boundary delineation
Visual interpretation and automatic classification are generally used to extract glacier boundaries. Compared to visual interpretation, automatic classification is suitable to obtain glacier boundaries over large areas. However, for glaciers covered by debris, manual validation and correction are needed (Xiang and others, Reference Xiang2013). Considering additional modifications and adjustments are required after automatic classification, we extracted glacier boundaries manually for the 2016 inventory. Images with snow and cloud coverage <10% were selected preferentially. Images with >10% overall cloud coverage but no cloud in glacierized areas were also considered. In regions where glaciers were persistently influenced by snow and cloud in 2016, we chose high-quality images from adjacent years. All images were projected onto the WGS84 UTM (Zone 43N) coordinate system and resampled to 15 m resolution, and SWIR1, Green, and Blue were used for band synthesis (Paul and others, Reference Paul and Svoboda2009). The boundaries of debris-covered glaciers were identified by using a combination of many features, such as the water ponds and ice cliffs on glacier surface, elevation profiles of valley bottoms and the visible drainage system at glacier termini. All those features could be verified in Google Earth (see detailed extraction examples in Supplementary Fig. S1). After boundary delineation, the area of each glacier was calculated on the WGS84 spheroid and Albers equal-area map projection centered on 95°E, 30°N with standard parallels at 15°N and 65°N (Ye and others, Reference Ye2017). Because of the limitations of the image resolution, glaciers with areas of <0.01 km2 cannot be reliably discriminated from Landsat images (Paul and others, Reference Paul and Svoboda2009; Leigh and others, Reference Leigh2019) and were excluded from this study (<36 pixels). Furthermore, glaciers in CGI-2, which were used to examine glacier surface elevation changes between 2000 and 2013, were modified with 2000 Landsat 7 ETM+ images by visual interpretation.
3.2.2 Geodetic method
Multi-temporal DEMs have been extensively used to calculate the geodetic glacier mass balance. In this study, SRTM DEM and TanDEM-X elevation data were used to estimate the geodetic glacier mass balance in the Kunlun Mountains. All the DEMs were projected to the coordinate system of WGS84 UTM zone 43 N/EGM96 and were resampled bilinearly to 30 m to reduce the impact of different resolutions. TanDEM-X DEMs were translated to the orthometric height using global EGM96 geoid model (https://cddis.nasa.gov/926/egm96/egm96.html) before co-registration. In order to obtain optimal results, the TanDEM-X DEMs were co-registered vertically to the SRTM DEM on stable ground using an iterative approach (Nuth and Kääb, Reference Nuth and Kääb2011). Stable regions were selected in the non-glacier areas with <15° slope (results from SRTM DEM), without vegetation and water.
Statistical methods were utilized to minimize errors from spatial matching and data resolution. The non-glacierized regions were generally considered stable over decades, which could help to assess the matching bias of spatial de-correlation (Berthier and others, Reference Berthier2007). The shift distance and std dev. of the TanDEM-X DEM are shown in Supplementary Table S4. In order to minimize error in the elevation change estimation, values above ±100 m (usually around data gaps and near DEM edges) were removed assuming these values were outliers (Berthier and others, Reference Berthier, Schiefer, Clarke, Menounos and Rémy2010). Gardelle and others (Reference Gardelle, Berthier and Arnaud2012) demonstrated the correlation between elevation difference of different DEMs and maximum surface curvature; the relationship between glacier and non-glacier regions showed a robust consistency, which was used to correct the residual error of elevation difference in our study. Detailed formulas and fitted equations were presented in the Supplementary material and Figure S2.
The elevation difference between two DEMs represented the volume change of the glacier surface, assuming that the density profile remained constant and only ice was lost or gained (Paterson, Reference Paterson1994; Zemp and others, Reference Zemp2010). To transform volume change into mass change, we assume a density of 850 kg m−3 and an additional uncertainty of 60 kg m−3 on account of the lack of measured data (Zemp and others, Reference Zemp2010). Mass balance was estimated by a density conversion formula (1), which could be expressed as follows:
where B is mass balance, ρ is transition density, S g is glacier area, n is the number of pixels in the glacier region, Δh i represents the surface elevation difference of i pixel in different times corresponding to the relevant DEMs and S p represents the area of a pixel.
3.2.3 Meteorological data processing
Monthly air temperature and precipitation observations from meteorological stations (Fig. 1) near to the Kunlun Mountains are used to calculate the variations of temperature and precipitation from 1960 to 2015. The observation data are derived from the National Meteorological Information Center, China Meteorological Administration. Considering the scarceness of meteorological stations around the Kunlun Mountains, especially around the west Kunlun, ERA5 reanalysis data are used to analyze the temporal and spatial distribution of temperature and precipitation variations. The data are the fifth-generation meteorological dataset of the European Centre for Medium-Range Weather Forecasts (ECMWF), providing data on the global weather and climate for the past 40–70 years, which has much higher temporal resolution and more accurate atmospheric conditions than the previous ERA versions and perform well in terms of precipitation in high altitude compared with some other reanalysis products (Hersbach and others, Reference Hersbach2020; Huai and others, Reference Huai, Wang, Sun, Wang and Zhang2021). ERA5-Land monthly averaged data from 1980 to 2020 are used to calculate the summer and winter precipitation changes. The average elevation of glaciers in the Kunlun Mountains is ~5000 m, so the ERA5 monthly averaged data on 500 hPa level (corresponding to ~5500 m a.s.l.) are selected for analyzing upper air temperature change in the glacierized area. Before using ERA5 reanalysis data, we evaluated it with meteorological station observations (see section 4.3).
3.3 Accuracy assessment
3.3.1 Glacier area
The error of co-registration (E M) can be ignored as the Landsat 8 images have already been orthorectified by the USGS with high accuracy (Guo and others, Reference Guo, Liu and Wei2013). Hence, the error of boundary identification can be related to image resolution, which was estimated by creating a buffer of half-pixel for the Landsat images and 15 m for the CGI-1 as the discriminant error (Paul and others, Reference Paul2013). The error of the glacier area can be calculated as formula (2):
where N is the number of pixels of glacier boundary and S is the resolution of the images.
The final uncertainty of glacier area change was estimated using the methods by Guo and others (Reference Guo, Liu and Wei2013). The errors (E C) are regarded as follows according to the error propagation (3):
where E 1 and E 2 represent the error of glacier area in each period.
3.3.2 Mass balance
The relative uncertainties of the DEMs before and after the adjustments are calculated on the basis of glacier-free terrain (derived from SRTM DEM). We calculated the uncertainty of the glacier elevation changes (U Δh) using the std dev. of the averaged elevation change (STDΔh) (Wan and others, Reference Wan, Wang, Hou, Huai and Liu2020), the number of values (after removing spatial autocorrelation) in the non-glacier area (Nuth and Kääb, Reference Nuth and Kääb2011; Paul and others, Reference Paul2015) and the averaged absolute difference (AAD) between the median elevation change on the glacier and off-glacier terrain (Berthier and Brun, Reference Berthier and Brun2019) in the formula (4):
where N eff is the number of pixels calculated. Due to the strong spatial autocorrelation, the influence of this should be removed. N eff is estimated by the formula (5) as following:
where N is the total number of pixels in the stable regions, P is the pixel resolution and d is the distance of spatial autocorrelation. We choose a de-correlation distance of 600 m (20 pixels) for the TanDEM-X DEM to minimize the effect of auto-correlation (Nuth and Kääb, Reference Nuth and Kääb2011).
For the uncertainty arising from glacier boundary, a previous study (Paul and others, Reference Paul2013) suggested a 3% error for alpine regions and a ratio between the glacier perimeter and area (R P/A) of 5.03 km−1. Braun and others (Reference Braun2019) observed R P/A values of individual regions range from 0.67 to 10.81 km−1, and they scaled the error and applied a larger area measurement error in valley glaciers. In this study, the uncertainty of glacier extent was calculated by the formula (6) proposed by Braun and others (Reference Braun2019):
where U A is defined as the uncertainty of glacier area and R P/A the ratio between the perimeter and area.
Additionally, since the radar data of C-band and X-band are different from each other, it is necessary to compensate for the signal of different radar penetration depth in snow and ice (Neelmeijer and others, Reference Neelmeijer, Motagh and Bookhagen2017). Previous studies (Lin and others, Reference Lin, Li, Cuo, Hooper and Ye2017; Li and others, Reference Li, Jiang, Liu and Wang2021) showed that the spatial variability of C-band radar penetration was large, and the radar penetration of X-band in high altitude regions was too large which limits their use for deriving annual geodetic mass balances. Millan and others (Reference Millan, Dehecq, Trouvé, Gourmelen and Berthier2015) observed the penetration depth of TanDEM-X could reach 5–6 m in region with a higher elevation (>3500 m) in October, and a maximum value (7 m) above 2500 m in February due to the dry snow condition.
Thus, we corrected the C-band penetration according to the results from Li and others (Reference Li, Jiang, Liu and Wang2021): the penetration depth of each SRTM scene was set to the maximum at the head of the glacier, and a correction for C-band radar penetration was applied to each altitude bin with a linear model. For the X-band penetration, we used the same correcting method as C-band and choose 7 m as the maximum value. This ensured methodological consistency in the study and the results could be an upper bound for the elevation change. Std dev. was used to estimate the uncertainty (U pen) of penetration depth.
When using elevation differencing to calculate glacier mass change, a final possible source of error is the timing of each acquisition during the mass-balance year. Huintjes and others (Reference Huintjes, Neckel, Hochschild and Schneider2015) explained that glacier mass balance is dominated by snowfall, sublimation, refreezing and snow/ice melt in general. For a glacier of summer-accumulation type in the Tibetan Plateau, snowfall as well as melt occur mainly in the ablation season (Maussion and others, Reference Maussion2014). In addition, during the dry and cold accumulation season (both winter and spring), the slight mass gain caused by snowfall might be largely consumed by the mass loss related to sublimation, while the small amounts of melting snow can be offset by the re-freezing process (Huintjes and others, Reference Huintjes, Neckel, Hochschild and Schneider2015; Li and others, Reference Li, Yao, Yang, Yu and Zhu2018). Liu and others (Reference Liu2019) considered that the seasonal variation between late winter and spring could be ignored when calculating annual mass balance. Hence, we did not take into account the impact of seasonal differences over a 13-year study period because the uncertainty assessment of mass balance should capture this impact.
The total uncertainty of mass balance in this study was estimated using formula (7) as Braun and others (Reference Braun2019). We considered the error from glacier elevation change, glacier boundaries, uncertainty from the conversion density and radar signal penetration:
where ΔM/Δt represents the mass-balance estimate, A the glacier area and ρ the density for volume-to-mass conversion.
4. Results and discussion
4.1 Glacier area changes
Glacier area changes of the different regions in the Kunlun Mountains in the past 50 years are summarized in the Supplementary materials (Tables S2, S5). The year of the previous inventories was determined by the date of images used, and glaciers marked high uncertainty in the CGI-1 were removed from the study. After removing glaciers marked as high-uncertainty in CGI-1 from the analysis, 7232 glaciers in the Kunlun Mountains, with a total area of 11 558 ± 375 km2 in the 1970s, were considered. From 1970s to 2016, a total glacier area change of −593 ± 424 km2 (−5.41 ± 4%) was observed over the Kunlun Mountains. The maximum absolute glacier area change occurred in the west Kunlun Mountains (−455 ± 313 km2), and the minimum occurred in the east (−60 ± 12 km2). These spatial patterns are partly related to the glacier size, as average glacier areas in the West Kunlun Mountains are larger than those in the east. Comparing the relative rate of glacier area change in the sub-regions of the Kunlun Mountains before and after CGI-2 (Fig. 2; Supplementary Table S2), we observe consistent and strongly negative rates of area loss in the eastern Kunlun Mountains (−0.43 ± 0.13% a−1 between 1970s and 2016), large increases in the rate of area loss in the central Kunlun Mountains (−0.01% a−1 prior CGI-2, and −0.24% a−1 after CGI-2), and consistent but relatively lower rates of area loss in the western Kunlun Mountains (−0.12% a−1 overall).
Fig. 2. Distribution of glacier area change rate in the Kunlun Mountains between CGI-1 to -2 (1970 to 2005–10) and CGI-2 to our inventory (2005–10 to 2016), the sub-regions are divided according to the basin and distribution of CGI-1 glacier to avoid the influence from position bias. Red values represent glacier area change rate between CGI-1 to -2, and blue ones are in the period of CGI-2-2016.
Glacier area changes aggregated over a $0.1 \times 0.1^\circ$
grid indicated that between CGI-2 and our 2016 inventory glaciers in the Kunlun Mountains decreased overall but there was spatial heterogeneity (Fig. 3 and Supplementary Fig. S3). The shrinkage of glaciers in the west was small (most <0.40% a−1), while in the east it was larger (>1.20% a−1), and the shrinkage of small glaciers in the west (most <0.5% a−1) was much smaller than it is in the east (>1.50% a−1). From 83°E to the east, the rate of glacier area shrinkage began to increase after the CGI-2 inventory (~2006), and the slope of glacier area shrinkage rate reached the maximum near 85°E (Fig. 4).
Fig. 3. Glacier area change from CGI-2 to our inventory (2005–10 to 2016) in Kunlun Mountains, aggregated over 0.1° grid cells.
Fig. 4. Longitudinal distribution of glacier area change rate before and after CGI-2 in 1° grids, black dash line represents the transition line. Red and blue dash lines represent the piecewise linear fitting of the CGI-1 to -2 and CGI-2 to 2016, respectively.
In order to explore the similarities of glacier changes between the Kunlun Mountains and the Karakoram and Pamir regions, we compared the glacier shrinkage rate in the Kunlun Mountains with the results reported for other regions on the Tibet Plateau over the past several decades (Supplementary Table S6). Glacier area decreased slightly in the western Karakoram by 0.03% a−1 from 1973 to 2014 (Qureshi and others, Reference Qureshi, Yi, Xu and Li2017). The rates of glacier shrinkage in the west Kunlun Mountains and central Kunlun Mountains are lower than those in the eastern Pamir and central Himalaya, and the value of eastern Pamir is slightly higher than that in the west Kunlun Mountains. The smaller glacier area reduction in these regions is different from those in other ranges of the Tibet Plateau, which is considered as the ‘Anomaly’ in the HMA.
4.2 Glacier surface elevation and mass-balance changes
The average surface elevation change of glacier in the Kunlun Mountains was 0.08 ± 0.14 m a−1 between 2000 and 2013. Assuming a surface density of 850 ± 60 kg m−3, the corresponding average mass change was 0.07 ± 0.12 m⋅w.e. a−1 during the period, which is consistent with the results (0.06 ± 0.17 m w.e. a−1) from Shean and others (Reference Shean2020). Our results confirm the mass-balance anomaly observed in the Kunlun Mountains and identify the longitude where surface elevation and mass changes transition from positive to negative values, moving west to east. The most positive elevation change of 0.15 ± 0.35 m a−1 (0.13 ± 0.30 m w.e. a−1) was observed in the west Kunlun Mountains (Fig. 5; Table 1), and was similar to the mass balance (0.14 ± 0.08 m w.e. a−1) estimated by Brun and others (Reference Brun, Berthier, Wagnon, Kääb and Treichler2017). The most negative change was found in the eastern Kunlun Mountains in Aemye Ma-chhen Range of −0.51 ± 0.18 m a−1, which is similar to the value of −0.68 ± 0.29 m a−1 calculated by Shean and others (Reference Shean2020) (Supplementary Table S7).
Fig. 5. Glacier elevation changes from 2000 to 2013. The sub-regions (red dash lines) segmentation is according to the glacierized area. The spatial distribution of glacier elevation change in 0.1° grid, and the spatial distribution of averaged glacier elevation change in the numbered sub-regions.
Table 1. The surface elevation change and mass balance of the sub-regional glaciers in the Kunlun Mountains over the period of 2000–2013
The number of each sub-region refers to Figure 6.
Glacier surface elevation change in four glacierized regions reflected typical characteristics of surge-type glaciers, namely lowering of the reservoir zone and thickening of the receiving zone (Neckel and others, Reference Neckel, Kropáček, Bolch and Hochschild2014) (Fig. 6). There were 38 advancing glaciers in the Kunlun Mountains from CGI-2 to 2016. Figure 5a shows that most glaciers were in a state of mass gain in the accumulation zone, while four glaciers surged from 2000 to 2013. In the Ulugh Muztagh (Fig. 6b), there was one surging glacier and most glaciers were in a state of mass loss. In the Bukatage Mountains (Fig. 6c) glaciers on northern slopes showed mass loss, but several glaciers on the southern slopes showed mass gains in their reservoir zones. The Monuomaha Glacier (GLIMS-id G091032E36060N) surged over 1 km from the 2009 (CGI-2) to 2016 (this inventory; Supplementary Table S8). Glaciers in the Aemye Ma-chhen Range (Fig. 6d) showed clear mass loss, and no surging glaciers were detected.
Fig. 6. Glacier surface elevation change in the west Kunlun peak (a), Ulugh Muztagh (b), Bukatage Mountains (c) and Aemye Ma-chhen Range (d) in the Kunlun Mountains over the period of 2000–2013.
For further analysis of the spatial distribution of surface elevation change, we provided the results averaged over 0.1 × 0.1° grids (Fig. 5). The surface elevation change was positive (~0.15 ± 0.35 m a−1) in the west but negative (~−0.51 ± 0.18 m a−1) in the east, with a transition in the western part of the central Kunlun Mountains near 85°E (Fig. 7). Glaciers in the west Kunlun Mountains (Sub-region 1) were in a state of mass gain (0.13 ± 0.30 m w.e. a−1, Table 1), consistent with previous values (Supplementary Table S9). A mass balance of 0.12 ± 0.17 m w.e. a−1 was found in the central Kunlun Mountains (Sub-region 2), with clear mass loss (−0.09 ± 0.05 m w.e. a−1) in the eastern Kunlun (Sub-region 3). Notably, glaciers in the Bukatage Mountains (Sub-region 4) were on average stationary (−0.00 ± 0.07 m w.e. a−1), and this may be linked to the surging glaciers in this region. Glaciers in the eastern-most sub-region (Aemye Ma-chhen Range) experience substantial mass loss (−0.43 ± 0.15 m w.e. a−1).
Fig. 7. Longitudinal distribution of the surface elevation change of glaciers with different sizes in the Kunlun Mountains over the period of 2000–2013.
Nevertheless, Zhou and others (Reference Zhou, Li, Li, Zhao and Ding2018) estimated mass balance in the Tibet Plateau and its surroundings areas from the mid-1970s to 2000 using Hexagon KH-9 and SRTM DEMs, and found mass balances of 0.02 ± 0.10 and −0.06 ± 0.12 m w.e. a−1 in the west Kunlun Mountains and Ulugh Muztagh, respectively. Our results indicated that mass gain in the west Kunlun Mountains increased between the mid-1970s to 2000 and 2000 to 2013, and mass loss in the Ulugh Muztagh accelerated. These can explain the different spatial distribution of the area shrinkage rate in the Kunlun Mountains, i.e. slightly decelerating in the west, slightly accelerating in the central and largely reduced in the east (Fig. 2). The area of glaciers (7187.15 km2) used to calculate the surface elevation change accounts for nearly 70% of the total glacier area (10 965.07 km2) in the whole Kunlun Mountains. Therefore, the state of mass balance should represent the variation in the entire mountains to a large extent. We inferred that the west Kunlun Mountains belong to the domain of the Karakoram Anomaly (Hewitt, Reference Hewitt2005) and the eastern limit of the ‘Kunlun-Pamir-Karakoram Anomaly’ was near 85°E.
4.3 Effect of air temperature and precipitation variations
Global warming has led to a general shrinkage of glacier area, but there is a clear spatial heterogeneity in the mass balance in different regions. Glacier change is a reflection of climate fluctuations: summer temperature and solid precipitation are the most important drivers of glacier mass balance in the HMA (Wang and others, Reference Wang, He, Pu, Jiang and Jing2010; Bao and others, Reference Bao, Liu, Wei and Guo2015; Kang and others, Reference Kang, Guo, Zhong and Xu2020; Bhattacharya and others, Reference Bhattacharya2021). Analysis of measurements of meteorological stations indicates that mean summer temperature shows a significant upward trend from 1960 to 1999. However, the upward trends are not significant over the period 2000–2015. Furthermore, the upward trends of annual precipitation from 2000 to 2015 are more significant than the value during 1960–1999 (Supplementary Fig. S4), which might cause the recent anomaly mass balance in Kunlun Mountains. However, the meteorological stations near Kunlun Mountains are far away from the glacier area, especially in the west Kunlun Mountains. Therefore, ERA5 reanalysis data were used to analyze the spatial distribution of climate change in the glacier area. Air temperature and precipitation from Wudaoliang station with an elevation near 5000 m a.s.l. near the Kunlun Mountains are used to evaluate the adaptability of ERA5 reanalysis data and the results show that their trends are consistent over the period of 1980–2015 (see in Supplementary Figs S5, S6). So, we use the ERA5 reanalysis data to draw a picture of climate change in our study area and analyze the possible driving factors of the spatial heterogeneity of glacier mass balance. There is no trend in winter precipitation, but an overall increasing trend in summer precipitation. Generally, more than 70% of annual precipitation occurs in the summer in the Kunlun Mountains. The summer temperature increased strongly in the east Kunlun Mountains while very slightly in the west Kunlun Mountains (Fig. 8), which is similar to the recent anomalous cooling of summer air temperature in Karakoram (Forsythe and others, Reference Forsythe, Fowler, Li, Blenkinsop and Pritchard2017). Furthermore, mass balance of glaciers in the east part is more sensitive to temperature than those in the west part (Sakai and Fujita, Reference Sakai and Fujita2017). Hence, the large mass loss in the east Kunlun Mountains is mainly caused by the obviously increasing summer temperature. Glacier mass balance in the dry, cold regions is mostly influenced by the solid precipitation (Bhattacharya and others, Reference Bhattacharya2021). Therefore, anomaly mass gain in the west Kunlun Mountains is probably influenced by the increasing trend of summer precipitation along with the very slight change in summer temperature.
Fig. 8. Trend rates of summer upper air (500 hPa) temperature, winter and summer precipitation from ERA5 data during 1980–2020. Black dot indicates the trends are significant at 99% confidence levels, and blue dot denotes the cell with glacier. The temperature in west Kunlun Mountains rises slightly in summer and gradually increases eastward, and there is no trend in winter precipitation, but an overall increasing trend in summer precipitation.
5. Conclusion
The Kunlun Mountains located in the northern Tibet Plateau, are connected to the Karakoram Mountains in the west and stretch to the far eastern part of the Tibet Plateau. Studying glacier change in this region can not only reveal the eastern limit of Kunlun-Pamir-Karakoram Anomaly, but also help us to understand the patterns in glacier change from west to east in the northern Tibet Plateau. Based on the glacier area and surface elevation variations, we found that the eastern limit of the positive mass-balance anomaly was located near 85°E. In the western part of this limit, the rate of glacier area change was −0.13 ± 0.10% a−1 before 2009 (CGI-2), and was −0.10 ± 0.57% a−1 after 2009 (CGI-2), the surface elevation increased slightly (~0.15 ± 0.35 m a−1) between 2000 and 2013. However, in the eastern part, the glacier area decreased significantly (−0.43 ± 0.13% a−1) and showed accelerating shrinkage. In the Aemye Ma-chhen Range the glacier surface elevation was lowered the most (−0.51 ± 0.18 m a−1). The spatial patterns of glacier change in the west and the east Kunlun Mountains appear to be linked to summer temperature, which increased strongly in the east Kunlun Mountains and very little in the west Kunlun Mountains. Increases in summer precipitation in the west Kunlun Mountains may also contribute to the mass-balance anomaly that extends to 85°E.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2022.30.
Data availability
The Chinese Glacier Inventory data are available at http://westdc.westgis.ac.cn/glacier. We downloaded the Landsat images and SRTM from the USGS (https://earthexplorer.usgs.gov/). TanDEM-X 90 m DEMs are available at https://geoservice.dlr.de/web/dataguide/tdm90/provided by German Aerospace Center (DLR). The python scripts for DEM coregistration provided by Shean and others (Reference Shean, Bhushan, Lilien and Meyer2019) are available on Github, and helped us a lot. The ERA-5 monthly averaged data on pressure levels from 1979 to present are available at https://www.copernicus.eu/en.
Acknowledgements
We thank the Scientific Editor Shea Joseph, the Associate Chief Editor Hester Jiskoot, and the anonymous reviewers for their constructive comments. Here, we thank all the institutions, groups and China Meteorological Administration for sharing data and methods. This work is supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant No. XDA20060201), the National Natural Science Foundation of China (42130516) and the Second Tibetan Plateau Scientific Expedition and Research Program (2019QZKK020102).
Author contributions
All authors contributed to the study conception and design. Ninglian Wang designed the study. Qian Liang, Anan Chen, Xuewen Yang, Ting Hua and Zhijie Li were involved in the discussions. Qian Liang compiled the data, performed the analysis and generated the figures. Qian Liang and Ninglian Wang wrote the manuscript with help from all other authors. Daqing Yang helped revise the manuscript.
Open access
