Glacier Changes and Their Linkage to the Climate-Topographic Context in the Borohoro Mountains, Tian Shan 1977–2018

Against the backdrop of climate change and socio-ecological sustainability, studying glacier changes provides essential knowledge to the basic water needs and security for regions and populations under such threats, such as Central Asia. Little attention has focused on glaciers in the northern periphery of the Chinese Tian Shan. This study aims to map a recent glacier inventory and examine the glacier area shrinkage and surface elevation change for the central massif of the Borohoro Mountains in the past 41 years. Using declassified Hexagon images (1977), Landsat 5 TM (1994 and 2007), Sentinel 2A (2018) and altimetry data from the Ice, Cloud and Land Elevation Satellite (ICESat) over 2003–2009 with the 30-m Shuttle Radar Topography Mission (SRTM) digital elevation model, multi-temporal glacier fluctuations and the influence of topographic and climatic factors were investigated. Results show that the glacier area decreased from 287.5 ± 8.2 km2 in 1977 to 215.8 ± 4.1 km2 in 2018, at a rate of 0.61 ± 0.01% year−1. Glacier disintegration has led to a gradual increase in the number of glaciers and reached 224 glaciers in 2018. The shrinkage was at the highest rate during the 1994–2007 period and the smallest during 1977–1994. Glacier size, hypsometry, and median, maximum, and range of elevation are the most significantly correlated parameters with the relative area change. The surface elevation changes from two of the largest glaciers revealed a stronger thinning on the southern slope compared to the northern slope. These observations of glacier loss are primarily driven by the marked warming trend since the 1970s and confirmed with the overall pattern of glacier retreat in the Tian Shan from previous studies.


Introduction
As reported in the Intergovernmental Panel on Climate Change (IPCC) special report, mountain surface air temperatures show warming at an average rate of 0.3 ± 0.2 • C/decade, outpacing the global average warming rate of 0.2 ± 0.1 • C/decade [1]. This enhanced warming has induced changes to elements and processes in the atmosphere, hydrosphere and biosphere in the mountain environment, and perhaps most prominently, in the cryosphere (glaciers, snow cover, and permafrost) that has responded sensitively to the changing climate. Scientific evidence has reached a coherent conclusion that global warming has caused an overall mountain glacier recession since the Little Ice Age (~1850AD), with an accelerated rate entering the 21st century [2][3][4]. Among major mountain areas in the world, the High Mountain Asia (HMA) represents the highest concentration of glaciers outside of the polar regions and exemplifies heterogeneous glacier mass budget with a contrasting pattern within the region [5]. The HMA contains the headwater of many large rivers such as Brahmaputra, Indus, Syr Darya, Amu Darya, Mekong, Ganges, Yellow and Yangtze Rivers. Knowing the glacier fluctuations in this key region is essential because of their importance to the water supply for over

Study Area
The Tian Shan, located in the hinterland of Eurasia, is a gigantic mountain range that originated from one of the largest intra-continental orogens in the world in the Paleozoic (~360-205 Ma) [33]. The modern landscapes were built by episodic faulting and uplift during the Cenozoic (~65 Ma) but have continued to be modified up until today [34]. To its south lie the Tibet-Pamir Plateau and the Tarim Basin, and to the north are the Junggar Basin and the Altai Mountains. The mountain ranges across the Tian Shan are characterized by diverse geology, topography, climate, and vegetation in its intervening mountains and valleys. The Borohoro Mountains are located at the northwestern periphery of the eastern Tian Shan (also the Chinese Tian Shan), with greener forests/steppes in the Ili Valley to the south and the drier oasis-desert ecosystem to the north (Figure 1). The bedrock of the mountain is mainly made of Late Carboniferous (305-300 Ma) granites and metamorphics, from the Borohoro pluton [33]. The glacier-capped mountain elevations range from 3000 m to 4700 m in the studied region. Ridgelines of the Borohoro divides two major river basins in Xinjiang, China: the Ili River Basin to the south and the Junggar Inland Basin to the north [18]. Meltwater from glaciers contributes significantly to the runoff of hundreds of tributaries in these basins. The Jinghe River, a glacier-fed river adjacent to the Borohoro, was estimated to receive 74% of annual runoff that occurs in summer, a time when glacier ablation is high [35]. Indeed, glaciers comprise a key component in sustaining lives and the economy for this water-deficient region of Central Asia.
Water 2020, 12, x FOR PEER REVIEW 3 of 23 periphery of the eastern Tian Shan (also the Chinese Tian Shan), with greener forests/steppes in the Ili Valley to the south and the drier oasis-desert ecosystem to the north (Figure 1). The bedrock of the mountain is mainly made of Late Carboniferous (305-300 Ma) granites and metamorphics, from the Borohoro pluton [33]. The glacier-capped mountain elevations range from 3000 m to 4700 m in the studied region. Ridgelines of the Borohoro divides two major river basins in Xinjiang, China: the Ili River Basin to the south and the Junggar Inland Basin to the north [18]. Meltwater from glaciers contributes significantly to the runoff of hundreds of tributaries in these basins. The Jinghe River, a glacier-fed river adjacent to the Borohoro, was estimated to receive 74% of annual runoff that occurs in summer, a time when glacier ablation is high [35]. Indeed, glaciers comprise a key component in sustaining lives and the economy for this water-deficient region of Central Asia.  Due to its central location on the Eurasia landmass, a typical arid continental climate dominates the Borohoro. But, like most places in the Tian Shan, the striking topographic relief gives a variety of vertical climatic zones and associated landscape belts that ultimately transition into nival belts at the highest reach [36]. The prevailing wind of the mid-latitude westerlies delivers moisture from the Caspian Sea to the mountains. The southern flank of the Borohoro receives relatively high precipitation because the westward opening of the Ili Valley facilitates the entrance of moisture-laden air masses. The northern flank connects to a rain shadow region formed by the Dzungarian Alatau mountains in the west, and thus has less precipitation. In this range of the Tian Shan, precipitation occurs mostly in summer caused by convective events and has been increasing over the past decades [8,11]. The annual temperature variation is pronounced, creating strong seasonality. Under the control of the Siberian High, winter is featured with cold-dry climate. According to the closest meteorological station, the Jinghe station, which is situated at 320 m a.s.l. and approximately 100 km northwest of the studied glaciers (Figure 1b), the mean annual precipitation and average air temperature are 251.9 mm and 7.8 • C, respectively [37]. However, detailed climatological data are sparse at higher elevations.

