Recent Accelerating Glacier Mass Loss of the Geladandong Mountain, Inner Tibetan Plateau, Estimated from ZiYuan-3 and TanDEM-X Measurements

The headwaters of many Asian rivers are at mountain glaciers of the Tibetan Plateau. Glacier melt-water is a non-negligible contributor of river runoff, especially for a drought year. However, the observation of mass glacier changes was scarce in recent years. Here, we estimated the recent glacier mass change of the Geladandong mountain, by differencing the digital elevation models (DEMs) produced from ZiYuan-3 images and TanDEM-X data. Moreover, we compared the SRTM-C DEM with TanDEM-X DEMs to retrieve glacier mass balances since 2000. The annual mass loss rates of −0.11 ± 0.03 and −0.47 ± 0.09 m w.e. yr−1 were derived in 2000–2012 and 2012−2018, respectively. This result revealed an accelerating rate of negative glacier mass changes during recent years, which is mainly caused by the significant increase of mass loss over non-surge glaciers, rather than surge-type glaciers, which held a slight increase of mass loss. In addition, we found a pronounced discrepancy of glacier mass change between non-surge and surge-type glaciers during 2012−2018, and suggested that this difference may be caused by the heterogeneous responses of surge-type glaciers to climate variations, because of the different timing and type of surge events.


Introduction
The Tibetan Plateau (TP), which is one of the most glacierized areas outside the Arctic and Antarctic, is the water tower of Asia [1]. In general, the headwaters of many great Asian rivers (e.g., Yangtze River) are located over mountain glaciers in the TP [2,3]. Glacier melt-water is a non-negligible contributor of these rivers' runoff [4][5][6][7][8]. Particularly in a drought year, runoff is mainly originated from glacier melt-water, rather than precipitation, over some river basins [9]. Therefore, the observation of glacier mass balances is necessary in the TP for investigating the variations of regional water resources.

Study Area
The GLDDM which is located in the inner plateau is the highest peak of the Tanggula mountain range [30]. Mountain glaciers over the GLDDM are summer-accumulation type because snowfall mainly occurs between June and August [45]. Glacier melt-water of the GLDDM mainly flows into the source region of the Yangtze River (SRYR) by the Tuotuo River and the Garang River and the endorheic basin of the TP (EBTP) (Figure 1a). A hydrological station (Tuotuohe station) has been installed downstream of the Tuotuo River to record river runoff (Figure 1a). The glacierized regions of the GLDDM are composed of more than 150 glaciers (Figure 1b). The regional climate of the GLDDM is mainly controlled by continental climatic conditions, rather than the Indian summer monsoon and the mid-latitude westerlies [13]. The records of ice core on the upper area of Guoqu Glacier revealed a continuous rise of air temperature for the past 500 years, and an accelerated climate warming in the 20th century [46]. Moreover, a significant increase in the mean annual temperature since the 1960s has been recorded by meteorological stations near the GLDDM as well [31].

Glacier Surface Topographic Data
We employed three space-borne DEMs to detect glacier mass balances of the GLDDM. The most recent DEM was generated from ZY3 tri-stereo optical images. The first satellite of the ZY3 series (ZY3-01), which was launched on 9 January 2012, equipped with three panchromatic cameras working in time-delay integration mode for linear push-broom imaging [48]. The ZY3-02 satellite, which is the second satellite of the ZY3 series, was launched on 30 May 2016 in the same orbit of the ZY3-01 satellite with 180˚ phase apart [49]. In general, the main sensor parameters of the ZY3-02 are similar to those of the ZY3-01. However, the spatial resolution of the forward and backward cameras of the ZY3-02 is higher than that of the ZY3-01 (from 3.5 m to 2.5 m) [50]. Here, considering spatial resolution and cloud cover, we employed three pairs of ZY3-02 tri-stereo images acquired in December 2016 and December 2017 to reconstruct the glacier surface topography of the GLDDM. The basic parameters of the used ZY3-02 optical scenes are summarized in Table 1. The ZY3 DEMs were produced by the Space Data Processor (SDP) software from the Land Satellite Remote Sensing Application Center (LSRSAC), Ministry of Natural Resources of the People's Republic of China. Detailed descriptions of the generation of ZY3 DEMs can be found in Tang et al. (2018) [51].
We extracted the glacier surface DEMs in the early 2010s from three pairs of TanDEM-X InSAR data obtained in March and November 2012 ( Table 1). The TanDEM-X mission, which is the first satellite bistatic interferometric configuration, is comprised of two identical satellites of TerraSAR-X and TanDEM-X [44]. These two satellites were launched by the German Aerospace Center in June 2007 and June 2010 [44]. Here, the used three pairs of TanDEM-X InSAR data were obtained in the CoSSC format. The spatial resolution of these TanDEM-X InSAR data is ~3 m in both azimuth and range directions. Moreover, we used the bistatic SAR interferometric method [37] to generate three TanDEM-X DEMs in March and November 2012.

