Recent 50-Year Glacier Mass Balance Changes over the Yellow River Source Region, Determined by Remote Sensing

Abstract

The A’nyêmaqên Mountains have the largest concentration of glaciers in the Yellow River basin, which play a crucial role in regulating the runoff regime of the Yellow River. Thus, the quantification of glacier mass balance and its effects on river runoff is greatly required. However, current studies mainly focus on mass changes since 2000. Here, we report for the first time region-wide glacier elevation and mass changes, which were derived from digital elevation models (DEMs) produced from historical topographic maps (TOPO), SRTM retrievals, and ASTER L1A stereo imagery spanning the past 50 years. The results indicated a negative mass balance (−0.24 ± 0.05 m w.e. a⁻¹) of all glaciers for the 1966–2018 timespan. The mass loss rapidly accelerated from −0.16 ± 0.09 m w.e. a⁻¹ in 1966–2000 to −0.36 ± 0.06 m w.e. a⁻¹ during the period from 2000–2018. The rise in mass loss rate from 2000 onwards was mainly associated with the rapidly increased summer warming.

1. Introduction

The Yellow River is the sixth longest river on the Earth and flows from its origin in the Bayan Har Mountains to the Bohai Sea, passing through nine provinces in China. It provides the water supply for the lives of more than 100 million people and also for agricultural production. The source region of the Yellow River, which is located in the northeastern Tibetan Plateau, is regarded as the water tower of the entire Yellow River basin due to its large contribution (~35%) to the total annual runoff [1]. In particular, glaciers in this area play a key role in regulating tributaries and runoff [2]. Against the background of climatic warming and increased human activities, the hydrological and ecological environment in this region has become more and more vulnerable, and as a result, a lot of environmental problems [3,4,5,6], such as glacier ablation, ice avalanche hazards, permafrost degradation, and shrinkage of wetlands and lakes have arisen. To address these, the estimation of glacier variations in the Yellow River source zone is highly needed, and also essential for the implementation of the major national strategies of China for ecological protection and high-quality development in the Yellow River basin.

The A’nyêmaqên Mountains are the largest glaciated area in the Yellow River source region, and glacial changes exert a significant impact on the local water resources and animal husbandries [4]. Since 2004, three ice avalanche events have happened in the Maqên Gangri Glacier, which resulted in debris flow, the formation of a glacier lake, and then outburst floods. Despite no human and livestock deaths, large amounts of grasslands were seriously damaged. Several previous studies [2,7] have investigated glacier area variability in this region based on satellite images and have examined their response to climate changes. Liu et al. [2] and Yang et al. [8] reported that the glacier area shrank rapidly between the Little Ice Age (LIA) and the year 2000. Previous studies [7,9] have indicated that the glacier shrinkage rate declined from 2000 onwards. Limited attempts have also been made to evaluate glacier elevation and mass balance in this region. Jiang et al. [4] reported the mass loss of three glaciers from 2000 to 2013 by means of SRTM DEM and TanDEM-X. Recently, Li et al. [10] reconstructed the annual mass balance of two glaciers in this area based on the gradient between snowline altitude and mass balance from 1990 to 2020. Despite the significant amount of large-scale, satellite-based glacier mass balance assessments in the Tibetan Plateau [11], High Mountain Asian (HMA) [12] and even the whole Earth [13] since 2000, the glaciers in the A’nyêmaqên Mountains are usually not included. In brief, we found it to be interesting that very little is known about glacier mass balance in this region due largely to the unavailability of long-term altimetry data.

In this study, the major objective was to utilize multi-source DEM data (Topographic maps, SRTM DEM, and ASTER L1A) to investigate glacier mass balance changes in the Yellow River source region during the past 50 years, and then to provide a reference for the estimation of water resources and natural hazards in the Yellow River basin.

2. Study Area

