Recent Glacier Mass Balance and Area Changes from DEMs and Landsat Images in Upper Reach of Shule River Basin, Northeastern Edge of Tibetan Plateau during 2000 to 2015

: Glacier changes in the Upper Reach of the Shule River Basin (URSRB) serve as a good indicator of climate change in the western part of the Qilian Mountains, located on the northeastern edge of the Tibetan Plateau. However, information on recent glacier changes in the URSRB is limited. In this study, the changes in ice surface elevation were determined using geodetic methods based on digital elevation models (DEMs) derived from the Shuttle Radar Topography Mission (SRTM) (2000), and from pairs of Third Resources Satellite (ZY-3) of China (taken around 2013). In addition, glacier area changes from 2000–2015, were derived from Landsat TM/ETM+/OLI images. The results suggest that 478 glaciers with an area of 375.1 ± 2.68 km 2 remained in the URSRB in 2015. Ice cover diminished by 57.5 ± 2.68 km 2 (11.9 ± 0.60%), or 0.79 ± 0.04% a − 1 and 35 small glaciers disappeared from 2000 to 2015 in the URSRB. The most pronounced glacier shrinkage occurred during 2004 to 2009. The average ice surface elevation of the URSRB from 1999 to 2013 reduced by about 4.98 ± 0.6 m, which is equal to a mass loss of 0.383 ± 0.046 m · a − 1 . This reduction indicates that the ice storage loss has accelerated since 1999, compared to a mass loss of 0.21 ± 0.04 m · a − 1 around Tuanjiefeng from 1966 to 1999, as reported by Xu et al. (2013).


Introduction
Mountain glaciers are a critical freshwater resource for populations inhabiting the downstream region or valleys at the base of mountains [1,2]. Worldwide glacier changes and the associated change in future runoff raise major concerns for the sustainability of global water resources [3]. The meltwater from glaciers on the Tibetan Plateau (TP), which is called "Asia's Water Tower", is extremely important to the vitality of large rivers, such as the Indus and Brahmaputra Rivers [1]. Therefore, reduction in meltwater can be a major threat to the food security of an estimated 60 million people who live in the downstream regions [1,4]. The glacier mass balance is a very useful indicator of climate variability [5], as it provides information about glacier changes in high Asia as an important issue in cryospheric and climate change studies [6].
Due to recent climate warming [7], many mountain glaciers in Asia have shrunk in area and lost mass, such as glaciers of the Tianshan Mountains [8][9][10][11], Qiangtang Plateau [12], Aksu-Tarim catchment [13], central Himalaya [14,15], and Qilian Mountains [16]. However, some stabilization or 38.2 • -40.0 • N) and has an area of about 1.14 × 10 4 km 2 . The URSRB is characterized by the continental climate of the plateau and an elevation between 2081-5764 m a.s.l. (Figure 1). The yearly annual average air temperature is about −4 • C, and the average annual precipitation is about 378.4 mm. Precipitation mainly occurs during the glacier melting season (May-September), accounting for more than 90% of the annual precipitation [34].