Glacier Surface Topographic Data
We employed three space-borne DEMs to detect glacier mass balances of the GLDDM. The most recent DEM was generated from ZY3 tri-stereo optical images. The first satellite of the ZY3 series (ZY3-01), which was launched on 9 January 2012, equipped with three panchromatic cameras working in time-delay integration mode for linear push-broom imaging [48]. The ZY3-02 satellite, which is the second satellite of the ZY3 series, was launched on 30 May 2016 in the same orbit of the ZY3-01 satellite with 180 • phase apart [49]. In general, the main sensor parameters of the ZY3-02 are similar to those of the ZY3-01. However, the spatial resolution of the forward and backward cameras of the ZY3-02 is higher than that of the ZY3-01 (from 3.5 m to 2.5 m) [50]. Here, considering spatial resolution and cloud cover, we employed three pairs of ZY3-02 tri-stereo images acquired in December 2016 and December 2017 to reconstruct the glacier surface topography of the GLDDM. The basic parameters of the used ZY3-02 optical scenes are summarized in Table 1. The ZY3 DEMs were produced by the Space Data Processor (SDP) software from the Land Satellite Remote Sensing Application Center (LSRSAC), Ministry of Natural Resources of the People's Republic of China. Detailed descriptions of the generation of ZY3 DEMs can be found in Tang et al. (2018) [51].
We extracted the glacier surface DEMs in the early 2010s from three pairs of TanDEM-X InSAR data obtained in March and November 2012 ( Table 1). The TanDEM-X mission, which is the first satellite bistatic interferometric configuration, is comprised of two identical satellites of TerraSAR-X and TanDEM-X [44]. These two satellites were launched by the German Aerospace Center in June 2007 and June 2010 [44]. Here, the used three pairs of TanDEM-X InSAR data were obtained in the CoSSC format. The spatial resolution of these TanDEM-X InSAR data is~3 m in both azimuth and range directions. Moreover, we used the bistatic SAR interferometric method [37] to generate three TanDEM-X DEMs in March and November 2012.
We derived the glacier surface DEM in 2000 from the Shuttle Radar Topography Mission (SRTM), which is the first spaceborne single-pass InSAR system [52]. Surface topographic data for most of the Remote Sens. 2020, 12, 472 4 of 18 continental regions (56ºS-60ºN) have been reconstructed by the SRTM mission [52]. In general, this topography mission was equipped with two InSAR systems of C-band and X-band. The vertical and horizontal accuracy of the X-band topographic product (the SRTM-X DEM) is better than that of the C-band topographic product (the SRTM-C DEM) [53]. However, only a few parts (less than 10%) of the glacierized regions of the GLDDM are covered by the SRTM-X DEM. Here, the SRTM-C DEM with a spatial resolution of 30 m, which was provided by the United States Geological Survey, was employed for better coverage of the glacierized regions over the GLDDM.

Glacier Outlines
Glacier boundaries, which were derived from the Second Chinese Glacier Inventory (SCGI), were employed to identify the glacierized regions of the GLDDM. In this study area, a Landsat Thematic Mapper (TM) image which was acquired in May 2007 was used to delineate the glacier outlines of the SCGI [47]. Therefore, in order to match with the SRTM-C DEM obtained in February 2000, the terminus locations of glaciers with apparent retreats during the studied periods were manually corrected with optical satellite images. Here, we employed one cloud-free Landsat TM image obtained in 1999 (Figure 1b), because the Landsat optical scenes in 2000 suffered from snow or cloud cover. Moreover, the recent terminus locations of some glaciers with surging events in 2000-2018 were updated with Landsat scenes corresponding to the acquisition dates of the used TanDEM-X InSAR data and ZY3-02 tri-stereo images.