The A’nyêmaqên Mountains are situated between latitude 34°20′–35°N and longitude 99°10′–101°30′E of the northeastern Tibetan Plateau and belong to the eastern branch of the Kunlun Mountains (Figure 1). The mountains stretch in a northwest to southeast direction, and their average length and width are 120 km and 80 km, respectively. They consist of 13 peaks, of which the highest peak is Maqên Gangri, with an elevation of 6282 m a.s.l. Affected by the Asian monsoon and westerlies, the climate is cold, moist, and windy, and precipitation largely happens in the summer. The lowest temperature of the year is in January, and the temperature at the mountain peaks can reach below −50 °C. Based on the first Chinese Glacier Inventory, there are 61 glaciers in total, covering 124 km² in the area [2]. Among them, the Halong Glacier, Weigeledangxiong Glacier, and Yahelong Glacier are three glaciers covering more than 10 km² in area. Valley-type glaciers are situated on the northeastern slopes, whereas hanging and cirque glaciers are mainly distributed over the southwestern slopes. The snowline altitude of the glaciers on the eastern slopes ranges from 4900 to 5190 m. The glaciers on the western slopes have an average snowline altitude of 5160 m.

3. Data and Methods

3.1. Weather Station Data

In order to investigate the climatic forcing on glacier variations in this study area, we analyzed variations in the monthly air temperature and precipitation records from four meteorological stations from 1966 to 2018: Xinghai (35°34′N, 99°58′E, 3324.3 m a.s.l.), Maduo (34°55′N, 98°13′E, 4273.3 m a.s.l.), Guoluo (34°28′N, 100°15′E, 3720.0 m a.s.l.), and Dari (33°45′N, 99°39′E, 3968.5 m a.s.l.), which are close to the glacier fields (see Figure 1). The meteorological station data used in this study were obtained from the daily meteorological dataset of basic meteorological elements from the China National Surface Weather Station (V3.0), which were downloaded from the China Meteorological Data Service Center (CMDC, https://data.cma.cn/). After controlling for data quality, the quality and integrity of the data for each meteorological parameter was significantly improved and compared with similar ground data products released in the past, and the integrity of the data was close to 100%. Missing meteorological data in December 2001 were provided using the records from the original version of the dataset from the China National Surface Weather Station (V2.0).

3.2. Mapping Glacier Boundaries

Based on the ArcGIS platform, the glacier boundaries were delineated by visual interpretation using topographic maps and Landsat images. In order to ensure its accuracy, we adjusted the topographic map and Landsat images to the same resolution (30 m), and two Chinese Glacier inventories (CGIs) were used to verify the mapping of the glacier boundaries. We used three 1:50,000 scale topographic maps (TOPO), which were produced from aerial stereo pairs in November 1966 using the Chinese Military Geodetic Service (CMGS). Their horizontal and vertical references were the Beijing Geodetic Coordinate System 1954 (BJ54) and the Yellow Sea 1956 datum (the mean sea level at the Qingdao Tidal Observatory in 1956), respectively. We re-projected the coordinate system into the World Geodetic System 1984 (WGS1984) and Earth Gravity Model 1996 (EGM96). When re-projecting, a seven-parameter transformation method was used due to the high accuracy of this method (<0.5 m; [14,15]).

The Landsat TM sensor consisted of 7 spectral bands and all bands except band 6 (thermal infrared radiation) had the horizontal resolution of 30 m. The ETM+ images added a new 15 m spatial resolution panchromatic band, and improved the spatial resolution of the thermal infrared band to 60 m. Compared to the ETM+ sensor, the OLI sensor launched onboard Landsat 8 in 2013 enhanced the radiometric signal-to-noise ratio, resolution, and scanning systems, which was composed of 9 spectral bands with a spatial resolution of 30 m for bands 1–7 and band 9, and 15 m for the panchromatic band 8 [4]. Here, we used one TM and one ETM image in 2000, and three OLI images obtained in 2015, 2016, and 2018, with minimal cloud and snow (

Table 1. List of data utilized for the determination of glacier changes over the A’nyêmaqên Mountains.

Data	Date	Pixel Size/Scale	Purpose	Resource
Topographic maps	1966	1:50,000	DEM/base image	CMGS
SRTMGL1	February 2000	30 m	Reference DEM	USGS
SRTM	February 2000	90 m	DEM	CGIAR-CSI
ASTER L1A	October 2018	30 m	DEM	EARTH DATA
HMA DEM	July 2015	8 m	Reference DEM	NSIDC
Landsat TM	September 2000	30 m	Base image	GLCF
Landsat ETM+	August 2000	30 m	Base image	GLCF
Landsat 8/OLI	July 2015	30 m	Base image	USGS
Landsat 8/OLI	June 2016	30 m	Base image	USGS
Landsat 8/OLI	June 2018	30 m	Base image	USGS

). Debris-covered and debris-free areas were distinguished by means of the false color composite of bands from these Landsat images. Before extraction of the glacier outlines, all Landsat images were co-registered to the topographic maps using >20 ground control points (GCPs), which were distributed on stable topography of the remote sensing images and topographical maps, such as bedrock, rock outcrops, and moraine structures.

3.3. Surface Elevation

Fifty-year glacier surface elevation and mass variations were determined by DEMs from three different sources including TOPO DEM (1966) generated from the topographic maps, SRTM DEM (2000), and the ASTER L1A stereo imagery (2018). Generally, they were produced by different methods and have different uncertainties.

The 20 m contours and elevation points from the topographical maps were digitized to generate a triangulated irregular network (TIN) and the TIN was then used to establish the DEM (TOPO DEM) at 30 m resolution by linear interpolation. According to the General Administration of Quality Supervision Inspection and Quarantine (2008) [16], its maximum error was present at the slope of >25°, and the nominal accuracy of these regions on the topographic map was ±14 m in verticality.

The SRTM DEMs were constructed by using C-band and X-band Spaceborne Imaging Radar operating in interferometric mode aboard the space shuttle Endeavour from 11 to 22 February 2000, which acquired radar data between 60°N and 56°S covering more than 80% of the global land surface. The near-global DEMs were temporally consistent with vertical uncertainty of ±10 m root mean square error (RMSE), and have been widely applied for geodetic glacier mass balance estimation, e.g., [17,18,19]. Due to the narrower swath of X-band radar, its spatial coverage is rather incomplete. Thus, we estimated glacier mass changes by means of C-band DEMs, which refers to the WGS84 in the horizon, and the EGM96 geoid in verticality. The SRTM Version 3.0 Global 1 arc second dataset (SRTMGL1) was released by the United States Geological Survey (USGS) in 2015, and data voids were filled with Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model 2 (ASTER GDEM2), combined with Global Multi-Resolution Terrain Elevation Data (GMTED) 2010, and National Elevation Dataset (NED) from the USGS [20]. The SRTM GL1 covering the A’nyêmaqên area was selected for glacier DEM co-registrations. We also used non-void-filled SRTM DEM to calculate the glacier elevation changes.

The Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) [21] is an advanced optical sensor on the TERRA satellite, which contains 14 spectral channels of Visible and Near Infrared (VNIR), Shortwave Infrared (SWIR), and Thermal Infrared (TIR). The AST_L1A data are a relatively early product of ASTER, and retain the original spatial resolution, with a temporal coverage from 2000 to the present. We selected the AST_L1A stereo imagery acquired in October 2018 with good quality. Taking the SRTM dataset as a reference, the NASA Ames Stereo Pipeline (ASP) tool was used to generate DEMs from the ASTER image based on the stereo pair composed of band 3N (nadir-viewing) and band 3B (backward-viewing). In the process, most errors caused by the cloud were effectively eliminated. The elevation values of the DEMs were provided in meters with respect to the EGM96 geoid, and their projection and datum were Universal Transverse Mercator (UTM) and WGS84, respectively. The absolute accuracy estimates were high in horizon and verticality. In addition, we used the stable off-glacier regions of HMA DEM to verify the altitude accuracy of ASTER DEM. The resulting root mean square error of elevation for the stable terrain area was ~7.4 m (See the Supplementary File).

