Temporal and Spatial Assessment of Glacier Elevation Change in the Kangri Karpo Region Using ASTER Data from 2000 to 2024

Abstract

Temperate glaciers in the Kangri Karpo region of the southeastern Qinghai–Tibet Plateau (QTP) have experienced significant ablation in recent decades, increasing the risk of glacier-related hazards and impacting regional water resources. However, the spatial and temporal pattern of mass loss in these glaciers remains inadequately quantified. In this study, we used ASTER L1A stereo images to construct a high-resolution elevation time series and provide a comprehensive spatial–temporal assessment of glacier elevation change from 2000 to 2024. The results indicate that almost all glaciers have experienced rapid ablation, with an average surface elevation decrease of −18.35 ± 5.13 m, corresponding to a rate of −0.76 ± 0.21 m yr⁻¹. Glaciers in the region were divided into the northern and southern basins, with average rates of −0.79 ± 0.17 m yr⁻¹ and −0.72 ± 0.13 m yr⁻¹, respectively. A notable difference in acceleration trends between the two basins was observed, with the elevation rate increasing from −0.78 ± 0.17m yr⁻¹ to −1.04 ± 0.17 m yr⁻¹ and from −0.52 ± 0.13 m yr⁻¹ to −0.92 ± 0.13 m yr⁻¹, respectively. The seasonal cycle was identified in glacier surface elevation change, with an accumulation period from November to March followed by a prolonged ablation period. The seasonal amplitude decreased with elevation, with higher elevations exhibiting longer accumulation periods and less ablation. Correlation analysis with meteorological data indicated that higher summer temperatures and increased summer rainfall intensify elevation loss, while increased spring snowfall may reduce ablation. Our analysis highlights distinct variations in glacier elevation changes across different locations, elevations, and climatic conditions in the Kangri Karpo region, providing valuable insights into glacier responses to environmental changes on the Tibetan Plateau.

1. Introduction

According to the IPCC’s sixth assessment report, the global air temperature was 1.1 °C higher in 2011–2020 than in 1850–1900, leading to irreversible hydrological changes driven by glacier melting [1]. The Qinghai–Tibet Plateau (QTP), often referred to as the “Water Tower of Asia”, holds the largest collection of glaciers outside the polar regions, supplying essential freshwater to some of the most densely populated areas on Earth [2,3,4]. Over the last half-century, temperatures on the QTP have increased at a rate of 0.3–0.4 °C (10 yr)⁻¹, approximately twice the global rate. Under global warming, it is shown that glaciers in the QTP experience continuous and accelerating mass loss but with distinct spatial–temporal variability across its subregions. In particular, the glaciers in the southeastern QTP have suffered from severe retreat [5,6,7,8,9].

Glaciers in the southeastern QTP are temperate, with annual rainfall ranging between 1000 and 3000 mm, primarily from the Indian monsoon [10,11,12,13]. It is reported that the summer air temperature around the equilibrium line ranges from 1 to 5 °C, while that of the glacier surface varies between 0 and 4 °C. This leads to the relatively large glacier ablation in summer [10,14]. In recent decades, the observed mass loss rate in the southeastern QTP has surpassed the global average, about three times the average of the QTP [6]. The Kangri Karpo region in the southeastern QTP is one of the most humid areas in the plateau and hosts a dense concentration of temperate glaciers [15]. Studies have shown that many glaciers in this region have been retreating since the early 20th century, raising concerns about long-term water availability and increasing the risk of glacier-related disasters [16,17]. Notably, two catastrophic glacier collapse events have been recorded in the vicinity of the Kangri Karpo [18]. Furthermore, a recent study highlights the increasing outburst flood risk of Guangxieco Proglacial Lake due to rapid glacier changes in the region [19]. However, due to the high altitude and complex terrain, it is challenging to conduct fieldwork, and there are limited reports from in situ observations in this region.

With the development of remote sensing technology, remote sensing data have been applied to estimate the mass balance of glaciers. Three main approaches were used: satellite altimetry, satellite gravimetry, and digital elevation model (DEM) differencing. Satellite altimetry usually has a higher accuracy and was originally designed to measure sea levels and ice sheets. But due to the sparse data and complex terrain, it is challenging to use repeat along-track and crossover analysis methods for alpine glaciers. It is often more suitable to select relatively flat glacier tongue areas or use the bin cluster and other statistical approaches to estimate elevation changes [9,20,21,22,23,24]. Satellite gravimetry measurements or Gravity Recovery and Climate Experiment (GRACE) measurements have better temporal continuity but typically very coarse spatial resolution. Moreover, it can be leaked by other signals, such as groundwater, lakes, and other water storage [25,26]. DEM differencing, which usually derives elevation changes from the difference of two DEMs at different time epochs, provides high-spatial-resolution results [27,28,29,30,31]. However, its temporal resolution is limited. To overcome this problem, Hugonnet et al. [32] developed a method to construct continuous elevation time series by filtering outliers and interpolating from all available DEMs. Some recent studies also aimed to achieve better spatial–temporal results by combining multiple approaches [26,33,34,35].

