Accelerated glacier mass loss in the largest river and lake source regions of the Tibetan Plateau and its links with local water balance over 1976–2017

Abstract

Variations in glacier meltwater in the source regions of the Tibetan Plateau's (TP) largest lake (Selin Co) and China's longest river (Yangtze River) regulate the local and downstream water balance under the warming climate. However, the magnitude of their variations over the past four decades is still unknown. Here, we examine long-term changes in glacier mass balance over 1976–2017 using KH-9 and CoSSC-TanDEM-X data. We find that the mean rate of glacier mass loss (GML) has accelerated from −0.21 ± 0.11 m a^–1 over 1976–2000, to −0.28 ± 0.14 m a^–1 over 2000–11, and subsequently to −0.48 ± 0.10 m a^–1 over 2011–17. Changes in temperature and precipitation are the major causes of GML. Over 1976–2017, the contribution of decadal GML to Tuotuohe sub-basin runoff ranges from 4.3 to 8.0%, while its contribution to increased lake volume change in Selin Co and Chibzhang Co-Dorsoidong Co ranges from 3.5 to 16.3% and 19.2 to 21.4%, respectively. The GML of source regions made relatively small contributions to river runoff and lake supply, but plays a vital role when precipitation decreases. The quantitative evaluation of the water balance for the sources of great rivers and lakes over the TP is therefore important for water resource management and hydrological cycle studies.

1. Introduction

The Tibetan Plateau (TP) (mean elevation ~4000 m a.s.l.) has the greatest abundance of glaciers outside the polar regions, with a total area of ~10 × 10⁴ km² (RGI Consortium, 2017) and a total ice volume of ~7 × 10³ km³ (Farinotti and others, 2019). Hence, together with Hindu Kush-Karakoram Himalayan region, they are often referred to as the Third Pole of the Earth (Yao and others, 2012a). Most large Asian rivers originate in this region and so it is also known as the water towers of Asia. Over the past few decades, the climate in the TP has become wetter and warmer (Li and others, 2018a), resulting in glacier mass loss (GML) (Yao and others, 2012b; Brun and others, 2017; Zhou and others, 2018), lake expansion (Lei and others, 2014; Zhang and others, 2019a) and an increase in river runoff (Lutz and others, 2014). These changes are altering the balance of the cryosphere and have significant impacts on the local hydro-climate-ecological cycle, and water supplies (Immerzeel and others, 2010; Lutz and others, 2014; Milner and others, 2017).

The Geladandong Mountains are located in the Tanggula mountain range in the central TP, an area influenced by both continental climate and the Indian summer monsoon (Tian and others, 2001; Yao and others, 2013). Glaciers are widely distributed throughout this region due to the low temperature and sufficient precipitation (Ding and others, 1992). They feed the Yangtze River, and the two biggest lakes on the TP, Selin Co and Chibzhang Co-Dorsoidong Co (CC-DC) (Fig. 1a).

Fig. 1. Location of study region and data coverage. (a) Location of study area and three sub regions: CC-DC (labeled A), Selin Co (labeled B) and Yangtze River basin (labeled C). The glacier outlines are derived from the Randolph Glacier Inventory (RGI 6.0). Measured glaciers are all located within the dashed rectangular area zoomed in (b). The left-upper inset indicates the location of the study area in the TP. (b) The spatial and temporal coverage of KH-9 and TanDEM-X CoSSC pair datasets used. The background is the orthorectified KH-9 image on 7 January 1976. The TanDEM-X data, marked in blue and black, are used to calculate the glacier elevation change in the periods 2000–11 and 2000–17, respectively. It should be noted that only one TanDEM-X scene on 7 April 2014, covering the eastern part of the study area, is adopted due to the lack of available data. The TanDEM-X data marked in red are used to calculate the seasonal glacier changes.

Due to the TP's extreme environment, almost no glaciological field measurements have been carried out, but several studies using Landsat images have suggested an accelerated shrinkage of glaciers in this region from 1969 to 2011 (Ye and others, 2006; Zhang and others, 2013b, 2013a). Geodetic glacier mass balances have also been monitored, as shown in Table 1. However, there are three main issues surrounding these studies. The first is due to the large uncertainties in the data or methods used. ICESat, with its sparse tracks, only covered a limited area of the glaciers and a limited time period (2003–2009), resulting in large discrepancies between studies (Neckel and others, 2014; Chao and others, 2017). Generated DEMs from optical images suffer from large gaps in accumulation areas such as Ziyuan-3 (Liu and others, 2020). Second, spatial and temporal data coverage is inconsistent, hindering comprehension of the glaciers' long-term response to environmental change (Table 1). Finally, synthetic aperture radar (SAR) signal penetration differences and the neglected seasonal glacier elevation change in these studies could have large effects on the results, especially over short time periods of <10 years (Liu and others, 2020).

Table 1. Glacier mass balance in the Geladandong region using different methods from previous studies

Location or covered area is shown in Figure 1a.

Over the period of 1976–2016, Selin Co expanded in area by 41% and CC-DC by 23%. Glaciers play a minor, but important role in this expansion (Tong and others, 2016; Zhang and others, 2017a). Hydrological modeling reveals that the combined contribution of snow and glacier melt runoff into Selin Co was 22.8% between 2003 and 2013 (Zhou and others, 2015), while 11.5% between 1981 and 2012 (Ding and others, 2018). A degree-day glacier-melt model quantifies the contribution of seasonal glacier meltwater to the lake as 5% over the period 1979–2013 (Tong and others, 2016). A quantitative assessment of lake water balance, also conducted for CC-DC (2003–14), but using glacier mass-balance data reconstructed from MODIS albedo products for 2000–16 (Qiao and others, 2019), suggests that the contribution of glacial meltwater to lake expansion is 19.3 ± 4.5%.

