The eastern limit of ‘Kunlun-Pamir-Karakoram Anomaly’ reflected by changes in glacier area and surface elevation

Abstract

Under global warming, many glaciers worldwide are receding. However, recent studies have suggested the extension of the Karakoram Anomaly, a region of anomalous glacier mass gain, into the western Kunlun and eastern Pamir mountains. However, the eastern limit of this anomaly in the Kunlun Mountains is unclear. This study, using changes in glacier area and surface elevation, estimates the eastern limit of the Kunlun-Pamir-Karakoram anomaly at ~85°E. Over the past 50 years, glaciers west of 85°E in the Kunlun Mountains decreased in area from 8401 to 7945 km² at a rate of −0.12 ± 0.07% a⁻¹, showed a reduction in the rate of retreat through time and have recently gained mass, with surface elevation changes of 0.15 ± 0.35 m a⁻¹ over the period of 2000–2013. Glaciers east of 85°E have experienced greater rates of area change (−61 ± 12 km² and −0.43 ± 0.13% a⁻¹) over the past 50 years, accelerated area loss in recent years and elevation change rate of −0.51 ± 0.18 m a⁻¹ between 2000 and 2013. These patterns of elevation and area change are consistent with regional increases in summer temperature in the eastern Kunlun Mountains and slight cooling in the western Kunlun Mountains.

1. Introduction

Since the beginning of this century, most glaciers have retreated worldwide at an increasing rate (Dehecq and others, 2019). 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, 2005). Gardelle and others (2013) 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⁻¹) and western Pamir (0.14 ± 0.13 m w.e. a⁻¹). Recent studies have shown that anomalous glacier behavior probably occurs in the nearby west Kunlun Mountains and Pamir (Kääb and others, 2015; Brun and others, 2017; Berthier and Brun, 2019; Farinotti and others, 2020), the extent of mass gain could even reach to the central of the Tibet Plateau from 2003 to 2009 (Neckel and others, 2014). 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, 2019). This accelerating ice velocity contrasted with the other parts of HMA, where the glacier flow decreased. Farinotti and others (2020) suggested more mass gain in the west Kunlun Mountains compared to the results in the Karakoram. Negi and others (2021) also reported that there was no glacier anomaly in the eastern Karakoram (−0.10 ± 0.07 m w.e. a⁻¹) 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, 2015). 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 (2007) revealed that 278 glaciers (with a total area of 2711.57 km²) 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 (2015) observed the glacier area (283 with a total area of 2822.04 km²) 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, 2014), and these glaciers gained mass during 2003–2009 (Bao and others, 2015) in the west Kunlun Mountains. Recently, Wang and others (2018), using multiple DEM sources, estimated average glacier mass change of −0.06 ± 0.13 m w.e. a⁻¹ in the west Kunlun Mountains between 1970s and 1999. However, Zhou and others (2018) 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, 2015; Chudley and Willis, 2018), 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, 2002; Tian and others, 2007; Yao and others, 2013; Ma and others, 2017; Thompson and others, 2018). 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, 1999). 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, 2012), there are 8452 glaciers with a total area of 11 099 ± 400 km² 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, 2015).

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, 2000). 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, 2020). 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, 2019).

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, 2013). 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, 2009). 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, 2017). Because of the limitations of the image resolution, glaciers with areas of <0.01 km² cannot be reliably discriminated from Landsat images (Paul and others, 2009; Leigh and others, 2019) 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, 2011). 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, 2007). 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, 2010). Gardelle and others (2012) 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, 1994; Zemp and others, 2010). To transform volume change into mass change, we assume a density of 850 kg m⁻³ and an additional uncertainty of 60 kg m⁻³ on account of the lack of measured data (Zemp and others, 2010). Mass balance was estimated by a density conversion formula (1), which could be expressed as follows:

(1)$$B = \displaystyle{\rho \over {S_{\rm g}}} \times \mathop \sum \limits_{i = 1}^n \Delta h_iS_{\rm p}, \;$$

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, 2020; Huai and others, 2021). 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, 2013). 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, 2013). The error of the glacier area can be calculated as formula (2):

(2)$$E_{\rm A} = N \times \displaystyle{1 \over 2}S, \;$$

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 (2013). The errors (E _C) are regarded as follows according to the error propagation (3):

(3)$$E_{\rm C} = \sqrt {E_1^2 + E_2^2 } , \;$$

where E ₁ and E ₂ 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, 2020), the number of values (after removing spatial autocorrelation) in the non-glacier area (Nuth and Kääb, 2011; Paul and others, 2015) and the averaged absolute difference (AAD) between the median elevation change on the glacier and off-glacier terrain (Berthier and Brun, 2019) in the formula (4):

(4)$$U_{\Delta {\rm h}} = \displaystyle{{{\rm ST}{\rm D}_{\Delta {\rm h}}} \over {\sqrt {N_{{\rm eff}}} }} + {\rm AAD, \;}$$

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:

(5)$$N_{{\rm eff}} = \displaystyle{{P \times N} \over {2d}}, \;$$

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, 2011).

For the uncertainty arising from glacier boundary, a previous study (Paul and others, 2013) suggested a 3% error for alpine regions and a ratio between the glacier perimeter and area (R _P/A) of 5.03 km⁻¹. Braun and others (2019) observed R _P/A values of individual regions range from 0.67 to 10.81 km⁻¹, 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 (2019):

(6)$$U_{\rm A} = \displaystyle{{R_{{\rm P/A}}} \over {R_{{\rm P/A\;Paul\;and\;others}}}} \times 0.3, \;$$

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, 2017). Previous studies (Lin and others, 2017; Li and others, 2021) 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 (2015) 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 (2021): 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 (2015) 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, 2014). 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, 2015; Li and others, 2018). Liu and others (2019) 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 (2019). We considered the error from glacier elevation change, glacier boundaries, uncertainty from the conversion density and radar signal penetration:

(7)$$ \! \! \! \! \! \! \! \! \! \!\!\!\!\!\!\!\!\! \! U_{\rm M} = \sqrt { \!\!\!{\left({\displaystyle{{\Delta M} \over {\Delta t}}} \right)}^{\!\!\!2} \left( \! \! \!{{\left({\displaystyle{{U_{\Delta {\rm h}}} \over {\Delta h}}} \right)}^{\!\!\!2} + {\left({\displaystyle{{U_{\rm A}} \over A}} \right)}^{\!\!\!2} + {\left({\displaystyle{{U_\rho } \over \rho }} \right)}^ {\!\!\!2}} \right) + {\left({\displaystyle{{U_{{\rm pen}}} \over {\Delta t}} \times \rho } \right)}^{\!\!\!2},} $$

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 km² in the 1970s, were considered. From 1970s to 2016, a total glacier area change of −593 ± 424 km² (−5.41 ± 4%) was observed over the Kunlun Mountains. The maximum absolute glacier area change occurred in the west Kunlun Mountains (−455 ± 313 km²), and the minimum occurred in the east (−60 ± 12 km²). 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⁻¹ between 1970s and 2016), large increases in the rate of area loss in the central Kunlun Mountains (−0.01% a⁻¹ prior CGI-2, and −0.24% a⁻¹ after CGI-2), and consistent but relatively lower rates of area loss in the western Kunlun Mountains (−0.12% a⁻¹ 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⁻¹), while in the east it was larger (>1.20% a⁻¹), and the shrinkage of small glaciers in the west (most <0.5% a⁻¹) was much smaller than it is in the east (>1.50% a⁻¹). 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⁻¹ from 1973 to 2014 (Qureshi and others, 2017). 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⁻¹ between 2000 and 2013. Assuming a surface density of 850 ± 60 kg m⁻³, the corresponding average mass change was 0.07 ± 0.12 m⋅w.e. a⁻¹ during the period, which is consistent with the results (0.06 ± 0.17 m w.e. a⁻¹) from Shean and others (2020). 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⁻¹ (0.13 ± 0.30 m w.e. a⁻¹) 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⁻¹) estimated by Brun and others (2017). The most negative change was found in the eastern Kunlun Mountains in Aemye Ma-chhen Range of −0.51 ± 0.18 m a⁻¹, which is similar to the value of −0.68 ± 0.29 m a⁻¹ calculated by Shean and others (2020) (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, 2014) (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⁻¹) in the west but negative (~−0.51 ± 0.18 m a⁻¹) 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⁻¹, Table 1), consistent with previous values (Supplementary Table S9). A mass balance of 0.12 ± 0.17 m w.e. a⁻¹ was found in the central Kunlun Mountains (Sub-region 2), with clear mass loss (−0.09 ± 0.05 m w.e. a⁻¹) 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⁻¹), 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⁻¹).

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 (2018) 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⁻¹ 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 km²) used to calculate the surface elevation change accounts for nearly 70% of the total glacier area (10 965.07 km²) 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, 2005) 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, 2010; Bao and others, 2015; Kang and others, 2020; Bhattacharya and others, 2021). 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, 2017). Furthermore, mass balance of glaciers in the east part is more sensitive to temperature than those in the west part (Sakai and Fujita, 2017). 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, 2021). 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⁻¹ before 2009 (CGI-2), and was −0.10 ± 0.57% a⁻¹ after 2009 (CGI-2), the surface elevation increased slightly (~0.15 ± 0.35 m a⁻¹) between 2000 and 2013. However, in the eastern part, the glacier area decreased significantly (−0.43 ± 0.13% a⁻¹) and showed accelerating shrinkage. In the Aemye Ma-chhen Range the glacier surface elevation was lowered the most (−0.51 ± 0.18 m a⁻¹). 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.