3.4. DEM Co-Registration and Corrections of Radar Penetration and Terrain Curvature

All DEMs were co-registered to C-band SRTM DEM based on the method by Nuth et al. [22], which applies the elevation difference as a function of aspect and slope over a stable off-glacier region to correct the horizontal biases of the DEMs. Before co-registration, all DEM data were re-projected to the same reference system (map projection: UTM Zone 47 North projection and WGS84; vertical datum: EGM96 geoid), and the pixels of DEM difference falling outside of the 5–95% quantile range were omitted to eliminate outliers over the DEM edges and near data gaps [23]. The corrected results are shown in

Table 2. The offsets in X, Y, and Z directions in the DEM dataset and the uncertainty in DEMs before and after co-registration.

	Offsets in X, Y, and Z Directions			Before Co-Registration		After Co-Registration
	X (m)	Y (m)	Z (m)	$\overline{\mathbf{dh}}$	SD (m)	$\overline{\mathbf{dh}}$	SD (m)
TOPO-DEM	−0.5	5.6	0.6	0.3	11.1	0.3	11.1
SRTM DEM 2000	−4.5	−23	0	−0.2	20.3	0.0	12.8
ASTER DEM 2018	−1.5	−2.8	0.9	0.4	10.8	0.6	10.6

. The DEMs with low resolution tend to underestimate the height of peaks or ridges with high terrain curvature. In order to correct the vertical bias due to the different spatial resolutions of the DEMs, following Gardelle et al. [24] we applied the relationship between the elevation difference and the maximum curvature over stable off-glacier regions.

