Decadal Changes in Glacier Area, Surface Elevation and Mass Balance for 2000–2020 in the Eastern Tanggula Mountains Using Optical Images and TanDEM-X Radar Data

Abstract

The response of lake-terminating glaciers to climate change is complex, and their rapid changes are often closely linked to glacial-lake outburst floods. However, the eastern Tanggula Mountains, which are the only area where lake-terminating glaciers are found within the Tibetan Plateau, have received little attention to date. In this study, to address this gap, we generated updated glacier boundaries and estimated the interdecadal area changes for 2000–2020 based on the interpretation of Landsat-5/8 and Sentinel-2 images. In addition, based on the method of digital elevation model (DEM) differencing, we quantified the changes in glacier thickness and mass balance using TanDEM-X radar data and SRTM DEM over almost the same periods. The final results show that the glaciers in the eastern Tanggula Mountains, as a whole, have experienced accelerated area shrinkage (with a rate of area loss increasing from −0.34 ± 0.83 km² a⁻¹ to −0.93 ± 0.81 km² a⁻¹ for 2000–2013 and 2013–2020, respectively) and accelerated ice thinning (changing from −0.19 ± 0.05 m a⁻¹ and −0.53 ± 0.08 m a⁻¹ for 2000−2012 and 2012–2020, respectively). Furthermore, the region-wide glacier mass balance was −0.16 ± 0.04 m w.e. a⁻¹ and −0.45 ± 0.07 m w.e. a⁻¹ for these two sub-periods, corresponding to a 1.8 times acceleration of mass loss rate. The average mass balance during 2000–2020 was −0.23 ± 0.04 m w.e. a⁻¹, which is equivalent to a rate of mass loss of −0.04 Gt a⁻¹. More specifically, within the region, the lake-terminating glaciers have exhibited more significant acceleration of area loss and mass loss, compared to the land-terminating glaciers. However, interestingly, the average thinning rate of the lake-terminating glaciers is always lower than that of the land-terminating glaciers over all study periods, which is in contrast with previous findings in other high mountain areas (e.g., the Himalaya Mountains). Field study and proglacial lakes monitoring suggest that the local topography plays a vital role in the evolution of the glacial lakes in this region, which further affects the glacier changes. Furthermore, the present status of the glacier changes in this region can be attributed to the long-term increase in air temperature. Our findings provide a comprehensive overview of the current state of glacier changes across the eastern Tanggula Mountains and will help to improve the understanding of the heterogeneous response of glaciers to climate change.

1. Introduction

As the atmosphere warms, glaciers in the Tibetan Plateau and its surroundings have experienced serious mass wastage, especially in the Himalaya and Nyainqentanglha ranges [1,2,3,4]. Meanwhile, glacier retreat, given appropriate local topography, leads to the gradual formation and expansion of the proglacial lakes at the terminus of some glaciers. Since the 1990s, the area of proglacial lakes has increased to 455 km² in the Himalayas, corresponding to an expansion of ~14% [5]. The formation and development of proglacial lakes makes the response of the connected glaciers to climate change more complicated, as this can lead to more mass loss through the processes of the frontal calving and subglacial ablation, and can further lead to an increase in ice flux at the glacier front, causing so-called dynamic thinning [6,7,8,9]. In addition, glacier variations and the expansion of proglacial lakes can sometimes induce an outburst event or mountain flood, especially under the influence of certain extreme events (e.g., earthquake, ice avalanche, landslide, heavy rainfall, severe frontal calving, and so on), posing a high risk to the safety of downstream residents, infrastructure and engineering projects [10,11,12]. Some studies have predicted that the glaciers of the TPs will continue to thin and retreat [13], and that the threats from glacial lakes outburst flood will intensify across the TPs [14]. Hence, at present, the changes of lake-terminating glaciers have attracted more and more attention, and thus conducing relevant studies is of significance.

For lake-terminating glaciers of the Tibetan Plateau area, researchers focusing on glacier mass change or area change have detected enhanced thinning rates or greater area loss, compared to land-terminating glaciers in the central-eastern Himalaya, from 2000 to the 2010s [15,16]. An explicit comparative analysis for these two types of glaciers is therefore necessary and important, as it could improve the understanding of the main driving factors of mass deficit. In recent years, some studies have quantitatively assessed the difference of the mass change for the two types of glaciers at a regional scale, and found that the rate of mass loss for lake-terminating glaciers from 2000 to the 2010s, on average, was 32% and 78% greater than the mass loss for land-terminating glaciers in the Everest [17] and eastern Nyainqentanglha [18] areas, respectively. For individual glaciers, the difference of the mass balance between the different glacier types was found to be more than three times in the Bhutan Himalaya [6]. These findings indicate that the magnitude of the impact of proglacial lakes on glacier mass change is heterogeneous in different regions. More accurately, the development stages of proglacial lakes are distinct in different regions [19]. Furthermore, benefiting from the utilization of historical stereo images (i.e., KH-7, KH-9), across the whole of the Himalayas, comprehensive evaluations for a long-term period further confirmed that the above phenomenon has been active since the 1970s, and found that the effect of exacerbated mass loss from proglacial lakes has been more significant since 2000 [20,21]. Overall, the variation of lake-terminating glaciers accounts for 32% of the region-wide mass wastage in the Himalaya, despite occupying only a small part of the total glaciers. Based on these findings, it is speculated that the sustained expansion of glacial lakes and the formation of new lakes will further enhance ice mass loss from this region in the coming decades [20].