The mass balance of glaciers was estimated in the Kangri Karpo region using different approaches, covering varying time spans. Through in situ investigation, Yang et al. [36] assessed the mass balance of the six glaciers in the Kangri Karpo and its surroundings, and negative mass balance was observed in all six glaciers between 2005/2006 and 2007/2008. Based on topographic maps, Shuttle Radar Topography Mission (SRTM), and TerraSAR-X/TanDEM-X data, Wu et al. [27] found thinning rates of −0.27 ± 0.18 and −0.79 ± 0.11 m yr⁻¹ from 1980 to 2000 and from 2000 to 2014, indicating an acceleration in glacier mass loss. Using ZY-3 satellite images and SRTM data, Ren et al. [28,37] determined the thinning rate of −0.66 ± 0.24 m yr⁻¹ between 2000 and 2017. Integrating satellite altimetry, satellite gravimetry, and DEMs, Zhao et al. [35] found a thinning rate of −1.16 ± 0.29 m yr⁻¹ from 2000 to 2019 in the Eastern Bomi area.

In this study, we utilize updated ASTER L1A data and the method developed by Hugonnet et al. [32] to provide better spatial–temporal results of mass balance in the Kangri Karpo region from 2000 to 2024. We first constructed a continuous DEM time series and then derived the elevation change rates, seasonal patterns, and acceleration trends of glacier surface elevation across the region. We further showed the spatial–temporal variation of glacier elevation change in this region and analyzed its relationship with climate factors. This study provides us with more details on the characteristics of glacier elevation change in the Kangri Karpo region and also shows the response of glaciers to topography and climate change.

2. Study Area and Data

2.1. Study Area

The Kangri Karpo mountain range is located in the southeastern QTP, approximately 280 km long in a northwest–southeast direction. It is bounded to the north by the Purlung Tsangpo River and to the south by the Gongri-Gabo River. The Purlung Tsangpo River separates the range from the towering Nyainqentanglha Mountains, while to the south of the Gongri-Gabo River, the elevation decreases to 2000–3000 m This unique topography forms a critical corridor into the QTP for the warm and humid southwest Indian monsoon, which travels northward through the Yarlung Tsangpo and Chayu River valleys and delivers substantial precipitation to the Kangri Karpo region. The annual average precipitation in the region is about 2500–3000 mm. It is reported that summer temperatures near the glacier equilibrium line vary around 1 °C, and many glaciers are at the pressure melting point, resulting in significant surface melting and rapid glacier velocity [15].

The study area is shown in Figure 1, between 29.0° N and 29.5° N in latitude and between 96.3° E and 97.1° E in longitude. According to the Randolph Glacier Inventory 6.0 (RGI 6.0) [38], the region contains 127 glaciers, covering a total area of 719.01 km², and 17 of them are larger than 10 km². Of these, the largest is the Yalung Glacier, labeled as RGI60-15.11909 in RGI 6.0, with an area of 179.59 km² and situated on the northern slope. The Azha Glacier, labeled as RGI60-15.12613, on the southern slope, has a long history of scientific investigation [15,17,39], with the lowest terminus elevation around 2450 m. The Midui Glacier, or Gongzo Glacier, labeled as RGI60-15.11888 and situated on the north, is also known as one of the most renowned glacier tourism destinations in the QTP [40,41].

Significant topographical differences exist between the northern and southern slopes of the Kangri Karpo region. Based on the second Chinese glacier inventory (CGI-2), we divided the glaciers of this region into two parts: the northern and the southern, as shown in Figure 1. Glaciers in the northern basin labeled 5O282B have an area of 430.912 km² with an average elevation of 5131.76 m, while the southern counterpart 5O291B has a smaller area of 288.102 km² with an average elevation of 4880.50 m.

2.2. Data

2.2.1. ASTER L1A V003