Increasing glacier meltwater has also led to an increase in river runoff across the TP (Lutz and others, 2014). In the source region of the Yangtze River, the meltwater from glaciers has resulted in an increased river discharge (Zhang and others, 2007a; Chao and others, 2017). Liu and others (2009) modeled the runoff and reported a contribution of ~11% (1961–2000) from glacier meltwater upstream of Zhimenda hydrological station, compared with a ~17% contribution (2016–18) determined from a study of environmental isotopes (Li and others, 2020) and ~17.5% (1986–2009) from mass-balance data using an area-scaling method (Yao and others, 2014). Recently, accelerated mass loss from the adjacent Puruogangri Glacier (shown in Fig. 1a) over 2011–16 has been reported by Liu and others (2019). However, it is not known if such accelerated mass loss occurred in a similar way in our study area over the past four decades. Given these differences and uncertainties in current estimates of glacier mass change, a consistent and long-term estimate of glacier mass balance from the Geladandong Mountains will help to improve our understanding of the response of glaciers to climate change and its hydrological consequences.

Here, we investigate three periods of glacier mass change from 1976 to 2017 using the KH-9 and TanDEM-X CoSSC datasets. The aim is to determine how the glaciers have changed on regional, basin and individual glacier scales in response to the changing climate. The influences of the glacier mass changes on surrounding lakes and downstream river runoff on the basin scale are quantitatively assessed in the corresponding periods.

2. Study area

The Geladandong Mountains region (33°03′–33°37′ N, 90°43′E–91°40′ E) is the source of the Yangtze River and some large lakes on the TP. Geladandong Peak, with an elevation of ~6621 m, is the highest mountain in the range. The regional topography varies from high in the northwest to low in the southeast. The climate is cold and dry, with mean annual precipitation and temperature, observed at the nearby weather stations, Tuotuohe and Anduo, of 381 mm and −3.1°C, respectively, over the period of 1980–2017. The precipitation is mainly concentrated in the summer months of June to September, with this period accounting for ~80% of the annual amount. Geladandong Mountains region comprises the largest concentration area of glaciers in the Tanggula Mountains. There are ~252 glaciers in the region, covering an area of 924 km², according to the RGI 6.0 (RGI Consortium, 2017). Thereinto, debris-covered glaciers account for 2% in total area (Herreid and Pellicciotti, 2020).

As shown in Figure 1a, there are three main lake/river basins around the glaciers: CC-DC Basin (A); Selin Co Basin (B) and the Yangtze Basin (C). Both Basins A and B are endorheic, and have areas of 1.35 × 10⁴ and 2.91 × 10⁴ km², respectively. Glaciers account for 1.79% of the area of A and 0.47% of B. These two basins contain two large lakes: CC-DC lake, which had an area of ~1077 km² in 2017; and Selin Co, the largest lake in Tibet, which had an area of 2391 km² in 2017. Basin C is exorheic, with four rivers converging into one. Tuotuohe basin is regarded as the source of Yangtze River (Zhang and others 2007a). Xiaodongkemadi and Qiangtang No. 1 glaciers, close to our study area, have in situ measurements of glacier mass balance available.

3. Data and methodology

3.1 Glacier and lake mapping