DEM Co-Registration and Differencing
In this study, most of ZY3-02 and TanDEM-X images were acquired in December 2017 and March 2012, respectively. For detecting the annual rate of glacier elevation changes, we neglected the possible inter-annual variation and assumed that the observed results could represent the mean variations of the following time periods: 2000-2012, 2012-2018, and 2000-2018. Specifically, we subtracted the SRTM-C DEM from TanDEM-X DEMs and ZY3 DEMs to detect glacier elevation changes in 2000-2012 and 2000-2018, respectively. Moreover, we also compared TanDEM-X DEMs with ZY3 DEMs to observe recent glacier elevation changes (2012-2018). Before the process of DEM differencing, a pair of DEMs were first accurately co-registered with the universal co-registration method [54]. Moreover, the bias related to different spatial resolutions of the used DEMs were removed by applying a polynomial function between elevation difference and altitude with pixels over ice-free regions. Similar to that in Liu et al. (2016) [29], we excluded the glacier surface pixels with slope angles of larger than 25 • . The outliers were also identified for elevation differences larger than ±15 m yr −1 for surge-type glaciers and ±10 m yr −1 for other glacierized regions. These thresholds were determined by investigating previous results in this study site (e.g., Chao et al. 2017 [31]) and our observed glacier elevation changes. The void areas in the elevation change map which are related to these outliers or data gaps in ZY3 DEMs due to snow or cloud cover were not filled by interpolation.
According to the assumption of similar elevation differences of glacier pixels at an altitude band [55], we accurately calculated the mean glacier elevation change of the GLDDM for 50 m altitude bands. Moreover, at a given altitude interval, three times standard deviation were used to select the Remote Sens. 2020, 12, 472 5 of 18 valid pixels for computing mean value, in order to suppress the impact of random error. Note, the feature of glacier elevation changes with an altitude over surge-type glaciers is generally different to that of non-surge glaciers [22]. Therefore, for each surge-type glacier, we calculated its mean glacier elevation change separately. In addition, the mean glacier elevation change of the entire GLDDM was computed with the area-weighted average of the observed results for non-surge glaciers and all surge-type glaciers.

Corrections of Systematic Biases
In this study, radar signal penetration depth and seasonal glacier elevation change are the two possible systematic biases for estimating glacier mass balance. In general, when InSAR DEM was applied to observe glacier elevation change, radar signal penetration depth is a non-negligible systematic bias [56][57][58]. Here, we employed the X-band radar penetration depth of the Puruogangri ice field which was estimated by comparing the TanDEM-X DEMs in January and April 2012 [37]. We computed the correction values of X-band radar penetration depth for 50 m altitude bands and applied it to each altitude interval separately. Moreover, the C-band radar signal penetration depth was estimated by adding the X-band radar penetration depth to the difference in penetration depth between the C-band and X-band radar signals. For the limited coverage (less than 10%) of the SRTM-X DEM over the GLDDM, we estimated the C/X-band radar penetration difference by comparing the SRTM-C DEM and the SRTM-X DEM on the Puruogangri ice field.
Seasonal glacier elevation change is also a possible systematic bias because the used DEMs were generated from space-borne images obtained in different months of the year. As listed in Table 1, possible seasonal variations needed to be corrected in winter and early spring. However, snowfall and melt mainly occur in the summer season, because glaciers over the GLDDM belong to the summer-accumulation type [45]. This indicates that the GLDDM experiences a slight glacier elevation change in the months of winter and early spring. Thus, no correction of seasonal glacier elevation change was used in this study.