The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) is a multispectral imaging sensor launched aboard NASA’s Terra satellite in December 1999 [42]. The primary science objective of the mission is to improve the understanding of surface–atmosphere interaction processes occurring on or near the earth’s surface and lower atmosphere at the local and regional scale [43]. It has three separate optical subsystems: the visible and near-infrared (VNIR), shortwave-infrared (SWIR), and thermal infrared (TIR) radiometer. The VNIR obtains the data with the 3N band in nadir and the 3B band in backward and provides the along-track stereo coverage at a spatial resolution of 15 m. This stereo coverage was further used to generate digital elevation models (DEMs) [44,45]. Up to now, the ASTER has provided continuous observations over two decades and is still ongoing. In this study, ASTER L1A V003 data, openly accessible on NASA Earthdata Search (https://search.earthdata.nasa.gov/search, accessed on 19 September 2024), were used to generate DEMs and build elevation time series.

2.2.2. TanDEM-X DEM

The TerraSAR-X add-on for Digital Elevation Measurements (TanDEM-X) mission produces high-precision and high-resolution DEMs globally [46]. The TanDEM-X DEM product has a typical absolute vertical accuracy of 3.5 m and few void pixels, even in complex terrain [46,47]. In this study, the 90 m TanDEM-X DEM was used as a reference DEM to co-register, remove along-track jitter, and filter outliers from the DEMs generated from ASTER images.

2.2.3. Sentinel-2

The Sentinel-2 mission comprises two identical satellites, Sentinel-2A and Sentinel-2B, equipped with multispectral imagers of 13 spectral bands, ranging from visible to shortwave infrared wavelengths. Its spatial resolution varies from 10 to 60 m, depending on the spectral band [48]. In addition, the 5-day repeat cycle makes it ideal for land monitoring applications. In this study, we used it to derive the Normalized Difference Snow Index (NDSI) for mapping snow cover areas.

2.2.4. Auxiliary Data

The Ice, Cloud and Land Elevation Satellite-2 (ICESat-2) mission was launched in 2018. It employed laser altimetry to measure Earth’s surface elevation with a narrow footprint of 14.5 m and a vertical precision of 3.0 cm [49,50]. In this study, Land Ice Height data labeled as ATLAS/ICESat-2 L3A ATL08 were used to validate our elevation time series. RGI 6.0 was also used to identify glacier and non-glacier areas [38].

3. Method

To derive the monthly DEM time series and further analyze each glacier in the Kangri Karpo region, we developed a workflow as shown in Figure 2. The workflow mainly includes the generation and co-registration of ASTER DEMs, followed by the establishment of the elevation time series with ASTER DEMs.

3.1. Generation and Co-Registration of ASTER DEMs

We first retrieved all available 1119 ASTER L1A V003 scenes with cloud coverage below 99% over the Kangri Karpo region (28.6° N–29.9° N, 95.9° E–97.3° E) from 2000 to March 2024. Before generating the ASTER DEMs, we deleted the data with no stereo bands, the Band 3N and the Band 3B, and 757 stereo pairs were kept

In this study, we employed an algorithm called MicMac ASTER (MMASTER) [51] to generate ASTER DEMs from 757 stereo pairs. Compared to NASA’s standard AST14DMO product, the DEMs generated by this method have fewer mismatched areas, less overall noise and higher accuracy. To increase the range of off-glacier terrain available for the co-registration, consecutive images from the same track were mosaicked into longer strip DEMs with a maximum of three images per strip, and 334 striped ASTER DEMs were generated.

To correct residual along-track jitter and align all DEMs to the same elevation system, it is necessary to co-register all striped DEMs to the TanDEM-X DEM. To exclude the unstable areas such as glaciers, snow cover and other areas for co-registration, Sentinel-2 images from January 2023 were selected with less than 10% cloud cover, and the NDSI was calculated to distinguish snow-covered areas. The NDSI is the Normalized Difference Index (NDI) between the green and shortwave (SWIR) infrared bands, calculated as follows [52]:

$$NDSI={\displaystyle \frac{\left(Green-SWIR\right)}{\left(Green+SWIR\right)}}$$

(1)

Sentinel-2 images have 13 spectral bands, and the green and SWIR infrared bands are those of B3 and B11 [48]. The NDSI threshold value was selected as 0.4 to extract the snow-covered area [53], and RGI 6.0 was combined to create the mask for unstable terrain.

After the co-registration, the root-mean-square error (RMSE) of the elevation difference between co-registered DEMs and TanDEM-X DEM on stable terrain was calculated, and it was kept when it was less than 20 m. As a result, 257 co-registered ASTER DEMs were selected.

3.2. Establishment of the Elevation Time Series

The method of Hugonnet et al. was used to derive a continuous DEM time series [32]. Firstly, we stacked all 257 co-registered DEMs in temporal order, yielding a three-dimensional array h (t, x, y) of elevations for a pixel (x, y) as the function of time t and spatial x and y. As array h (t, x, y) may contain outliers due to clouds, photogrammetric blunders, and other reasons, we first used TanDEM-X DEM as a reference to remove extreme outliers with a spatial filter and a temporal filter defined by Hugonnet et al. [32]. To further filter outliers of each pixel, we used weighted least squares (WLS) to fit a linear elevation trend and remove outliers outside the 99% confidence interval. Then, we used WLS fit again and applied an adaptive local percentile filter to remove the remaining outliers.

To further filter and interpolate array h (t, x, y), a non-parametric and empirically based Gaussian process (GP) method was applied to elevation time series, and the kernel function was expressed as [32]

$$\begin{array}{c}{\mathsf{\sigma }}_{\mathrm{h}}{\left(\mathrm{x},\mathrm{y},\mathsf{\Delta }t\right)}^2=\mathrm{P}\mathrm{L}\left(\mathrm{x},\mathrm{y},\mathsf{\Delta }t\right)+\mathrm{E}\mathrm{S}\mathrm{S}\left({\mathsf{\varphi }}_p,{{\mathsf{\sigma }}_p}^2,\mathsf{\Delta }t\right)+\mathrm{R}\mathrm{B}\mathrm{F}\left(\mathsf{\Delta }{t}_l,{{\mathsf{\sigma }}_l}^2,\mathsf{\Delta }t\right)+\\ \mathrm{R}\mathrm{Q}\left(\mathsf{\Delta }{t}_{nl},{\mathsf{\sigma }}_{nl}^2,{\mathsf{\alpha }}_{nl},\mathsf{\Delta }t\right)·\mathrm{P}\mathrm{L}\left(\mathrm{x},\mathrm{y},\mathsf{\Delta }t\right)+{\mathsf{\sigma }}_{\mathrm{h}}{\left(\mathrm{t},\mathrm{x},\mathrm{y}\right)}^2\end{array}$$

(2)

where ${\mathsf{\sigma }}_{\mathrm{h}}$ is the elevation error and $\mathsf{\Delta }t$ is the time lag between observations. The pairwise linear (PL) kernel shows the long-term trend, the exponential sine-squared (ESS) kernel shows seasonal variation, the local radial basis function (RBF) kernel shows the temporal correlation between observations, and the rational quadratic (RQ) kernel shows a nonlinear elevation change trend. The process was shown in Hugonnet et al. [32], where ESS, RBF, and RQ parameters were set the same for all pixels, but PL varied with each pixel. RGI 6.0 was used to identify glacier and non-glacier areas and apply the complete set of kernels to glacier areas, while only the PL and ESS kernels were applied to non-glacier areas.

For each pixel, GP regressions were iteratively applied five times, and the threshold value was set as 20$\mathsf{\sigma }$, 12$\mathsf{\sigma }$, 9$\mathsf{\sigma }$, 6$\mathsf{\sigma }$, and 4$\mathsf{\sigma }$, respectively. Finally, the GP regression was applied to interpolate the remaining time series for each pixel, and we derived the interpolated monthly elevation time series from January 2000 to January 2024.

4. Results

4.1. Accuracy Analysis with ICESat-2

The vertical accuracy of ASTER DEM standard data products ranges from 10 to 25 m [54]. Here, ICESat-2 ATL08 data shown in Figure 1 were employed to assess the accuracy of the elevation time series, covering the glacier and non-glacier regions during the summer (JJA) between 2019 and 2023. Each ICESat-2 ATL08 dataset we selected contains more than 1000 valid observation points within the study area. The elevation time series were linearly interpolated according to the time epoch of the ICESat-2 data. As shown in Figure 3a, the ICESat-2 data and corresponding elevation time series were consistent. We further calculated the root-mean-square error (RMSE) between them, and the overall RMSE was 17.45 m, matching the ASTER DEM performance. However, our result is larger than 11.47 m observed between ICESat-2 ATL08 and ASTER GDEM V3 in the Qinghai–Tibet Plateau [55]. Hence, the comparisons including extra glacier area may induce extra error. We further the point–point comparison for all data, as shown in Figure 3b, and the coefficient of determination (R²) between them was 1.00 with p < 0.0001. Based on the analysis, the accuracy of our elevation time series is sufficient for further glacier mass balance calculation.

4.2. Temporal Change of Glacier Elevation from 2000 to 2024

Based on the monthly elevation time series derived from ASTER, we extracted the surface elevation change of each glacier in the Kangri Karpo region, as shown in Appendix A. Almost all glaciers in the Kangri Karpo region experienced a notable elevation decrease between 2000 and 2024. Using the area weight, we derived the elevation time series from 2000 to 2024, as shown in Figure 4a. Obvious seasonality and elevation decreases were observed. Moreover, the average glacier surface elevation decrease during 24 years was −18.35 ± 5.10 m, corresponding to an average rate of −0.76 ± 0.21 m yr⁻¹.

Fitting a quadratic trend of the elevation change, we found the acceleration was statistically significant. To further show this acceleration, we calculated the elevation change rate, which was derived from 12-month samples of each year, and the time series is shown in Figure 4b. The time series shows the rate increased from −0.68 m yr⁻¹ to −0.99 m yr⁻¹, showing the acceleration. Moreover, there exist three distinct phases in the Kangri Karpo region from the figure. Before 2008, the elevation change rate was relatively stable, even with a slight deceleration trend. Between 2008 and 2018, the elevation change rate became more negative, showing the acceleration. Furthermore, it fluctuated with relatively minimal values in 2008, 2009, and 2014 and a maximum in 2016. After 2018, the elevation change rate became more negative.

We further investigate the elevation change of the northern 5O282B basin and the southern 5O291B basin. Similar to the whole region, obvious elevation decrease and seasonal variation were also observed for the two parts, as shown in Figure 4a. Compared with that of Kangri Karpo, 5O282B had a larger elevation decrease rate of −0.79 ± 0.16 m yr⁻¹, while 5O291B showed a smaller −0.72 ± 0.13 m yr⁻¹. The acceleration was also statistically significant for the two parts. We also derived the elevation change rate shown in Figure 4b. Despite the larger elevation decrease rate, a slower acceleration was observed in 5O282B, which increased from −0.78 m yr⁻¹ in 2000 to −1.04 m yr⁻¹ in 2023. Similar to Kangri Karpo, there also exist three distinct phases in 5O282B from the figure. From 2008 to 2018, it fluctuated with relative minimum values in 2008, 2009, and 2014 and with a maximum in 2016. However, there is no obvious fluctuation in 24 years for 5O291B. 5O291B has experienced continuous acceleration since 2000, with a larger acceleration from −0.52 m yr⁻¹ in 2000 to −0.92 m yr⁻¹ in 2023. This demonstrates that the variation of Kangri Karpo mainly depends on 5O282B, which has a larger area.

We also investigated the seasonal variation of the Kangri Karpo region and its two parts, 5O282B and 5O291B.To capture the absolute elevation changes of this seasonal variation, the mean elevation changes for each month, compared to the beginning of the year, are shown in Figure 5a. Distinct seasonal variation was observed in the Kangri Karpo region and 5O282B, but not very obvious in 5O291B. The accumulation starts in about November and reaches the peak in about next March for the former two. Then, the elevation decreases until next November with the minimum elevation. For 5O291B, a shorter accumulation period was observed from December to next January, and a longer ablation from next February to next November. This seasonal pattern is similar to those reported in previous studies on southeastern Tibet [21,22]. Overall, the Kangri Karpo region has a gain of 0.17 ± 0.21 m during the accumulation period, while a decrease of −0.93 ± 0.21 m during the ablation. Notably, our analysis revealed some changes in seasonal variation patterns, characterized by a delayed onset of the accumulation period from October to November and an earlier termination from March to February. Additionally, the magnitude of elevation loss during the ablation period showed an increasing trend, while the elevation gains during the accumulation period showed a slight decreasing trend. Again, the seasonal variation of Kangri Karpo mainly depended on 5O282B.

There exist different seasonal variations between 5O282B and 5O291B, implying the elevation change may depend on the altitude. To further investigate this relationship, we further analyzed the seasonal elevation change variation with altitude. We divided the region by 500 m spaced elevation bins and extracted the seasonal change of each bin, as shown in Figure 5b. From the figure, it can be seen that the seasonal amplitude decreases with the increasing altitude. Moreover, the accumulation period and the ablation period vary with the altitude. The seasonal variation below 3500 m is similar to that of 5O291B, only with the 2-month accumulation period from December to next January. Furthermore, the seasonal variation between 3500 and 6500 m is similar to that of the Kangri Karpo region and 5O282B, with a longer accumulation period from October or November to next March or April. For the regions with an altitude between 6500 and 7000 m, there is no obvious seasonal variation.

4.3. Spatial Variation of Glacier Surface Elevation

We further calculated the elevation change of each pixel between January 2024 and January 2000, as shown in Figure 6a, with significant spatial variation across the Kangri Karpo region. The overall elevation decrease was −18.35 ± 5.10 m, and almost all glaciers experienced elevation decreases. Large elevation decreases were observed at some terminal zones, especially in Gongzo, Yalong, and Azha Glacier, exceeding −96 m in 24 years with a corresponding rate of −4 m yr⁻¹. Meanwhile, elevation increases were observed in some areas, such as in areas above 5800 m and some glacier termini. To systematically analyze these elevation changes, we further calculated the mean elevation change for each 200 m elevation bin between January 2000 and January 2024, as shown in Figure 6b. The elevation bins started from 2400 m and ended at 6800 m. The elevation change decreased quickly from −24.18 m in the 2500 m bin to −120.14 m in the 2700 m bin. Then, it increased quickly to −46.05 m in the 3300 m bin, decreased gradually to −59.43 m in the 4300 m bin, increased quickly again to near 0 m in the 5900 m bin, and then remained stable.

Similar elevation decreases over 24 years were also found in 5O282B and 5O291B, with −19.00 ± 3.98 m and −17.35 ± 3.23 m, respectively. The mean elevation changes of the 200 m bins for 5O282B and 5O291B are also shown in Figure 6b. Their minimum elevation bin is different, where the former is 2500 m and the latter 3900 m. Hence, the elevation change in 5O291B is the same as that of the Kangri Karpo region below 3900 m bin. From 3900 to 5900 m, the elevation change increased quickly to near 0 m for both 5O282B and 5O291B. Then, they remained stable. Moreover, the glaciers in 5O291B experienced less ablation than those in 5O282B between 3900 and 5900 m, and this may be due to the higher warming trend in the northern area [27].

5. Discussion

5.1. Variations in Elevation Change Across Different Glaciers

To elucidate the spatiotemporal characteristics of glacier mass balance in the Kangri Karpo region, we selected several typical glaciers, including Yalong Glacier and Gongzo Glacier on the northern slope and Azha Glacier and RGI60-15.12587 Glacier on the southern slope. As shown in Figure 6a, there exists spatial variation across the Yalong Glacier. We specifically selected the upper region of Yalong Glacier with obvious mass gain labeled as Area A to show the spatiotemporal variation across the glacier.

We further showed the elevation time series of selected glaciers in Figure 7a, showing seasonal changes. For the northern 5O282B, Yalong Glacier and Gongzo Glacier have an accumulation period from November to April and an ablation period from May to October, which is consistent with the results shown in Figure 5a. For the southern 5O291B, Azha Glacier and RGI60-15.12587 Glacier have a shorter accumulation period and a longer ablation period, which is also shown in Figure 5a. Moreover, the amplitude is different for the northern and southern parts. For example, it is 1.09 and 0.85 m for the northern Yalong Glacier and Gongzo Glacier, and these values are much larger than those of 0.07 and 0.31 m for the southern RGI60-15.12587 and Azha Glacier. Area A showed a large amplitude of 1.56 m. This may be due to different topography and climate conditions.

As shown in Figure 7a, all selected glaciers experienced mass loss. During the last 24 years, the northern Yalong Glacier and Gongzo Glacier exhibited close elevation decreases, which were −17.84 m and −14.59 m, corresponding to average decrease rates of −0.74 m yr⁻¹ and −0.61 m yr⁻¹, respectively. However, there exists an obvious difference between the southern Azha Glacier and RGI60-15.12587 Glacier. During the last 24 years, the Azha Glacier showed the largest elevation decrease in the region, −58.31 m, corresponding to an average rate of −2.43 m yr⁻¹. However, the latter only exhibited a very small elevation decrease of −4.52 m. Different from the four glaciers, Area A showed an elevation gain of 41.17m.

To better show the elevation change of RGI60-15.12587 Glacier, we selected its terminal Area B with a slope less than 5°, as shown in Figure 7b. The elevation time series of Area B exhibits a small fluctuation, as shown in Figure 7c. Moreover, there was an overall elevation increase of 1.99 m during the 24 years, and this is consistent with Wu et al.’s study using only the data from 1999 and 2014 [27]. Furthermore, the elevation trend change was identified before and after 2008. Before 2008, the elevation continued to decrease for −2.15 m, at a rate of −0.27 m yr⁻¹, which is less than the other glaciers. Then, the elevation increased by 4.14 m at a rate of 0.26 m yr⁻¹. The underlying causes of these changes remain to be investigated further.

5.2. Influences of Climate Factors in Kangri Karpo

To analyze the impact of climate factors on the mass balance of the Kangri Karpo glaciers, we extracted monthly temperature, total precipitation (TP), and snowfall (SF) data from the ERA5-Land monthly meteorological dataset, covering the area between 29.0° N and 29.5° N in latitude and between 96.3° E and 97.1° E in longitude. From 2000 to 2024, the mean temperature in the Kangri Karpo region was −3.98 °C, with annual mean temperatures ranging between −4.89 °C and −2.96 °C. Warming was observed with a rate of 0.39 °C (10 yr)⁻¹ derived with simple linear fitting, which is about twice the global warming rate of 0.2 °C (10 yr)⁻¹ during recent decades [1]. We also estimate the annual total TP and SF between 2000 and 2024. The former was 1809 mm, ranging between 877 and 2632 mm, while the latter was 1142 mm, varying between 475 and 1992 mm. Moreover, a similar decreasing trend was observed for both TP and SF, with a rate of −127.80 and −142.30 mm (10 yr)⁻¹, respectively. The decline in SF was more significant, suggesting that the reduction in TP is primarily driven by the decrease in SF. The pronounced warming trend and decreasing precipitation/snowfall likely led to accelerated glacier ablation and reduced accumulation in the Kangri Karpo region.

To further analyze the influence of climate factors on elevation change, we remove the linear trend, apply z-score normalization, and then calculate the Pearson correlation coefficients between elevation and three parameters. Their correlation coefficients between elevation and temperature, TP, and SF were estimated as −0.12, −0.11, and 0.26, respectively. These results suggest no significant correlation at the annual scale. We further the correlation analyses at a seasonal scale, and here we divide a year into four seasons labeled as winter (DJF), spring (MAM), summer (JJA), and autumn (SON). Among the results, significant correlation coefficients were observed between elevation and temperature in JJA, TP in JJA, and SF in MAM, shown in Figure 8. Specifically, temperature in JJA showed a significant negative correlation coefficient of −0.53 (p < 0.05) with elevation in JJA, indicating that higher summer temperature is associated with more pronounced elevation decreases. Moreover, TP in JJA also showed a negative correlation of −0.35 (p < 0.1), and this also suggests that increasing summer rainfall may accelerate glacier ablation and is consistent with previous results [56]. Apart from these results, SF in MAM exhibited a positive correlation coefficient of 0.32 (p < 0.13), implying that higher snowfall during spring may increase the elevation.

Similar to Figure 1, we divided the area into northeastern and southwestern parts. The temperature increase rate in JJA was estimated as 0.25 °C (10 yr)⁻¹ in the south, while a higher rate of 0.30 °C (10 yr)⁻¹ was observed in the north. This temperature pattern corresponds with the smaller decrease in elevation observed in the south. In terms of TP during JJA, a slightly decreasing trend was found in both regions: −20.07 mm (10 yr)⁻¹ in the north and −9.90 mm (10 yr)⁻¹ in the south. With more summer rainfall in the south, this could be related to the larger acceleration of elevation decrease in the south. Furthermore, SF in MAM showed a decreasing trend of −67.49 mm (10 yr)⁻¹ in the north, and a larger −110.38 mm (10 yr)⁻¹ in the south, indicating a greater reduction in accumulation in the southern area. As shown in Figure 6a, there is an obvious spatial variation in elevation change across the Kangri Karpo region; further research is needed to analyze this difference at the glacier or pixel scale, incorporating influencing factors such as topography, ice lakes, and debris cover. Additionally, exploring the feedback mechanisms between glacier dynamics and local climate systems could yield deeper insights into the cascading effects of climate change on glaciers. Future studies should also focus on the development of predictive models that integrate observed trends and underlying mechanisms with projected climate scenarios, which would enable more accurate assessments of glacier vulnerability and improve regional water resource management.

5.3. Comparison with Previous Mass Balance Estimates in Kangri Karpo

Studies have shown that glaciers in the Kangri Karpo region and the surrounding areas are experiencing significant mass loss primarily derived from remote sensing data. However, the results vary due to differences in methods, data sources, and study periods Here, we compare our results with previous studies shown in

Table 1. Comparison of our results with previous studies on Kangri Karpo and surrounding areas.

Study	Region	Period	Elevation Change Rate (m yr−1)	Elevation Change Rate of This Study Mean (m yr−1)
Wu et al. [27]	Kangri Karpo	2000–2014	−0.79 ± 0.11	−0.71 ± 0.33
Ren et al. [28]	Kangri Karpo	2000–2017	−0.66 ± 0.24	−0.71 ± 0.28
Zhao et al. [35]	Eastern Bomi	2000–2019	−1.16 ± 0.29	−0.74 ± 0.25
		2011–2020	−1.14 ± 0.28	−0.79 ± 0.50
Hugonnet et al. [32]	Kangri Karpo	2000–2020	−0.88 ± 0.35	−0.73 ± 0.25
Kaab et al. [57]	Eastern Nyainqentanglha range	2003–2008	−1.34 ± 0.29	−0.69 ± 0.84
Neckel et al. [33]	Eastern Nyainqentanglha range and Hengduan Mountains	2003–2009	−0.81 ± 0.32	−0.70 ± 0.84
Gardner et al. [8]	Hengduan Shan	2003–2009	−0.40 ± 0.41	−0.70 ± 0.84
Brun et al. [6]	Nyainqentanglha	2000–2016	−0.73 ± 0.27	−0.71 ± 0.32
Shean et al. [7]	Nyainqentanglha	2000–2018	−0.59 ± 0.18	−0.72 ± 0.28

[6,7,8,27,28,32,33,35,57]. At the regional scale, all studies reported the glaciers in the Kangri Karpo region and southeastern QTP were experiencing mass loss, with the elevation change rate ranging from −0.57 m yr⁻¹ to −1.16 m yr⁻¹. However, most studies calculated an average result derived only over two epochs. In this study, we derive the elevation time series from 2000 to 2024 shown in Figure 4, and hence, we can derive the elevation change rate of any time span. So, we derived the corresponding results from the time series as shown in Table 1. Consistent results were observed, especially with better agreements for studies utilizing similar data and analysis methods [6,27,28,35]. We also conducted compared the elevation change rates of each year with the results reported by Hugonnet et al. [32]. Our results have the same fluctuation patterns and trends, with only minor numerical differences. These slight variations can be attributed to differences in data sources, temporal coverage, and subtle methodological distinctions.

In situ measurements are the most reliable observation. We further the comparison at the glacier scale using multi-year in situ measurement data for Parlung Glacier No. 10, No. 12 of the Kangri Karpo region, and No. 94 next to the region, which were obtained from the World Glacier Monitoring Service (WGMS), as shown in Figure 9 [58,59]. Among these three glaciers, Parlung Glacier No. 94 has a longer dataset from 2005 to 2021, as shown in Figure 9a. Field measurements may yield large variations in results due to the limited density of measurement points. A moving 3-year average was applied to the time series. The mean and standard deviation of their difference were −0.03 and 0.40 m, showing consistency. However, our results did not capture the accelerated trend observed in the latter period of the in situ measurements. This may be attributed to the lack of observational data, since this glacier is located outside the Kangri Karpo region defined in this study, at the edge of our data coverage. For Parlung Glacier No. 10 and No. 12, we directly calculated the mean and standard deviation of their differences, as 0.02 m and 0.48 m for the former, −0.52 m and 0.46 m for the latter. In these shorter periods, our results exhibited slightly larger deviations, which may be attributed to the limitations of our method in accuracy or capturing short-term abrupt events or fluctuations in the in situ measurements.

6. Conclusions

This study provides a comprehensive analysis of glacier mass balance in the Kangri Karpo region from 2000 to 2024, revealing significant and accelerating glacier mass loss. We utilized ASTER images to construct a high-resolution elevation time series and validated the results with ICESat-2 data, with an overall RMSE of 17.45 m. The results showed that the mean glacier surface elevation has decreased at an average rate of −0.76 ± 0.15 m yr⁻¹ over the past 24 years. The 5O282B basin has a higher surface elevation but experienced a larger average ablation rate of −0.79 ± 0.17 m yr⁻¹, showing an acceleration from −0.78 ± 0.17 m yr⁻¹ in 2000 to −1.04 ± 0.17 m yr⁻¹ in 2023. In contrast, the 5O291B basin had a lower average rate of −0.72 ± 0.13 m yr⁻¹, with a larger acceleration from −0.52 ± 0.13 m yr⁻¹ in 2000 to −0.92 ± 0.13 m yr⁻¹ in 2023. Distinct seasonal variations were observed. On average, the region exhibited a gain of 0.17 ± 0.21 m during the accumulation period of November to March and a loss of −0.93 ± 0.21 m during the ablation period of April to October. A clear altitude-dependent pattern was found: as elevation increased, the seasonal amplitude decreased, the accumulation period extended, and ablation reduced.

The overall elevation decrease for the Kangri Karpo region was −18.35 ± 5.10 m, with the 5O282B and 5O291B basins showing decreases of −19.00 ± 3.96 m and −17.35 ± 3.22 m, respectively. Elevation bin analysis revealed a decrease in ablation with increasing elevation, with the ELA around 5900 m. Notably, glaciers in 5O291B experienced less ablation at the same elevation compared to those in 5O282B. We also selected four glaciers: Yalong and Gongzo glaciers from 5O282B and Azha and RGI60-15.12587 glaciers from 5O291B. Yalong and Gongzo glaciers exhibited elevation decreases of −17.84 m and −14.59 m, respectively, and mass gain was observed in the upper area of Yalong Glacier. In contrast, Azha Glacier showed the largest decrease of −58.31 m, while RGI60-15.12587 Glacier had a slight decrease of −4.52 m but an increase of 1.99 m in its terminal area.

Meteorological data reveal significant trends of increasing temperature and decreasing precipitation and snowfall. Correlation analysis suggests that higher temperature in JJA and increased TP of JJA could accelerate elevation loss, while more SF during MAM could reduce ablation. These climatic factors likely contribute to the spatial variations observed between the northern and southern parts of the region. However, further research is needed to fully understand the underlying causes of these spatial discrepancies in elevation change across the Kangri Karpo region.