When estimating mass balance changes based on C-band SRTM DEM, the impact of the penetration should be taken into account due to the penetration of C-band radar into snow and ice [25,26]. Previous studies usually used the simultaneously acquired X-band DEM with the C-band DEM to diminish the influence of penetration. However, there was no available X-band DEM for this study area because of its narrower data coverage. Here, we used the significant linear positive correlation between elevation and the penetration depth difference (PDD) of the C/X-band radar signals on glaciers over the whole HMA region as revealed by Li et al. [27] to estimate the penetration depth of the C-band. In this method, the whole HMA region was divided into several 1° × 1° grids. The A’nyêmaqên Mountains were exactly in the 1° × 1° grid (N34E099). The penetration correction on glaciers was performed for each pixel using the linear correlation between elevation and PDD. The resulting average penetration depth in the study region was about 4.5 m (the accuracy verification of the method for penetration depth correction is in the Supplementary File).

3.5. Glacier Elevation Change and Mass Balance, and Their Uncertainty

After differencing the DEMs to yield surface elevation changes, it was necessary to exclude outlier values to reduce the uncertainty in the elevation change calculation. We first excluded DEM pixels with absolute elevation changes of more than 100 m, possibly resulting from stereo-matching errors [28]. Then, a slope threshold of 45° was selected to serve as the lower bound for excluding any rather steep parts in the accumulation regions. Furthermore, we removed the pixels with absolute elevation variability of more than 3 standard deviations in each 100 m elevation band. Following Holzer et al. [29], pixels outside of the 31.7–68.3% quantiles of elevation changes in each glacier accumulation zone were excluded as outliers. Lastly, we filled gaps in DEM difference data caused by outlier filtering using the average value of elevation changes from the appropriate 100 m elevation [30].

The elevation difference in glacier surface was obtained by comparing DEMs over the different time phases and then it was converted into water equivalent to estimate the glacier mass balance. The combined ice and snow density should be considered. We used a constant density of 850 ± 60 kg·m⁻³ to convert the glacier elevation changes to glacier mass balance [31].

Uncertainty of each co-registered DEM was quantified using the quadratic sum of the standard deviation of the DEM difference in relation to SRTM DEM over the stable of the median elevation changes on and off glaciers [32], which can be calculated by:

$$U_{\Delta h}=\frac{S{D}_h}{\sqrt{n_z}}$$

(1)

where n_z is the effective number of observations and is estimated by the total number of elevation change measurements (n_b), DEM resolution (r), and distance of spatial autocorrelation (d) [25]. The formula is following:

$$n_z=\frac{n_b\ast r^2 }{\pi \ast d^2}$$

(2)

The uncertainty in glacier area was calculated based on standard principle of error propagation developed by Braun et al. [33]:

$$U_a=\frac{P/A}{P/A_{\left(Paul et al.\right)}}\ast 0.03$$

(3)

where U_a represents the error of the glacier area and P/A is the perimeter–area ratio; P/A was a constant value of 5.03 km⁻¹ as suggested by the reference [34]. Furthermore, when calculating the final uncertainty of glacier volumetric mass balance, we also accounted for the ice density uncertainty (U_d) (± 60 kg·m⁻³) and the penetration depth error (p_e) (the uncertainty of the SRTM radar penetration, ±1.4 m [35]), as suggested by Kääb et al. [36]. The error of gap filling by extrapolation (U_e) was equivalent to the maximum of the SDs of glacier elevation change in each 100 m elevation bin [29,34]. The final mass balance uncertainty (δ_m) was calculated, also considering systematic and random uncertainty in glacier elevation changes (U_h), the extrapolation uncertainty (U_e), the area uncertainty (U_a), and the ice density error (U_d):