Glacier Mass Balance Estimation and Error Analysis
We employed the conversion factor of 850 kg m −3 [59] to estimate glacier mass balance from the observed mean glacier elevation change. In summary, the glacier mass balance of the GLDDM (∆MB) was calculated with the following equation.
where S total is the total glacier area of the entire GLDDM. ∆MB i and S i are glacier mass change and glacier area over an altitude interval, respectively.
where ∆h i and P i are the mean glacier elevation change and radar signal penetration depth over an altitude interval, respectively. Note, when subtracting the SRTM-C DEM from TanDEM-X DEMs, P i is the penetration difference of the C/X-band radar signals. T year is the integer number of years. ρ and ρ w are conversion factors (850 kg m −3 ) and the density of water (1000 kg m −3 ), respectively. We employed the standard law of error propagation to evaluate the uncertainty of the estimated glacier mass balance. Therefore, the error of glacier mass balance (σ ∆MB ) can be computed with Equation (3). where σ S total is the error of the glacier area of the entire GLDDM. In this study, we assumed an uncertainty of ±5% for the total glacier area [60], because glacier terminus locations were corrected with satellite optical data. σ ∆MB i is the error of our observed glacier mass change for an altitude band.
where σ ρ is the uncertainty of the conversion factor. Huss (2013) suggested that σ ρ (±60 kg m −3 ) is approximately ±7% of the used conversion factor [59]. σ ∆h i and σ P i are the uncertainties of mean glacier elevation change and radar signal penetration depth over an altitude band, respectively. Here, these two uncertainties were estimated with Equation (5), because mean values were calculated for 50 m altitude bands.
where ST i and Ne i are standard deviation and the number of glacier pixels with independent observations over an altitude bin, respectively.
For an altitude band, Nt i is the number of all glacier pixels. SR is the spatial resolution of topographic data. Here, D sa (spatial autocorrelation distance) was calculated by applying the semivariogram model [37].

Glacier Surface Elevation Changes
Annual glacier elevation changes of the GLDDM are given in Figure 2 for the periods of 2000-2012, 2012-2018, and 2000-2018. In general, we found a clear spatial variation of glacier elevation changes for the three time intervals. Specifically, an apparent glacier surface lowering was basically detected over the terminus zones of the GLDDM, whereas the upper zones experienced a slight glacier surface thinning or even surface thickening. However, a reverse spatial variation of surface elevation changes was derived for several surge-type glaciers (Figure 2), because of surging events in 2000-2018. For glaciers with surging events in the study period, we detected a significant increase of ice thickness at the terminus and a pronounced surface lowering at the upper regions.  changes for the three time intervals. Specifically, an apparent glacier surface lowering was basically detected over the terminus zones of the GLDDM, whereas the upper zones experienced a slight glacier surface thinning or even surface thickening. However, a reverse spatial variation of surface elevation changes was derived for several surge-type glaciers (Figure 2), because of surging events in 2000-2018. For glaciers with surging events in the study period, we detected a significant increase of ice thickness at the terminus and a pronounced surface lowering at the upper regions. The altitudinal distribution of the observed annual elevation changes shows that glacier elevation changes over the GLDDM were affected by altitude in the studied periods ( Figure 3). This altitude-dependent result is basically consistent with the spatial pattern of glacier elevation changes ( Figure 2). Furthermore, we found that the values of glacier surface thinning in 2012-2018 were greater than that in 2000-2012 for almost all altitude bands over glacier tongue regions (lower than 5650 m). In addition, over the upper regions, the GLDDM experienced an apparent glacier surface thickening for all altitude bands in 2000-2012, whereas glacier surface thinning or a slight increase The altitudinal distribution of the observed annual elevation changes shows that glacier elevation changes over the GLDDM were affected by altitude in the studied periods ( Figure 3). This altitude-dependent result is basically consistent with the spatial pattern of glacier elevation changes ( Figure 2). Furthermore, we found that the values of glacier surface thinning in 2012-2018 were greater than that in 2000-2012 for almost all altitude bands over glacier tongue regions (lower than 5650 m). In addition, over the upper regions, the GLDDM experienced an apparent glacier surface thickening for all altitude bands in 2000-2012, whereas glacier surface thinning or a slight increase of ice thickness was detected for most of the altitude bands between 2012 and 2018. This difference of altitudinal distribution of glacier surface elevation changes in 2000-2012 and 2012-2018 generally supports the accelerating rate of glacier lowering over the GLDDM during recent years. It is noteworthy that the altitude band of 5400-5450 m experienced serious surface thinning, which is probably caused by glacier termination. Specifically, the terminus altitudes of some glaciers over the GLDDM are between 5400 and 5450 m.  (Table 2). This result shows that the negative glacier mass change in 2012-2018 is much more than that in 20002012, and thus reveals a probably accelerated glacier mass loss during recent years. Furthermore, -for non-surge glaciers, we detected a similar temporal variation of glacier mass balance as that over the entire GLDDM. The glacier mass loss of non-surge glaciers in 2012-2018 (−0.53 ± 0.12 m w.e. yr −1 ) was about four times larger than that in 2000-2012 (−0.12 ± 0.04 m w.e. yr −1 ). However, the surge-type glaciers of the GLDDM experienced only a slight increase of glacier mass loss in recent years, from −0.08 ± 0.03 m w.e. yr −1 in 2000-2012 to −0.15 ± 0.05 m w.e. yr −1 in 2012-2018. Consequently, the observed accelerating rate of glacier mass loss over the GLDDM is mainly attributed to the significant increase of negative mass changes of non-surge glaciers since the early 2010s.