However, to date, the eastern Tanggula Mountains, which are the only area where lake-terminating glaciers are found in the interior of the Tibetan Plateau, have attracted little attention. Furthermore, to better understand the changing status of lake-terminating glaciers, time-series monitoring of the glacier evolution process is essential, although this is presently extremely difficult, due to the challenges of data acquisition. However, an investigation at an appropriate timescale (e.g., the decadal change) is an effective complementary way to understand the evolution of lake-terminating glaciers. Since 2000, only a few local regions, such as the Geladandong and Nyainqentanglha mountains, have been investigated using advanced stereo images (e.g., ZY-3) and the TanDEM-X DEM [22,23,24]. Despite this, knowledge of interdecadal changes is still lacking.

Aiming at the current issue (i.e., how the glaciers of the eastern Tanggula Mountains have evolved over the first 20 years of the 21st century, especially at an interdecadal time scale), the main purpose of this study was: (1) to obtain updated glacier boundaries for 2000, 2013, and 2020 and estimate the rates of area change; and (2) to quantify the decadal changes in glacier thickness and mass balance for 2000–2020 based on the use of the geodetic method with TanDEM-X data. We hope that the findings of this study will provide a comprehensive assessment of the present status of the glaciers in the eastern Tanggula Mountains.

2. Study Region

In the eastern Tanggula Mountains, benefiting from the relatively high altitude and abundant precipitation, lots of mountain glaciers form. According to the records of the Randolph Glacier Inventory (RGI) inventory (version 6.0), there are about 120 glaciers with a total of area of ~163 km² in the region. The main large-scale glaciers are basically distributed on north- and south-facing slopes within the region, and are typically connected with glacial lakes (Figure 1). Moreover, the glacier termini on the south-facing slopes extend to the lower altitudes (i.e., 4200 m a.s.l.) than the glacier termini on the north-facing slopes (above 4600 m a.s.l.). Analyzing the ERA5-Land reanalysis data for this region shows that the annual average precipitation was about 940 mm for 1981–2020, of which the summer precipitation was 485 mm, accounting for 52% of the total precipitation. Moreover, the mean annual air temperature (MAAT) was −7.72 °C, and the average temperature was 5.58 °C and −20.24 °C in summer and winter, respectively. In addition, based on the median elevation information of the RGI inventory, the region-wide equilibrium line altitude (ELA) was about 5400 ± 120 m.

3. Data and Methods

In this study, we collected two sets of TanDEM-X CoSSC data from 14 November 2012 and 02 January 2020, to monitor the changes in glacier thickness and mass balance at the decadal timescale. The C-band SRTM DEM acquired in February 2000 was used as the benchmark data for the estimation of the glacier mass balance. Specifically, the one arc-second global SRTM DEM was employed for calculating the glacier thickness change. The newly released NASA DEM, which is an improved version combining multiple-source topographic data, was used to assist with the generation of the TanDEM-X DEM. The X-band SRTM DEM was also used to estimate the penetration depth difference between the C-band and X-band radar in the glacierized areas. In addition, multiple optical satellite images, including Sentinel-2 and Landsat 5/8 images, were applied to retrieve the boundaries of the glaciers and glacial lakes. We also downloaded the ERA5-Land reanalysis data to analyze the long-term climate change in the study region. More detailed information with regard to each type of data is provided in

Table 1. The remote sensing data sets used in this study.

Sensor	Path/Row	Date of Acquisition	Spatial Resolution	Purpose
SRTM C-band DEM		2000/02	30.0 m 90.0 m	Reference DEM Elevation change
SRTM X-band DEM		2000/02	25 m	Penetration depth estimation
TanDEM-X		2012/11/14 2020/01/02	~2.0 m	Elevation change
Landsat-5	136/38	2000/07/22 2001/07/09	30 m	Glacier boundaries
Landsat-8	136/38	2013/09/28 2014/11/18	15 m	Glacier boundaries
Sentinel-2		2020/10/30	10 m	Glacier boundaries
ERA5-Land		1981–2020	0.1° × 0.1°	Climate analysis

.

TanDEM-X DEM generation: We followed the commonly-used differential interferometry method used in previous studies to process the TanDEM-X data and generate the corresponding DEMs [25,26,27]. The first step was to directly generate the interferogram, with a multi-look factor of 4 and 5 for the range and azimuth directions, respectively. The NASA DEM was then used for modeling the topography-dependent phase, and was further subtracted from the original interferogram to generate the differential phase. The second critical step was to unwrap the differential phase using the minimum cost flow method, with a coherence threshold of 0.7. Notably, during the data processing, we found that phase jumps clearly appeared in local areas for the TanDEM-X data in 2020, despite the relatively flat slope (Figure 2a). To address this problem, we manually delineated the boundary of the area with phase jumps, and then added 2pi to this area. The corrected differential phase is shown in Figure 2b. Next, the unwrapped differential phase was converted to the height difference using the phase-to-height transformation formula [25]. Finally, the TanDEM-X DEMs with a spatial resolution of 10 m were obtained by adding the height difference to the referenced DEM.

Data co-registration and bias correction: After obtaining the two TanDEM-X DEMs, we co-registered these two DEMs to the reference SRTM DEM. In order to avoid the impact of spatial trend bias on the data co-registration as much as possible, we firstly checked and corrected the trend bias using the approach of iterative polynomial fitting [2]. In particular, to accurately model the trend bias during this step, we identified the layover and shadow areas using the incidence angle information for the TanDEM-X data from two different orbits [28], and excluded pixels in these potential distorted areas. Subsequently, the detrended TanDEM-X DEMs were co-registered to the referenced SRTM DEM using the 3-D analytical method [29]. The final elevation difference maps for the ice-free areas are shown in Figure 3, which also includes the corresponding statistical information for the three-periods (2000–2012, 2012–2020, 2000–2020).

