Characterizing 4 decades of accelerated glacial mass loss in the west Nyainqentanglha Range of the Tibetan Plateau

Glacier retreat is altering the water regime of the Tibetan Plateau (TP) as the region's climate changes, but there remain substantial gaps in our knowledge of recent glacier loss in this region due to the difficulty of making direct high-mountain observations, and this limits our ability to predict the future of this important water resource. Here, we assessed 44 years of glacier area and volume changes in the major west Nyainqentanglha Range (WNR) that supplies meltwater to the densely populated Lhasa River basin and Nam Co, the second largest endorheic lake on the TP. Between the two periods 1976–2000 and 2000–2020, we found that the glacier areal retreat rate more than doubled (from $-\mathrm{0.54}\pm \mathrm{0.21}$ % a⁻¹ to $-\mathrm{1.17}\pm \mathrm{0.30}$ % a⁻¹), and surface lowering also accelerated (from $-\mathrm{0.26}\pm \mathrm{0.09}$ to $-\mathrm{0.37}\pm \mathrm{0.15}$ m w.e a⁻¹) with particularly intense melting after 2014. This acceleration is similar in both timing and magnitude to that observed for Himalayan glaciers farther south. Besides, the areal retreat rate and mass loss rate of most glaciers in the WNR were not synchronized. To understand the sensitivity of WNR glaciers to climate forcing, we examined the effects of topography, debris cover and the presence of proglacial lakes on our observed changes. We found consistently faster areal retreat but slower thinning rates on steeper slopes and an inconsistent relationship with aspect. We concluded that our observed spatial and temporal patterns of glacier change were dominated by observed local variations in temperature and precipitation, the melt-reducing role of supraglacial debris, and the increasing influence of ice-marginal lakes on glacier ablation.

1 Introduction

The Tibetan Plateau (TP), known as the “Water Tower of Asia”, is the source of several of Asia's major rivers (Bolch et al., 2010), and glacial melt on the TP plays an important role in water supply for downstream populations, agriculture and industries along these rivers (Pritchard, 2019; Viviroli et al., 2007). Climate change over recent decades has boosted river discharge by increasing runoff from shrinking glaciers (Lin et al., 2020; Yao et al., 2007; S. Q. Zhang et al., 2011), but this boost will eventually decrease as glacier area declines (Zhao et al., 2019). The sensitivity of ice loss to climate change is variable, however, and often poorly known, being a function of glacier size, hypsometry, aspect, debris cover, and the presence of proglacial lakes and ice cliffs, for example. Combined with uncertainties in ice thickness and future climate scenarios, the timing of peak water runoff and the rate of its subsequent decline remain key unknowns (Maurer et al., 2019; Nie et al., 2021; Su et al., 2016; Zhao et al., 2019). It is therefore critical to monitor and analyze glacier change to improve our understanding of its climate drivers and to assess its impacts on glacier-fed river basins.

Compared with the interpolation of sparse in situ measurements, satellites can observe glacier change over much larger areas of remote terrain (Q. Wang et al., 2021). In recent years, our understanding of the state of TP glaciers has been greatly improved by the increasing coverage and accuracy of multisource remote sensing observations of glacier area, volume and mass change from KH-9 (Hexagon military satellites), Landsat, ASTER, ICESat altimetry, and other digital elevation models (DEMs) constructed by geodetic techniques and from GRACE gravimetry (Guo et al., 2015; Kääb et al., 2012; Q. Wang et al., 2021; Zhou et al., 2018). Based on the KH-9 images and SRTM, for example, Zhou et al. (2018) found that from the mid-1970s to 2000 glaciers in the northwest TP thinned at $-\mathrm{0.11}\pm \mathrm{0.13}$ to 0.02±0.10 m w.e. a⁻¹, while those in the southeast part thinned faster at $-\mathrm{0.30}\pm \mathrm{0.12}$ to $-\mathrm{0.11}\pm \mathrm{0.14}$ m w.e. a⁻¹. Brun et al. (2018a) employed ASTER DEMs from 2000–2006 and showed that glacier mass balance in High Mountain Asia varied from $-\mathrm{0.62}\pm \mathrm{0.23}$ m w.e. a⁻¹ in eastern Nyainqentanglha to $+\mathrm{0.14}\pm \mathrm{0.08}$ m w.e. a⁻¹ in the Kunlun Mountains and averaged $-\mathrm{0.14}\pm \mathrm{0.14}$ m w.e. a⁻¹ over the large inner TP that includes the west Nyainqentanglha Range (WNR). Maurer et al. (2019) found a doubling of the Himalayan average loss rate between the periods 1976–2000 ($-\mathrm{0.22}\pm \mathrm{0.13}$ m w.e. a⁻¹) and during 2000–2016 ($-\mathrm{0.43}\pm \mathrm{0.14}$ m w.e. a⁻¹) using KH-9 and ASTER DEMs. These studies showed that glacier changes on and around the TP have marked spatial and temporal heterogeneity, likely associated in part with variable glacial sensitivity to climate change (Yao et al., 2012).

The drivers of regional glacier loss include, for example, a jump in mean annual temperature and precipitation in the Yarlung Zangbo River basin around 1997 (J. R. Wang et al., 2021) and an accelerating warming trend over the TP between the periods 1980–1997 and 1998–2013 (from 0.21 to 0.25 ^∘C decade⁻¹) (Duan and Xiao, 2015). Modulating the effect of these climatic drivers are local factors including glacier topography, debris cover and glacial lakes (Brun et al., 2018b, 2019; Ke et al., 2020; Maurer et al., 2019; Pandey et al., 2017; Yao et al., 2012). Some studies have analyzed the melt-inhibiting effect of debris cover and melt-promoting effect of proglacial lakes on glacier ablation since 2000 (Ke et al., 2020; Vincent et al., 2016), but with the potential for KH-9 in 1976, SRTM in 2000, ASTER in 2000–2020 and aerial mapping Landsat in 1976–2020 through time, we are now able to assess glacier area and mass change in the WNR in relation to both regional climatic drivers and local modulating factors.