Comparison with Previous Studies
Until now, the field measurement of glacier mass balances was not conducted over the GLDDM. In this study, the measured glacier mass changes were compared to the published geodetic estimates and the modeled annual glacier mass changes (Table 3)  Note: a Geod., geodetic method; Albed., albedo-based regression method.
Our geodetic estimates revealed an accelerating rate of negative mass changes over the GLDDM since the early 2010s. This recent accelerated glacier mass loss was also detected over several other glacierized regions in the ITP (e.g., the Xiaodongkemadi Glacier [61] and the Puruogangri ice field [37]), or even some glaciers in the Pamir mountains [24] and the southeast TP [11]. However, we cannot simply infer that mountain glaciers in the TP and its surroundings generally experienced an accelerating rate of negative mass changes, because of the heterogeneous responses of glaciers to regional climate variations [17,64].

The Influence of X-band Radar Penetration in Geodetic Estimates
When InSAR data was employed in the geodetic method, a bias caused by radar signal penetration depth is needed to correct for glacier mass balance estimation [56][57][58]. For the nearly exhaustive coverage of the glacierized areas in the TP and its surroundings, the SRTM-C DEM has been widely used as a reference topographic data to estimate glacier mass balance [22,23,35]. Moreover, previous works have mainly focused on investigating the penetration depth of the C-band radar signal [56,57]. The penetration depth of the X-band radar signal is generally neglected for glaciers in these glacierized areas [35,56]. However, the X-band radar penetration depth in the winter season was found to be non-negligible over the Fedchenko Glacier, Pamir mountains [24], and the Puruogangri ice field, ITP [37].
In this study, we extracted the glacier surface DEMs in the early 2010s from three pairs of X-band InSAR data. The correction values of the penetration depth of the X-band radar signal, which were computed for 50 m altitude bands over the Puruogangri ice field, were applied to each altitude interval of the GLDDM separately. A mean value of 0.53 ± 0.05 m was detected for the X-band radar penetration depth of the entire GLDDM. Generally, this X-band radar penetration caused a bias of 0.07 ± 0.01 m w.e. yr −1 for the result of glacier mass change during 2012−2018. This indicates that X-band radar penetration cannot be simply neglected for glacier mass change estimation on the GLDDM, or even in the ITP. Basically, the impact of radar penetration depth in geodetic estimates is related to the integer number of years for the observation period ( Figure 4). Consequently, for a short study period of less than ten years (especially for 1−5 years), accurate correction of X-band radar penetration depth is necessary for geodetic estimates of glacier mass balance obtained from the X-band InSAR DEMs.
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 19 Figure 4. The impact of X-band radar penetration in mass balance estimation versus integer number of years for the study period. For visualization, the error bars indicate five times the estimated uncertainties.