$${\delta }_m=\Delta M\ast \sqrt{{\left(\frac{U_{\Delta h}}{\Delta h}\right)}^2+{\left(\frac{U_e}{\sqrt{N_{eff}}}\right)}^2/\Delta h^2+p_e^2+{\left(U_a/a\right)}^2+{\left(U_d/{\rho }_i\right)}^2}$$

(4)

where ΔM denotes the mass balance estimate, Δh is the estimated elevation change, a is the glacier area, and ρ_i is the glacier ice density uncertainty.

4. Results

4.1. Climatic Changes

The average annual air temperature at the four stations from 1966 to 2018 was 0.61 °C: 9.18 °C for the summer average and −11.33 °C for the winter average. The time series of air temperature anomalies averaged from the four stations showed significantly increasing trends in each season, with warming rates of 0.23 °C decade⁻¹ (p < 0.05) in the spring, 0.34 °C decade⁻¹ (p < 0.05) in the summer, 0.40 °C decade⁻¹ (p < 0.05) in the autumn, and 0.49 °C decade⁻¹ (p < 0.05) in the winter, respectively, for the 1966–2018 period (Figure 2). After 2000, the overall air temperature showed a warming trend, and the spring and summer warming rate significantly increased, from 0.10 °C decade⁻¹ (p > 0.05) and 0.19 °C decade⁻¹ (p < 0.05) before 2000 to 0.63 °C decade⁻¹ (p < 0.05) and 0.55 °C decade⁻¹ (p > 0.05), respectively. Obviously, the summer warming rate during the 2000–2018 period was threefold higher than that during the 1966–2000 period (Figure 2). Similarly, compared with the trend of annual mean air temperature before 2000 (0.21 °C decade⁻¹, p < 0.05), an increased trend appeared in the average annual temperature over the 2000–2018 period (0.49 °C decade⁻¹, p > 0.05).

The annual precipitation in this region was 445.15 mm from 1966 to 2018, and for the same period, the annual precipitation showed an insignificantly increasing trend (0.98 mm decade⁻¹, p > 0.05). The precipitation in all four seasons showed increasing trends from 1966 to 2018, with rates of 1.47 mm decade⁻¹ (p > 0.05) in the spring, 1.98 mm decade⁻¹ (p < 0.05) in the summer, 0.22 mm decade⁻¹ (p > 0.05) in the autumn, and 0.25 mm decade⁻¹ (p > 0.05) in the winter, respectively. As shown in Figure 3, the annual precipitation from 2000 to 2018 showed a significantly increasing trend (5.03 mm decade⁻¹, p < 0.05), while before 2000 the annual precipitation exhibited a slowly decreasing trend (−0.14 mm decade⁻¹, p > 0.05). The upward trends in the spring and summer precipitation over the 2000–2018 period (7.71 mm decade⁻¹, p < 0.05) in the spring and 10.05 mm decade⁻¹, (p < 0.05) in the summer) were higher than for the 1966–2000 period (0.23 mm decade⁻¹ in the spring and 0.28 mm decade⁻¹ in the summer). The autumn precipitation exhibited a slowly decreasing trend during the 1966–2000 period (−1.64 mm decade⁻¹, p < 0.05), but an insignificantly increasing trend between 2000 and 2018 (1.71 mm decade⁻¹, p > 0.05). Although the winter precipitation increased for the 1966–2000 and 2000–2018 periods, the increasing rate during the former timespans (0.67 mm decade⁻¹, p > 0.05) was slightly lower than the latter (0.59 mm decade⁻¹, p > 0.05).

4.2. Glacier Area Changes

The area covered by glaciers in the A’nyêmaqên was 124.01 km² in 1966. Dominated by small glaciers, the glaciers with an area of less than 1 km² accounted for 59% of the total number of glaciers, but their area was only 11.3% of the total area. There are 3 glaciers covering more than 10 km², accounting for 47% of the total area. These glaciers shrank to 105.63 ± 3.04 km² in 2000, and to 96.57 ± 2.78 km² in 2018. For the 1966–2000 and 2000–2018 periods, the glacier area shrinkage rate was 0.53 km² a⁻¹ and 0.48 km² a⁻¹, while the glacier area decreased by 0.42 % a⁻¹ and 0.45% a⁻¹, respectively (

Table 3. Changes in glacier areas over the A’nyêmaqên Mountains.

Period	1966–2000	2000–2018	1966–2018
Δarea (km2 a−1)	−0.53	−0.48	−0.52
Δarea (% a−1)	−0.42	−0.45	−0.42

). The results indicated a slight reduction in mean glacier shrinkage rates after 2000 in relation to the period 1966–2000, but the shrinkage percentage increased.