Satellite Data
Three criteria were applied when selecting satellite images for the study area: (1) the acquisition time around the end of the ablation season (late August and September), (2) minimal seasonal snow, and (3) minimal cloud cover. Only a limited number of images are qualified. Four images were selected from the declassified Hexagon imagery, Landsat 5 TM, and Sentinel 2A collections, covering a temporal range of four decades and spaced with similar time intervals (1977-1994-2007-2018) (Table 1). They were obtained from the USGS web server (earthexplorer.usgs.gov). The 30 m Shuttle Radar Topography Mission (SRTM) digital elevation model (DEM) was also downloaded from the USGS for orthorectification and extraction of topographic parameters (Table 1). The Sentinel 2A image of 20 August 2018 shows the best quality at the latest time and thus was used for producing the most recent inventory of glaciers in the Borohoro. It was also used as a base image for co-registration of other images. The short repeat interval (10 days at the equator) and high resolution (10 m in the visible and near-infrared bands, 20 m in the shortwave infrared bands) of the scenes has shown advanced glacier extraction since its launch in 2015 [38,39]. Two Landsat scenes (30 m resolution) are both from Landsat 5 TM Level-1TP collection, acquired on 22 August 1994 and 10 August 2007, respectively. These radiometrically and geometrically corrected images, with minimum influence from snow and cloud and in the similar time of the ablation season, can greatly help to reduce uncertainties in the extraction of glaciers. The co-registration root-mean-square errors (RMSE) to the base image for 1994 and 2007 Landsat 5 images, respectively, are 15.6 m and 14.5 m using 32 and 35 ground control points (GCPs). The Hexagon KH-9 image is from the declassified archive of the U.S. reconnaissance satellite KH-9 operated from 1973 to 1980, and some images are previously scanned and available to download free of charge by USGS. The selected film was acquired on 26 September 1977, extending the temporal coverage to the 1970s comparable to the time when the first Chinese Glacier Inventory [17] was compiled. The size of the film is 9 × 18 inch, covering a footprint at 125 × 250 km 2 with a resolution of~8 m. Because of no terrain correction, the Hexagon image needs to be orthorectified and co-registered. A subset of the film was generated for the Borohoro study area, and following a commonly used two-steps approach [26,40,41], the subset was applied with the first step of projective transformation using GCPs and the SRTM DEM and then the second step of a spline adjustment in ArcGIS. The RMSE of the co-registration using 60 GCPs is calculated as 1 m. All of the above images were cropped to the study area and projected into the World Geodetic System (WGS) 1984 Universal Transfer Mercator (UTM) Zone 44N in ArcGIS Pro.

Glacier Mapping
A semi-automated approach was applied to extract glacier boundaries and assisted with the use of band combinations (11-8-4 for Sentinel 2A MSI and 7-4-3 for Landsat 5 TM) and Google Earth imagery. As the 20 August 2018 Sentinel 2A image is at a higher resolution (10 m) and contains almost no impact from snow and cloud, it was used to generate a baseline inventory of glaciers in the Borohoro at the first step. After resampling the 20 m shortwave infrared (SWIR) band (Band 11) to 10 m using bilinear interpolation, the MIS4/MIS11 (red/SWIR) band ratio was applied to effectively and efficiently distinguish glaciers from other surface features [38]. A threshold value of 1.8 for the ratio was determined to be the best after visually inspecting the results of multiple testing values. Then a median filter of 3 × 3 kernel size was used to smooth out the isolated pixels and fuzzy boundaries generating a binary image that was converted to a vector layer. The above image processing steps were performed in ERDAS Imagine ® (Hexagon AB, Stockholm, Sweden). Manual adjustments are necessary to remove or modify the misclassified polygons and were conducted with the false color composite image underlaid in ArcGIS Pro. Major sources of misclassification come from proglacial lakes, shadows and glacial fronts where debris-cover exists occasionally. To separate polygon vectors into individual glacier entities, the mountain ridgelines were first derived using the drainage basin function with a reversed SRTM DEM, and then they were manually corrected especially at the upper part of glaciers to ensure the consistency with the outlines of the glacier inventory data (e.g., GLIMS www.glims.org/MapsAndDocs/guides.html). Glacier outlines of 1977, 1994, and 2007 were manually delineated from the Hexagon and Landsat 5 TM scenes, assuming the upper glacier divides did not change much during the investigated period. Major challenges in manual mapping occurred with the shadow problems and over-exposure (in Hexagon film), and were visually checked back and forth in Google Earth and the Sentinel 2A imagery.
Topographic information of each glacier entity was extracted by applying analysis tools in ArcGIS Pro. Glaciers of an area <0.02 km 2 were eliminated to reduce the potential error of misinterpreted snow patches. In cases of a separated glacier during the investigated time period, multiple features were treated as one joint feature in correspondence with the pre-fragmentation extent. Parameters of each glacier, including area, minimum, maximum, and median elevation, elevation range, mean slope, mean aspect, and hypsometric integral (HI), were calculated based on the 2018 glacier outlines obtained and the 30 m SRTM DEM. In particular, the median elevation is the altitude dividing the glacier area (map area) into two equal areas; the mean aspect is calculated as the directional mean of the aspect Water 2020, 12, 1502 6 of 22 from each raster cell within the glacier boundary. Descriptive statistics and quantitative relationships between the change of the glacier area and these parameters were explored.