Difference in Mass Balances of Surge-Type Glaciers
Glacier surging events which are characterized by the quick movement of ice masses from upper to lower zones [65,66] may lead to geological hazards [67,68]. Many glaciers in the Pamir and Karakoram mountains have surged during recent decades [57,69]. Moreover, in this mountain range, the mass balance of surge-type glaciers was almost equivalent to non-surge glaciers during 1999−2011 [57]. In this study, six glaciers of the GLDDM experienced surging events between 2000 and 2018 ( Figure 2). Furthermore, we observed a similar glacier mass change for non-surge (-0.12 ± 0.04 m w.e. yr -1 ) and surge-type (-0.08 ± 0.03 m w.e. yr -1 ) glaciers in 2000−2012. However, an apparent difference of glacier mass balance was detected between non-surge (-0.53 ± 0.12 m w.e. yr -1 ) and surge-type (-0.15 ± 0.05 m w.e. yr -1 ) glaciers over the GLDDM in 2012-2018 (Table 2). This discrepancy may be partly caused by the spatially heterogeneous response of surge-type glaciers to regional climatic variations in recent years ( Figure 5).
In general, we detected a relatively balanced glacier mass change (between -0.22 ± 0.08 and 0.09 ± 0.03 m w.e. yr -1 ) over all the six surge-type glaciers of the GLDDM in 2000−2012 ( Figure 5). However, the difference in the observed glacier mass balances of these surge-type glaciers was pronounced between 2012 and 2018 ( Figure 5). Glacier mass gains (up to 0.47 ± 0.12 m w.e. yr -1 ) were detected over three surge-type glaciers of 5K451F0030, 5Z213A0007, and 5Z213B0003, whereas the other three surge-type glaciers experienced an apparent negative mass change (up to -0.51 ± 0.13 m w.e. yr -1 ). Furthermore, we found that the surge-type glaciers with mass gains in 2012−2018 have surged in both time periods of 2000−2012 and 2012−2018 ( Figure 6). In contrast, the glaciers of 5K451F0036 and 5K451F0008 only surged between 2000 and 2012 (Figure 7). This means that mass losses of these two glaciers in 2012−2018 are possibly attributed to rapid ablation in the post-surge phase. In addition, by analyzing the changes in the terminus location and surface elevation, we found that the recent surge

Difference in Mass Balances of Surge-Type Glaciers
Glacier surging events which are characterized by the quick movement of ice masses from upper to lower zones [65,66] may lead to geological hazards [67,68]. Many glaciers in the Pamir and Karakoram mountains have surged during recent decades [57,69]. Moreover, in this mountain range, the mass balance of surge-type glaciers was almost equivalent to non-surge glaciers during 1999−2011 [57]. In this study, six glaciers of the GLDDM experienced surging events between 2000 and 2018 ( Figure 2). Furthermore, we observed a similar glacier mass change for non-surge (−0.12 ± 0.04 m w.e. yr −1 ) and surge-type (−0.08 ± 0.03 m w.e. yr −1 ) glaciers in 2000−2012. However, an apparent difference of glacier mass balance was detected between non-surge (−0.53 ± 0.12 m w.e. yr −1 ) and surge-type (−0.15 ± 0.05 m w.e. yr −1 ) glaciers over the GLDDM in 2012-2018 (Table 2). This discrepancy may be partly caused by the spatially heterogeneous response of surge-type glaciers to regional climatic variations in recent years ( Figure 5). events of the other five glaciers belong to Svalbard-type (active phase of more than 3 years). Consequently, the different responses of these six surge-type glaciers over the GLDDM to recent climate variations may be caused by the different timing and type of surge.  In general, we detected a relatively balanced glacier mass change (between −0.22 ± 0.08 and 0.09 ± 0.03 m w.e. yr −1 ) over all the six surge-type glaciers of the GLDDM in 2000−2012 ( Figure 5). However, the difference in the observed glacier mass balances of these surge-type glaciers was pronounced between 2012 and 2018 ( Figure 5). Glacier mass gains (up to 0.47 ± 0.12 m w.e. yr −1 ) were detected over three surge-type glaciers of 5K451F0030, 5Z213A0007, and 5Z213B0003, whereas the other three surge-type glaciers experienced an apparent negative mass change (up to −0.51 ± 0.13 m w.e. yr −1 ). Furthermore, we found that the surge-type glaciers with mass gains in 2012−2018 have surged in both time periods of 2000-2012 and 2012−2018 ( Figure 6). In contrast, the glaciers of 5K451F0036 and 5K451F0008 only surged between 2000 and 2012 ( Figure 7). This means that mass losses of these two glaciers in 2012-2018 are possibly attributed to rapid ablation in the post-surge phase. In addition, by analyzing the changes in the terminus location and surface elevation, we found that the recent surge timing of Glacier 5K444B0064 is 2014-2015. We thus inferred that the Alaskan-type surge (active phase of 1-3 years) may be the reason for recent mass loss over this surge-type glacier. The surge events of the other five glaciers belong to Svalbard-type (active phase of more than 3 years). Consequently, the different responses of these six surge-type glaciers over the GLDDM to recent climate variations may be caused by the different timing and type of surge. Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 19