4.3. Surface Elevation and Geodetic Mass Balance

A remarkable surface thinning was observed for most of the investigated glaciers during the 1966–2018 period. Over the 1966–2000 period, we calculated a mean elevation change of −0.19 ± 0.12 m a⁻¹, while for the 2000–2018 period, the mean elevation change was −0.42 ± 0.07 m a⁻¹. The mean geodetic mass balance for the entire A’nyêmaqên glaciers was estimated to be −0.16 ± 0.09 m w.e. a⁻¹ during the 1966–2000 period, and −0.36 ± 0.06 m w.e. a⁻¹ during the 2000–2018 period, respectively. This indicates that the mass loss accelerated during the 2000–2018 period, compared with that during the 1966–2000 period.

As shown in Figure 4, an obvious decline in surface elevation was dominant at the ablation zone, whereas a slight increase in elevation occurred in part of the accumulation zone (>5800 m) of the glaciers for the 1966–2000 and 2000–2018 periods. Over the 1966–2000 period, the gains in glacier surface elevation occurred at the terminus of glaciers A, D, and E (Figure 8A,D,E), whereas the surface lowered upstream in these glaciers. For the eastern branch of glacier C, the surface distinctly rose in the middle of this branch, whereas the surface elevation significantly reduced over the upstream. These findings suggest that surging events happened at glaciers A, C, D, and E during the 1966–2000 period. Over the 2000–2018 period, a very large elevation decrease occurred in the middle of glacier C, but the terminus of glacier C showed a strong accumulation. For the same period, the elevation increase was shown upstream of glaciers A and E, and elevation reduction occurred downstream of the two glaciers.

The mean elevation thinning rate of each 100 m altitude bin during the three periods (Figure 5) indicated that the surface thinning was observed mainly in the low altitudes and surface thickening occurred in the high altitudes. The thinning rate of the glaciers gradually decreased with increased altitude. Compared with the trend for the 1966–2000 period, the glacier elevation changes were more negative for each elevation bin of <5400 m elevation over the 2000–2018 period. For both periods, the gains in glacier elevation were only present above 5800 m elevation.

Figure 6 shows the mass losses of each glacier aspect during the three different timespans (1966–2000, 2000–2018, and 1966–2018). Strong mass losses were observed in the south, southwest, and west aspects for the three periods, with the highest loss in the southwest aspect. During the past 50 years, the mean mass loss rates in the south, southwest, and west aspects were estimated to be −0.48 ± 0.05 m w.e. a⁻¹, −0.52 ± 0.05 m w.e. a⁻¹, −0.38 ± 0.05 m w.e. a⁻¹, respectively. The mass loss rate in all aspects including the northeast and western aspects during the 1966–2000 period was smaller than those during the 2000–2018 period. The larger mass loss in the southwestern aspects may be related to less precipitation on the eastern slopes than on the western slopes, mainly caused by the East Asian monsoon [9], whereas the mass loss on the western slopes after 2000 may be due to rapid warming. Furthermore, the glaciers facing south generally receive more solar radiation than those facing north. The slopes of the mountains in the south, west, and southwest are much steeper than those in the north, northwest, and east.

The debris-covered areas have been concentrated in the A, B, and C glaciers of the A’nyêmaqên Mountains (Figure 7). In this study, the mass balance changes of debris-covered glaciers (A, B, and C glaciers) and clean glaciers from 1966 to 2018 were calculated, respectively, and the results showed that the average ablation rate of debris-covered glaciers (−0.15 ± 0.05 m w.e. a⁻¹) was much slower than that of the clean glaciers (−0.33 ± 0.05 m w.e. a⁻¹) (

Table 4. Glacier mass balance changes in debris-covered and clean glaciers.

	Debris-Covered Glaciers (m w.e. a−1)	Clean Glaciers (m w.e. a−1)
1966–2000	–0.05 ± 0.09	–0.26 ± 0.09
2000–2018	–0.32 ± 0.06	–0.39 ± 0.06
1966–2018	–0.15 ± 0.05	–0.33 ± 0.05

). Among them, the average ablation rate of the debris-covered glaciers in the 1966–2000 and 2000–2018 periods were −0.05 ± 0.09 m w.e. a⁻¹ and −0.32 ± 0.06 m w.e. a⁻¹, respectively, which was also slower than that of the clean glaciers (−0.26 ± 0.09 m w.e. a⁻¹ and −0.39 ± 0.06 m w.e. a⁻¹).