Data
Satellite images from Landsat 5, 7, and 8 were used to investigate the changes in the glacier distribution and area. Landsat 5 Thematic Mapper (TM) images contained seven bands, of which the spatial resolution of band 1-5 and band 7 is 30 m. Landsat 7 Enhanced Thematic Mapper (ETM+) images included 8 bands, and the spatial resolutions of bands 1-5, 7, and 8 were 30, 30, and 15 m, respectively. Landsat 8 Operational Land Imager (OLI) images consisted of nine bands, including 8 bands with a spatial resolution of 30 m and one pan band of 15 m. All satellite images were derived from the geospatial data cloud platform (http://www.gscloud.cn/) of the Computer Network Information Center (CNIC) of the Chinese Academy of Sciences. All products were Level L1T and were calibrated by radiation correction, geometry correction by ground control points, and terrain correction by the digital elevation model (DEM). All selected Landsat images had cloud cover less than 15%. Additional information about the images are listed in Table 1. All Landsat data on glacier area changes from 2000 to 2015 were extracted between July and September, mainly in August, to reduce the effect of seasonal snow cover in the images.
SRTM was launched at the Kennedy Space Center on 11 February 2000, which used microwaves to build a global DEM. The average elevation accuracy of the SRTM DEM is better than 16 m (90% confidence interval), but the accuracy is closely related to topographic factors, such as surface slope and slope direction [35]. The unfilled version SRTM DEM (containing data voids) was obtained from interferometry of the C-band of a microwave sensor and post-processing with 90 m resolution. Information from the SRTM DEM can be found on the website http://glcf.umd.edu/data/srtm. The X-band SAR system carried with the SRTM C-band SAR system had a swath width of 45 km, leaving large data gaps in the derived X-band DEM [36]. SRTM X-band DEMs only cover 21% of URSRB glaciers, which can be downloaded from DLR multi-mission web portal https://geoservice. dlr.de/egp/. The SRTM DEM and X-band DEM were resampled with 30 m resolution by a bilinear interpolation method. The ZY-3, launched on 9 January 2012, is China's first civilian high-resolution stereo mapping satellite [37]. ZY-3 contains three panchromatic cameras that are organized in a triple-linear array imaging geometry, and positioned front-facing at 22 • , ground-facing, and rear-facing at 22 • . The ground pixel resolutions of the front-facing and rear-facing cameras are about 3.5 m, and that of the ground-facing camera is about 2.1 m. We selected ZY-3 images that showed the least cloud cover and that were acquired in the same season. The information on the 5 selected ZY-3 images is list in Table 2. A small part of the URSRB was not imaged properly by the ZY-3 stereo pair cameras, due to the influence of cloud cover. ZY-3 DEMs were resampled to 30 m in comparison to SRTM DEMs.  (7), which correspond with the elevation changes in the glacier sub-regions.
There are 12 national meteorological stations in and around the URSRB, which have obtained observation data from 1970 to 2013. Some precipitation stations with T200B and total rain gauges, are set up in alpine regions, which have acquired observed precipitation data after 2008. All precipitation data observed by 40 total rain gauges at elevations ranging from 1139 to 4156 m a.s.l. were used to obtain the precipitation gradient. The monthly precipitation gradient of 12.0 mm/100 m and air temperature lapse rate were calculated using linear regression between meteorological and elevation data [38,39], and the regional averages of the annual average air temperature and annual precipitation were used to analyze the glacier response to long-term climate changes in the URSRB.

Glacier Area Extraction
A semi-automated approach, using the band ratio between different Landsat bands and corrected by visual interpretation, was applied to outline the glaciers. The band ratio of the Landsat TM to Landsat ETM+ with TM3/TM5 band was used. The extraction threshold of the glacier was 2.1 [40]. Due to the Band 1 of Landsat 8 OLI image is new, all other bands were shifted compared to the older Landsat series, and the band ratio between Red and SWIR1 was shifted to TM4/TM6 for Landsat 8 OLI with a glacier extraction threshold of 2 [41]. To eliminate ice patches less than 0.01 km 2 in area from the images, a 3 × 3 median filter was first applied to eliminate isolated pixels [42], and then "Eliminate" tool with the area 0.01 km 2 in ARCGIS software was used.
The derived glacier polygons were checked manually against images from adjacent years with little or no snow and cloud coverage, to discriminate supraglacial lakes, seasonal snow, surface boulders, and debris-covered ice. The glacier outlines from all Landsat images were combined to produce the glacier polygons of five periods (2000/2001, 2004/2006, 2008/2009, 2012/2013, 2014/2015).
The distribution of mountain ridge was derived from the elevation of DEMs according the method of Guo et al. [43], and then, the boundary of each single glacier was determined by dividing the glacier polygon by ridge. Finally, the outlines of each glacier were manually compared using pseudo-color composed of bands 4/3/2 from Landsat images and glacier boundaries from the Second Chinese Glacier Inventory(SCGI), which contains glacier outlines from 2006 to 2010 [43].
The identification number (ID) of each glacier was obtained from the Chinese Glacier Inventory (CGI). The area changes of each single glacier were calculated by comparing the area in the five periods, and CGI and SCGI with the same ID.

Glacier Ice Surface Elevation Generation
The stereo pair images composed by rear-facing and ground-facing sensors were used to generate DEMs, and the stereo pair images composed by front-facing and ground-facing sensors were taken as assist images. In the process of extracting DEMs from the ZY-3 stereo pair images, the rational polynomial coefficients (RPCs) from the image files were automatically identified by ENVI software. The same homologous points on both images were manually identified and combined, which included about 80 to 100 homologous points on each pair of images, then the selection stopped when the largest y-parallax was less than 2 • . Both the RPCs and selected homologous points were used for map orthorectification. The average maximum y-parallax of the five ZY-3 image pairs was about 1.2 degrees, suggesting that the orthorectification has high precision. Then, the absolute elevation was obtained by adding the measured control points in the URSRB.
The absolute orientation parameters of high-resolution satellite stereo image pairs, such as ZY-3, can be calculated from precise ephemeris data and sensor model's parameters without ground control points, to obtain a high-precision DEM [44]. The dynamic exterior systematic errors (including camera installation errors) and static interior distortion of ZY-3 images were eliminated using 1:2000 digital orthophotomaps and DEMs of the Dengfeng (Henan) and Tianjin areas of China as control data. The plan and height accuracy of ZY-3 DEMs without ground control points were found to be better than 3 and 2 m over the Anping area in Hebei Province and Zhaodong area in Heilongjiang Province, respectively [45]. Our study also proves that the ZY-3 can obtain good quality DEMs.
For changes in glacier surface elevation, differences between the ZY-3 DEM and SRTM DEM were used to construct a map of the ice surface elevation change. Relative horizontal and vertical distortions between the two datasets were corrected with statistical approaches based on the relationship between elevation difference, slope, and aspect [46]. Elevation differences in off-glacier regions were used to analyze the consistency of the ZY-3 DEM and SRTM DEM.
We first compared the spatial patterns of elevation differences with that of a greyscale hill shadow map by visual interpretation, where consistency in the off-glacier regions indicates that there is no obvious horizontal spatial mismatch between DEMs [47]. We determined that the relationship between the elevation difference, which divides the tangent of the slope of terrain and the cosinusoidal of the terrain aspect [47], was far from ideal. In addition, a close relationship was found between the elevation difference and terrain aspect ( Figure 2). Since the 5 ZY-3 satellites have different orbits, the relationships between the elevation difference and terrain aspect varied for different images, and the following corrections, Equations (1) to (5), were used for the elevation differences between the 5 ZY-3 DEM of a-e in Table 2 and SRTM DEM, respectively.
where y is the elevation difference between SRTM DEM and ZY-3 DEM in off-glacier regions, and x is terrain aspect derived from SRTM DEM. After co-registration, the histogram statistics showed that median elevation difference in offglacier regions changed from −64.7, −33.6, −53.6, −22.9, and −28.4 m (left of Figure 3) before coregistration to 0.00, 0.00, 0.00, 0.00, and 0.13 m after co-registration (right of Figure 3) for a, b, c, d, and e images, respectively, and the elevation difference was concentrated between −5 and 5 m. It is concluded that elevation difference in off-glacier regions have been stabilized, and that the preprocessed DEMs were suitable for the estimation of changes in glaciers mass balance. Then, Equations (3) to (7) were applied to the glacierized regions in each image.
Considering the penetrating effect of C band of radar microwaves on snow, SRTM DEM is thought to reflect the glacier ice surface elevation at the end of melting season in October 1999 [33]. The penetration depth of the SRTM C-band radar beam into snow and ice must be considered when assessing glacier elevation changes [13,48,49]. The penetration depth can range from 0 to 10 m depending on snow temperature, density, and water content [48].  Figure 3) before co-registration to 0.00, 0.00, 0.00, 0.00, and 0.13 m after co-registration (right of Figure 3) for a, b, c, d, and e images, respectively, and the elevation difference was concentrated between −5 and 5 m. It is concluded that elevation difference in off-glacier regions have been stabilized, and that the pre-processed DEMs were suitable for the estimation of changes in glaciers mass balance. Then, Equations (3) to (7) were applied to the glacierized regions in each image. Considering the penetrating effect of C band of radar microwaves on snow, SRTM DEM is thought to reflect the glacier ice surface elevation at the end of melting season in October 1999 [33]. The penetration depth of the SRTM C-band radar beam into snow and ice must be considered when assessing glacier elevation changes [13,48,49]. The penetration depth can range from 0 to 10 m depending on snow temperature, density, and water content [48].
Kaab et al. [50] developed a first-order method to correct regionally-based elevation change measurements for the average penetration depth of the C-band signal into ice and snow. In this work, we used an average penetrate depth of 2.4 m in the eastern Karakoram region [50], which has a similar climate to URSRB due to the influence of Westerly winds.