Contribution of Glacier Mass Loss to River Runoff
As shown in Figure 1a, the glacier melt-water of the GLDDM flows into the SRYR by the two rivers of Tuotuo and Garang. The runoff of the Tuotuo River has been recorded by a hydrological station downstream. Here, we detected that the glacier mass balances of the GLDDM in the SRYR were −0.029 ± 0.011, −0.187 ± 0.043, and −0.078 ± 0.032 Gt yr −1 in the periods of 2000-2012, 2012-2018, and 2000-2018, respectively (Table 4). Moreover, over the catchment of the Tuotuo River, we derived the glacier mass changes of −0.025 ± 0.010 (2000−2012) and −0.102 ± 0.031 Gt yr −1 (2012−2018). Considering the annual river runoff of 1.35 Gt yr −1 recorded by the Tuotuohe station in the 2000s [70], the contribution of glacier mass loss to the runoff of the Tuotuo River was only~2% in the first decade of the 21st century. However, for the entire SRYR, glacier mass loss contributed to 17.5% of the total river runoff between 1986 and 2009 [6]. This difference may be related to the pronounced mass losses of other glaciers in the SRYR. For example, Xiaodongkemadi Glacier experienced a mass loss of −0.36 m w.e. yr −1 in 2000−2010 [13], which is about three times that over the GLDDM in 2000-2012.
Overall, during the past decades, glacier melt-water was not a major contributor to river runoff in the SRYR.

Conclusions
In this study, glacier mass changes of the GLDDM since 2000 were investigated by employing the SRTM-C DEM and glacier surface topographic data generated from TanDEM-X InSAR images and ZY3 tri-stereo scenes. Glacier mass balances of −0.11 ± 0.03, −0.47 ± 0.09, and −0.24 ± 0.07 m w.e. yr −1 were derived for the periods of 2000-2012, 2012-2018, and 2000-2018, respectively. These estimates revealed an accelerating rate of glacier mass loss for the GLDDM during the past few years. This finding is basically supported by the modeled annual glacier mass balances between 2000 and 2016 [61]. Moreover, we found that this accelerated mass loss is mainly attributed to a significant increase in negative glacier mass changes of non-surge glaciers in 2012-2018.
Glacier mass change of surge-type glaciers was generally similar to that of non-surge glaciers over the GLDDM in 2000−2012. It is basically consistent with the similarity of mass balances of non-surge and surge-type glaciers in the Pamir and Karakoram between 1999 and 2011. However, in the period 2012−2018, we detected an apparent difference in glacier mass balance for the two types of glaciers, which may be partly caused by the heterogeneous responses of these surge-type glaciers to climate variations in recent years due to the different timing and type of surge.
Glacier mass loss only contributed to~2% of the runoff of the Tuotuo River during the first decade of the 21st century. In general, glacier mass loss was not the major contributor to river runoff in the SRYR in the past decades, although glacier mass loss provides 17.5% of the river runoff in the entire SRYR between 1986 and 2009. A significant increase of glacier melt-water was possibly detected for several recent years due to accelerated glacier mass loss over some glacierized regions such as the GLDDM and the Xiaodongkemadi Glacier. In the future, continuous measurement of glacier mass balance is needed in the SRYR to investigate whether the recent accelerating rate of glacier mass loss extends to a long-term trend.