The WNR, in the southeastern TP (Fig. 1), is located in the transition zone between the two large-scale atmospheric circulation patterns characterized by dominant westerlies and the Indian summer monsoon. It holds an abundance of glaciers and glacial-fed lakes, notably Nam Co (Fig. 1), whose rising water levels indicate a water imbalance primarily due to recently intensified glacier melting (Bolch et al., 2010), as supported by mass balance data from the Zhadang glacier and other hydrological observations from 2007 to 2011 (Zhou et al., 2013). The number and area of supraglacial lakes (of > 0.0036 km²) in the WNR also increased between 1976 and 2018 by 56 % and 35 %, respectively, due to the increase in glacial meltwater (Luo et al., 2020). In the relatively densely populated Lhasa Basin to the southeast of the WNR, Lin et al. (2020) found that a water imbalance also existed using the first and second Chinese Glacier Inventory in 1960 and 2009. Despite these extensive changes and the large affected population, logistical constraints have meant that in situ glacier mass balance records are limited to a few low-lying, small glaciers that are unlikely to be representative of the broader region (Kääb et al., 2012; Li and Lin, 2017; Yao et al., 2012). Similarly, glacier volumes in the Chinese Glacier Inventory were primarily calculated indirectly by area–volume scaling, and limited direct observations mean that these volumes have larger uncertainty (Bahr et al., 1997, 2015). Detailed investigation of the WNR glacier area change and mass balance on a longer timescale is therefore a high priority.

Figure 1 (a) Overview of the study area. (b) Glacier distribution. Label I in the large, dotted red rectangle represents the SW section of the WNR, and Label II in the small, dotted dark red rectangle represents the NE section.

Investigations of the WNR glacier area have so far focused on the period before 2014 (Bolch et al., 2010; Wu et al., 2016). For glacier mass balance, most studies focus on after 2000 (Li and Lin, 2017; Neckel et al., 2014; Ren et al., 2020; Zhang and Zhang, 2017). There are limited discussions on local glacier changes in the WNR region from before 2000, although Zhou et al. (2018) included this area in his study of glacier mass balance on the TP and its surroundings from the mid-1970s to 2000 and did not present the characteristics of glacier changes in detail. Furthermore, the warming rate of the TP is heterogeneous both spatially and temporally in recent decades (Duan and Xiao, 2015; Wu et al., 2015). Under a changing runoff regime (Lin et al., 2020), the lack of a detailed survey of glacier changes over a long timescale is a major impediment to water resource management and decision-making (Lutz et al., 2014).

The key purpose of this study is therefore to provide an internally consistent dataset of glacier area and mass change in the WNR over the past 44 years and comparative analysis of the impact of topography, debris and proglacial lakes on glacier change during 1976–2000 and 2000–2020. Although the areas of both debris cover and lake-terminating glaciers are relatively small, the characteristics of their influence on individual glaciers can be used as a reference for glacier changes in other regions. We have compiled a complete glacier inventory in the years 1976, 2000, 2014 and 2020 with the Landsat and KH-9 images and quantified the geodetic glacier mass balance during 1976–2000 and 2000–2020 with DEMs derived from KH-9, ASTER and SRTM3.0. We report area and mass changes for the periods 1976 to 2000 and 2000 to 2020, and we examine and compare the influence of topographic, climatic and glaciological factors on these changes during 1976–2000 and 2000–2020.

2 Materials and methods

2.1 Study area

The WNR has a mean slope of 15^∘, and its elevation spans 4150–7125 m, with an average of 4930 m in the whole region. Its primary mountain ridge runs 230 km in a southwest–northeast direction, bounded by the Nam Co Basin to the north and the Lhasa River basin to the south (Yao et al., 2010). Nam Co, the second largest lake after Selin Co on the TP, is mainly fed by glacier meltwater (Luo et al., 2020; Zhang et al., 2017). The Lhasa River basin to the southeast is a major branch of the Yarlung Zangbo River and forms part of the route taken by the warm and humid monsoon airflow into the plateau, making it warmer and wetter than the Nam Co Basin (Luo et al., 2020). The annual air temperature and precipitation in the WNR range from −0.6 to 2.8^∘ and from 37 to 500 mm, respectively (Yu et al., 2013).

Being in a climatic transition zone, the glaciers in this area range from the maritime-influenced glaciers of Gangrigabu (southeast TP) to the subcontinental- and continental-type glaciers of the Tanggula Mountains (Li and Lin, 2017). There are 845 glaciers covering 675.85 km² and 15 debris-covered glaciers with a total area of 71.74 km² (10.61 % of the total glacier area in the WNR) (RGI 6.0) (RGI Consortium, 2017). Of these, only the small, ∼3 km² Zhadang and Gurenhekou glaciers (red polygon in Fig. 1) have in situ observations available to validate the satellite observations, and these run from 2005 to 2010 (Yao et al., 2012).

2.2 Methods and data

2.2.1 Glacier outlines