Ice Mass Change
The change in glacier ice volume can be obtained according to the elevation change in glacier surface, which considers glacier area and elevation change in ice surface, and is calculated by where ∆V represents the ice volume change, S is the glacier area, and ∆H is the change in the elevation of the ice.

Glacier Area
The uncertainty in the glacier area extraction is composed of a system error and random error [51]. The system error, such as identifying snow and mountain shadow, can be reduced by selecting high quality remote sensing images and glacier expert interpretation, but the system error is difficult to quantitatively estimate. In the process of glacier outlier extraction, the random error is inevitably affected by the image spatial resolution. The uncertainty in the glacier outlier is considered to be equal to the perimeter of the grid number through the glacier boundary. The glacier outlier from raster to vector includes the echelon form, straight line, and tooth shape forms. The majority of the glacier's boundary is in echelon form, while only a few are in straight line form. The tooth shape boundary is usually at the end of the glacier, and the influence of this type boundary on the number of total glacier grids can be ignored [33]. Therefore, when comparing the total glacier area, the area error σ S caused by the remote sensing image resolution is expressed as where Perimeter is the perimeter of the stepped glacier outlier; Perimeter2 represents the perimeter of the straight line of the glacier outlier; and resolution represents the spatial resolution of the remote sensing image.

Ice Surface Elevation
Generally, the best way to evaluate the elevation error is to use a more precise ground control point. Since ZY-3 has a very high resolution, we could not achieve a higher accuracy of the ground control point, so we used the residual error of the off-glacier region to estimate the error. The error of the sinusoidal model was calculated according to the standard error of each parameter [33], as follows: where k represents the slope of the linear model, σ is the standard deviation of the regression model, and the last two terms on the right give the standard error of the sinusoidal part. The relative errors of the five regions were calculated using Equation (9).
The average error of the elevation difference σ H in the basin was calculated as where n is the number of images, and σ i is the standard error of the sinusoidal curve for the ith image. The ZY-3 image displayed snow cover and mountain shadows in parts of the image, which can cause high elevation overestimation in areas with thick snow, and underestimation in areas in the mountain's shadow. Therefore, these anomalies need to be removed before analyzing the changes. In this study, the threshold scheme was adopted [52,53], and the thread difference between DEMs was ±100 m.