5. Discussion

5.1. Glacier Changes and Glacier Dynamics

Liu et al. [2] and Yang et al. [8] reported an accelerating shrinkage of glacier extent in the A’nyêmaqên glaciers since the LIA. In particular, the area of the A’nyêmaqên glaciers shrunk by about 17% during 1966–2000, consistent with our glacier area change estimation for the same period. Another study showed that the glacier area was reduced by 5.87% in this region from 1992 to 2001 [9]. However, after 2000 the glaciers have been shrinking at a significantly slower rate [9]. The total glacier area loss over the A’nyêmaqên glaciers was estimated to be 3.37 % from 2001 to 2010 by Wang et al. [9]. These results were broadly consistent with our glacier area change estimation, i.e., about 14.82% shrinkage of the A’nyêmaqên glaciers from 1966 to 2000, and a slight reduction rate after 2000. The small differences between the different studies may be related to the different calculation periods and the quality of the satellite images.

For changes in the surface mass balance of the A’nyêmaqên glaciers, previous studies have suggested mass loss of the glaciers in this region [4,10]. Based on SRTM and Tandem-X, Jiang et al. [4] reported the elevation changes of the A’nyêmaqên glaciers during the 2000–2013 period in more detail and found that the elevation of the Yehelong Glacier decreased most significantly by about −0.93 ± 0.26 m a⁻¹. This was larger than what we calculated for this glacier between 2000 and 2018 (−0.69 ± 0.07 m a⁻¹), which may be associated with no correction of SRTM C band penetration by Jiang et al. [4]. Our result (−0.42 ± 0.07 m a⁻¹ in the A’nyêmaqên glaciers over the 2000–2018 period) agrees comparatively well with the previously published mean change results in this region over the 2000–2019 period from Hugonnet et al. [37].

Glacier surging, ice avalanches, and debris cover could also affect glacier changes.

Glacier surging is caused by internal ice instability, which is related to the glacier entropy balance, length, area, and slope [38,39,40]. Surge-type glaciers are generally characterized by terminus advance, and obvious thickening and thinning of the ablation and accumulation zones, respectively [41,42]. From Figure 4, it is clearly found that there were accumulation trends in the terminus of A, D, and E glaciers (see Figure 8) for the 1966–2000 period, suggesting the occurrence of surge events, which was in accordance with the previously identified results by [4,43]. Surge events also happened over the B glacier during the 1966–1981 period as observed by Wang et al. [44]. Only part of the debris-covered tongue of the C glacier appeared to surge during the 2000–2018 period [4], and the glacier tongue with high thickening (Figure 4) also hinted at a past surge event. To examine the impact of surging on the A’nyêmaqên glacier changes, we calculated the mass balance of the surging glaciers and non-surging glaciers. With the exception of the five surging glaciers, the mass balance of the A’nyêmaqên glaciers from 1966 to 2018 was −0.39 ± 0.05 m w.e. a⁻¹, which was faster than the overall mass balance of these glaciers containing surging glaciers (−0.28 ± 0.05 m w.e. a⁻¹). However, although surging glaciers presented numerical mass loss (−0.08 ± 0.09 m w.e. a⁻¹ during the 1966–2000 period; −0.39 ± 0.06 m w.e. a⁻¹ during the 2000–2018 period), the glacier surging did not affect the ablation trends of the whole glaciers before and after 2000.

Ice avalanches are one of the predominant forms of glacier ablation [43]. According to their report [45], three large-scale ice avalanches have occurred in the Maqên Gangri Glacier since 2004, affected by climate warming. Ice avalanche activity has resulted in the significant mass loss (Figure 4) and terminus retreat (Figure 8) in the Maqên Gangri Glacier. The elevation of this glacier decreased by −0.40 ± 0.07 m a⁻¹ during the 2000–2018 period, higher than the elevation decline (−0.31 ± 0.12 m a⁻¹) during the 1966–2000 period. A positive impact of ice avalanches on the reduction in the overall glacier elevation was exhibited. Additionally, influenced by rapid melting of the glacier, incidental debris flow occurred easily in this area. As Figure 9a–c shows, debris flow suddenly happened on 8th August 2022, largely induced by the melting of the Maqên Gangri Glacier. With the further warming climate, the snowline will continue to rise, and there is still the possibility of another ice avalanche event in the region.