Cloud-free Landsat MSS images (Table 2) were used to map the glacier boundaries in 1976 in our study area by using the ISODATA-based unsupervised classification method in ENVI software. The boundary obtained was then manually revised according to the original Landsat imagery. The RGI 6.0 Glacier inventory was adopted to define the basin boundary at which to split the glacier into different parts. Lake outlines between 1976 and 2017 were derived from Landsat MSS/TM/ETM+/OLI data (https://glovis.usgs.gov). Cloud-free images from climatically relatively stable seasons (around October) are chosen to avoid potential seasonal variability (Zhang and others, 2019b). An automatic water classification algorithm, based on global–local threshold segmentation for normalized difference water index grayscale image, was used to distinguish surface water from non-water features (Zhang and others, 2017b). The final lake dataset was obtained by visually editing and excluding other water bodies (such as rivers and small streams) from the automatically extracted water bodies. The lake polygons used in this study are adopted from Zhang and others (2019a).

Table 2. Datasets used to calculate glacier surface elevation change

The spatial coverage of the KH-9 and TanDEM-X datasets for the corresponding dates is shown in Figure 1b.

3.2 Decadal glacier surface elevation change

3.2.1 DEMs generation and differencing

The TerraSAR-X (TSX) satellite, launched in 2007, and the TanDEM-X (TerraSAR-X add-on for Digital Elevation Measurement, TDX), launched in 2010, are coupled into one satellite constellation with a slanting incline to obtain high-resolution terrain measurements of the Earth's surface. Their spatial coverage is 30 × 50 km², with a ground resolution of 1.7–3.5 m. According to previous studies (Neckel and others, 2014; Liu and others, 2019), TDX DEM was generated using the InSAR TanDEM-X bistatic DEM Workflow in ENVI. This workflow mainly consists of interferogram generation, adaptive filter, coherence generation, phase unwrapping, phase to height conversion and geocoding. The TanDEM-X CoSSC datasets were used to generate DEMs for 2011 and 2017 (Table 2). A pair of images from 2014 (Fig. 1b) was selected to cover the eastern part of our study area because of limited data availability, which covers 8% of glacier area only. The DEMs produced for 2011 and 2017 are from different tracks and do not have a good spatial overlap (Fig. 1b), so we calculated the elevation differences between 2000 and 2011 and between 2000 and 2017, respectively; the elevation change between 2011 and 2017 is then calculated by subtracting the elevation change over 2000–11 from that over 2000–17.

A benefit of the declassification of the Hexagon KH-9 mission and its overlap of more than 70%, is that DEMs can be produced back to the 1970s, from which glacier mass changes can be estimated (e.g. Zhou and others, 2018). KH-9 has a high spatial resolution of 6–9 m. One pair of images covering our study area, acquired on 7 January 1976, is available. The KH-9 images (accessed from https://earthexplorer.usgs.gov) have been scanned from archived films by NASA, and could be distorted (Surazakov and Aizen, 2010). A total of 1081 fiducial marks in the KH-9 image automatically recognized by the MATLAB code were utilized to correct any such distortion (Surazakov and Aizen, 2010). The nonmetric camera model was then selected, due to the lack of fiducial information. Twenty-three ground control points, such as mountain ridge lines, and river and road intersections from Google Earth, were used for exterior orientation. Combining these with the other 89 tie points generated by the ERDAS LPS (Figure S1), enabled the triangulation calculation to be performed. The RMSE of total image unit-weight is 0.61 pixel. The generated DTM was then resampled to a resolution of 30 m to match the SRTM-C DTM. The elevation difference between different DTMs was calculated. The orbit-induced systematic spatial trend, elevation-, aspect-, slope- and curvature-dependent differences in non-glacier regions were used to calculate the correction value (Nuth and Kääb, 2011; Li and others, 2017).

The C-band and X-band SAR aboard the space shuttle Endeavor in the Shuttle Radar Topography mission (SRTM) (11–22 February 2000) simultaneously acquired the global digital terrain model. The C-band SRTM generated by NASA (https://earthexplorer.usgs.gov) covered over 80% of the Earth's land surface between 60° N and 56° S. We use these data as a reference in the decadal glacier surface change process. The X-band SRTM generated by DLR (https://download.geoservice.dlr.de/SRTM_XSAR) is only available along ~50 km wide orbital ground traces. Therefore, the available dataset with a grid-like coverage, only covering southern part of study area, was used to determine the SAR penetration value by DEM differencing.

3.2.2 Outlier removal

For KH-9, brightness saturation in accumulation areas may cause poor texture in optical images. For SAR images, although TanDEM-X CoSSC performs better than KH-9 in the accumulation areas, layover and shadow can introduce uncertainty in elevation estimates. These issues can result in outliers in elevation differences for both SAR and KH-9 data. The maximal absolute magnitude of glacier thinning is calculated according to assumed data distribution in each 50-m elevation bin (Pieczonka and Bolch, 2015; Zhou and others, 2017). Here, an elevation difference beyond the range of this magnitude is defined as an outlier. Then for each 50 m elevation bin, a 1.5 normalized median absolute deviation filter is used to remove such outliers (Brun and others, 2017). All data gaps are not interpolated, as such a procedure may result in larger uncertainties. The elevation changes of pixels in the void area are assumed to have the same mean value in the same 50 m altitude bin (Gardelle and others, 2012). A glacier density of 850 ± 60 kg m⁻³ is used to convert glacier volume change to mass balance (Huss, 2013).

3.3 Seasonal change and penetration depth corrections

The seasonal glacier mass balance was examined using eight pairs of TanDEM-X CoSSC data from 2016 and 2017. In general, seasonal changes in glacier surface elevation are small compared with their annual changes. Hence, small seasonal errors could lead to significantly abnormal results. To maximize the coverage of the data used, TanDEM-X CoSSC data pairs acquired on 8 January 2017 were selected to generate the reference DEM with a resolution of 5 m. This high-resolution TDX DEM was then used for the geocoding process of other images to avoid the potential misalignment among images (Fig. 2). In the following error removal process (Li and others, 2017), almost no shifts are found among these generated TDX DEMs. All mean errors of the non-glacier area are close to zero (Table S1).

Fig. 2. Flowchart illustrating the methods used to derive seasonal glacier surface elevation change.

The TDX data for estimating the decadal glacier surface elevation changes are obtained from January to October (Table 2). To remove the seasonal signal, we need the correction values of months (Table 2) corresponding to the acquisition date of the base SRTM data. As the months of the data used for decadal and seasonal glacier surface elevation changes do not coincide (Table 2), a linear relationship between two adjacent dates is assumed (Lin and others, 2017a, 2017b), and the correction value for the corresponding date is estimated from this linear regression.

The penetration depth of X- and C-bands varies for different types of glaciers. Temperature, snow depth, water content and debris cover can all affect penetration depth (Rignot and others, 2001; Gardelle and others, 2012), therefore it is not reasonable to ascertain a general penetration value of C- or X-band for all glaciers. The penetration depth of X-band is usually smaller than that of C-band. SRTM-X (DLR), acquired in the same period as SRTM-C (NASA), offers one way to estimate the approximate C-band penetration over glaciers (Gardelle and others, 2012; Gardelle and others, 2013). The difference between X- and C-band SRTM was used to estimate the difference between C-band and X-band signal penetration (Gardelle and others, 2012). Penetration depths of 2.1–4.7 m have been measured for the radar signal at 10 GHz in the Antarctic (Davis and Poznyak, 1993). Here, to describe the influence of the penetration depth of X-band, an uncertainty of 150% of the difference between the SRTM-C and SRTM-X elevation measurements was added.

3.4 Lake volume change

There are few studies on the lake levels of CC-DC and Selin Co using altimetry satellites prior to the 1990s (Crétaux and others, 2016), therefore the water levels of these lakes in 1976 were reconstructed by using a regression between lake area and lake elevation (Zhang and others, 2017a). The water levels from satellite-altimetry data are obtained from the Hydroweb (http://hydroweb.theia-land.fr/?lang=en). The lake volume changes are estimated according to Taube and others (2000). An additional complication is that the twin lakes, DC and CC, only became connected in June 2005 (Jiang and others, 2017). For a better estimate of the contribution of GML to lake volume increase, the lake volume change over the period 1976–2000 is based on CC. After June 2005, until 2007, CC started to discharge into DC (Jiang and others, 2017). Therefore, the water volume changes for CC and DC are divided into three stages. The first period is the water volume change of CC between 2000 and 2005. The second period, 2005–07, covers the DC volume change, and we assume that this change in volume is supplied by CC. The third period is 2007–11 when CC and DC change in volume synchronously, so we calculated a single volume change for the two lakes together.

3.5 Other data

Monthly and yearly climate data from the National Meteorological Information Center (CMA Meteorological Data Center) (http://data.cma.cn/) are used. Anduo and Tuotuohe stations are the two nearest stations to our study area. Anduo station is located at 32.35° N, 91.10° E with an elevation of 4800 m, while Tuotuohe station is located at 34.22° N, 92.43° E with an elevation of 4533 m. Annual discharge of Tuotuohe runoff data from 1976 to 2012 is adopted from Chao and others (2017).

3.6 Uncertainty analysis

The uncertainty of glacier surface elevation change includes mean elevation change error E _Δh in each elevation bin, error in the off-glacier area with slope <15° (E _{Δoff-glacier}), seasonal error (E _Δs) and penetration depth error (E _ΔP). E _Δh is calculated according to Gardelle and others (2013):

(1)$$E_{\Delta h} = \displaystyle{{E_{\Delta hi}} \over {\sqrt {N_{{\rm eff}}} }}{\rm , \;\ }\quad N_{{\rm eff}} = \displaystyle{{N_{{\rm tot}}PS} \over {2d}}$$

where E _Δhi is the std dev. and N _eff is the number of independent measurements in each elevation bin. N _tot is the total number of measurable pixels. PS is the pixel size (30 m), d is the spatial autocorrelation distance, determined by covariance modeling in the ArcGIS as 262 ± 34 m for KH-9 and 635 ± 55 m for TDX data. The initial uncertainty of mass balance (E _ΔM0) introduced by E _Δh is estimated according to Eqn (2) (Zhou and others, 2017):

(2)$${E_{\Delta M}}_0 = \sqrt {{{\left( {{E_{\Delta V}}{{{\rho _{\rm{i}}}} \over {{\rho _{\rm{w}}}A}}} \right)}^2} + {{\left( {{E_{\Delta {\rho _{\rm{i}}}}}{V \over {{\rho _{\rm{w}}}A}}} \right)}^2} + {{\left( {{E_{\Delta A}}{{{\rho _{\rm{i}}}V} \over {{\rho _{\rm{w}}}{A^2}}}} \right)}^2}} $$

where V, A, ρ _i, and ρ _w are glacier volume change, glacier area, ice density (850 ± 60 kg m⁻³) and water density (1.0 × 10³ kg m⁻³), respectively; E _ΔV, $E_{\Delta \rho _{\rm i}}$, and E _ΔA are the uncertainties of glacier volume change, glacier density and glacier area, respectively. E _ΔV = E _Δh × A, E _ΔA = E _ΔA = (1/2)resolution × perimeter. Ideally, the mean value of the off-glacier area should be close to 0. Standard value E _off-glacier is also included in uncertainty analysis. The seasonal error E _△s is calculated using ${\rm } \scale80% \sqrt {E_{\Delta h}^2 + E_{\Delta {\rm off}\hbox{-}{\rm glacier}}^2 }$. The uncertainty in penetration depth difference between SRTM-C and SRTM-X is 0.97 m. Based on this value, the uncertainty of the C-band penetration is 1.5 m, which is 150% of 0.97 m. Taking all these uncertainties into account, the final uncertainty in the mass balance (E _ΔM) is given by Eqn (3):

(3)$$E_{\Delta M} = \sqrt {E_{\Delta M_0}^2 + E_{\Delta P}^2 + E_{\Delta S}^2 + E_{\Delta {\rm off}\hbox{-}{\rm glacier}}^2 } $$

4. Results

4.1 Seasonal change and penetration depth correction

The glaciers underwent ablation from July to January of the following year, but accumulated mass from January to July. Therefore, the glaciers in this region are characterized as the summer accumulation type. Between 8 January 2017 and 30 January 2017, the glacier surface elevation presented a slight positive increase of 0.08 m. From 8 January 2017 to 17 April 2017, the glaciers had a net accumulation of 0.24 m, which is similar to the value during the period 17 April 2017 to the 9 May 2017 and indicates a mass-balance gain in this region. During these two periods, the accumulation occurs in the lowest part of the glaciers and shows little variation in each elevation bin (Fig. 4). A more obvious accumulation started in June and July of 2017. This is especially evident from 11 June 2017 to 3 July 2017, when the surface elevation increased by 2.51 m, accounting for 77% of the preceding total accumulation (3.26 m). These accumulations mainly occur in the upper area and there is only slight melting in the terminus (Figs 3f, 4a). From 3 July 2017 to 7 September 2017, the Geladandong glaciers had a thinning of −1.53 m (Figs 3a, 3g, 4). An abnormal ablation of 2.15 m occurred between 9 September 2016 and 8 January 2017. This observation may be connected to the differences in SAR penetration and is discussed further in Section 5. Overall, the glaciers show a negative mass balance, with a surface elevation lowering of −0.42 m in the 2016/2017 mass-balance year. Seasonal glacier change results are used for the seasonal corrections in decadal glacier elevation change. The correction values on the dates of 8 and 30 January, 8 March, 7 and 28 April, 13 May, 4 June and 3 October were 0.08, 0, −0.16, −0.28, −0.37, −0.59, −0.80 and −1.70 m, respectively.

Fig. 3. Seasonal glacier surface elevation changes in 2016/17 referenced to 8 January 2017. Data from two orbits are utilized to acquire monthly change. The outlines in red in (h) are the coverage of (a, c, e, g). The overlapping region of red and blue outlines in (h) are the coverage of (b, d, f). The dates are in the format YYYYMMDD.

Fig. 4. (a) Changes in glacier surface elevation in each elevation bin. Coverage of data over each period is shown in Figure 3. The histogram at the bottom shows the glacier area (covered by outlines in red in Fig. 3h) percentage in each 50-m elevation bin. (b) Time series of mean glacier surface elevation change over eight time periods.

Figure S2 shows the differences between SRTM-X and SRTM-C in glacier and off-glacier regions after error removal. A penetration correction of 1.59 ± 0.97 m is adopted over 2000–11 and 2011–17, which is comparable with the value of 2.1 m used by Liu and others (2017a, 2017b). Many previous studies have adopted different penetration depths for the C-band. For example, a value of ~3 m is assumed by Chao and others (2017) and Xu and others (2018). Chen and others (2018) used a value of 2.5 ± 0.5 m in East Nepal, according to Kaab and others (2012). Zhou and others (2018) used a 1.9 m penetration depth for the C-band, but did not consider the penetration depth of the X-band. However, little is known about the penetration of C- or X-band in this study area. To describe the effect of X-band penetration depth, we added a 150% uncertainty of the 1.59 ± 0.97 m, and assumed the penetration of the C-band to be 2.39 ± 1.46 m for correction over 1976–2000, compared with 2.10 ± 4.55 m estimated from ICESat data (Chen and others, 2018).

4.2 Decadal glacier elevation change

GML accelerated in the Geladandong region from 1976 to 2017, especially from 2011 to 2017 (Fig. 5). The elevation change rate of the off-glacier area is shown in Figure S3. Mean elevation change rates in the three periods in the off-glacier area are all close to zero (Table S1). The glaciers thinned with mean rates of −0.21 ± 0.11 m a^–1 over 1976–2000 and −0.28 ± 0.14 m a^–1 over 2000–11. During 2011–17, the thinning rate almost doubled from the 1976–2000 value to −0.48 ± 0.10 m a^–1, similar to the rate in 2016/17 (−0.42 m a^–1) (Fig. 4b).

Fig. 5. Region- and basin-scale changes in glacier surface elevation from 1976 to 2017: (a) 1976–2000, (b) 2000–11 and (c) 2000–17. The rates of glacier elevation change for the whole region are labeled in black at the top of each figure. The rates of glacier elevation change in the three subregions are labeled in dark blue text.

Figure 1b shows that there are three major continuous ice bodies, i.e. West, Middle and East, distributed along the longitude direction. The area-averaged change rates of these three parts in 1976–2000, 2000–11 and 2011–17, were −0.21, −0.19 and −0.25 m a^–1 (West), −0.14, −0.28 and −0.60 m a^–1 (Middle) and −0.61, −0.49 and −0.09 m a^–1 (East), respectively. GML was greater in the East than that in the West, except for the East part in 2011–17 (Fig. 6c). The mass loss rate of glaciers in the Middle part significantly increases along the longitude direction (Fig. 6). This trend also appears in the altitude direction (Figure S4), but not in the latitude direction

Fig. 6. Changes of glacier surface elevation for different types with longitude: (a) 1976–2000, (b) 2000–11 and (c) 2011–17. The red box indicates an outlier area, where the 2014 TanDEM-X images shown in Figure 1b were used. The size of the circles denotes the glacier area. A 1.5 std dev. range is adopted to denote the error range of each glacier.

Differences in basin-scale thinning rates are apparent (Fig. 5). The thinning rate of glaciers in Basins A and B accelerated during 1976–2017. Specifically, the glaciers in Basin A had a mean rate of −0.06 ± 0.10 m a^–1 over 1976–2000, which increased to −0.27 ± 0.13 m a^–1 over 2000–11, followed by a higher rate of −0.68 ± 0.10 m a^–1 over 2011–17. The mass loss accelerated sharply over 2011–17 compared with the two preceding periods. The glaciers in Basin B present a similar temporal ablation pattern as Basin A. The GML rate in Basin C was relatively stable in 1976–2011. The glaciers had a mean thinning rate of −0.26 ± 0.10 m a^–1 over 1976–2000, −0.25 ± 0.10 m a^–1 over 2000–11, followed by an increase to −0.48 ± 0.12 m a^–1 over 2011–17.

Glaciers of different sizes exhibit different behaviors (Fig. 7). Glaciers of size 0–1 km² account for 53.1% of the total number of glaciers, but only 5.7% of the total glacier area (Fig. 7d). These glaciers fluctuate with no consistent temporal change characteristics, perhaps because of uncertainties in delineating the glacier outlines and the fact that DEM differences could be greater than the magnitude of mass changes for glaciers in the size of 0–1 km² (Maurer and others, 2019). The mean change rates of glaciers in the category of 1–10 km² are consistent, with an initial decrease from the period of 1976–2000 to 2000–11, followed by an acceleration in 2011–17. Glaciers larger than 10 km² account for a small percentage of the total number of glaciers, but make up a large fraction of the area. These glaciers present an obvious accelerating loss, consistent with the mass change of the whole study region.

Fig. 7. Rates of glacier surface elevation changes with sizes in 1976–2000, 2000–11 and 2011–17 in Basins A (a), B (b), C (c) and whole study area (d). The coverage of the three basins is shown in Figure 1. The Black circles denote the area averaged glacier elevation change rate. The three columns in each size category denote the periods 1976–2000, 2000–11 and 2011–17.

4.3 Contribution of glacier excess meltwater to lake volume and river discharge

The areas of lake CC-DC and Selin Co expanded by large amounts over 1976–2017 (Figure S5). Lake CC, with its glacier meltwater supply, expanded in area by 5.6% (27 km²) and had a water storage gain of ~1.1 Gt between 1976 and 2000. During the same period, DC expanded in area by 3.7% (14 km²). After 2000, both CC and DC underwent rapid expansion; they connected in 2005, and have changed synchronously since 2007. During the period 2000–11, Lake CC had an area expansion of 2.4%, significantly smaller than the 24.5% expansion of DC. In total, CC-DC expanded by 14.4% and gained water mass of 5.3 ± 0.8 Gt over 2000–11. Over the period 2011–17, CC-DC had an area increase of 8.4% and water mass gain of ~2.4 Gt. Overall, CC and DC were relatively stable before 2000, but expanded rapidly in 2000–11. Their expansion slowed over 2011–17. The area of Selin Co expanded by 12.5% over 1976–2000, 17.0% over 2000–11 and 7.3% over 2011–17. The lake level change over 1976–2000, constructed using the area-elevation relationship, is 4.63 m, compared to an observed change of 4.55 m over the period 1979–2000 (Tong and others, 2016). Accordingly, the lake had water gains of 8.3 ± 0.5, 14.0 ± 0.1 and 0.5 ± 0.1 Gt in the three corresponding periods.

Table 3 presents estimates of the contribution of GML to lake volume and river discharge changes. Here, we assumed that the glacier excess meltwater is totally discharged into the lakes/rivers without evaporation or other loss as with many other studies such as Zhou and others (2019). In the three basins, glacier excess meltwater plays a minor, but important, role in the water balance and varies over time. The relative contribution of glaciers becomes greater as the precipitation decreases. For Basin A, the total lake volume increased by 5.4 ± 0.1 Gt over 2000–17 compared with 4.4 Gt over 2003–14 (Qiao and others, 2019). CC gains water at the rate of 0.06 Gt a^–1, while the glaciers in this basin lost mass at a rate of −0.01 Gt a^–1, contributing 20.0% to the lake water gain over 2000–17. Zhang and others (2017b) reveal that increasing precipitation was the main driver of lake expansion (74%) in the inner TP. In our study area, the precipitation (taken as an average of observations from the Tuotuohe and Anduo stations) increased by 21% over 2000–11 and decreased by 32% over 2011–17. In particular, over 2000–11, the glacier loss rate accelerated, and the relative contribution also decreased. This observation also supports the dominant role of precipitation. Over 2000–17, the contribution of glacier excess meltwater to lake volume change was 20.2%, which should be slightly larger as the Puruogangri Glacier (20% of total glacier area in CC-DC basin) is excluded. Using mass-balance data from 2000 to 2016 from Zhang and others (2018), Qiao and others (2019) suggest that the contribution over 2003–14 is 19.3%.

Table 3. Changes in lake volume and river discharge corresponding to glacier mass balance over three periods

LVC is lake volume change.

Compared to Basin A, the proportion of glaciers area relative to basin area in Basin B is smaller. The glacier excess meltwater contributing to the lake volume change is smaller, but it still exhibits a change pattern similar to that in Basin A. The relative contribution changed counter to the precipitation change. However, considering the effect of the long distance between Selin Co and Geladandong mountains, the glacier runoff will be consumed partly by evaporation and infiltration, the contribution of glacier excess meltwater to the lake increase should in fact be smaller.

Unlike Basins A and B, Basin C is exorheic. Only discharge from the Tuotuohe sub-basin has been measured since 1961. Therefore, the excess meltwater from glaciers located in Tuotuohe is estimated (Table 3). The GML, relative to the river discharge change, was <10%, falling to 4.3% over the period 2000–11. Glacier meltwater plays a modest role in this region (Immerzeel and others, 2010), making a much smaller contribution than that from seasonal glacier melt (Zhang and others, 2007a; Liu and others, 2009; Li and others, 2020).

5. Discussion

5.1 Impact of penetration depth

SAR penetration depth is an important source of errors in measuring seasonal glacier variations. Strong ablation occurred in the accumulation area and increased with elevation from September 2016 to Jan 2017 (Fig. 4). Sublimation, the dominant ablation mechanism in very cold environments, is one cause, but it is only of the order of centimeters per year (Cuffey and Paterson, 2010). Differences of penetration depth in snow, firn and ice in September and January probably contribute to this anomalous ablation (Rignot and others, 2001). Temperature and snow wetness can both influence the penetration depth. Temperature difference between September and January can reach more than 20°C inferred from Figure 8. Gay and Laurent (2016) reported that the difference for X-band for −20 and −5°C in dry snow can reach 2.9 m due to the absorption loss, with low temperature snow having a deeper penetration depth (Gusmeroli and others, 2017). In addition, penetration depth decreases as snow wetness increases (Arslan and others, 2005). With higher above zero temperatures, snow wetness values will be larger in September than in January (Fig. 8). Additionally, drier air in January than in September (Fig. 8), could also potentially affect the snow wetness, since, if the air is more humid, more heat is transferred to the ice from the air, causing the snow to melt. The increased grain size of snow with time can also increase the penetration depth (Cuffey and Paterson, 2010; Gay and Laurent, 2016).

Fig. 8. Seasonal temperature and precipitation at the Anduo and Tuotuohe weather stations. (a, b) Monthly mean and maximum temperature and precipitation from 1980 to 2017. (c, d) Monthly mean and maximum temperature and precipitation from January 2016 to September 2017. Note that several months of temperature data from the Anduo station are not available.

The effect of variations in penetration depth may also exist in other periods, when snow conditions changed under changing precipitation and increased temperature. Therefore, the seasonal elevation change detected here is a mixed process of mass balance, snow compaction and TDX radar penetration signal. The difference between 8 and 30 January 2017 is close to zero, illustrating the feasibility of the method. If we were to assume, following Liu and others (2019), that there is no snowfall from January to May (as snow/ice meltwater in May can greatly reduce the penetration), the X-band penetration would be 0.48 m, similar to the value of 0.61 m found by Liu and others (2019). The abnormal accumulation on 3 January 2017 would then be reduced to 2.78 m. However, we could not follow that approach, because the snowfall, detected by MODIS in the spring, changed the condition of the glacier surface (Figure S6). The final seasonal correction values we used will not have much influence on our decadal glacier change result. This seasonal correction value is a result of mixed process, it could also compensate the effect of neglected difference in X-band penetration depth and snow compaction in different months in estimating decadal glacier change.

5.2 Comparison with other studies

The seasonal glacier mass balance has not previously been determined in this region. The accumulation and ablation periods detected in our study (Fig. 3) are supported by observations on Xiaodongkemadi Glacier (XDG) (Seko and others, 1994) and Qiangtang No. 1 (Li and others, 2018a, 2018b) (locations shown in Fig. 1a). For these two glaciers, major melting occurs in July and August, with accumulation occurring from January to June.

Chemical analysis of ice core records provides evidence that the glacier mass balances in Geladandong have been negative since the 1980s (Zhang and others, 2007b; Kang and others, 2015). This finding is also supported by observations of glacier area change (Ye and others, 2006). In particular, glaciers retreated continuously from 1973 to 2002 (Ye and others, 2006), which is consistent with our glacier elevation change results. The mass balance of Xiaodongkemadi was −0.36 m w.e. a^–1 over 2000–10 (Yao and others, 2012b), while the mass balance of Qiangtang No. 1 around the equilibrium-line altitude was −0.36 m w.e. a^–1 over 2012–16 (Li and others, 2018a, 2018b). Here, we found the glacier mass balance was −0.24 ± 0.12 m w.e. a^–1 over 2000–11 and −0.41 ± 0.10 m w.e. a^–1 over 2011–17.

Heterogeneous behavior of glacier surface elevation change can appear due to discontinuities in time series and spatial coverage (Neckel and others, 2014; Chao and others, 2017; Chen and others, 2018) (Table 1). However, our results are consistent with the two studies covering most parts of our study region. The elevation change rate is −0.22 ± 0.12 m w.e. a^–1 by KH-9 over 1976–2000 (Zhou and others, 2018), and −0.27 m w.e. a^–1 over 2000–16 (Brun and others, 2017), comparable with our results of −0.21 ± 0.11 over 1976–2000 and −0.35 ± 0.09 m w.e. a^–1 over 2000–17. In addition, the accelerated mass loss trends detected in our study are also confirmed by nearby glacier mass-balance observations (Yao and others, 2012b) and DEM-differencing (Chen and others, 2018; Liu and others, 2019) (Table 1).

5.3 Factors influencing glacier elevation change

5.3.1 Glacier size

Small glaciers are more sensitive to climate variations (Bahr and others, 1998; Huss and Fischer, 2016). They usually respond to climate change on shorter timescales as small glaciers have less area to preserve enough precipitation (Li and others, 2017). Here, the glaciers in the 1–10 km² size range show the same change trend in three separate basins, concurrent with the variation of precipitation (Fig. 9) (r = 0.68 and r = 0.99 in Anduo and Tuotuohe station, respectively). The mean annual precipitation at Anduo and Tuotuohe stations increased by 10 and 40% between 1980–2000 and 2000–2011, and then decreased by 38 and 24% between 2000–11 and 2011–17, respectively. The mass loss rates of glaciers in the 1–10 km² category in the three basins all decreased over 2000–11, compared with the rates over 1976–2000, and then increased over 2011–17.

Fig. 9. Long-term temperature and precipitation at the Anduo and Tuotuohe weather stations. (a) Annual summer temperature (June–August). (b) Annual precipitation. (c) Annual summer precipitation. The numbers at the top of (b) and (c) are mean precipitation of the two stations during the three periods.

The mean change rates of the glaciers in the 10–100 km² category (62% of the total glacier area) (Fig. 7d), show the same trend as the basin-scale glacier change rate (Fig. 5). In contrast to the glaciers in the 1–10 km² range, the glacier change rates in the 10–100 km² range, in all three basins, accelerated concurrently with the variation of summer temperature (r = 0.99 for Anduo station and r = 0.98 for Tuotuohe station) (Figs 6, 9a). The average summer temperatures recorded at Anduo and Tuotuohe stations increased by 0.49 and 0.29°C between 1980–2000 and 2000–11, and 0.62 and 0.40°C between 2000–11 and 2011–2017, respectively (Fig. 9a). Melt rates increase rapidly as air temperatures rise past the melting point (Cuffey and Paterson, 2010). Also, the increased summer temperature results in an increase of the upper limit of precipitation, a reduction in the accumulation of snow on glaciers and the albedo, and faster glacier melting (Gardelle and others, 2013). Although under the same warming background, the glacier elevation change rate in Basin C only shifted slightly from −0.26 to −0.25 m a^–1 between 1976–2000 and 2000–11, not as significantly as the changes in Basin A (−0.06 to −0.27 m a^–1) or Basin B (−0.23 to −0.42 m a^–1) (Fig. 5). Braithwaite and Zhang (2000) suggest that a 25% increase in precipitation can partly offset a 1°C temperature increase. Comparing the increase of precipitation (40% in Basin C versus 10% in Basin B between 1980–2000 and 2000–11), it can be seen that the substantially increased precipitation in Basin C may moderate the effect of increased temperature (0.29°C) (Zhu and others, 2017; Liang and others, 2018; Liang and others, 2019), preventing the glaciers from melting as much as in the other two basins.

Above all, it can be concluded that despite the originally increased, but then reduced, precipitation over 1976–2017, the trend of enhanced mass loss at the basin scale has changed little. However, the accelerated glacier melt was synchronized with the continuously rising summer temperature. The precipitation, which is neglected by Ye and others (2006), is also crucial, especially for small glaciers. The basin-scale glacier mass balance is regulated by the large glaciers. However, the glaciers under investigation in Tibet are mostly smaller than 10 km² (Yao and others, 2012b) ‒ a fact which should be noted when analyzing glaciers' response to climate change over long timescales.

5.3.2 Terrain factors

Terrain factors, such as elevation, slope and aspect are also responsible for spatial differences in glacier change. In our study area, the temperature decreases northwestward with the increase of elevation of glaciers (Ding and others, 1992). The precipitation also decreases with the increasing elevation in this region (Sun and others, 2019), being larger in the south and east than in the north and west parts of the Tanggula Mountains (Ding and others, 1992; Chao and others, 2017). However, the glacier surface elevation change in the west (lower temperature and less precipitation) is smaller than that in the east (higher temperature and more precipitation) (Fig. 6). Therefore, differences of location and elevation, resulting in the differences of temperature and precipitation, could affect the patterns of glacier surface elevation change with longitude (Fig. 6) and elevation (Figure S4).

Slope and aspect influence glacier mass balance by creating differences in energy receipt, snow accumulation and glacier discharge (Arnold and others, 2006; DeBeer and Sharp, 2017). A linear trend was found between glacier elevation change and slope (r = 0.57) (Fig. 10). The bigger the slope, the smaller the glacier melt. This suggests that slope has an effect, as it determines not only the sensitivity of the changes in ELA, but also influences the glacier advance by dynamics (Venkatesh and others, 2012). A simple solar radiation index (SRI) was used to describe the solar radiation in different aspects (Keating and others, 2007). The solar radiation is strongest from the south direction and decreases toward the north (Fig. 11). The rate of glacier melt increased with SRI (r = 0.64). Accordingly, the glaciers mainly concentrate in the north and west aspect directions, where solar radiation is weaker (Table 4).

Fig. 10. Relationship between glacier surface elevation change and glacier mean slope. The 95% prediction bounds are shown.

Fig. 11. Relationship between glacier surface elevation change and mean glacier aspect. (a) A simple SRI is used to define the solar radiation in different aspects. (b) The relationship between glacier elevation change and SRI. The glacier elevation change rates for each 0.05 SRI interval are shown. The trend line is calculated based on the median value of glacier elevation change rate in each interval.

Table 4. Glacier area percentage and SRI in different directions

5.3.3 Glacier surging and surface roughness

Surge-type glaciers exhibit quiescent and active phases, leading to increases in flow velocity, terminus advances, crevassing and surface elevation changes (Meier and Post, 1969). Several surge-type glaciers on the TP have previously been identified (Zhang, 1992; Guo and others, 2013; Sevestre and Benn, 2015; Yasuda and Furuya, 2015). Here, we identify six active surge glaciers over 1976–2000 and six active surge glaciers over 2000–17. The advancing glaciers observed by Ye and others (2006) are all surge glaciers. Surge glaciers can lead to serious mass loss as glacier mass at high elevations is released to the lower parts where there are warmer temperatures (Jiskoot, 2011). The glaciers surging over 1976–2000 were all followed by heavy mass loss over 2000–17 (Table 5). Their mean elevation change over 1976–2000 was −0.39 m a^–1 which almost doubled over 2000–17. Conversely, the surging glaciers over 2000–17 exhibited a higher ablation rate over 1976–2000. Singh and others (2020) reported that surge-type glaciers exhibit an acceleration in motion before the onset of the surge, which may be related to this higher ablation rate.

Table 5. Comparisons of glacier elevation change rates before and after glacier surging

The No. relates to the labels in Figure 5.

^a Denotes the surging period.

Glacier surface roughness also facilitates glacier melt. Varying glacier surface roughness influences radiative incidence angles, surface albedo and the glacier energy balance via its effect on the turbulent fluxes of sensible and latent heat (Smith and others, 2016; Irvine-Fynn and others, 2017). It could have a more significant contribution to GML in a warming climate (Braithwaite and Olesen, 1980). Owing to the high-resolution of KH-9 (6–9 m), we can visually identify glaciers with a rough surface in the ablation area in contrast to others (Fig. 5, examples shown in Figure S7). Except for the surge-type glaciers, these types of glaciers have all melted strongly in the terminus throughout the three periods (Fig. 5), changing more drastically than other types of glaciers nearby (Fig. 6), with values of −0.41 against −0.29 m a^–1 in 1976–2000, and −0.50 against −0.37 m a^–1 in 2000–17.

6. Conclusions

For the first time, the seasonal glacier mass balance over 2016–17 was estimated using a geodetic method to determine the accumulation and ablation characteristics of glaciers in the Geladandong Mountains region. The glaciers in this region are characterized as the summer accumulation type. In 2016/17, the glaciers had a negative mass-balance equivalent to −0.42 m in thickness. The major accumulation and ablation events occur in the warm season. Furthermore, we also determine decadal glacier elevation changes over 1976–2000 (−0.21 ± 0.11 m a^–1), 2000–01 (−0.28 ± 0.14 m a^–1) and 2011–17 (−0.48 ± 0.10 m a^–1), to provide a long-term geodetic glacier mass balance. Our results indicate that the glaciers have undergone accelerated mass loss in the past four decades at the basin scale. Glaciers of <10 km² in size are more sensitive to short-term climatic fluctuations, such as precipitation change. Large glaciers dominate the mass balance at the basin scale. Glaciers larger than 10 km² exhibit accelerated melting over the entire study period, controlled by the continuous increase in temperature. Furthermore, lake volume changes and in situ river discharge are estimated, to explain the water balance from decadal GML. The contributions of GML to increased CC-DC and Selin Co lake volume and Tuotuohe river discharge changes decreased from 20.0, 8.6 and 8.0% over 1976–2000 to 19.2, 3.5 and 4.3 over the wetter period of 2000–11, respectively; but then they increased to 21.4 and 16.3% in CC-DC and Selin Co with precipitation decline over 2011–2017, respectively. Overall, the contribution is small, but changes significantly with the variations of seasonal water supply, and functions as an important component in regulating lake expansion and runoff downstream.