We identified glacier boundaries mainly from Landsat MSS, ETM+ and OLI scenes from various years (Table S1 in the Supplement), orthorectified automatically by the USGS using the level 1T SRTM3 DEM (from http://glovis.usgs.gov/, last access: 5 January 2021). We selected high-quality images with minimal cloud and snow coverage between June and November and used a semi-automated approach with a TM3 $/$ TM5 band ratio (2.0±0.2) to produce glacier outlines. This method is widely used and appropriate for glacier mapping over large study areas (Guo et al., 2015; Ye et al., 2017). We used a 3×3 median filter to eliminate isolated pixels and those likely to have been misclassified due to debris or boulders on the glacier (Bolch et al., 2010). We manually checked and edited the glacier outlines, including the debris-covered glaciers, with height change maps and a coherence map formed by Sentinel-1 images observed on 5 and 29 August 2016, to help distinguish debris-covered ice from ice-free areas. Finally, referring to the second glacier inventory, we assigned contiguous ice masses to drainage basins in order to obtain a glacier inventory.

2.2.2 Glacier elevation change

We used the KH-9 images and SRTM DEM (version 3) to estimate glacial elevation changes for the period 1976 to 2000 and ASTER DEMs for the period 2000–2020.

DEM data

The declassified KH-9 images were obtained by the Hexagon mission from 1971 to 1986, with a ground resolution of 6 to 9 m (Surazakov and Aizen, 2010). We downloaded images from 7 January 1976 via the EarthExplorer user interface (https://earthexplorer.usgs.gov, last access: 5 January 2021) and adopted the Hexagon imagery automated pipeline method to generate digital elevation models. This method is coded in MATLAB and uses the OpenCV library for Oriented FAST and Rotated BRIEF (ORB) feature matching, uncalibrated stereo rectification, and semiglobal block matching algorithms (Maurer and Rupper, 2015).

The SRTM mission carried out in February 2000 produced two types of DEM datasets, the C-band DEM with a coverage range of 60^∘ N–60^∘ S and the X-band DEM with a smaller coverage. We used version 3 of the C-band SRTM DEM (https://earthexplorer.usgs.gov/, last access: 5 January 2021) at 1 arcsec resolution (about 30 m) in our primary processing and masked out areas with gaps in the unfilled SRTM3 version 2.1 DEM at 3 arcsec resolution (about 90 m).

The ASTER instrument was launched on the Terra satellite in December 1999, and a single DEM covers approximately 3600 km². We downloaded 250 “Data1. l3a.demzs” ASTER DEMs at 30 m resolution in GeoTIFF format with cloud coverage of less than 40 % from the METI AIST Data Archive System (MADAS) satellite data retrieval system (https://gbank.gsj.jp/madas, last access: 5 January 2021). After cloud and outlier removal we fitted a linear regression through the time series of coregistered ASTER DEMs and set the minimum stack interval per pixel to 15 years to estimate the rate of elevation change for each 30 m pixel (Maurer et al., 2019).

Coregistration and bias correction of DEMs

All DEMs were coregistered to the SRTM master DEM using a standard elevation–aspect optimization procedure (Nuth and Kääb, 2011). Then, the elevation correlation deviation of all the DEMs was corrected by a third-order polynomial. In addition, we used a 2 km buffer zone around the union of glacier boundaries to define stable (unchanging) terrain for DEM alignment, bias correction and uncertainty calculation. Supplement Fig. S1a shows the coverage of the KH-9 images and the number of valid ASTER DEM grids after removal of clouds and outliers in the buffer. The glacier area covered by the dataset from 1976 to 2000 and from 2000 to 2020 accounted for 70.85 % and 81.94 % of the total glacier area, respectively (shown in Fig. S1b, c). We used only the area common to both of these datasets to measure elevation change between the 1976–2000 and 2000–2020 periods. After correction for alignment and elevation-related deviation, apparent elevation changes over stable terrain (masked glaciers and lakes in square buffer zone) had no change trend with elevation, slope and aspect, as shown in Fig. S2.

SRTM penetration depth correction

Over the WNR, the average penetration depth of C-band SRTM is 1.67±0.53 m calculated using X-band SRTM DEM as the reference (Li and Lin, 2017). Linear regression between the glacier elevation and penetration showed that the penetration depth varies from 1.29 to 2 m at altitudes of 5550 and 6250 m, respectively (Li and Lin, 2017). We used this more accurate linear altitude-dependent correction, and the result is similar to several other study regions on the TP (Gardelle et al., 2013; Kääb et al., 2012; Li and Lin, 2017).

Glacier mass change

Estimation of average glacier thickness changes based on elevation difference maps involves noise filtering and glacier-hypsometry-weighted averages in an approach widely employed to calculate regional glacier mass balance where glacier thinning is highly dependent on altitude. Firstly, we subjected the thickness change maps to outlier removal using a 5 m a⁻¹ threshold. We then masked slopes > 45^∘, where uncertainties are large, before visually inspecting the final thickness change maps. We additionally masked out any remaining anomalous pixels, which occurred almost exclusively in low-contrast, snow-covered accumulation zones. Finally, we separated thickness changes into 50 m elevation bins by referring to the SRTM at different spatial scales, i.e., the whole glacierized area, subregions, different glacier types and individual glaciers of area > 2 km². In each altitude bin, we filtered out any height-change values that differed by more than 3 standard deviations from the median and removed any bins with less than 100 pixels. For elevation bins with no observations (mostly over the low- and high-elevation limits), we assumed zero mean elevation changes. We calculated the mean glacier thickness changes for the spatial unit/group (dh) as a hypsometric average:

$$\begin{array}{}\text{(1)}& \mathrm{d}h=\sum _{i=\mathrm{1}}^n{\displaystyle \frac{S_i}{S}}\cdot \stackrel{\mathrm{‾}}{\mathrm{d}{h}_i},\end{array}$$

where i and n denote the ith 50 m elevation bin and the number of total bins, respectively; S_i is the glacier area of the ith elevation bin; S is the total glacier area; and $\stackrel{\mathrm{‾}}{\mathrm{d}{h}_i}$ is the mean dh in the bin.

We calculated the final geodetic mass balance (B) using Eq. (2).

$$\begin{array}{}\text{(2)}& B=\mathrm{d}{h}_i\times {\displaystyle \frac{{\mathit{\rho }}_{\mathrm{ice}}}{{\mathit{\rho }}_{\mathrm{water}}}}\end{array}$$

We translated glacier thickness changes into mass balance by the ratio of column-averaged glacier density ρ_ice (850 kg m⁻³) to water density (ρ_water, 1000 kg m⁻³).

2.2.3 Uncertainty

Uncertainty of the glacier area

Similar to previous studies (Wu et al., 2016; Ye et al., 2017), we obtained the uncertainty of the glacier area (δ_s) using Eq. (3):

$$\begin{array}{}\text{(3)}& {\mathit{\delta }}_s=L_{\mathrm{c}}{E}_{\mathrm{pc}}+L_{\mathrm{d}}{E}_{\mathrm{pd}},\end{array}$$

where L_c and L_d represent the lengths of the clean-ice and debris-covered glacier outlines, and E_pc and E_pd denote the positional accuracies. We calculated the uncertainty in glacier area change (δ_sc) by combining the area uncertainties using Eq. (4).

$$\begin{array}{}\text{(4)}& {\mathit{\delta }}_{s\mathrm{c}}=\sqrt{({\mathit{\delta }}_{s\mathrm{1}}{)}^{\mathrm{2}}+({\mathit{\delta }}_{s\mathrm{2}}{)}^{\mathrm{2}}}\end{array}$$

Guo et al. (2015) compared glacier outlines derived from Landsat images with real-time kinematic differential GPS (RTK-DGPS) measurements and found an average difference of ±11 and ±30 m for the delineation of clean and debris-covered ice. Using a buffer size of 10 m for areas from the Hexagon images (Bolch et al., 2010), our combined uncertainty in the glacier area is 3.9 %, 5.1 %, 5.1 % and 5.9 % in 1976, 2000, 2014 and 2020, respectively.

Uncertainty of glacier thickness change

The uncertainty in surface-elevation change derived from ASTER DEMs can be estimated using the point elevation error (E_pt) and extrapolation error (E_ext) (Nuth and Kääb, 2011; Maurer et al., 2016).

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{\delta }}_{hi}=\sqrt{{\left({\displaystyle \frac{E_{\mathrm{pt}}}{\sqrt{n_i}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{E_{\mathrm{ext}}}{\sqrt{n_i}}}\right)}^{\mathrm{2}}}\text{(6)}& {\displaystyle }{n}_i={\displaystyle \frac{n_{ib}\cdot r^{\mathrm{2}}}{\mathit{\pi }\cdot d^{\mathrm{2}}}}\text{(7)}& {\displaystyle }{\mathit{\delta }}_h=\sqrt{\sum _{i=\mathrm{1}}^{i=n}({\mathit{\delta }}_{hi}\ast {\displaystyle \frac{S_i}{S}}{)}^{\mathrm{2}}}.\end{array}$$

Here, E_pt refers to the standard deviations of the relative elevation change over the off-glacier areas; E_ext is the standard deviations of glacial elevation change within each 50 m bin; n_ib and n_i represent the total number of pixels and the number of independent measurements of pixels, respectively; r is the DEM spatial resolution (30 m in our study); and d is the autocorrelation length. We used an autocorrelation length of 500 m, which is a conservative value based on semivariogram analyses of mountainous regions in previous studies (Brun et al., 2018a; Maurer et al., 2019). We combined the uncertainty of surface-elevation changes derived from the KH-9 DEM with the SRTM penetration uncertainty, estimated as ±0.53 m (Li and Lin, 2017). This study ignored the errors caused by seasonal changes in glacier thickness due to the lack of observations of such seasonal changes.

We estimated the overall uncertainty in the total glacier mass change (δ_m, in kg) by including the uncertainty in the assumed ice, firn and/or snow density (δρ= 60 kg m⁻³, which is 7.1 % of ρ_ice= 850 kg m⁻³), errors in the glacier area (δ_S, m²) and glacier elevation change (δ_h, m), using Eq. (8).

$$\begin{array}{}\text{(8)}& {\mathit{\delta }}_m=\sqrt{(S\cdot \mathrm{d}h\cdot {\mathit{\delta }}_{\mathit{\rho }}{)}^{\mathrm{2}}+({\mathit{\delta }}_s\cdot \mathrm{d}h\cdot {\mathit{\rho }}_{\mathrm{ice}}{)}^{\mathrm{2}}+(S\cdot {\mathit{\delta }}_h\cdot {\mathit{\rho }}_{\mathrm{ice}}{)}^{\mathrm{2}}}\end{array}$$

2.2.4 Lake data

We identified glacier-marginal lakes as those lying within 50 m of a glacier boundary, using lake data from the 1970s to 2018 (Luo et al., 2020; http://data.tpdc.ac.cn, last access: 10 May 2021).

2.2.5 Meteorological data

There are three meteorological stations adjacent to the WNR, at Bange (31^∘23^′ N, 90^∘01^′ E; elevation of 4700 m), Lhasa (29^∘40^′ N, 91^∘08^′ E; elevation of 3648 m) and Damxung (30^∘29^′ N, 91^∘06^′ E; elevation of 4200 m). We obtained air temperature and precipitation data during 1976–2020 from the Climatic Data Center, National Meteorological Information Center, of the China Meteorological Administration.

We also obtained gridded data of precipitation and temperature with a spatial resolution of 0.1^∘ × 0.1^∘ and a 3 h time interval for 1979–2018 from the China Meteorological Forcing Dataset (Yang et al., 2019; http://data.tpdc.ac.cn, last access: 5 June 2021), which has been widely utilized in land processes, hydrological modeling and other studies (Qiao et al., 2021; J. R. Wang et al., 2021). This dataset is made by fusing the conventional meteorological observation of the China Meteorological Administration based on the Princeton reanalysis data, GLDAS data, GEWEX-SRB radiation data and TRMM precipitation data as the background field (He et al., 2020; Yang et al., 2010).

2.2.6 Hydrological data

In order to assess hydrological changes under glacier retreat, we collected runoff data of the Lhasa River station during 1976–2013 and the Yangbajain station during 1979–2013 from the Tibet Autonomous Region's Hydrology and Water Resources Survey Bureau.

We calculated the ratio of total glacier mass change to runoff in Lhasa River basin (R_r, %) and the total lake water storage change of Nam Co (R_l, %) as follows:

$$\begin{array}{}\text{(9)}& {\displaystyle }{R}_{\mathrm{r}}={\displaystyle \frac{\mathrm{\Delta }M\cdot S_{\mathrm{g}}}{S_{\mathrm{r}}{R}_{\mathrm{a}}}}\text{(10)}& {\displaystyle }{R}_{\mathrm{l}}={\displaystyle \frac{\mathrm{\Delta }M\cdot S_{\mathrm{g}}}{\mathrm{\Delta }V}},\end{array}$$

where ΔM, S_g, S_r, R_a and ΔV represent the average annual glacier mass balance, glacier area, area of the Lhasa River basin, average runoff depth and lake water storage increase.

3 Results

3.1 Glacier area change

There were 921 glaciers with a total area of 589.17±31.72 km² in 2020 in the WNR (Fig. 2a). Small glaciers dominated the number (those ≤1 km² occupy 83.17 % of the total number) and a large proportion of the area (those ≤1 km² occupy 30.42 % of the total area). Glaciers larger than 5 km² accounted for 21.39 % of the total area and only 1.63 % of the total number. Glaciers were mainly distributed in the eastern-oriented zone with an altitude of 5600–6100 m and a slope of 5–40^∘ (Fig. 2b, c, d).

Figure 2 Glacier distribution in the WNR in 1976, 2000, 2009, 2014 and 2020. (a) Number and area of glaciers by size category. (b) Distribution of glacier area with altitude. (c) Distribution of glacier area with slope. (d) Distribution of glacier area with aspect. Data in 2009 came from the Chinese Glacier Inventory II.

Glaciers in the WNR experienced significant retreat from 1976 to 2020, and altitude, slope and aspect all appear to have influenced this retreat (Figs. 3 and 4). The glacier area decreased from 884.90±29.71 km² in 1976 to 589.17±31.72 km² in 2020, with an average annual decrease of $-\mathrm{0.76}\pm \mathrm{0.11}$ % a⁻¹. The retreat rate of the glacier area in 2000–2020 ($-\mathrm{1.17}\pm \mathrm{0.30}$ % a⁻¹) was more than twice as fast as in 1976–2000 ($-\mathrm{0.54}\pm \mathrm{0.21}$ % a⁻¹) (Figs. 3 and 4, Table 1). The glacier area declined faster in the northeast and southwest but slower in the middle, except for a few small glaciers with an area of less than 2 km² during 1976–2000. During 2000–2020 the glacier area receded faster in the whole region except for a few small glaciers with an area of less than 1 km² (Fig. 3). Retreat was greatest in the area classes of 1–3 and 3–5 km², and glaciers with significant areal retreat were mainly distributed below 6000 m altitude. Glaciers in the Nam Co Basin retreated slightly faster than those outside this basin between 2000 and 2020. Retreat was particularly rapid at lower altitudes and decreased at higher elevations. As for the effect of slope and aspect, glaciers retreated more rapidly with increasing slope between 5 and 40^∘, but the retreat rate decreased as slope increased between 0–5 and 40–60^∘, where relatively few glaciers are distributed. In each elevation band, we found a positive correlation between areal retreat rates and slope (faster retreat with steeper slope) for most elevation bands and in both time periods (Fig. 9a and b). The only areas where this relationship differed were on flat or shallow slopes at lower altitudes (slopes below about 5^∘ at elevations below about 5500 m, e.g., blue lines in Fig. 9a), which also experienced relatively rapid retreat. During both 1976–2000 and 2000–2020, the retreat rate was slowest on the north-facing slopes. During 1976–2000, retreat was most rapid in the southeast quadrant, while from 2000 to 2020, rapid retreat occurred at similar rates in all aspects other than in north and southeast, i.e., the effect of an aspect on glacier area retreat varied in space and time.

Figure 3 The distribution of glacier area change in the WNR from (a) 1976 to 2000, (b) 2000 to 2020 and (c) 1976 to 2020.

Figure 4 Glacier area changes with (a) time, (b) elevation, (c) slope and (d) aspect. The short lines on either side of the point indicate the margin of error in panels (a), (b) and (c).

Table 1 Glacier area changes over the WNR from 1976 to 2020.

Besides, we noticed that the area decreased, but glacier numbers increased in the WNR (Fig. 2a). The reason for this is that intact glaciers break down into several smaller glaciers in the process of glacier ablation, e.g., a large glacier in 1976 may become several smaller glaciers in 2020 (shown in Fig. S3).

3.2 Geodetic mass balance

Glacier height changes for the past 44 years are shown in Fig. 5. Substantial and near-ubiquitous thinning occurred in the WNR during 1976–2020 (with mean surface lowering of $-\mathrm{0.37}\pm \mathrm{0.13}$ m a⁻¹, a water-equivalent loss rate of $-\mathrm{0.31}\pm \mathrm{0.12}$ m w.e. a⁻¹ or a mass loss rate of $-\mathrm{0.26}\pm \mathrm{0.09}$ Gt a⁻¹), with a widespread increase in the most recent decades. From 1976 to 2000, glaciers experienced a mean thinning rate of $-\mathrm{0.31}\pm \mathrm{0.10}$ m a⁻¹ ($-\mathrm{0.26}\pm \mathrm{0.09}$ m w.e. a⁻¹) or $-\mathrm{0.24}\pm \mathrm{0.08}$ Gt a⁻¹. From 2000 to 2020, the mean elevation rate was $-\mathrm{0.44}\pm \mathrm{0.13}$ m a⁻¹, (0.37±0.12 m w.e. a⁻¹) or $-\mathrm{0.29}\pm \mathrm{0.09}$ Gt a⁻¹. Several glacier tongues have suffered severe thinning, exceeding −1.5 m a⁻¹ from 1976 to 2000, notably several long, debris-free glaciers on the southwestern slope. From 2000 to 2020, the range of glacier tongues with losses exceeding 1.5 m a⁻¹ expanded, and losses were greater in the central WNR (see the red rectangular box in Fig. 5). In both 1976–2000 and 2000–2020, the glacier thinning rate was slightly higher inside the Nam Co drainage basin than outside it (Table 2, Fig. 6), though these rates do not differ by more than their combined uncertainties.

Figure 5 Mean annual glacier surface elevation changes in the WNR from (a) 1976 to 2000, (b) 2000 to 2020 and (c) 1976 to 2020. Label I in (a), (b) and (c) represents the SW section and label II in (b) represents the NE section of the WNR (on the same scale). The red rectangular box in (b) shows an area of the central WNR referred to in the paper.

Figure 6 Glacier elevation change with altitude (m a.s.l) in the whole WNR, inside the Nam Co drainage basin and outside the Nam Co drainage basin from (a) 1976 to 2000 and (b) 2000 to 2020. The dots represent the mean elevation change in each 50 m elevation bin, and shaded regions in the altitudinal distributions indicate the uncertainty.

Table 2 Glacier elevation change, mass balance and total mass change over the WNR from 1976 to 2020.

For glaciers with an area of more than 2 km², we found high loss rates in the northeast, followed by the southwest, and moderate in the middle during 2000–2020, but there was no obvious spatial varied trend of mass loss during 1976–2000 (Fig. 7). Mass loss was substantially more intense in 2000–2020 with no glaciers in a state of positive balance (Fig. 7, blue dots), and loss from some glaciers in the northeast exceeded −0.6 m w.e. a⁻¹. Moreover, we found that glacier area retreat and mass loss were not synchronized between the two periods 1976–2000 and 2000–2020. The glacier with the fastest area retreat did not correspond to the glacier with the fastest mass decrease, and the spatial varied trend of glacier area retreat rate was inconsistent with that of mass loss rate (Figs. 3 and 7).

Figure 7 The distribution of glacier-wide mass balance for individual glaciers (> 2 km²) in the WNR from (a) 1976 to 2000 and (b) from 2000 to 2020. Label I represents the SW section, and label II represents the NE section of the WNR (on the same map scale).

Figure 8 Glacier elevation change from 1976 to 2000 and from 2000 to 2020 with (a) elevation, (b) slope and (c) aspect. The dots in panel (a) represent the mean elevation change in each 50 m bin, and shaded regions in (a) indicate the uncertainty in the altitudinal distributions. Panel (b) is the boxplot of elevation change in 2^∘ slope bins, and four lines from the bottom to top for one box represent minimum value, 25th percentile, 75th percentile and maximum value, respectively, and dots in panel (c) represent the mean elevation change in each 2^∘ slope bin. Panel (c) represents the mean elevation change in each 45^∘ aspect bin.

Figure 9 Glacier area changes with slope during 1976–2000 (a) and 2000–2020 (b) and glacier elevation changes with slope during 1976–2000 (c) and 2000–2020 (d).

Finally, glacier elevation change as a function of elevation slope and aspect is shown in Fig. 8. Elevation is inversely correlated with thickness change, while slope and aspect appear to have a weak relationship with thickness change. In both 1976–2000 and 2000–2020, the elevation change rate was the largest at lower altitudes and gradually decreased with the increasing of altitude. The thinning rate also exhibited a weak inverse relationship with slope, becoming somewhat stronger in the 2000–2020 period. A very similar but inverse relationship between slope and glacier area change rate in each elevation band is that thinning rates were highest on shallow slopes and decreased over steeper slopes, except for flat or shallow slopes at lower altitudes where thinning rates were relatively low (Fig. 9c and d). For the impact of aspect, thinning for 1976–2000 was most rapid in the southwest and northwest quadrants, but by 2000–2016 high thinning rates were affecting all aspects, i.e., the effect of aspect on thinning rates also varied through time.

3.3 The effect of debris cover and proglacial lakes on glacier mass changes

In our WNR study area, there are five debris-covered glaciers, covering 55.42±1.25 km² in 1976 and 51.59±1.77 km² in 2000. Lake-terminating glaciers occupied a similar proportion, with an area of 70.29±2.69 km² in 1976 and 49.60±1.82 km² in 2000. Only one glacier was both covered by debris and terminated in a proglacial lake.

The thinning rate of different types of glaciers varied somewhat, though with greater uncertainty given the relatively small sample (Fig. 10, Table 3). During 1976–2000, the lake-terminating glaciers thinned more rapidly, followed by the regular and debris-covered glacier types. From 2000 to 2020, the ablation rate of debris-covered glaciers was slightly lower than that of regular glaciers at low altitudes but progressively greater at higher altitudes, leading to a slightly more negative total mass balance for the debris-covered type. In the same period, the thinning of lake-terminating glaciers continued to exceed that of regular glacier. Our results suggest that debris cover in the WNR suppressed glacier thinning to some extent and enabled the debris-covered ice to survive at lower elevations than adjacent clean-ice glaciers. In contrast, a glacial lake at the end of a glacier accelerated its retreat, and this behavior was more pronounced at lower elevations.

Figure 10 Rate of glacier elevation change with elevation of different glacier types during (a) 1976–2000 and (b) 2000–2020. Plots represent the mean values of glacier elevation change in each 50 m elevation bin, and shaded regions indicate the uncertainty in the altitudinal distributions.

Table 3 Statistics of area, quantity and mass balance of different types of glaciers.

4 Discussion

4.1 Comparison to previous studies

4.1.1 Glacier area change

Based on space-borne imagery, we found that glacier area in the west Nyainqentanglha Range (WNR) has changed by $-\mathrm{12.98}\pm \mathrm{4.91}$ % during 1976–2000 and $-\mathrm{23.49}\pm \mathrm{5.99}$ % during 2000–2020. The comparison between the Chinese Glacier Inventory (CGI) I and CGI II over the WNR and previous studies is shown in Tables S2 and S3. The CGI II of the WNR in 2009 is in good agreement with the areal retreat trend in our study (also shown in Fig. 2). The CGI I of the WNR in 1970 is slightly smaller than the glacier area in 1976 in our study, but it is within the margin of uncertainty. The CGI I was mapped based on the Chinese topographic maps, while glacier area in our study was mapped based on Landsat images. The difference between them might come from this difference in data source used to extract the glacier outlines. Besides, Frauenfelder and Kääb (2009) reported that there are georeferencing errors in the areas in the CGI I. Differences between studies may have arisen from the georeferencing errors in the areas for 1970 used by Shangguan et al. (2008) and Wu and Zhu (2008) which came from the CGI I. Discrepancies may also have arisen from differences in the methods used to distinguish glaciers from seasonal snow and debris-covered glaciers from neighboring moraine or rock slopes (Bolch et al., 2010). The deviation between our results and those from Wu et al. (2016) and Wang et al. (2012) over the whole WNR and the southwest WNR is within the margin of uncertainties. In addition, the 789.15 km² area reported for the WNR by RGI V4.0, which used Landsat images obtained on 6 December 2001, agrees with our result.

4.1.2 Glacier mass balance

Field measurements of mass balance are available from the small Zhadang glacier for 2005–2008 and Gurenhekou glacier for 2005–2010 on the southeastern slope of the WNR (Table S4). Although the period of our study is longer and provides a much larger sample size, the mass balance results are similar to these field measurements.

Previous studies have also reported region-averaged glacier mass balance over a similar spatial extent to ours, obtained from DEMs using various sensors (Table S5). Our results during 2000–2020 are more negative than those of Neckel et al. (2014), Li and Lin (2017), and Zhang and Zhang (2017) but agree within the uncertainties over comparable time periods, even though these studies differ in data processing, glacier mask, penetration correction and data coverage. For comparison, we calculated the change for the 2000–2014 period from ASTER DEMs (Fig. 11). Our estimated mass balance in this area ($-\mathrm{0.28}\pm \mathrm{0.15}$ m w.e. a⁻¹) is very similar to the other studies (Table S5). It is also similar to that of 1976–2000, suggesting that the more strongly negative average for the longer 2000 to 2020 period ($-\mathrm{0.37}\pm \mathrm{0.12}$ m w.e. a⁻¹) is mainly due to the intensified glacier ablation after 2014, although cloud-free ASTER data are insufficient for direct calculation of the mass balance from 2014–2000. The glacier area retreat during 2014–2020 (1.53±1.14 % a⁻¹) is also faster than that during 2000–2014 (1.13±0.42 % a⁻¹), though the change is within the uncertainties. This interpretation is supported by Ren et al. (2020) who also calculated a higher 2013–2020 thinning rate ($-\mathrm{0.43}\pm \mathrm{0.06}$ m w.e. a⁻¹), twice as negative as in 2000–2013. Though the difference in rate is within the combined uncertainties for these periods, this apparent acceleration in thinning in the WNR (from $-\mathrm{0.26}\pm \mathrm{0.06}$ m w.e. a⁻¹ in 1976–2000 to $-\mathrm{0.37}\pm \mathrm{0.15}$ m w.e. a⁻¹ in 2000–2020) is similar to the broader regional pattern of accelerating loss across the Himalayas and Kangri Karpo Mountains (Maurer et al., 2019; Wu et al., 2018, 2019).

Figure 11 (a) Glacier elevation change in the WNR during 2000–2014. (b) Glacier elevation changes with altitude in the WNR, inside the Nam Co drainage basin and outside the Nam Co drainage basin from 2000 to 2014. The dots represent the mean elevation change in each 50 m elevation bin, and shaded regions indicate the uncertainty in the altitudinal distributions. (c) Total area of glaciers and that area covered by the datasets during 1976–2000 and 2000–2014.

4.2 The influences of debris cover and proglacial lakes on glacier mass changes

Debris can inhibit or enhance glacial ablation depending on its thickness (Maurer et al., 2016). A shallow layer of debris usually enhances melt rates due to its low surface albedo, while thicker layers could suppress melt rates through thermal insulation (Reid et al., 2012). Our results (Table 3) suggest that the debris-covered glaciers in our study thinned more slowly than the regular, debris-free glaciers in the 1976–2000 period, though the difference is not statistically significant, and the small sample size (5) of the debris-covered glaciers compared to regular glaciers (> 600) limits our ability to compare these classes. In the 2000–2020 period, the thinning rate of the debris-covered glaciers increased significantly, to double its previous rate, though it remains indistinguishable from the thinning rate for regular glaciers at that time. While several previous studies indicated that on the glacier scale debris-covered glaciers thin more slowly than debris-free glaciers (Nicholson and Benn, 2006; Scherler et al., 2011; Vincent et al., 2016), large-scale geodetic studies reported no significant differences in the thinning rates between debris-covered and clean glaciers on timescales more than a decade after 2000 (Brun et al., 2019; Ke et al., 2020; Maurer et al., 2019), a finding that is supported by this study, and Banerjee (2017) suggested that the thinning rate of a debris-covered glacier is initially slower than that of a similar clean glacier at the early stage of warming but subsequently matches and then overtakes the clean counterpart. In this theory, the time required for their respective melting rates to cross is controlled by the rate of warming, with little difference between their thinning rates at low rates of warming (Banerjee, 2017). The large difference in the 1976–2000 mean melt rates of the regular versus debris-covered glaciers in our study provides some support for this theory, but a larger sample with lower uncertainty is needed to verify this.

For debris-covered glaciers, the area of debris cover actually increased from 6.60±1.15 km² in 1976 to 7.37±1.49 km² in 2020 in our study (Table S6), and we note that this is not necessarily inconsistent with an overall glacier retreat. This is because increased melt rates that lead to surface lowering drive retreat of the glacier front, while also promoting a greater concentration of debris on the wider surface of the glacier ablation area, as more debris melts out from ice below. A spatial expansion of the debris layer has, for example, been observed on different debris-covered glaciers during retreat and sustained mass loss (Kirkbride and Deline, 2013; Stokes et al., 2007; Tielidze et al., 2020; Xie et al., 2020). Unfortunately, no data are available relating to changes of the thickness of the debris cover itself, and we assume that all glacier thickness changes resulted from loss of ice, without considering the thickness change of the debris-covered layer. We think that this is reasonable, because in most areas debris layers are typically thin (order of 1 m or less) compared to elevation changes we map, and because most debris cover in the ablation area emerges from englacial transport rather than direct deposition by new, local rockfall (e.g., McCarthy et al., 2017), so changes in the debris-layer thickness represent a redistribution of existing glacier volume, not a change in volume.

Glaciers with proglacial lakes can experience relatively high mass loss through calving and thermal undercutting, (Maurer et al., 2016; Thompson et al., 2012) and the expansion of such lakes can cause dynamic thinning to propagate glaciers (Ke et al., 2020). Glaciers terminating in proglacial lakes in our study area had the highest mean thinning rates of all the classes in both time periods and more negative mass balance compared to both regular and debris-covered glaciers during 2000–2020.

The area of debris-covered glaciers and lake-terminating glaciers decreased, while surface lowering also accelerated, mainly driven by the continuous increase in temperature in the WNR region during 1976–2000, especially after 2014 (Figs. 12 and 13), which was discussed in Sect. 4.3. In terms of the number and area of lake termination, we identified glacier-marginal lakes as those lying within 50 m of a glacier boundary. As glaciers retreat, the distance between the end of the glacier and their proglacial lake increased, and some of the lake-terminating glaciers in 1976 no longer belonged to the lake-terminating class in 2000. This helps to explain the area decrease for this glacier class, as shown in Table 3.

Figure 12 Temperature and precipitation changes for the study area at Damxung, Lhasa and Bange stations from 1976 to 2020. Annual average temperature and precipitation (a, b), ablation season (June to September) average temperature and precipitation (c, d), and accumulation season (January to May and October to December) average temperature and precipitation (e, f).

Figure 13 Gridded temperature and precipitation change during specific time periods.

4.3 Topographic and climatic controls of varying glacier mass loss

If climate is the driving force behind glacier change, topographical parameters can modulate this change (Pandey et al., 2017). Controls on glacier thickness and areal change are complicated with additional factors including local variations in climate, glacier thickness, morphology, the presence of proglacial and/or supraglacial lakes and debris cover, and latitude and longitude (Brun et al., 2018b, 2019; Ke et al., 2020; Maurer et al., 2019).

We found that both glacier areal retreat rate (Fig. 4b) and thinning rate (Figs. 8 and 9) generally decreased with increasing altitude, agreeing with previous studies (Li and Lin, 2017; Wu et al., 2016; Ye et al., 2017; Zhou et al., 2019). However, the effect of slope and aspect on glacier thickness has rarely been studied. We found that in the slope range of 8–40^∘, where the glaciers were mainly distributed, the rate of areal retreat increased as slope increased (Fig. 4c), but the thinning rate decreased (Fig. 8b). This occurred because steep slopes are associated with thinner ice (Linsbauer et al., 2012), which means that any given thinning rate will tend to drive more rapid areal retreat on steeper slopes as the thinner ice there is depleted first, explaining the broadly positive correlation between retreat and slope. Besides, steeper slopes are biased towards higher elevations, where the colder climate leads to slower thinning rates (dh), explaining the broadly negative correlation between slope and thinning rate. The somewhat different behavior of the low-elevation flat areas (relatively rapid retreat, relatively slow thinning, Fig. 9) may in part reflect the modulating effects of proglacial lakes (quicker retreat) and thicker debris cover (slower thinning) near the terminus. Overall, the relationship between aspect and both areal retreat and thinning was spatially inconsistent and varied in time (Figs. 4d and 8c).

Mean glacier mass thinning and retreat rates were consistently higher in the Nam Co Basin than in the Lhasa River basin (Tables 1 and 2), in agreement with Bolch et al. (2010) and Li and Lin (2017), and the glaciers in the central WNR showed particularly strong melting from 2000 to 2020. While the glacier distribution on the TP broadly follows the regional atmospheric circulation pattern (Yao et al., 2012), the variability in glacier loss within regions cannot always be fully explained by the changes in precipitation and temperature on this scale (Wu et al., 2018).

The increasingly negative mass balance through time is consistent with the temperature record from the three weather stations that shows a consistent warming trend (averaging 0.0485 ^∘C a⁻¹) (Fig. 12) and gridded temperature data showing a more rapid increase during 2000–2018 than 1979–2000 (Fig. 13), alongside precipitation that increased during 1979–2020 but decreased during 2000–2018. The accelerated warming from 2014 to 2018 (red rectangle in Fig. 13 (2014–2018)) corresponds geographically to the substantial central WNR glacier thinning highlighted in the red rectangle in Fig. 5. Precipitation also increased substantially in this region from 2014 to 2018, and glacier melting can be particularly intense under combined warm and wet conditions (Li et al., 2020; Oerlemans and Fortuin, 1992).

While the overall trends of temperature and precipitation in the ablation season (June to September) and accumulation season (October to December and January to May) were similar to annual changes, the temperature and precipitation data from 2014 to 2018 described above offer a compelling explanation for the main temporal and spatial variations in glacier change in the WNR, particularly the high rates of thinning from 2014–2018. They do not directly explain why the Nam Co glaciers thinned more rapidly than elsewhere, however. Other possible explanations include the difference in the impact of black carbon and dust in reducing surface albedo (Lau et al., 2010; Li and Lin, 2017; Ming et al., 2008), and Qu et al. (2014) observed a decrease in albedo at Zhadang glacier (Nam Co Basin) from 2001–2012.

4.4 Hydrological effect

The glacier melt contribution to streamflow decreases significantly from the glacier terminus to the lowlands, as it becomes diluted by other water sources (Kaser et al., 2010; Lutz et al., 2014; Pritchard, 2019), and this is reflected in our finding that the average annual glacier mass loss during 1976–2014 ($-\mathrm{0.26}\pm \mathrm{0.14}$ m w.e a⁻¹) equates to 8.5±4.6 % of the mean annual runoff depth for the Yangbajain Basin, in the upper reaches of Lhasa River (location shown in Fig. 1a), but only 1.6±1.0 % for the Lhasa River basin as a whole.

Through this period, the annual runoff in the Yangbajain Basin showed a significant increase trend of 1.22 mm a⁻¹ and the Lhasa River basin a non-significant increase trend of 0.84 mm a⁻¹ (Fig. S4). Increasing runoff may in part be explained by a coincident 1.36 mm a⁻¹ increase in precipitation observed over the Lhasa River basin (Fig. 12b), though the glacier ablation increase in the Lhasa River basin and Yangbajain Basin (4.63±2.49 and 23.52±12.67 mm a⁻¹, respectively) was substantially greater than the increase in precipitation, and evaporation losses from glacier melt water tend to be substantially smaller than those from evaporation of precipitation over the basin (Pritchard, 2019), suggesting that increased glacial meltwater primarily drove increased runoff. This is supported by Lin et al. (2020) who attributed increased streamflow at the Yangbajain station to accelerated glacier retreat and J. R. Wang et al. (2021) who argued that glacier melt has increased its contributions to the surface runoff by 12 %–43 % among the subbasins of the Yarlung Zangbo River basin (the mainstream of Lhasa River) after 1997.

Some components of basin hydrology remain poorly observed; however we note that the combined increase in precipitation and ablation detailed above was notably much greater than the observed increase in runoff especially in Yangbajain, a discrepancy that due to some combination of increased residential, industrial or agricultural water use (Pritchard, 2019) increased evaporation (Han et al., 2021) and possible deep seeps in the upper Lhasa River (Lin et al., 2020).

In the Nam Co Basin, increased glacier runoff also appears to have been important in controlling the level of Nam Co. The Nam Co Basin glacier mass balance (0.32±0.16 m w.e a⁻¹) that we found for 1976–2014 equated to 30.9±15.4 % of the reported increase in Nam Co water storage (B. Zhang et al., 2011). This glacier contribution is comparable to previous estimates of 52.9 % for the 1971–2004 period (Zhu et al., 2010), 28.7 % for 1999–2010 based on a mass balance of 0.59 m w.e a⁻¹ (Zhadang glacier) (Lei et al., 2013), 10.50±9.00 % for 2003–2009 based on a mass balance of $-\mathrm{0.27}\pm \mathrm{0.13}$ m w.e. a⁻¹ (Li and Lin, 2017) and 17.5±7.6 % for 2000–2014 based on a mass balance of −0.32 m w.e. a⁻¹ (Ke et al., 2022). Differences in these contributions of glaciers to increases in lake level reflect differences in the time periods studied and variability in the rate of change in the lake. For example, Ke et al. (2022) reported that their average lake-level change (0.26±0.04 m a⁻¹ for the 2000s to 2014) is substantially higher than 0.14±0.18 m a⁻¹ for 1994–2015 reported by Brun et al. (2020).

5 Conclusions

Based on KH-9, Landsat, SRTM and ASTER satellite data, we have quantified the changes of glacier area, surface elevation and mass balance in the WNR over the past 44 years and compared the effects of topography, debris cover and proglacial lakes on glacier change during 1976–2000 and 2000–2020. Our major conclusions are as follows:

1. Glaciers in the WNR retreated by 295.73±43.45 km², or 33.42±4.9 % of their area, from 1976–2020, at a mean rate of $-\mathrm{0.76}\pm \mathrm{0.11}$ % a⁻¹. Over this time, they lost a total of 11.56±0.12 Gt of ice. In addition, increased glacial meltwater in the WNR primarily drove increased runoff in the Lhasa River basin and appeared to have been important in controlling the level of Nam Co. The Nam Co Basin glacier mass balance (0.32±0.16 m w.e a⁻¹ during 1976–2014) equated to 30.9±15.4 % of the increase in Nam Co water storage.

2. The average retreat rate from 2000 to 2020 (1.17±0.30 % a⁻¹) was more than twice that from 1976 to 2000 (0.54±0.21 % a⁻¹). Similarly, the mean glacier mass balance from 2000 to 2020 ($-\mathrm{0.37}\pm \mathrm{0.12}$ m w.e a⁻¹) was more negative than that from 1976 to 2000 ($-\mathrm{0.26}\pm \mathrm{0.09}$ m w.e. a⁻¹) (though the change is within the uncertainties). The more rapid ice loss from 2000 to 2020 was mainly due to intensified glacier melting after 2014, which was likely associated with particularly strong warming of the region after that year. Besides, areal retreat rate and mass loss rate of most glaciers were not synchronized during 1976–2000 and 2000–2020.

3. In the WNR the spatial and temporal patterns of glacier loss can largely be explained by the observed patterns of regional climate change. Locally, the mass balance varied between different types of glaciers with proglacial lakes associated with the most rapid loss, particularly during 2000–2020. The mass balance of debris-covered glaciers was similar to debris-free glaciers during 2000–2020.

4. Topographic setting influenced retreat and thinning, with loss rates decreasing with increasing elevation. The rate of both glacier retreat and thinning decreased with elevation, but the relationship between the parameters of slope and aspect with thinning rates differed from their relationship with retreat rates, spatially and through time. For slopes of 8–40^∘ (which includes most glaciers), for example, the retreat rate increased with slope, while the thinning rate decreased.

In this study, we observed accelerated glacier loss in the WNR on multi-year timescales. However, factors such as precipitation, temperature and altitude could not yet fully explain the heterogeneity of glacier changes. Thus, more detailed data and glacier ablation models are needed to fully understand the mechanism of glacier change in the future.