Ice Storage Volume Change
Considering that the change in ice storage volume is affected by the differences between the elevation error and glacial area error, the error of ice volume change ε was estimated as where σ H is the error of elevation change, and σ S the error in the area estimation. The average density of 850 ± 60 kg·m −3 was used as the volume-mass conversion parameter [54], so that the ice storage change could be converted into water equivalent change.

Glacier Area Change
There were 516 glaciers with areas larger than 0.01 km 2 , with a total area of 481.5 km 2 according to the revised CGI. This is less than the 660 glaciers with a total area of 509.9 km 2 in the entire Shule River Basin [19]. The difference is because parts of some glaciers, including Laohugou glacier No. 12, were not above Changma hydrological station, so they were not included. A previous study [19] reported that the total glacier area in the Shule River Basin decreased by about 24.5 km 2 per decade, from the 1970s to 2010, according to the CGI (1970s) and SCGI (2006 to 2010). Comparatively, our results suggest that the glacier area decreased by about 21.9 km 2 per decade, from the 1970s to 2008, which is consist with Sun's [19] results.
Between   It is interesting to find that although the 6 glaciers with areas greater than 8 km 2 remained unchanged, the median glacier area obvious decreased ( Figure 5). The number of medium-size glaciers with areas between 3 and 8 km 2 decreased. The glacier areas between 1 and 3 km 2 displayed similar patterns with those between 3 to 8 km 2 . The number of glaciers with area less than 1 km 2  It is interesting to find that although the 6 glaciers with areas greater than 8 km 2 remained unchanged, the median glacier area obvious decreased ( Figure 5). The number of medium-size glaciers with areas between 3 and 8 km 2 decreased. The glacier areas between 1 and 3 km 2 displayed similar patterns with those between 3 to 8 km 2 . The number of glaciers with area less than 1 km 2 decreased from 417 to 393 from 2000/2001 to 2014/2015. These results indicate that the change in large glaciers is mainly displayed as area shrinkage, where 35 glaciers with an area less than 1 km 2 disappeared. The average size of the 35 glaciers that disappeared was ~0.24 km 2 in 2000/2001. The elevation of these glaciers ranged from 4422 to 5526 m a.s.l., indicating that the small glaciers are affected by multiple factors, excluding elevation.