Estimation of Elevation Change
To detect the interannual changes of the glacier surface elevation, data from the Ice, Cloud and Land Elevation Satellite (ICESat) were used and obtained from the U.S. National Snow and Ice Data Center (NSIDC; https://nsidc.org/data/icesat/data.html) ( Table 1). The Geoscience Laser Altimeter System (GLAS) on board measured time series of surface elevation along its repeat-track orbits two or three times a year for the period of 2003-2009 [42]. Laser pulses were sent to sample elevations within a diameter of 70 m footprint spaced 172 m along the track [43]. In HMA, studies have been applied to GLAS/ICESat data to understand the elevation/mass/volume changes of glaciers [15,44,45]. In the studied Borohoro region, one track passed over two of the largest glaciers, as shown in Figure 1. But since the repeated track patterns did not overlap exactly, the comparison of ICESat elevations can only be achieved by using a reference elevation. In particular, the GLAS/ICESat Global Land Surface Altimetry HDF5 data (GLAH14, release 33) were used to compare with the 1 arcsec (~30 m) SRTM DEM, which represents elevations acquired in February 2000. Before differencing the two elevations, the ICESat values, which are referred to the TOPEX/POSEIDON Ellipsoid, need to be converted to the same ellipsoid of the SRTM data, which is the WGS84 Ellipsoid based on EGM96 Geoid, using the following formula proposed by [46]: The 0.7 m is the vertical offset estimated between the two ellipsoids. When using the calculated ICESat elevation to subtract the SRTM elevation, dh can be obtained to represent the surface elevation change over time between the DEM source date (2000) and the ICESat acquisition time, with negative indicates thinning, positive indicates thickening and zero means no change. Elevation differences that exceed a threshold of 100 m (|dh| > 100 m) were excluded as outliers, as suggested in previous studies [44,47]. The filtered on-glacier data were further divided into two zones, the accumulation zone and the ablation zone, roughly separated by the transient snowline using the end of ablation season image of the 2018 Sentinel 2A [48].

Climate Data
To assess the impact of changing climate on glaciers, the 0.5 • × 0.5 • gridded, monthly University of Delaware air temperature and precipitation V5.01 dataset (provided by esrl.noaa.gov/psd) [49] was used because no station-based climate data are available in the study area. The latitude/longitude grid that overlaps most with the study area is centered at 43.75 • N, 83.75 • E (Table 1). Time series of monthly and annual average air temperature (T) and total precipitation (P) from 1900 to 2017 were retrieved for examining the climatic trends concerning the glacier changes.

Uncertainty Assessment
Sources of uncertainty and error in mapping glacier outlines are embedded in the mapping methods (either automated or manual) and issues such as co-registration of images and the use of the various satellite datasets with different resolutions and different clouds, seasonal snow, shadow, and debris cover conditions. It is crucial to include the error when reporting results [50,51]. Previous studies have suggested different methods to perform error analysis, including both qualitative and quantitative approaches [50]. The lack of higher resolution data in this study area prevented a qualitative comparison and the identification of misinterpretation in the derived outlines. But a statistical extrapolation method, the buffer method [41,52], was adopted here to estimate the uncertainty for each glacier. A buffer size of 5 m (half of a pixel) was applied to the 2018 Sentinel 2A glaciers; half of the RMSE values (15.6 m and 14.5 m) from the co-registration of Landsat 5 images to the Sentinel 2A base image were used for 1994 and 2007 outlines, respectively; and an uncertainty value of 10 m from the literature [41] was used for the Hexagon results. After calculation, the precision was determined as~1.9%,~2.6%, 2.5%, and~2.9% for 2018 Sentinel 2A, 2007 Landsat 5, 1994 Landsat 5, and 1977 Hexagon glaciers, respectively. When reporting multi-temporal changes of glacier areas, the error propagation (e) as in Equation (2) was applied to estimate the uncertainty from different sources (e 1 and e 2 ): Uncertainties in the ICESat altimetry data primarily come from the system inconsistency between campaigns, the atmospheric interference, and possible terrain roughness within footprints [50]. Studies have shown that elevation biases are of centimeters to decimeter level during the lifetime of ICESat, which can be neglected [47]. By filtering the dh with a cutoff threshold of ±100 m, it helps to reduce the outliers of those footprints on clouds. Moreover, the surface heights were sampled on the large glaciers, as in this case, they may not be as variable as over the small glaciers because of the less rugged terrain.

Glacier Characteristics
Based on the mapping from the 2018 Sentinel 2A base image, I identified 224 glaciers larger than 0.02 km 2 , which had a total area of 215.8 ± 4.1 km 2 in the study area; 70% of the glaciers are smaller than 0.5 km 2 , while only 10 out of 224 glaciers are larger than 5 km 2 , comprising 47.5% of the total area ( Table 2; Table S1). The largest one reached 17.8 km 2 , and along with the next three largest glaciers (15.5 km 2 , 14.3 km 2 , and 13.8 km 2 ), they form the central massif of the area. The average glacier size is approximately 1.0 km 2 , which is similar to other published results in regions across the Tian Shan, e.g., 0.8 km 2 in the Jinghe River Basin in 2004 [37], 0.8 km 2 in Ala Archa in 2003 [27], 1.0 km 2 in the eastern Terskey Ala-Too in 2003 [6], and 1.5 km 2 in the Sary-Jaz River Basin in 2010 [53]. Small glaciers with a size <1 km 2 outnumbered other size classes in this mid-latitude mountain range ( Figure 2).
The lowest altitude that a glacier tongue reached is 2345 m a.s.l., while the highest peak of a glacier is 4680 m a.s.l. The hypsometric distribution (Figure 2a inset) suggests that most glacierized areas are concentrated in elevations ranging from 3600 m-3800 m a.s.l. As the median elevation is commonly used to estimate the balanced-budget equilibrium line altitude (ELA) [54], it is adopted in describing the altitudinal location of individual glaciers and later used in the topographic analysis, instead of using the mean elevation. The median elevation ranges from 3316 m to 4527 m a.s.l., with an average of 3797 m a.s.l. (Table 2). Across different glacier size classes, the smallest glaciers (0-0.5 km 2 ) have a wide spread of their median elevations but a narrow elevation range, while the elevation range of the large glaciers (>5 km 2 ) is much wider but their median elevation is similarly distributed (Figure 2a). The lowest altitude that a glacier tongue reached is 2345 m a.s.l., while the highest peak of a glacier is 4680 m a.s.l. The hypsometric distribution (Figure 2a inset) suggests that most glacierized areas are concentrated in elevations ranging from 3600 m-3800 m a.s.l. As the median elevation is commonly used to estimate the balanced-budget equilibrium line altitude (ELA) [54], it is adopted in describing the altitudinal location of individual glaciers and later used in the topographic analysis, Other topographic parameters also reflect the characteristics of glaciers in the region. The mean aspect represents the orientation direction of each glacier, expressed using the azimuth angle. When the aspect is aggregated into eight compass directions, it is clear that most glaciers are generally facing north, including northwest (315 • ), north (0 • ), and northeast (45 • ) ( Figure 2b). However, the largest glaciers are generally oriented toward the east (91.7 • ) ( Table 2). The mean slope of all glaciers is 25.2 • , but almost half of them are steeper than 30 • (Figure 2c). The slope becomes gentler as the glacier size increases (Table 2), similar to other studies [40,55]. The hypsometric integral (HI) ranges 0-1 and indicates the shape of the hypsometric curve, with large values for a convex shape and small values for a concave shape. The HI values are normally distributed with a mean of 0.5, and the average HI is nearly the same across size classes ( Figure 2d and Table 2).

Temporal Changes of Glacier Areas
The temporal sequence of glacier area changes has been investigated for 1977, 1994, 2007, and 2018, and heterogeneous extent and rate of change are found during different periods of these past 41 years ( Figure 3). From 1977 to 2018, the total area of glaciers decreased from 287.5 ± 8.2 km 2 to 215.8 ± 4.1 km 2 , by 71.6 ± 9.2 km 2 (24.9 ± 3.2%). During three sub-periods, it decreased 25.0 ± 10.5 km 2 (8.7 ± 3.7%) from 1977 to 1994, then 30.8 ± 8.9 km 2 (11.7 ± 3.4%) from 1994 to 2007, and 15.8 ± 7.3 km 2 (6.8 ± 3.1%) from 2007 to 2018. The rate of recession is the highest (0.90 ± 0.07% year −1 ) for 1994-2007 and slowest (0.51 ± 0.03% year −1 ) for 1977-1994. The number of observed glaciers increased from 200 in 1977 to 224 in 2018, primarily due to the disintegration from large ice bodies to smaller entities ( Table 3). As the glacial extent has been shrinking in the time series, new developments of land features, such as exposed rocks, water channels, debris cover, and proglacial lakes, can be observed with a comparison of satellite images most commonly in areas adjacent to glacier boundary (Figure 4a). The sequence of images also revealed that at some places small glaciers gradually disappeared in the past decades (Figure 4b).
Water 2020, 12    The absolute and relative area changes were further examined by glacier size classes as defined in 1977 ( Figure 5). In cases of disintegration, the sum of the fragmented glacier areas in the later period was used to compare with the area of the glacier in an earlier year. It is shown that the large glaciers (>5 km 2 ) had the greatest absolute area shrinkage (Figure 5a), while their relative percentage of change was the smallest among all size classes during any investigated time period (Figure 5b). The small glaciers (<1 km 2 ) collectively contributed little to the total area loss, but the relative change was averaged to more than 50% between 1977 and 2018. They also expressed high variability in relative area loss, with a few with 100% loss in the smallest class (<0.  The absolute and relative area changes were further examined by glacier size classes as defined in 1977 ( Figure 5). In cases of disintegration, the sum of the fragmented glacier areas in the later period was used to compare with the area of the glacier in an earlier year. It is shown that the large glaciers (>5 km 2 ) had the greatest absolute area shrinkage (Figure 5a), while their relative percentage of change was the smallest among all size classes during any investigated time period (Figure 5b). The small glaciers (<1 km 2 ) collectively contributed little to the total area loss, but the relative change was averaged to more than 50% between 1977 and 2018. They also expressed high variability in relative area loss, with a few with 100% loss in the smallest class (<0.

Glacier Change and Topographic Factors
Topographic characteristics of individual glaciers, such as size, elevation, slope, aspect, and hypsometry, may influence the magnitude of glacier changes, e.g., [56][57][58][59]. As the change was defined  The absolute and relative area changes were further examined by glacier size classes as defined in 1977 ( Figure 5). In cases of disintegration, the sum of the fragmented glacier areas in the later period was used to compare with the area of the glacier in an earlier year. It is shown that the large glaciers (>5 km 2 ) had the greatest absolute area shrinkage (Figure 5a), while their relative percentage of change was the smallest among all size classes during any investigated time period (Figure 5b). The small glaciers (<1 km 2 ) collectively contributed little to the total area loss, but the relative change was averaged to more than 50% between 1977 and 2018. They also expressed high variability in relative area loss, with a few with 100% loss in the smallest class (<0.

Glacier Change and Topographic Factors
Topographic characteristics of individual glaciers, such as size, elevation, slope, aspect, and hypsometry, may influence the magnitude of glacier changes, e.g., [56][57][58][59]. As the change was defined

Glacier Change and Topographic Factors
Topographic characteristics of individual glaciers, such as size, elevation, slope, aspect, and hypsometry, may influence the magnitude of glacier changes, e.g., [56][57][58][59]. As the change was defined as the loss from 1977 glaciers to the 2018 glaciers, the total number of 200 glaciers in 1977 was used to extract the topographic factors, including size, minimum, maximum, and median elevation, elevation range, mean slope, hypsometric integral, and the sine and cosine terms of the aspect [57,60]. It should be noted that for cases of disintegrated glaciers from 1977 to 2018, a merge rule of summing was applied to obtain the correspondent area in 2018 (a merge rule of averaging for other parameters); while for cases of disappeared glaciers from 1977 to 2018, a 100% relative area change was assigned.
For all glaciers in the study area, the relative reduction in area has a significant negative correlation with the size, indicating that small glaciers tend to lose more than the large ones (Figure 6a). The relative area change is also negatively correlated with the median elevation, elevation range, and hypsometry (Figure 6b,c,e). The relationship between mean slope and the percent of change is slightly negative, but not significant (Figure 6d). East-facing glaciers experienced less proportional shrinkage compared to the glaciers with other orientations (Figure 6f). Pearson correlation analysis was performed on the relative area change (∆A), the natural logarithm transformed area (ln(A), to normalize the distribution), and eight other factors as mentioned above. Results show that ln(A), median elevation, maximum elevation, elevation range, and hypsometry were negatively correlated with ∆A and significant at the 99% confidence interval (p < 0.01), while other factors are not significant (p > 0.1) (Figure 7). When examining the interaction between factor pairs, some strong collinearities are shown. For example, the area was significantly positively correlated with the elevation range (r = 0.75, p < 0.01) and maximum elevation (r = 0.28, p < 0.01), and the median elevation was strongly correlated with other elevation-dependent variables (min., max., and range), slope, and hypsometry (p < 0.01) (Figure 7). was applied to obtain the correspondent area in 2018 (a merge rule of averaging for other parameters); while for cases of disappeared glaciers from 1977 to 2018, a 100% relative area change was assigned.
For all glaciers in the study area, the relative reduction in area has a significant negative correlation with the size, indicating that small glaciers tend to lose more than the large ones ( Figure  6a). The relative area change is also negatively correlated with the median elevation, elevation range, and hypsometry (Figure 6b,c,e). The relationship between mean slope and the percent of change is slightly negative, but not significant (Figure 6d). East-facing glaciers experienced less proportional shrinkage compared to the glaciers with other orientations (Figure 6f). Pearson correlation analysis was performed on the relative area change (∆A), the natural logarithm transformed area (ln(A), to normalize the distribution), and eight other factors as mentioned above. Results show that ln(A), median elevation, maximum elevation, elevation range, and hypsometry were negatively correlated with ∆A and significant at the 99% confidence interval (p < 0.01), while other factors are not significant (p > 0.1) (Figure 7). When examining the interaction between factor pairs, some strong collinearities are shown. For example, the area was significantly positively correlated with the elevation range (r = 0.75, p < 0.01) and maximum elevation (r = 0.28, p < 0.01), and the median elevation was strongly correlated with other elevation-dependent variables (min., max., and range), slope, and hypsometry (p < 0.01) (Figure 7).       The gradient in the elevation change rate is also demonstrated by differentiating the ablation zone at the lower elevations and the accumulation zone at the higher elevations, separated by the transient snow line as an approximation of the ELA (Figure 8). In general, the dh values were averaged as positives in the accumulation zone, while they are mostly negative in the ablation zone for both glacier ID-27 and glacier ID-16 during the period of 2003-2009. The overall trends through the time series went downslope for both zones of the glaciers, basically translating into more loss of thickness in the ablation zone and less accumulation in the accumulation zone as time progressed (Figure 8). It should also be noted that the south-facing glacier ID-16 showed a more significant absolute thickness loss compared to the north-facing glacier ID-27, at both ablation and accumulation zones. Considering the inter-annual elevation trends derived from all dh values with a linear fitting, the accumulation zone of the ID-27 experienced a higher thinning rate (−0.783 m year −1 ) than that of the ID-16 (−0.031 m year −1 ), while the ablation zone of the ID-16 is reducing its thickness at a higher rate (−1.064 m year −1 ) than the ID-27 (−0.742 m year −1 ). However, it would be difficult and farfetched to compare the mass budget changes of the two glaciers in the period because the datasets of the ICESat sampling are not balanced between the accumulation and the ablation parts and not extensive enough to cover the full scope of glaciers.

Climatic Trends
Because no weather station data are available in the study area, the most recent version of the University of Delaware dataset of air temperature and precipitation at the grid (43.75° N, 83.75° E) that covers the target region was used to reflect the trends of climate between 1900 and 2017 [49]. As shown in Figure 9a, the mean annual temperature (MAT) was the coldest at the beginning of the twentieth century and then fluctuated with positive and negative trends throughout the century. The overall trend of MAT is increasing, and the 1977-2017 period, which overlaps closely with the investigated glacier change period, showed a significant positive trend coefficient (0.03 °C year −1 ) using linear regression (Figure 9a). Regarding the total annual precipitation (TAP), it can be observed that the variability was higher for the first half of the century than the second half (Figure 9b). The 10-year smoothed average helped to identify a gently increasing trend of TAP since the 1970s, and the trend analysis using linear regression for 1977-2017 resulted in an increasing rate of precipitation of approximately 1.5 mm per year (Figure 9b).

Climatic Trends
Because no weather station data are available in the study area, the most recent version of the University of Delaware dataset of air temperature and precipitation at the grid (43.75 • N, 83.75 • E) that covers the target region was used to reflect the trends of climate between 1900 and 2017 [49]. As shown in Figure 9a, the mean annual temperature (MAT) was the coldest at the beginning of the twentieth century and then fluctuated with positive and negative trends throughout the century. The overall trend of MAT is increasing, and the 1977-2017 period, which overlaps closely with the investigated glacier change period, showed a significant positive trend coefficient (0.03 • C year −1 ) using linear regression (Figure 9a). Regarding the total annual precipitation (TAP), it can be observed that the variability was higher for the first half of the century than the second half (Figure 9b).

Comparison with Glacier Inventory Data
The 2018 inventory of glaciers at the central massif of the Borohoro Mountains was generated using a reproducible, time-efficient semi-automated procedure based on the Sentinel 2A imagery [38]. To compare and cross-check with the most recent, publicly available glacier inventory data, the glacier outlines of the new 2018 inventory were overlaid with outlines of the second Chinese Glacier Inventory (2nd-CGI) [18] (Figure 10a). The 2nd-CGI was released in December 2014 by the Chinese Academy of Sciences and was an updated version of the first CGI completed in 2002 [17]. The primary data sources for the 2nd-CGI glaciers were Landsat TM and the Enhanced Thematic Mapper Plus (ETM+) images from 2006-2010. The Randolph Glacier Inventory version 6.0 (RGIv6.0) and the GLIMS database incorporated the 2nd-CGI data into their world glacier inventories, making it readily and freely accessible to a broader audience. A direct comparison of the 2018 results of this study with the counterparts extracted from the RGIv6.0 demonstrated a larger coverage of 236.1 km 2 in the RGIv6.0, by 8.6% (or 20.3 km 2 ) larger than the 2018 extent. When comparing the RGIv6.0 extent with the mapping from the 2007 Landsat TM image, the total area decreased similarly, with a minor overestimation (~5 km 2 ) in the RGIv6.0 dataset. However, the number of glaciers was found to be significantly different among datasets. The inventory from this study contains 224 glaciers in 2018 and 221 glaciers in 2007, while the RGIv6.0 possesses 176 entities. There may be two reasons for the variation in the number of glaciers. The first is that the methods applied to divide ice polygons may be different, as well as uncertainties in the use of different elevation data in deriving basin boundaries; the second possibility is that some small glaciers are missing from the compilation of the 2nd-CGI, as seen in Figure 10b,c. The disintegration also contributed to the increase of the count of ice bodies in recent years, but is not sufficient to explain a difference of 45 glaciers in the 2007 inventory vs. the RGIv0.6 data since they were sampled around the same time.

Comparison with Glacier Inventory Data
The 2018 inventory of glaciers at the central massif of the Borohoro Mountains was generated using a reproducible, time-efficient semi-automated procedure based on the Sentinel 2A imagery [38]. To compare and cross-check with the most recent, publicly available glacier inventory data, the glacier outlines of the new 2018 inventory were overlaid with outlines of the second Chinese Glacier Inventory (2nd-CGI) [18] (Figure 10a). The 2nd-CGI was released in December 2014 by the Chinese Academy of Sciences and was an updated version of the first CGI completed in 2002 [17]. The primary data sources for the 2nd-CGI glaciers were Landsat TM and the Enhanced Thematic Mapper Plus (ETM+) images from 2006-2010. The Randolph Glacier Inventory version 6.0 (RGIv6.0) and the GLIMS database incorporated the 2nd-CGI data into their world glacier inventories, making it readily and freely accessible to a broader audience. A direct comparison of the 2018 results of this study with the counterparts extracted from the RGIv6.0 demonstrated a larger coverage of 236.1 km 2 in the RGIv6.0, by 8.6% (or 20.3 km 2 ) larger than the 2018 extent. When comparing the RGIv6.0 extent with the mapping from the 2007 Landsat TM image, the total area decreased similarly, with a minor overestimation (~5 km 2 ) in the RGIv6.0 dataset. However, the number of glaciers was found to be significantly different among datasets. The inventory from this study contains 224 glaciers in 2018 and 221 glaciers in 2007, while the RGIv6.0 possesses 176 entities. There may be two reasons for the variation in the number of glaciers. The first is that the methods applied to divide ice polygons may be different, as well as uncertainties in the use of different elevation data in deriving basin boundaries; the second possibility is that some small glaciers are missing from the compilation of the 2nd-CGI, as seen in Figure 10b,c. The disintegration also contributed to the increase of the count of ice bodies in recent years, but is not sufficient to explain a difference of 45 glaciers in the 2007 inventory vs. the RGIv0.6 data since they were sampled around the same time. The advent of Sentinel 2A since 2015 has benefited glacier mapping with its short repeat interval and high-resolution imagery, e.g., [38][39][40]58], however, it has not yet been utilized often in creating glacier data across the Tian Shan. As the high temporal resolution gives a better opportunity of obtaining cloud-free and snow-free images, the finer spatial resolution should allow a more accurate mapping, especially for regions dominated by small glaciers, like the Borohoro Mountains. As illustrated in Figure 10, despite that glacier area loss during ~2007 (RGIv6.0 data acquisition) and 2018 (the current mapping), a better quality of the 2018 outlines shows a more precisely digitized boundary with details along the edge, nunataks, and in shadows. Also, the very small glaciers (e.g., <0.1 km 2 ) that might not be detectable when using Landsat TM/ETM+ can be recognized now, although they may not affect much the total mapped glaciated area.

Comparison with Other Regions across the Tian Shan
The results of the areal and surface elevation change of glaciers in the Borohoro Mountains are consistent with an expected and widely published trend of glacier retreat in the most recent decades across the Tian Shan range. The retreat of the 200 glaciers in 1977 was quantified with a loss of 24.9 ± 3.2% (0.61 ± 0.01% year −1 ) of the total area in 41 years onward to 2018. Although it is difficult to compare the rate with different time frames and different sizes of glaciers, this magnitude of change is among the highest reported in many parts of the Tian Shan (Table 5; Figure 11). Previous studies reported that the strongest annual area shrinkage (0.38-0.76% year −1 ) occurred in the outer ranges of the Tian Shan or in peripheral, low-elevation ranges near the densely populated forelands [8]. For examples, in northern and western Tian Shan, Narama et al. [22] examined four mountain regions and found the most dramatic change occurred in the outer regions, with 19% change during ~1970-2000 (~0.63% year −1 ); Bolch [20] used the Soviet Glacier Inventory and the Landsat scenes for the Zailiyskiy and Kungey Alatau valleys in the northern Tian Shan and found a decrease of >32% between 1955 and 1999 (~0.72% year −1 ). In the inner ranges where glaciers are larger and at higher elevations, studies found them stayed relatively stable with smaller area shrinkage rates (0.15-0.40% year −1 ), in regions like the Terskey Alatau, the Akshiirak and Ala Archa, and the Aksu River Basin [6,8,53,61,62]. The Borohoro is located at the northern periphery of the eastern Tian Shan and is in line with the broader pattern of where glacier area loss is faster. The advent of Sentinel 2A since 2015 has benefited glacier mapping with its short repeat interval and high-resolution imagery, e.g., [38][39][40]58], however, it has not yet been utilized often in creating glacier data across the Tian Shan. As the high temporal resolution gives a better opportunity of obtaining cloud-free and snow-free images, the finer spatial resolution should allow a more accurate mapping, especially for regions dominated by small glaciers, like the Borohoro Mountains. As illustrated in Figure 10, despite that glacier area loss during~2007 (RGIv6.0 data acquisition) and 2018 (the current mapping), a better quality of the 2018 outlines shows a more precisely digitized boundary with details along the edge, nunataks, and in shadows. Also, the very small glaciers (e.g., <0.1 km 2 ) that might not be detectable when using Landsat TM/ETM+ can be recognized now, although they may not affect much the total mapped glaciated area.

Comparison with Other Regions across the Tian Shan
The results of the areal and surface elevation change of glaciers in the Borohoro Mountains are consistent with an expected and widely published trend of glacier retreat in the most recent decades across the Tian Shan range. The retreat of the 200 glaciers in 1977 was quantified with a loss of 24.9 ± 3.2% (0.61 ± 0.01% year −1 ) of the total area in 41 years onward to 2018. Although it is difficult to compare the rate with different time frames and different sizes of glaciers, this magnitude of change is among the highest reported in many parts of the Tian Shan (Table 5; Figure 11). Previous studies reported that the strongest annual area shrinkage (0.38-0.76% year −1 ) occurred in the outer ranges of the Tian Shan or in peripheral, low-elevation ranges near the densely populated forelands [8]. For examples, in northern and western Tian Shan, Narama et al. [22] examined four mountain regions and found the most dramatic change occurred in the outer regions, with 19% change during~1970-2000 (~0.63% year −1 ); Bolch [20] used the Soviet Glacier Inventory and the Landsat scenes for the Zailiyskiy and Kungey Alatau valleys in the northern Tian Shan and found a decrease of >32% between 1955 and 1999 (~0.72% year −1 ). In the inner ranges where glaciers are larger and at higher elevations, studies found them stayed relatively stable with smaller area shrinkage rates (0.15-0.40% year −1 ), in regions like the Terskey Alatau, the Akshiirak and Ala Archa, and the Aksu River Basin [6,8,53,61,62]. The Borohoro is located at the northern periphery of the eastern Tian Shan and is in line with the broader pattern of where glacier area loss is faster.
With four time slices selected in this study, it revealed that the 1994-2007 period had the highest rate of change (0.90 ± 0.07% year −1 ), followed by 0.62 ± 0.06% year −1 during 2007-2018 and 0.51 ± 0.03% year −1 during 1977-1994. Meltwater from glaciers in the Borohoro Mountains flows into tributaries of the Jinghe River system within the Ebinur Lake Basin in Xinjiang, China. Previous studies estimated a shrinkage rate of glacier area of 0.37% year −1 in the Ebinur Lake Basin [32] and of 0.38% year −1 in the Jinghe river basin during 1964-2004 [37] (Table 5). The higher rate of glacier retreat in Borohoro during 1977-2018 could be a result of accelerated ablation trend after the 1970s [8].
Comparing the first and second CGIs, Che et al. [63] reported a change rate of 0.72% year −1 for the sub-region that contains the Borohoro Mountains, which is similar to my result, but the quality issues (e.g., distortion, location shift, and missing parts) in the first CGI might have sacrificed the accuracy of the estimation to some extent.
Water 2020, 12, x FOR PEER REVIEW 17 of 23 Figure 11. Spatial pattern of glacier recession rates of the Tian Shan in the referenced literature (in brackets) in correspondence with the data presented in Table 5.
With four time slices selected in this study, it revealed that the 1994-2007 period had the highest rate of change (0.90 ± 0.07% year −1 ), followed by 0.62 ± 0.06% year −1 during 2007-2018 and 0.51 ± 0.03% year −1 during 1977-1994. Meltwater from glaciers in the Borohoro Mountains flows into tributaries of the Jinghe River system within the Ebinur Lake Basin in Xinjiang, China. Previous studies estimated a shrinkage rate of glacier area of 0.37% year −1 in the Ebinur Lake Basin [32] and of 0.38% year −1 in the Jinghe river basin during 1964-2004 [37] (Table 5). The higher rate of glacier retreat in Borohoro during 1977-2018 could be a result of accelerated ablation trend after the 1970s [8]. Comparing the first and second CGIs, Che et al. [63] reported a change rate of 0.72% year −1 for the sub-region that contains the Borohoro Mountains, which is similar to my result, but the quality issues (e.g., distortion, location shift, and missing parts) in the first CGI might have sacrificed the accuracy of the estimation to some extent.

Location
Investigated Period

Number/Total Area (km 2 ) of Glaciers
Change Rate (% year −1 ) Figure 11. Spatial pattern of glacier recession rates of the Tian Shan in the referenced literature (in brackets) in correspondence with the data presented in Table 5.

Influence of Topo-Climatic Factors
Within the study area, the glacier area loss in the past 41 years ranges from a minimum of 4.9% to complete disappearances (100%). The vast variability in the degree of individual glacier recession was suggested to be largely dependent on the topographic and climatic conditions, e.g., [57,59,65,66]. According to the results of this study, the average loss of area in three small size classes (<0.5 km 2 , 0.5-1.0 km 2 , and 1.0-5.0 km 2 ) is drastically lower than that of the other two large size classes (5.0-10.0 km 2 and >10.0 km 2 ), but the relative rates of recession are much higher for small glaciers ( Figure 5). The four largest glaciers located in the central massif lost 12.0% on average from 1977 to 2018, which is converted to 0.29% per year. Such relatively stable recession can be attributed to the type and topographic characteristics of these glaciers. Large compound valley glaciers that originate from the extensive accumulation area, often joined by multiple tributaries, occupy the highest elevations in the range and flow downstream to a narrow tongue. Where the topography of their accumulation area is gentle and frequently replenished by snowfall, the surface extent might not shrink significantly, but the thickness could be thinning. Some studies have shown that sublimation is observed as a major mass loss for some glaciers in the semiarid and arid environments [67,68]. Reduction in the area is mostly found in the frontal areas, even though the occasional debris cover might slow the ablation [26,69]. Small glaciers (usually <1 km 2 ) resulted in a rapid recession, and half of them had lost more than 48.8% of its area from 1977 to 2018. The markedly retreated glaciers are generally located along the periphery of the Borohoro massif, and more are found on the southern slope than the northern slope. Existing publications have claimed that the mid-latitude small glaciers are subject to strong retreat and are prone to give a sensitive response to the warming climate, and in fact, many have already vanished in the last century, e.g., [13,70,71].
The morphology of a glacier has a significant influence on its responsiveness to the changing climate. The area and elevation-related factors (minimum, maximum, median, and range) of a glacier are interconnected because they both tie closely to the intrinsic type of the glacier: large compound glaciers are more likely to occupy a larger range of elevations and have a high peak and a low front; while small hanging glaciers range small altitudinal band but may exist as a survivor at various altitudinal locations, depending on the surrounding topography. The aspect of individual glaciers did not correlate well with their glacier area change, but the overall spatial pattern as shown in Figure 3 depicted a higher level of loss on the southern side of the major divide than on the northern side. The hypsometric integral was found to be a significant negative correlation (r = −0.43, p < 0.01) with the area change of glaciers, meaning glaciers that possess a convex profile (more areas distributed at higher elevations) in the cumulative elevation-area plot resulted in less retreat, while glaciers with a concave profile (more areas distributed at lower elevations) experienced a greater retreat. The hypsometry was suggested to be used for making a reasonable first-order estimation of the mass-balance sensitivity of glaciers without in situ measurements [72]. Overall, glacier topomorphology, such as area, elevation, and hypsometry, might be important factors to the retreat rates of glaciers in the Borohoro, but to determine a more accurate assessment on the pattern of glacier loss, field measurements or adequate models would be necessary.
Reductions in glacier area and lowering of the glacier surface elevation are the outcome of a negative glacier mass balance, which is ultimately driven by the change in the synoptic climatic pattern and local meteorological conditions. The regional climatic trend in Central Asia was confirmed with rising temperatures since the 1970s by a number of studies, e.g., [8,11,73,74]. The IPCC Fourth Assessment Report (AR4) reported decadal trend coefficients between +0.1 and +0.2 • C and a pronounced warming in winter season based on nearly all meteorological stations for Central Asia (35-50 • N, 40-75 • E) [75], but the inconsistency and lack of stations at high elevations may cause some difficulties in the estimations. Haag et al. [11] examined the climate of Central Asia (35)(36)(37)(38)(39)(40)(41)(42)(43)(44)(45) • N, 65-80 • E) using the Climatic Research Unit Timeseries (CRU TS) v. 4.01 dataset between 1950-2016 and found that the mean annual temperatures increased at a rate of 0.28 • C per decade, with the maximum rate of 0.32 • C per decade in winter. Five meteorological stations most adjacent to the study area revealed a rate of temperature increase of 0.31 • C/decade over 1959-2005, and the rate accelerated to 0.54 • C/decade during the last 15 years [32]. For the studied Borohoro region, an increase of 0.3 • C per decade was observed in the mean annual temperature from 1977 to 2017, very similar to the previously reported values. This uniform, continuous temperature increase in the last few decades can be considered as the major driver of the ongoing glacier recession. It has been suggested that a weakening Siberian anticyclone system might have caused such a change in temperature [76]. Precipitation often shows marked spatial heterogeneity and thus the trends in precipitation are quite different across ranges of the Tian Shan. The total annual precipitation showed an insignificant or decreased trend in the central and western Tian Shan, e.g., [8,53,73], whereas an increasing trend is commonly found in the eastern and northern ranges, e.g., [11,64,77]. The precipitation in Borohoro experienced a steady upward trend in 1977-2017 with a rate of about 150 mm/decade. A weakening and displaced Siberian anticyclone can strengthen the advection of the mid-latitude westerlies into the peripheral Tian Shan [73], thus causing increased precipitation. In the Jinghe River Basin, orographic thunderstorms are dominant in summer while winter is cold and dry [37]. A major proportion of rain versus snowfall precipitation, together with the rising temperature, leads to unfavorable conditions for the glacier mass budget and accelerates the ablation processes. In addition, it is also important to note the existence of the response time of glaciers to climate signals because, counting such delays, the level of glacier's retreat is likely to be retained in the near future even if the warming trend was stopped or decelerated.

Conclusions
In this work, glacier changes during the past 41 years in the Borohoro Mountains, Tian Shan, were examined with declassified Hexagon image, Landsat 5 TM, Sentinel 2A, SRTM DEM, and GLAS/ICESat remote-sensing data. A total of 200 glaciers in the study area which covered 287.5 ± 8.2 km 2 in 1977 shrank by 71.6 ± 9.2 km 2 (0.61 ± 0.01% year −1 ) by 2018, with an increase to a total of 224 glaciers due to fragmentation. Among three sub-periods of 1977-2018, the rate of glacier area change was found the fastest (0.90 ± 0.07% year −1 ) in 1994-2007. The rates are similar to the published results for the peripheral and eastern ranges of the Tian Shan, where glaciers experienced a more rapid recession than those in the inner ranges. Small glaciers (<1 km 2 ) take up 80% of the 2018 glacier inventory of the region. The inverse relationship between the glacier size and shrinkage rate suggests that small glaciers are more sensitive and vulnerable to a warming climate, in agreement with findings of many studies worldwide. An individual glacier's median, maximum, and range of elevation, and its hypsometry also showed significant correlations with the rate of retreat, but it could be difficult to determine the degree of their influence due to the collinearity among the factors. However, apparently the aspect of the cirques did not influence the rate of glacier shrinking in the studied period of time. The ICESat dataset that crossed two of the largest glaciers in the Borohoro helped to reveal the surface elevation change from 2003 to 2009, which showed thinning at both accumulation and ablation zones. The climatic conditions in the last 41 years are remarkably different from the early 20th century. Both increases in mean annual temperature and total annual precipitation were clearly shown in the gridded climate dataset for this area as well as in other research of nearby regions. However, most of the precipitation arrives in summer, which is unfavorable for glacier growth. While the mesoor synoptic scale of the weather system, i.e., the interplay of the Siberian High and the westerlies, may have caused the trends in the climatic variables, it is reasonable to interpret that the warming temperatures are the major driver of glacier recession and such recession is very likely to continue in forthcoming decades.