Estimation and correction of radar penetration depth: In this study, we employed the simultaneously acquired SRTM X-band and C-band DEM data and performed the differencing process to estimate the penetration depth difference between the above two radar bands [2,30,31]. The first step was to convert the ellipsoid height (X-band DEM) to orthometric height (C-band DEM) by subtracting the geoid height, which was estimated to be, on average, −38.4 m in this region. The second step was to perform data co-registration and bias correction, to estimate the penetration depth difference in the glacierized areas. The average penetration depth difference we obtained is 1.80 m for this region. As for the penetration depth correction, similar to previous studies, we applied the strategy of using the average penetration for each 100-m altitude band, as shown in Figure 4. From the figure, it is clear that the penetration depth difference shows an approximately linear trend with increasing altitude, and reaches ~3 m above 6000 m a.s.l.. It should be noted that, given that the two TanDEM DEMs of this study were acquired in the same season (i.e., 11 November and 2 January), we did not carry out the penetration correction for 2012–2020, assuming a negligible penetration difference within a season. This can be supported by the climatic conditions in this region, in that the precipitation (i.e., snowfall) in winter is not significant.

Calculation of glacier thickness change and mass balance: In this study, we selected some representative glaciers to analyze their glacier thickness changes over the accumulation and ablation zones. During this step, the ELA of each glacier was determined by the median elevation information from the RGI inventory [17], which represents an approximate long-term ELA. For the calculation of the glacier mass balance, the first step was to calculate the glacier area and thickness change for each 100-m altitude band, and to integrate them to obtain the overall volume change [16]. A density conversion factor of 850 ± 60 kg m⁻³ was then used to convert the volume change to the mass change [32].

Extraction of glacier boundaries: For the extraction of the glacier boundaries, taking the high-quality Landsat-5/8 and Sentinel-2 true-color images as the referenced base map, we manually edited and adjusted the boundaries from the RGI inventory (version 6.0) to generate new glacier boundaries for 2000, 2013, and 2020.

Climate change analysis: Monthly ERA5-Land reanalysis data for 1981–2020 with a grid resolution of 0.1° were used for analyzing the region-wide climate change, with the average of 30 years (1981–2010) used as the reference [33]. For the glacierized areas of this study, we selected the data centered by the glacierized areas with a window of 3 × 3 and carried out the climate change detection, considering that individual pixels may lack representation. Only two variables (total precipitation and 2 m air temperature) were used.

Uncertainty assessment: For the uncertainty of the glacier boundaries, this was evaluated by the product of the boundary perimeter and the image resolution (30 m, 15 m and 10 m for 2000, 2013, and 2020, respectively), based on an assumption of one pixel error for the glacier and lake boundaries [34]. Furthermore, the uncertainty of the area change (${\sigma }_S$) was calculated based on the principle of error propagation.

$${\sigma }_S=\sqrt{\Delta S_1+\Delta S_2}$$

(1)

where $\Delta S_1$ and $\Delta S_2$ represent the uncertainty of the glacier boundaries at the initial and final moments, respectively.

With regard to the uncertainty of the geodetic mass balance, it was estimated as a combination of the uncertainties related to the glacier thickness change, the glacier extent and the conversion factor, mainly based on the law of error propagation [15,35,36]. More specifically, the uncertainty of the glacier thickness changes (${\sigma }_{\Delta h}$) for each altitude band ($i$) in this study was further decided by the uncertainty of the initial glacier elevation difference (${\sigma }_{\Delta he}$) and the uncertainty of the penetration depth estimation (${\sigma }_{\Delta hp}$), assuming that they are independent of each other [35]. Furthermore, for the uncertainty of the above two elevation differences, we used the standard deviation of the elevation difference (${\sigma }_{\Delta hx}$) in ice-free areas as an approximate representation. Accordingly, the final uncertainty of the thickness variations can be calculated by Equations (2) and (3) based on the law of error propagation:

$${\sigma }_S=\sqrt{\Delta S_1+\Delta S_2}$$

(2)