Based on the debris-covered glaciers, it is generally believed that a thick debris cover will act as a protective layer for glaciers due to its insulation [46,47], and the thicker debris cover has a negative influence on ablation. According to the mass balance changes of the debris-covered glaciers and clean glaciers from 1966 to 2018 calculated in this study, it was found that the average ablation rate of debris-covered glaciers was much slower than that of clean glaciers (Table 4). Especially for the 1966–2000 period, the rate of debris-covered glacier melting was only −0.05 ± 0.09 m w.e. a⁻¹.

5.2. Climate Forcing

The high sensitivity of the glacier mass balance to air temperature is suitable for the whole HMA [48]. Air temperature changes stood out as a stronger control on glacier thinning in our study region than the other meteorological parameters [9]. Observations by Liu [49] showed that the climate on the A’nyêmaqên Mountains has changed to wet and warm during the past 60 years. Our study also confirmed that a warming trend occurred in this region from 1966 to 2018 (Figure 2). Bhattacharya et al. [12] found a high and significant correlation (r² = 0.60, p < 0.01) between summer temperatures measured at the weather stations and glacier mass balance. The increasing trend in the summer air temperature was significant for the 2000–2018 period, which was more than three times as high as the summer warming rate during the 1966–2000 period. The increase in the summer warming rate after 2000 accelerated the melting of the glaciers. The winter warming over the 2000–2018 period may have shortened the accumulation time of the glaciers.

In addition, the annual precipitation showed an insignificantly increasing trend in this region over the 1966–2018 period. From Figure 3, we observed a weakly upward trend in the winter precipitation over the 2000–2018 period, slightly higher than the 1966–2000 period, and even the increase in total annual precipitation in this region from 2000 to 2018 was larger than that during the 1966–2000 period. This could have resulted in slightly more snow accumulation on the glaciers over the 2000–2018 period than the 1966–2000 period. However, the increased recharge of the glaciers by precipitation cannot compensate for more intensification of melting due to the rapid warming since 2000. As a consequence, each 100 m elevation bin of the >6000 m elevation showed only slight gains in surface elevation for the 2000–2018 period, but more negative surface elevation changes occurred in all elevation bins of <5400 m elevation during the 2000–2018 period than the 1966–2000 period (Figure 5). The acceleration of summer warming ablated the glaciers much more than the accumulation from the insignificant increase in precipitation. Thus, we reasonably conclude that the accelerated summer warming and the insignificant increase in precipitation seemed to be the main reasons for the decelerated mass loss of the A’nyêmaqên glaciers after 2000.

6. Conclusions

This study estimated the glacier surface elevation and mass balance changes for all glaciers in the A’nyêmaqên Mountains over the 1966–2018 period using geodetic measurements. A negative mass balance has been obtained since 1966, and the glacier mass loss rate showed an obviously increasing trend after 2000. The larger mass loss since 2000 may be associated with the acceleration of summer warming, which ablation much more than the accumulation caused by the insignificant increase in precipitation. To better estimate the response of the A’nyêmaqên glacier to future climate scenarios, it is essential to use the energy balance models for the accurate quantification of glacier responses to changes in the air temperature and precipitation. More broadly, this is also of vital importance to assess the changes in the A’nyêmaqên glacier for the future water adjustment of the Yellow River source.

In addition to climatic forcing, morphologic factors, and glacier dynamics also contribute to glacier mass balance changes. The glacier mass balance decreased mainly in the south, southwest, and west directions during the 1966–2018 period, and the debris-covered glaciers (−0.15 ± 0.05 m w.e. a⁻¹) underwent thinning at a lower rate than the clean glaciers (−0.33 ± 0.05 m w.e. a⁻¹). Five surging glaciers (6.4% of the total) were found in this study, but glacier surging did not affect the mass balance change trends of the whole glaciers during our study period. We also observed ice avalanche activities occurring during the study period, which exerted an intensive impact on glacier mass loss, and resulted in hazards.