Spatial Pattern of Glacier Area Change
According to the spatial patterns of glacier change, the whole URSRB can be divided into (I) Northwest, (II) Middle, and (III) Southeast subregions ( Figure 6). The

Spatial Pattern of Glacier Area Change
According to the spatial patterns of glacier change, the whole URSRB can be divided into (I) Northwest, (II) Middle, and (III) Southeast subregions ( Figure 6).

Glacier Surface Elevation Change
Due to the fact that 5 ZY-3 image pairs were used to calculate the ice surface elevation changes in the study, and there are still some areas which were not covered by the 5 ZY-3 image pairs in the URSRB, we therefore use the 5 subregions (a-e) which are in line with each ZY-3 image pair to investigate the characteristics of glacier surface elevation change in the URSRB.
The glacier surface elevation in the URSRB decreased from 1999 to around 2013, and large differences were found between the a, b, c, d, and e subregions that correspond with each ZY-3 image (Figure 7, Table 3 It is pertinent to note that the change rates of glacier areas have not always been in line with those of ice surface elevation (Table 3). For instance, in the ice elevation change in d subregion was only −0.184 m·a −1 , where the middle has changed in glacier area by −0.87% a −1 . The large variations in ice surface elevation in each subregion reflect the complexity of glacier mass change. The total glacier area in region (II) was about 342.42 km 2 , accounting for 78.7% of the total glacier area in the URSRB. In this region, the glacier size was relatively large, among which, 6 glaciers displayed areas larger than 10 km 2 . The average elevation of glaciers in this region was 4983 m a.s.l. The glacier shrinkage mainly occurred at the terminus of the glacier in the elevation ranging from 4700 to 5100 m a.s.l., and the average area was reduced by about 5.

Glacier Surface Elevation Change
Due to the fact that 5 ZY-3 image pairs were used to calculate the ice surface elevation changes in the study, and there are still some areas which were not covered by the 5 ZY-3 image pairs in the URSRB, we therefore use the 5 subregions (a-e) which are in line with each ZY-3 image pair to investigate the characteristics of glacier surface elevation change in the URSRB.
The glacier surface elevation in the URSRB decreased from 1999 to around 2013, and large differences were found between the a, b, c, d, and e subregions that correspond with each ZY-3 image (Figure 7, Table 3    It is pertinent to note that the change rates of glacier areas have not always been in line with those of ice surface elevation (Table 3). For instance, in the ice elevation change in d subregion was only −0.184 m·a −1 , where the middle has changed in glacier area by −0.87% a −1 . The large variations in ice surface elevation in each subregion reflect the complexity of glacier mass change.

Uncertainty
Uncertainty in the delineation of glacier outlines is the result of positional and processing errors [55] in which seasonal snow, cloud, and debris cover complicate the precision of the mapping [56]. There are limited debris-covered glaciers in URSRB. The accuracy of the outlines in this study was assessed by comparisons with a previous study by Sun et al. [19], and the results in CGI and SCGI. The uncertainty due to the resolution of the remote sensing data was also evaluated.
The penetration depth of the SRTM C-band radar beam contributed to some uncertainties when the SRTM DEM was used for geodetic mass-balance calculations. A previous study [57] indicates that the penetration depth decreases as the temperature and water content of the surface snow increases. When the penetration depth was 1.4 m, the ice surface elevation decrease rate was 0.31 m·a −1 , which increased to 0.46 m·a −1 when the penetration depth increased to 3.4 m.
Another challenge faced in glaciers studies is data voids in SRTM DEM, especially in the glacier accumulation area. In this study, information was available for elevation change from 4300 to 5500 m a.s.l. Therefore, we compared the unfilled version SRTM DEM and SRTM DEM4.1, found that the area of data voids was too small to affect the ice surface elevation significantly (1.0% of the area above 5000 m a.s.l.), so this issue was neglected.

Changes in Glacier Area and Ice Surface Elevation
The glacier area in the URSRB decreased by 57.5 km 2 (11.91% relate to CGI) from 2000/2001 to 2014/2015, with an average annual shrinkage rate of 0.79% a −1 , which is about 1.7 times the shrinkage rate from CGI to SCGI. This result indicates that the glacier shrinkage accelerated after 2000, which is consistent with studies on Tuanjiefeng region near to the URSRB [33], where the glacier area decreased 0.20% a −1 between 2000-2006, and 0.45% a −1 during 2006-2010. Thist also agrees with the shrinkage of Laohugou glacier No. 12, where the area decreased by 0.008 km 2 ·a −1 from 1957 to 1994, 0.043 km 2 ·a −1 from 1994 to 2000, and 0.041 km 2 ·a −1 from 2000 to 2007 [29]. It is also similar to Yanglonghe glaciers No. 1 and No. 5 in the middle Qilian Mountain, where the glacier experienced accelerated retreat since 1999 [30]. A similar pattern was also found in Kangri Karpo region, where the glacier area decreased by 0.41% a −1 from 1980 to 2000, and decreased by 0.52% a −1 from 2000 to 2015 [58]. According to the spatial pattern of glacier change for the entire Qilian Mountain range [19], the glacier shrink rate in the URSRB is relatively lower than other parts of the Qilian Mountains, which indicates that accelerated glacier shrinkage is probably more pronounced in other regions of the Qilian Mountains.
In terms of ice surface elevation changes in the Qilian Mountains, an average decrease of 0.64 m·a −1 was observed for nine glaciers in the Ningchan and Shuiguan River regions in the eastern Qilian Mountains from 1972 to 2010 [41]. Average elevation decreases of 0.4 ± 0.22 m·a −1 and 0.33 ± 0.22 m·a −1 were measured in the ablation area of Yanglonghe glacier No. 1 and No. 5 in the eastern Qilian Mountains from 1956 to 2007 [30]. During 2010-2015, an elevation decrease of 0.9 ± 0.5 m a −1 (water equivalent) was observed for Ningchan glacier No. 1 [32]. These results demonstrate that the mass loss of both a single glacier and multiple glaciers in one region has accelerated in the 21st century.
Xu et al. [33] determined the average ice thinning rate of 0.21 ± 0.04 m·a −1 for the glacier near Tuanjiefeng based on SRTM DEM and topographic mapping during 1966-1999. Similarly, our study suggests that the average ice thinning on the URSRB was about 0.38 m ± 0.05 a −1 during 1999-2013, which also indicates that glacier thinning has accelerated. This trend was also observed in the Kangri Karpo Mountain region in southeastern TP, where glaciers thinned on average by 0.24 ± 0.16 m w.e.a −1 during 1980-2000 and 0.71 ± 0.10 m w.e.a −1 during 2000-2014 [59]. This pattern of accelerated loss matches the global record, as well [60].
The recent mass loss of glaciers in the URSRB is larger than the loss of 0.24 ± 0.22 m·a −1 in central Tianshan Mountain, based on SRTM DEM and TanDEM-X during 1999 to 2012 [11]. In addition, glacier mass loss in URSRB is less than the 0.71 ± 0.10 m·a −1 that occurred in the Kangri Karpo Mountain region in southeastern TP from 2000 to 2014 [58], which is consistent with the observed glacier area change pattern in western China [59].

Glacier Response to Climate Change
The annual average air temperature has significantly (p = 0.01) increased by 0.465 • C per decade between 1971 and 2013. The average air temperature during the snow and ice melting season (May to September) has significantly (p = 0.01) increased by 0.433 • C per decade. No significant (p = 0.05) trend was found for annual precipitation (Figure 8), suggesting that the recent glacier area loss and ice storage volume losses are mainly controlled by increased air temperature. The air temperature increase is more than twice than the global average of warming (0.12 • C per decade, 1951-2012) [7], which reveals that the URSRB has experienced pronounced warming, like many other mountain areas in the world [61].

Glacier Response to Climate Change
The annual average air temperature has significantly (p = 0.01) increased by 0.465 °C per decade between 1971 and 2013. The average air temperature during the snow and ice melting season (May to September) has significantly (p = 0.01) increased by 0.433 °C per decade. No significant (p = 0.05) trend was found for annual precipitation (Figure 8), suggesting that the recent glacier area loss and ice storage volume losses are mainly controlled by increased air temperature. The air temperature increase is more than twice than the global average of warming (0.12 °C per decade, 1951-2012) [7], which reveals that the URSRB has experienced pronounced warming, like many other mountain areas in the world [61].

Conclusions
The following conclusions can be made from the investigation of recent glacier changes in URSRB using multisource remote sensing data.

Conclusions
The following conclusions can be made from the investigation of recent glacier changes in URSRB using multisource remote sensing data.
Thirty-five glaciers disappeared from 2000 to 2015, with a total glacier area shrinkage of 57.5 ± 2.68 km 2 (about 11.9% of the area in the 1970s), which equals an annual decrease of 3.83 ± 0.18 km 2 ·a −1 (0.79 ± 0.04% a −1 ). From 2004/2006 to 2008/2009, glacier shrinkage was most pronounced in which it decreased by 9.21 km 2 ·a −1 (1.90% a −1 ). The glacier shrink rate in the URSRB is relatively lower than other parts of the Qilian Mountains, which indicates that glacier shrinkage is probably more pronounced in other regions of the Qilian Mountains. (2) The average glacier surface elevation of the URSRB from 1999 to 2013 was reduced by about 4.98 ± 0.6 m, which is equal to a mass loss of 0.383 ± 0.046 m·a −1 . The ice storage loss has accelerated since 1999. This pattern of accelerated loss is in agreement with the global records reported. The mass loss of glaciers in the URSRB is larger than that in Tianshan, and less than that in southeastern TP, which is consistent with the observed glacier area change pattern in China.
The pronounced glacier shrinkage and mass loss probably had a significant impact on the annual volume and seasonal distribution of river runoff in the URSRB, increasing the difficulty of water supply management of the URSRB and downstream of Shule River. The potential changes in glacier and the subsequent impact on river runoff in URSRB require further study.