$${\sigma }_{\Delta hx,i}=\{\begin{array}{c}{\sigma }_{\Delta hx,i}, S_i\le S_{cor,i}\\ {\sigma }_{\Delta hx,i}\sqrt{\frac{S_{cor,i}}{5S_i}},S_i>S_{cor,i}\end{array}$$

(3)

where x represents the type of elevation difference (i.e., “$e$” for the initial glacier elevation difference map, and “$p$” for the penetration depth difference map). $S$ is the glacier area for each altitude band, and $S_{cor}$ denotes the effective correlation area, which were computed as $\pi d^2$, where $d$ denotes the autocorrelation distance for ice-free areas [36]. In this study, the autocorrelation distances were assumed to be 330 m and 240 m for the glacier elevation difference map and the penetration depth difference map, respectively [35]. Subsequently, the uncertainties of the volume change (${\sigma }_{\Delta V}$) and mass balance (${\sigma }_{\Delta m}$) can be respectively estimated by Equation (4) and Equation (5):

$${\sigma }_{\Delta V}=\sqrt{{{\displaystyle \sum }}_{i=1}^n{\left({\sigma }_{\Delta h,i}\times S_i\right)}^2}$$

(4)

$${\sigma }_m=\sqrt{{\left(\frac{\mathsf{\Delta }V\times {\sigma }_f}{S_t}\right)}^2+{\left(\frac{{\sigma }_{\mathsf{\Delta }v}\times f}{S_t}\right)}^2+{\left(\frac{\mathsf{\Delta }V\times f\times {\sigma }_S}{S_t}\right)}^2}$$

(5)

where $\Delta V$ and $S_t$ represent the volume change and the total glacier area. $f$ and ${\sigma }_f$ represent the conversion factor (850 kg m⁻³) and its uncertainty (60 kg m⁻³), respectively [32].

4. Results

4.1. Glacier Area Change

For the changes in glacier area, as shown in Figure 5 and

Table 2. The changes in glacier areas and the change rates in the study region. The asterisk in the first column represents glaciers that are connected to proglacial lakes.

Glacier Name	Glacier Area (km2)			Area Change (km2)			Rate of Area Change (km2 a−1)
	2000	2013	2020	2000–2013	2013–2000	2000–2020	2000–2013	2013–2000	2000–2020
G59	4.61 ± 0.29	4.32 ± 0.15	4.10 ± 0.10	−0.29 ± 0.33	−0.22 ± 0.17	−0.51 ± 0.31	−0.02 ± 0.02	−0.03 ± 0.02	−0.02 ± 0.01
G88	1.28 ± 0.14	1.10 ± 0.07	1.07 ± 0.04	−0.18 ± 0.14	−0.03 ± 0.08	−0.21 ± 0.14	−0.01 ± 0.01	0.00 ± 0.01	−0.01 ± 0.01
G71 *	2.33 ± 0.27	2.30 ± 0.14	2.16 ± 0.09	−0.03 ± 0.30	−0.14 ± 0.16	−0.17 ± 0.28	0.00 ± 0.02	−0.02 ± 0.02	−0.01 ± 0.01
G22 *	2.17 ± 0.27	2.06 ± 0.14	1.89 ± 0.09	−0.11 ± 0.30	−0.17 ± 0.16	−0.28 ± 0.28	−0.01 ± 0.02	−0.02 ± 0.02	−0.01 ± 0.01
G21 *	5.39 ± 0.47	5.21 ± 0.23	4.89 ± 0.16	−0.18 ± 0.51	−0.32 ± 0.27	−0.50 ± 0.49	−0.01 ± 0.04	−0.04 ± 0.04	−0.02 ± 0.02
ZGMN *	15.81 ± 0.77	15.50 ± 0.39	15.02 ± 0.26	−0.31 ± 0.86	−0.48 ± 0.46	−0.79 ± 0.81	−0.03 ± 0.06	−0.06 ± 0.06	−0.04 ± 0.04
West Beijia *	8.59 ± 0.57	8.20 ± 0.28	7.85 ± 0.19	−0.39 ± 0.63	−0.35 ± 0.34	−0.74 ± 0.59	−0.03 ± 0.04	−0.05 ± 0.05	−0.04 ± 0.03
Beijia *	13.58 ± 0.57	13.13 ± 0.28	12.73 ± 0.19	−0.45 ± 0.63	−0.40 ± 0.34	−0.85 ± 0.59	−0.04 ± 0.04	−0.05 ± 0.05	−0.04 ± 0.03
East Beijia *	10.48 ± 0.60	10.15 ± 0.30	9.55 ± 0.20	−0.33 ± 0.66	−0.60 ± 0.35	−0.93 ± 0.62	−0.02 ± 0.05	−0.09 ± 0.05	−0.05 ± 0.03
Zuxuehui *	35.87 ± 1.23	35.15 ± 0.62	35.00 ± 0.41	−0.72 ± 1.37	−0.15 ± 0.73	−0.87 ± 1.29	−0.05 ± 0.10	−0.02 ± 0.10	−0.04 ± 0.06
Poge *	21.11 ± 1.20	20.91 ± 0.60	20.64 ± 0.40	−0.20 ± 1.33	−0.27 ± 0.71	−0.47 ± 1.25	−0.01 ± 0.10	−0.04 ± 0.10	−0.02 ± 0.06
All glaciers	143.23 ± 9.69	138.83 ± 4.73	132.35 ± 3.17	−4.40 ± 10.78	−6.48 ± 5.69	−10.88 ± 10.19	−0.34 ± 0.83	−0.93 ± 0.81	−0.54 ± 0.51

, overall, glaciers in the eastern Tanggula Mountains exhibited accelerated area loss, but the rate and evolution of glacier retreat is heterogenous within this region. Specifically, the average rate of the area loss has increased from −0.34 ± 0.83 km² a⁻¹ (for 2000–2013) to 0.93 ± 0.81 km² a⁻¹ (for 2013–2020) over the whole region. Among the different glaciers, the most prominent increase in rate of area loss occurred at the East Beijia Glacier, with a rate increasing from −0.02 ± 0.05 km² a⁻¹ to −0.09 ± 0.05 km² a⁻¹ for the two study periods mentioned above. For most other glaciers, the increase in the rate of area loss ranges from 0.01 km² a⁻¹ to 0.03 km² a⁻¹ (i.e., the difference between the rates of the two periods) (Table 2), and it seems that the acceleration at lake-terminating glaciers (0.03 km² a⁻¹), in general, is more significant than that at land-terminating glaciers (0.01 km² a⁻¹). In particular, a contrary pattern can be observed at the Zuxuehui Glacier (lake-terminating type), where the rate of area loss has clearly declined from −0.05 ± 0.10 km² a⁻¹ to −0.02 ± 0.10 km² a⁻¹ for the above two periods. For the whole study period (i.e., 2000–2020), the average rate of area loss is −0.54 ± 0.51 km² a⁻¹, and it can be found that the lake-terminating glaciers generally show a higher shrinkage rate (ranging from −0.04 to −0.05 km² a⁻¹) than the land-terminating glaciers (ranging from −0.01 to −0.02 km² a⁻¹). It is worth noting that Poge Glacier, which belongs to the lake-terminating type, experienced an insignificant area loss at a rate of −0.02 ± 0.06 km² a⁻¹ for 2000–2020, which is comparable to the changes of the land-terminating glaciers.

4.2. Glacier Thickness Change

With regard to glacier thickness change, as shown in Figure 6, enhanced surface thinning can be clearly observed at an approximately inter-decadal time scale, especially in the ablation zone. Correspondingly, the rate of glacier surface lowering of the ablation zone has increased from −0.73 ± 0.07 m a⁻¹ (during 2000–2012) to −0.99 ± 0.10 m a⁻¹ (during 2012–2020) (Figure 6d,

Table 3. Glacier thickness change and mass balance for the different periods in the eastern Tanggula region. The “Acc.” and “Abl.” are the abbreviations for accumulation and ablation, respectively. The asterisk in the first column represents the glaciers that are connected to proglacial lakes.

Glacier Name	ELA (m)	2000–2012			2012–2020			2000–2020
		Thickness Change Rate (m a−1)		Mass Balance (w.e. a−1)	Thickness Change Rate (m a−1)		Mass Balance (w.e. a−1)	Thickness Change Rate (m a−1)		Mass Balance (w.e. a−1)
		Acc. Zone	Abl. Zone		Acc. Zone	Abl. Zone		Acc. Zone	Abl.
G59	5400	−0.02 ± 0.44	−0.51 ± 0.38	−0.22 ± 0.21	−0.24 ± 0.60	−0.75 ± 0.48	−0.40 ± 0.24	0.01 ± 0.30	−0.53 ± 0.27	−0.21 ± 0.14
G61	5294	−0.39 ± 0.71	−0.83 ± 0.61	−0.56 ± 0.34	−0.15 ± 0.88	−1.01 ± 0.73	−0.54 ± 0.36	−0.41 ± 0.53	−0.93 ± 0.43	−0.58 ± 0.23
G88	5364	−0.27 ± 0.79	−1.10 ± 0.74	−0.65 ± 0.35	−0.37 ± 0.88	−0.88 ± 0.86	−0.59 ± 0.36	−0.21 ± 0.52	−1.01 ± 0.49	−0.58 ± 0.23
G71 *	5383	0.10 ± 0.80	−0.58 ± 0.55	−0.20 ± 0.23	−0.02 ± 0.88	−0.79 ± 0.65	−0.33 ± 0.25	0.27 ± 0.54	−0.57 ± 0.36	−0.12 ± 0.15
ZGMN *	5560	0.18 ± 0.23	−0.70 ± 0.20	−0.21 ± 0.12	−0.14 ± 0.28	−0.55 ± 0.25	−0.29 ± 0.14	0.10 ± 0.13	−0.54 ± 0.14	−0.18 ± 0.08
W-Beijia *	5587	0.29 ± 0.28	−0.63 ± 0.26	−0.14 ± 0.15	−0.21 ± 0.34	−1.13 ± 0.32	−0.53 ± 0.18	0.12 ± 0.19	−0.64 ± 0.18	−0.22 ± 0.10
Beijia *	5593	0.40 ± 0.24	−0.57 ± 0.22	−0.07 ± 0.13	−0.15 ± 0.28	−0.75 ± 0.26	−0.41 ± 0.15	0.20 ± 0.16	−0.44 ± 0.15	−0.11 ± 0.08
E-Beijia *	5345	−0.01 ± 0.27	−1.12 ± 0.25	−0.48 ± 0.16	−0.17 ± 0.34	−1.41 ± 0.29	−0.67 ± 0.18	−0.01 ± 0.19	−1.24 ± 0.16	−0.52 ± 0.12
G21 *	5615	0.30 ± 0.45	−0.61 ± 0.38	−0.13 ± 0.20	−0.41 ± 0.69	−1.22 ± 0.48	−0.66 ± 0.23	−0.10 ± 0.29	−0.80 ± 0.25	−0.39 ± 0.14
G22 *	5486	0.30 ± 0.79	−0.20 ± 0.62	0.10 ± 0.21	−0.48 ± 0.88	−1.77 ± 0.75	−0.91 ± 0.25	−0.16 ± 0.50	−0.73 ± 0.40	−0.35 ± 0.15
Zucxuehui *	5726	0.33 ± 0.14	−0.16 ± 0.14	0.07 ± 0.08	−0.27 ± 0.18	−0.76 ± 0.18	−0.42 ± 0.10	0.09 ± 0.09	−0.34 ± 0.09	−0.10 ± 0.05
Poge *	5890	0.40 ± 0.18	−0.35 ± 0.20	0.05 ± 0.10	−0.21 ± 0.22	−0.45 ± 0.23	−0.25 ± 0.12	0.18 ± 0.12	−0.26 ± 0.13	−0.01 ± 0.06
Lake-terminating	5576	0.33 ± 0.08	−0.57 ± 0.08	−0.08 ± 0.05	−0.21 ± 0.10	−0.86 ± 0.10	−0.43 ± 0.07	0.12 ± 0.05	−0.57 ± 0.05	−0.17 ± 0.04
Land-terminating	5390	−0.18 ± 0.15	−0.70 ± 0.11	−0.38 ± 0.08	−0.29 ± 0.27	−0.91 ± 0.17	−0.54 ± 0.10	−0.19 ± 0.13	−0.73 ± 0.10	−0.42 ± 0.07
All glaciers	5400	0.17 ± 0.06	−0.73 ± 0.07	−0.16 ± 0.04	−0.23 ± 0.09	−0.99 ± 0.10	−0.45 ± 0.07	0.06 ± 0.04	−0.76 ± 0.05	−0.23 ± 0.04

). Meanwhile, over the accumulation zone, there has been a prevailing shift from surface thickening (at a rate of 0.17 ± 0.06 m a⁻¹) to surface thinning (at a rate of −0.23 ± 0.09 m a⁻¹) for the above two decadal periods. Overall, the whole region experienced accelerated surface thinning during the first two decades of the 21st century, as the average rate of glacier thickness change has changed from −0.19 ± 0.05 m a⁻¹ and −0.53 ± 0.08 m a⁻¹ for the two sub-periods, corresponding to a 1.8 times acceleration. The mean thickness change of the whole study period (i.e., 2000–2020) is −0.27 ± 0.05 m a⁻¹, for which the rates of average thickness change are 0.06 ± 0.04 m a⁻¹ and −0.76 ± 0.05 m a⁻¹ over the accumulation and ablation zones, respectively, suggesting a nearly stable state in high altitudes and a remarkable surface lowering in lower altitudes (Figure 6).

To better understand the variations of the different types of glaciers within the study region, we separately calculated the thickness changes for the land-terminating glaciers (with an area of 47.41 km²) and lake-terminating glaciers (114 km²), as shown in Figure 7. Clearly, the fundamental trends of the thickness variations are almost identical for both types, which both showed accelerated surface thinning. However, the acceleration of the glacier changes (i.e., the ratio of the increase of surface thinning rate to the initial change rate) are different for the two types of glaciers. The land-terminating glaciers exhibited a slight acceleration (~45%) across all of the altitude bands (Figure 7a), as the glaciers thinned at a rate of −0.44 ± 0.09 m a⁻¹ and −0.64 ± 0.12 m a⁻¹ during the two decadal periods. In contrast, there was a 4.7 times acceleration for the lake-terminating glaciers, for which the thinning rates increased from −0.09 ± 0.06 m a⁻¹ to −0.51 ± 0.08 m a⁻¹ for the two sub-periods. From Figure 7b, it can be seen that there has been a remarkable increase in the thinning rates at the fronts of the lake-terminating glaciers, especially in areas with an altitude of less than 5000 m a.s.l. This may largely account for the corresponding significant acceleration for the lake-terminating glaciers. In addition, interestingly, the average change rates of the land-terminating glaciers are higher than those of the lake-terminating glaciers for all of the study periods (Table 3), e.g., −0.50 ± 0.08 m a⁻¹ vs. −0.20 ± 0.05 m a⁻¹ for the whole study period.

Moreover, from Figure 6, we can see that the pattern of glacier thickness change varies significantly between individual glaciers. To more clearly reveal the changing characteristics of individual glaciers, we further selected some representative glaciers, as shown in Figure 8 and Table 3. Clearly, the thickness change rate of the land-terminating glaciers (e.g., G59, G61, and G88) shows an approximately linear trend with the increase of altitude (Figure 8a–c). However, in terms of the change rates, except for the G59 Glacier with an accelerated surface lowering (from −0.26 ± 0.25 m a⁻¹ to −0.47 ± 0.28 m a⁻¹), the other two glaciers underwent continuous thinning at a nearly constant rate (Table 3). For the nine lake-terminating glaciers, clearly, the most remarkable thinning rate of −4 m a⁻¹ to −6 m a⁻¹ took place at the termini of the Zuxuehui Glacier (below 4500 m a.s.l.) and G22 Glacier between 2012 and 2020 (Figure 8e–f), which accounts for the aggravated thinning in lower altitudes for the whole region (Figure 6d and Figure 7b). Similar conditions can also be found at the West Beijia, Beijia, and East Beijia glaciers (Figure 6 and Figure 7j–l), whose fronts, at altitudes of approximately 4700 m to 5100 m, generally kept a higher rate of surface lowering than the other glaciers, and exhibited more notable acceleration. For the G22 Glacier in particular, all of the altitude bands experienced significantly accelerated thinning (Figure 8e), and the corresponding thickness change rates changed from 0.12 ± 0.25 m a⁻¹ to −1.07 ± 0.29 m a⁻¹ for the two sub-periods. However, unexpectedly, some lake-terminating glaciers did not exhibit such a change pattern. For example, at the Poge Glacier, ZGMN Glacier, and the G71 Glacier, it appears that the rate of glacier thickness change varies linearly with increasing altitude for the above two periods. Interestingly, this pattern is highly similar to that of the land-terminating glaciers (where there has been stable or accelerated thinning).

4.3. Glacier Mass Balance

With regard to the glacier mass balance at the interdecadal time scale, the change trend is consistent with that of the glacier thickness change described above. The region-wide glacier mass balance was estimated to be −0.16 ± 0.04 m w.e. a⁻¹ (2000–2012), −0.45 ± 0.07 m w.e. a⁻¹ (2012–2020), and −0.23 ± 0.04 m w.e. a⁻¹ (2000–2020), which corresponds to a mean rate of mass loss of −0.04 ± 0.01 Gt a⁻¹ over the past 20 years. In addition, the calculated mass balance of the lake-terminating glaciers was estimated to be −0.08 ± 0.05 m w.e. a⁻¹, −0.43 ± 0.07 m w.e. a⁻¹, and −0.17 ± 0.04 m w.e. a⁻¹ for the above three study periods, while the result for the land-terminating glaciers was estimated to be −0.38 ± 0.08 m w.e. a⁻¹, −0.54 ± 0.10 m w.e. a⁻¹, and −0.42 ± 0.07 m w.e. a⁻¹, respectively. For the individual glaciers we selected, the glacier mass balances we obtained ranged from −0.65 ± 0.35 m w.e. a⁻¹ to 0.10 ± 0.21 m w.e. a⁻¹ for 2000–2012, from −0.91 ± 0.25 m w.e. a⁻¹ to −0.25 ± 0.12 m w.e. a⁻¹ for 2012–2000, and from −0.58 ± 0.23 m w.e. a⁻¹ to −0.01 ± 0.06 m w.e. a⁻¹ for 2000–2020.

5. Discussion

5.1. Impact of the Evolution of Proglacial Lake on Glacier Change

With regard to the heterogeneous pattern of the thickness changes of lake-terminating glaciers in this study region, although the results of the glacier area change, to some extent, can account for the thickness change, we further analyzed the evolution of the proglacial lakes over the past 20 years, as shown in Figure 9. At the terminus of the Zuxuehui Glacier, the proglacial lake experienced a rapid expansion in ~2005/2006. This was caused by the collapse of the glacier tongue area, where the calving ice did not melt completely until 2013, which largely accounts for the notable thinning in areas closely to the proglacial lake (Figure 6). At the adjacent Poge Glacier, the extent of the proglacial lake did not show pronounced variation, despite being connected with the glacier terminus. Based on an on-site investigation in August 2021, we found that the local topography restricted the development of the proglacial lake (Figure 10). The elevation of the glacier tongue is clearly higher than the lake level, although the glacier terminus is connected to the glacial lake and there is ice calving into the lake (Figure 10). Accordingly, the influence of the proglacial lake on the mass change of the Poge Glacier is insignificant. As for the ZGMN Glacier, the proglacial lake did not expand prominently after about 2011, indicating the observed separation of the glacier from the lake. Similar conditions also occurred at the East Beijia Glacier, G22 Glacier, and G71 Glacier (Figure 9), where all three proglacial lakes experienced rapid expansion and then remained stable, due to the glacier terminus gradually becoming detached. This turns lake-terminating glaciers into land-terminating glaciers. Consequently, the rate of the thickness changes at their tongues is less negative than at the other land-terminating glaciers. In particular, at the West Beijia Glacier and Beijia Glacier, the two proglacial lakes have continued to expand over the last 20 years, corresponding to significant surface thinning and more rapid acceleration. Nevertheless, based on the field expedition we conducted, we speculate that it is likely that the Beijia Glacier will transform from lake-terminating type into land-terminating type, as the bedrock at the end of the glacier has been exposed, and the elevation of the glacier is evidently higher than the current lake level (Figure 10). To sum up, our findings further confirm that the evolution of the proglacial lakes does indeed have an important impact on the glacier mass loss, especially at the terminus [17,19,20]. Furthermore, we also found that the development of the glacial lake is mainly controlled by the local topography, which varies significantly between different regions. In addition, we speculate that most of the eastern Tanggula glaciers will become separated from the proglacial lakes in the near future.

5.2. Climate Analysis

From the perspective of regional climate change, as shown in Figure 11, we found that the air temperature has shown a statistically significant increase across the whole of the eastern Tanggula region and its surroundings (with a warming rate of ranging from 0.12 °C dec⁻¹ to 0.39 °C dec⁻¹) over the past four decades (1981–2020), while the long-term trend of the precipitation has been statistically insignificant in most areas. More specifically, for the eastern Tanggula region, the rate of increase of air temperature has reached 0.16 °C dec⁻¹ (statistically significant) (Figure 11c), but there has not been a statistically significant change in precipitation. These findings are basically consistent with the measurements of nearby meteorological stations [37]. Hence, we believe that the accelerated area loss and mass loss at the decadal time scale can be attributed to the increase in air temperature.

5.3. Comparison with Previous Studies

With regard to the glacier mass balance, given that there are no field-based measurements available in eastern Tanggula region, we made a comparison with adjacent regions, as summarized in

Table 4. Results for the geodetic mass balance since 2000 in areas adjacent to the eastern Tanggula region. Notes that this table includes the results for glacier mass balance only after 2000.

Region	Data	Periods	Glacier Mass Balance (m w.e. a−1)	Sources
Western Tanggula (Geladandong)	SRTM DEM, TanDEM-X, ZY-3	2000–2012 2012–2018 2000–2018	–0.11 ± 0.03 –0.47±0.09 –0.24 ± 0.07	Liu et al. [23]
	SRTM DEM, TanDEM-X	2000–2011 2011–2017 2000–2017	–0.24 ± 0.12 –0.41 ± 0.09 –0.30 ± 0.08	Chen et al. [22]
Western Nyainqentanglha	SRTM DEM, TanDEM-X	2000–2013	−0.24 ± 0.13	Li et al. [38]
	SRTM DEM, TanDEM-X	2000–2014	–0.25 ± 0.08	Luo et al. [39]
	SRTM DEM, ZY-3	2000–2013 2013–2017 2000–2017	−0.22 ± 0.23 −0.43 ± 0.06 −0.30 ± 0.19	Ren et al. [24]
Eastern Nyainqentanglha (Yigong Zangbo)	SRTM DEM, TanDEM-X	2000–2013	−0.60 ± 0.20 (all)	Wu et al. [40]
	SRTM DEM, TanDEM-X DEM, AW3D	2000–2010s	−0.66 ± 0.19 (all) −0.50 ± 0.32 (land) −0.89 ± 0.36 (lake)	Ke et al. [18]
Himalaya	SRTM DEM HMA DEM	2000–2015	–0.39 ± 0.12 (all) –0.37 ± 0.12 (land) –0.55 ± 0.12 (lake)	King et al. [20]
	ASTER	2000–2016	–0.43 ± 0.14 (all) –0.41 ± 0.08 (land) –0.56 ± 0.08 (lake)	Maurer et al. [21]

. Specifically, our results (–0.16 ± 0.04 m w.e. a⁻¹, –0.45 ± 0.07 m w.e. a⁻¹, and –0.23 ± 0.04 m w.e. a⁻¹ for 2000–2012, 2012–2020, and 2000–2020, respectively) are basically consistent with those for the western Tanggula Mountains (i.e., the Geladandong region, with a distance of 370 km) over similar study periods. For example, Liu et al. [23] reported that the mass balance was –0.11 ± 0.03 m w.e. a⁻¹, –0.47 ± 0.09 m w.e. a⁻¹, and –0.24 ± 0.07 m w.e. a⁻¹ for 2000–2012, 2012–2018 and 2000–2018, respectively, while the results of Chen et al. [22] were –0. 24 ± 0.12 m w.e. a⁻¹, –0.41 ± 0.09 m w.e. a⁻¹ and –0.30 ± 0.08 m w.e. a⁻¹ for the periods of 2000–2011, 2011–2017, and 2000–2017, respectively. In addition, the results we obtained are also comparable to the findings for the western Nyainqentanglha Mountains (400 km apart), for which the mass balance results range from −0.22 ± 0.23 m w.e. a⁻¹ to –0.25 ± 0.08 m w.e. a⁻¹ for 2000–2014 (Table 4) [24,38,39], and the rates of mass wastage for 2013–2017 and 2000–2017 are −0.43 ± 0.06 m w.e. a⁻¹ and −0.30 ± 0.19 m w.e. a⁻¹, respectively [24]. In contrast, from 2000 to the 2010s, our result for the eastern Tanggula area (–0.16 ± 0.04 m w.e. a⁻¹) is significantly less negative than that for the Yigong Zangbo Basin in the eastern Nyainqentanglha Mountains, despite the short distance (~140 km), where the mass balance was estimated to be −0.60 ± 0.20 m w.e. a⁻¹ [40] and −0.66 ± 0.19 m w.e. a⁻¹ [18]. The above comparisons indicate that the response of the eastern Tanggula glaciers to climate change is highly consistent with that of the glaciers in adjacent areas within the Tibetan Plateau, but is quite different from the response of the glaciers in southeast Tibet. This seems to imply that there may be a dividing line with regard to glacier mass change between southeastern Tibet (e.g., the Yigong Zangbo Basin) and the eastern Tanggula Mountains.

Moreover, with regard to the mass changes of the different glacier types (i.e., land-terminating and lake-terminating), the change mode of the eastern Tanggula region (i.e., the rate of mass loss of the land-terminating glaciers is higher than that of the lake-terminating glaciers during all study periods) is completely contrary to the general understanding that lake-terminating glaciers, in general, show higher rates of mass loss than land-terminating glaciers (e.g., –0.55 ± 0.12 m w.e. a⁻¹ vs. –0.37 ± 0.12 m w.e. a⁻¹ in the Himalayas [20] and –0.89 ± 0.36 m w.e. a⁻¹ vs. –0.50 ± 0.32 m w.e. a⁻¹ in the eastern Nyainqentanglha Mountains [18]). The main reasons for such a pattern are that the land-terminating glaciers are relatively few in number and small in size in this study region, and these glaciers are almost all located in the relatively lower altitudes (Figure 7a). Furthermore, with regard to the eastern Tanggula lake-terminating glaciers, for one thing, the accumulation zones exhibited a certain degree of thickening in the first decade of the 21st century, which compensates for the mass loss in the ablation zone and further leads to relatively little mass wastage for the whole region. Furthermore, the termini of some of the lake-terminating glaciers have gradually separated from the proglacial lakes, and they have further evolved into land-terminating glaciers (e.g., the ZGMN Glacier, East Beijia Glacier). Furthermore, the Zuxuehui Glacier in the south of this region is almost unaffected by the proglacial lake that is connected to it. These findings indicate that, in terms of proglacial lakes, there is a trend toward gradual weakening (or even disappearance) of their influence on glacier mass loss in this region.

6. Conclusions

Aiming at the issue of the lack of glacier change monitoring in the eastern Tanggula Mountains, which is the only region where lake-terminating glaciers are found within the Tibetan Plateau, and based on multiple-source optical satellite data (i.e., Landsat 5/8 and Sentinel-2), we generated updated glacier boundaries and estimated the decadal area changes for the first twenty years of the 21st century. Meanwhile, two TanDEM-X bistatic mode data sets, together with the SRTM DEM, were used to estimate the changes in glacier thickness and mass balance in the same study periods as the area change monitoring. The final results showed that the glaciers of the eastern Tanggula Mountains have experienced an accelerated rate of area shrinkage, with the rate of area loss changing from −0.34 ± 0.83 km² a⁻¹ to −0.93 ± 0.81 km² a⁻¹ for 2000–2013 and 2013–2020, respectively. Simultaneously, the glacier thickness of this region has shown an accelerated thinning trend, with the rate of surface lowering changing from −0.19 ± 0.05 m a⁻¹ and −0.53 ± 0.08 m a⁻¹ for 2000–2012 and 2012–2020, respectively, corresponding to a 1.8 times acceleration. Moreover, we found that the patterns of both the area change and the thickness change (or mass change) are closely related to the glacier types. Specifically, for the lake-terminating glaciers, the rate of area loss has been generally higher than that for the land-terminating glaciers, and the corresponding acceleration is also more significant than for the land-terminating glaciers. Similarly, with regard to glacier thickness variation, the acceleration of surface thinning for the lake-terminating glaciers is higher than for the land-terminating glaciers. It is worth noting, however, that the surface thinning rates of the lake-terminating glaciers for all of the study periods were always less significant than the surface thinning rates of the land-terminating glaciers. Consequently, the final glacier mass balance was −0.16 ± 0.04 m w.e. a⁻¹ and −0.45 ± 0.07 m w.e. a⁻¹ for the two sub-periods of 2000–2012 and 2012–2020, respectively, and the results for the whole study period (during 2000–2020) was −0.23 ± 0.04 m w.e. a⁻¹, which is equivalent to a rate of mass loss of −0.04 Gt a⁻¹. A further analysis based on our filed investigation and proglacial lakes monitoring showed that the local topography plays a vital role in the evolution of glacial lakes. We speculate that all of these glaciers will turn into land-terminating glaciers in the near future. Furthermore, the present status of the glacier changes in this region can be attributed to the long-term increase in air temperature. Our findings will provide an overview of the overall changes, and will improve the understanding of the response of glaciers to climate change in this region.