Topography-Related Glacier Area Changes in Central Tianshan from 1989 to 2015 Derived from Landsat Images and ASTER GDEM Data

Studies have investigated the glacier projected area (2D Area) on a horizontal plane, which is much smaller than the glacier topographic surface extent (3D Area) in steep terrains. This study maps the glacier outline in Central Tianshan using Landsat images from four periods, i.e., 1989, 2002, 2007 and 2015, by an object-based classification approach, and analyzes the glacier 2D and 3D area changes related to topographic factors based on Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global DEM data. This approach shows an accuracy of 90.8% for clean ice mapping. The derived clean ice outlines are in good agreement with the 2nd Chinese Glacier Inventory (CGI2) and the Global Land Ice Measurements from Space (GLIMS). The fields with a northern aspect receive the least surface solar radiation, leading to dominant existing glaciers. Glaciers are near evenly distributed in slope zones of 0◦ to 50◦ and have a mean slope angle of 28.8◦, resulting in a 30.3% larger 3D area than the 2D area in 2015 in Central Tianshan. The glacier 2D area decreased by 404 km2 (−8.1%) between 1989 and 2015, while the 3D area declined by 516 km2 (−7.9%). The glacier 2D area showed a reduction of −1.8% between 1989 and 2002, −3.8% between 2002 and 2007, and −2.7% between 2007 and 2015, and these retreating rates closely responded to the variations of regional mean air temperature and precipitation. Topographically, most reductions occurred in elevation bands of 3000–4000 m and in slope zones of 10–20◦ and 40–50◦, and in the eastern aspect fields. The northern Tekes River catchment had the largest shrinking rate of −17.0% (2D area), followed by the southern Karasu River (−14.2%) and Muzart River (−7.7%) catchments. In contrast, glaciers in the Kumerik/Aksu and Tailan River catchments in the Tuomuer region showed little change (−2%).


Introduction
Glaciers and snow in Central Tianshan are crucial water resources in this arid and semiarid region [1,2].Almost all oasis agriculture and residents in the surrounding areas are set up along river valleys and are dependent on water resources from the upper Tianshan Mountains, where glacier melt water contributes a significant share of stream inputs and is a buffer against severe drought [3].The contributions of glacier melt to runoff vary from 3.5% to 67.5% (mean 24.0%) among the twenty-four catchments in Tianshan Ranges, and the annual runoff stability increases with the glacier area ratio following a power function curve [4].A small change in the glacier area ratio may cause dramatic changes in hydrologic regimes.Glaciers are sensitive to the rising air temperature and are among the best natural indicators of climatic change [5,6].Accurate estimation of glacier variation is significant for water resource prediction and planning [7,8].Repetitive mapping of glacier areas in different time periods could better assess glacier-climate interactions and predict the future runoff variations in glacierized catchments [9,10].
The full glacier inventories in the Tianshan Mountains has been produced primarily by the First and Second Chinese Glacier Inventory (CGI1 and CGI2) data, and World/Randolph Glacier Inventory (WGI/RGI), Global Land Ice Measurements from Space (GLIMS), and Glacier Area Mapping for Discharge from the Asian Mountains (GAMDAM, or GGI) by manual mapping from Landsat images [11,12].CGI1 and CGI2 are the original data sources for WGI and GLIMS within China, GGI is the source data of RGI in High Mountains Asia, and RGI is a parallel inventory to GLIMS [11,[13][14][15].
CGI1 was compiled mainly from different aerial photographs acquired during the 1950s-1980s, and CGI2 was derived from Landsat images acquired mainly from 2006-2010 [16].Wang et al. (2011a) analyzed 3000 glaciers (total 9035) in the Chinese Tianshan Mountains based on CGI1 and CGI2 and found the glacier areas declined by −11.5% from the 1960s to 2009 [17].In contrast, by using the same CGI1 and CGI2 data, Xu et al. (2015) found that the glacier areas (1511 km 2 in 1970s) decreased by −22.0% in the Tekes River catchment of Central Tianshan [18].There are large uncertainties in evaluating glacier changes by direct comparison of CGI1 and CGI2 and other glacier inventories due to their different data sources and glacier classification methods applied.
Glacier change detection can be best estimated by comparison of glacier areas derived from a same classification approach and similar image sources.Numerous studies have identified different change rates of individual glaciers at a catchment scale in the Tianshan Mountains mainly by using airborne photos and satellite images.Seasonal snow cover and the hanging cloud on high mountains are the two dominant factors that constrain optical sensors to map glacier change in catchment studies and further limit knowledge of the glacier status at a regional level [5].Nonetheless, studies focusing on a catchment could identify detailed variations with high confidence.For example, in Central Tianshan, the glacier areas showed little changes during 1999-2007 in the Inylchek Glacier, Kyrgyzstan [19], while the glacier areas decreased by −12.6% (−0.29% a −1 ) between 1965 and 2003 in the eastern Terskey-Alatoo Range, Kyrgyzstan [5], and by −7.8% (−1.30% a −1 ) in the Muzart Glacier catchment, China during 2007-2013 [20].In Northern Tianshan, the glacier areas decreased from 91.3 km 2 in 1964 to 77.4 km 2 in 2004 (−0.38% a −1 ) in the Jinghe River Basin, China [21], while they dramatically decreased from 142.8 km 2 to 109.0 km 2 (−1.02% a −1 ) from 1989 to 2012 in the Karatal River Basin, Kazakhstan [22].The largest shrinkage rate was found to be −1.67% a −1 between 2007 and 2013 in the Mountain Karlik of Eastern Tianshan, China [23].In contrast, in Central Tianshan of China territory, Landsat images identified that Qinbingtan Glacier No. 74 advanced ~90 m and the glacier area increased by 6.4% between 2002 and 2010, and the West Qongterang Glacier had a significant volume gain in the central part of the glacier tongue although it lost mass in the upper part, a typical characteristic of glacier surge [24].Those studies demonstrate large spatial heterogeneity of glacier variations in the Tianshan Mountains, while there is a limited overall budget of the glacier areas at a regional scale based on the same approach in the Landsat images era, e.g., in Central Tianshan.
Moreover, glacier inventories and most previous investigations are based on the planar area in two horizontal dimensions (2D area) projected onto the ellipsoid Earth surface, which is defined within the glacier ice mass balance community [25], rather than the real topographic surface extent (3D area) projected onto the slope normal, which is practiced within the snow cover community [26].Accordingly, their thickness is defined as the vertical length measured parallel to the horizontal normal (glacier thickness) and slope normal (snow depth), separately [20].Geometrically, the 3D area is 5% larger than the 2D area as the slope angle is larger than 18 • , 15% larger for 30 • , and 41% larger for 45 • .For instance, the mean slope angle is 31.5 • in the Muzart Glacier Catchment, resulting in a 38.1% larger glacier 3D area than the 2D area [20].The mean glacier surface slope angle of CGI2 in China is 19.9 • , while over one third of the surface slope angles are greater than 30 • in steep terrains, such as the Central Tianshan, Pamir Plateau, Qilian Mountains and Altai Mountains [16].It is worth investigating the distribution of the differences between the glacier 3D area and 2D area and how they are related to glaciers retreating, especially in the high Asian mountain glaciers with steep slopes.
The primary objectives of this study are to (1) map and quantify the glacier 2D and 3D area distribution and difference in Central Tianshan by using topographic data (ASTER GDEM V2) and an object-based classification approach from Landsat images in four time periods : 1989, 2002, 2007 and 2015; and (2) investigate how glacier area changes are related to topography such as surface slope zones, elevation bands, aspect fields (orientation), and natural glacier catchments.Thus, a whole budget of the glacier area variations can be best assessed at a regional wide and in the five major glacier (river) catchments over the entire Central Tianshan, the water tower of Central Asia.

Study Area
Central Tianshan (also called Tien Shan in some literature) mainly lies in Northwest China (in Xinjiang Uygur Autonomous Region), Kazakhstan and Kyrgyzstan, extending from 41.5 • N to 43.0 • N and 79.0 • E to 82.5 • E (Figure 1).The mountain ranges run mainly from northeast to southwest, with the highest peaks of 7010 m (Tuomuer) and 7439 m (Jengish Chokusu).The study areas are constrained above the 2500-m elevation contour, where most perennial snow and ice occur.The total study areas are 23,000 km 2 , and about 20% of the study areas are covered by glaciers or perennial snow.According to the World Glacier Inventory (WGI) in 2012, there are a total of 2600 glaciers and over 6000 km 2 in this region, i.e., about 40% of the total glacier areas in the entire Tianshan ranges [5].Most of these glaciers are valley glaciers (~80%) or hanging glaciers within rugged and steep terrains.Consequently, the glacier 3D surface area is 38.1% larger than their 2D planar area in the Muzart Glacier Catchment [20].There are the two largest debris-covered glaciers in China around the Tuomuer Peak.The Tuomuer Glaciers had a debris-covered area of 63 km 2 (17% of its total area), and the Tugbelqi glaciers had a debris-covered area of 39 km 2 (13% of its total area) [16].
Central Tianshan provides critical water resources for the surrounding streams, oasis wetland, farmland, industry and urban areas via summer rainfall runoff, snow and glacier melting.The annual precipitation in the high mountains (400-600 mm) is two to three times as much as it is in the lowland valley [5].Over 70% of the streamflow in the Tarim River (south) comes from Central Tianshan via its main tributaries of Aksu River, Tailan River and Muzart River [27].In the Ili River Basin (north), the Tekes River collects surface runoff from the northern sides of Central Tianshan, contributes the primary streamflow for the downstream Ili River, and finally flows into the Lake Balkhash in Kazakhstan.The west-facing and V-shaped Ili River Valley is surrounded by the Tianshan Mountains on three sides, which form a humid and temperate climate in spite of being distant from oceans.This valley is the precipitation center of the Tianshan Mountains [18].Glacier melt is a good buffer against severe droughts, and glacier changes in the high mountains considerably affect the regional water supply in the surrounding lowland areas [3,18].

Data Sets
Multiple Landsat images are acquired to map the glacier area changes from 1989 to 2015, including Landsat 5 thematic mapper (TM), Landsat 7 enhanced thematic mapper plus (ETM+), and Landsat 8 operational land imager (OLI) (Table 1).All Landsat images are downloaded from the USGS Global Visualization Viewer (http://glovis.usgs.gov).They are orthorectified images and re-projected into 1984 World Geodetic System (WGS) Universal Transverse Mercator (UTM) zone N44 projection.Four tiles of Landsat images are needed to cover the entire study areas (Figure 1), and most of glaciers distribute within the image of 147/031 (path/row, Table 1).These images in each time period span two to three years in order to select the optimal ones with minimum cloud and seasonal snow cover.Both cloud and seasonal snow cover are the most important factors affecting the accuracy of glacier classifications [28].Meanwhile, a suitable high-resolution GeoEye01 (nominal 4-m spatial resolution)

Data Sets
Multiple Landsat images are acquired to map the glacier area changes from 1989 to 2015, including Landsat 5 thematic mapper (TM), Landsat 7 enhanced thematic mapper plus (ETM+), and Landsat 8 operational land imager (OLI) (Table 1).All Landsat images are downloaded from the USGS Global Visualization Viewer (http://glovis.usgs.gov).They are orthorectified images and re-projected into 1984 World Geodetic System (WGS) Universal Transverse Mercator (UTM) zone N44 projection.Four tiles of Landsat images are needed to cover the entire study areas (Figure 1), and most of glaciers distribute within the image of 147/031 (path/row, Table 1).These images in each time period span two to three years in order to select the optimal ones with minimum cloud and seasonal snow cover.Both cloud and seasonal snow cover are the most important factors affecting the accuracy of glacier classifications [28].Meanwhile, a suitable high-resolution GeoEye01 (nominal 4-m spatial resolution) acquired on 29 July 2013 is downloaded to validate the glacier outline derived from Landsat images acquired on 31 July 2013.The derived glacier areas in this study are also compared to those from WGI (CGI1) and GLIMS (CGI2).CGI1 and CGI2 are downloaded from the West Data Center for Glaciology and Geocryology(WDCGG), Lanzhou, China at http://wdcdgg.westgis.ac.cn/ [29].WGI and GLIMS) are downloaded from the US National Snow and Ice Data Center (NSIDC, http://nsidc.org).Table 2 summarizes the features of those products [13][14][15].CGI1 was compiled mainly from different aerial photographs acquired during the 1950s-1980s (simplified as 1970s hereafter), which spans a long period and cannot easily represent the contemporary glacier status in China.CGI2 was derived from 218 Landsat TM/ETM+ scenes (30 m of spatial resolution) acquired mainly from 2006 to 2010 by using the band ratio segmentation method [16].Those images of 2006-2007 used in this study are the same as those used in CGI2.Glacier positioning errors were about 10 m for clean-ice outlines and 30 m for debris-covered outlines, and area errors were 3.2% [16].Both CGI1 and CGI2 are the original data sources of WGI and GLIMS in China.
WGI contains over 130,000 glaciers and is based primarily on aerial photographs and maps like CGI1.Inventory parameters normally include geographic location, 2D area, length, orientation, elevation, and classification, whereas most original glaciers have one data entry and are viewed as a snapshot of the glacier distribution in the second half of the 20th century [15].GLIMS is an international joint effort with over 60 institutes world-wide to repeatedly survey of the world's 200,000 glaciers from space.The primary data sources used in GLIMS are collected by the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) instrument aboard the Terra satellite and the Landsat series of satellites, along with historical information derived from maps and aerial photographs (Table 2).The goal of GLIMS is to create a globally comprehensive inventory of land ice including measurements of glacier area, geometry, surface velocity, and snow line elevation [14].
The ASTER Global Digital Elevation Model (GDEM) Version 2 (V2) data were downloaded from Japan Space Systems (http://gdem.ersdac.jspacesystems.or.jp/), and had a 30-m spatial resolution, 17-m vertical accuracies and 71-m horizontal accuracies [30] (Frey and Paul, 2012).It is used to establish the Triangulated Irregular Network (TIN) 3D surface model, delineate the glacier catchment, compute slopes and aspects, and classify elevation bands in this study.
The daily precipitation and daily mean air temperature data were downloaded from the Global Historical Climatology Network (GHCN), http://www.ncdc.noaa.gov/oa/climate/.These GHCN stations are distributed at lower elevation ranges below the study area in Xinjiang Uygur Autonomous Region of China, and the two nearby stations are located in Paicheng County of the lower Muzart River Basin and in Aksu County of the lower Aksu River Basin (Figure 1).The time spans from 1961 to 2015.

Object-Based Image Classification
The actual terminus of supraglacial debris cover is difficult to delineate by general image classification because of their similar spectral signature with adjacent periglacier debris [9].Although there are large areas of debris-covered glaciers (729 km 2 in 2007) in the Central Tianshan (Southern Inylchek, Tomur, and Koxkar glacier) and their surface thinning is considerable, their areas were relatively stable and only shrank slightly by 1.7 ± 4.8% from 1975 to 2008 [27,31].Therefore, this study only focuses on clean ice mapping and change analysis from Landsat images using the object-based classification approach, which makes full use of the image reflectance, texture and topographic information and is efficient in manual correction after the initial classification [20].It is quite advanced to map the clean ice glaciers based on optical satellite images, such as the normalized difference of snow cover index (NDSI) used in the MODIS snow cover product [32], band ratios of b3/b5 used in CGI2 [16], and object-based classification [33][34][35][36].The object-based classification approach is applied in this study to delineate the glacier area in software eCognition 9.0 (Trimble Inc., Sunnyvale, CA, USA).The specific procedures are illustrated in Figure 2.
Although there are large areas of debris-covered glaciers (729 km 2 in 2007) in the Central Tianshan (Southern Inylchek, Tomur, and Koxkar glacier) and their surface thinning is considerable, their areas were relatively stable and only shrank slightly by 1.7 ± 4.8% from 1975 to 2008 [27,31].Therefore, this study only focuses on clean ice mapping and change analysis from Landsat images using the object-based classification approach, which makes full use of the image reflectance, texture and topographic information and is efficient in manual correction after the initial classification [20].
It is quite advanced to map the clean ice glaciers based on optical satellite images, such as the normalized difference of snow cover index (NDSI) used in the MODIS snow cover product [32], band ratios of b3/b5 used in CGI2 [16], and object-based classification [33][34][35][36].The object-based classification approach is applied in this study to delineate the glacier area in software eCognition 9.0 (Trimble Inc., Sunnyvale, CA, USA).The specific procedures are illustrated in Figure 2.  [20]).
The object-based image classification method mainly includes three steps [37].The first step is image segmentation with a multi-resolution segmentation method.It minimizes the mean heterogeneity of image objects for a given resolution and maximizes their respective homogeneity by a bottom-up segmentation algorithm.The original Landsat (GeoEye01) images are first segmented into different pixel patches (objects) using a scale parameter of 50, shape factor of 0.9 and compactness factor of 0.5.The scale parameter determines the size of image objects and is obtained by using try-and-error analysis.The value of the Shape factor determines the weight of the Shape and Color criteria, here Shape = 0.9 (allowed maximum value), and Color = 1-Shape = 0.1, which indicates that objects are optimized for spatial homogeneity or shape, instead of color.The The object-based image classification method mainly includes three steps [37].The first step is image segmentation with a multi-resolution segmentation method.It minimizes the mean heterogeneity of image objects for a given resolution and maximizes their respective homogeneity by a bottom-up segmentation algorithm.The original Landsat (GeoEye01) images are first segmented into different pixel patches (objects) using a scale parameter of 50, shape factor of 0.9 and compactness factor of 0.5.The scale parameter determines the size of image objects and is obtained by using try-and-error analysis.The value of the Shape factor determines the weight of the Shape and Color criteria, here Shape = 0.9 (allowed maximum value), and Color = 1-Shape = 0.1, which indicates that objects are optimized for spatial homogeneity or shape, instead of color.The value of the Compactness factor = 0.5, which indicates a neutral option between compact and non-compact objects.
The second step is to build the class hierarchy, which determines the classified classes' name and variable thresholds in the classification.All pixels within an object are classified into 6 classes.Background: value < 0, non classification area; Glaciers: NDSI > 0.8 and B3/B5 > 15, the criteria for class Glaciers; Vegetation: NDVI > 0.09, the criteria for class Vegetation; Other classes such as Bare rock, Cloud and River do not set classification criteria and will be determined using known features during the Classifier-training.
The third step is to build and apply the classifier.The first step is to train a classifier using the classified objects in the class hierarchy domain.Then, build the configuration: a, which is a self-built classification configuration using the class hierarchy and various features.The classifier-Train is built using the Support Vector Machine (SVM) with a linear kernel [38] based on features such as Brightness, Reflectance from all bands, Length/Width, NDSI, B3/B5, and NDVI of the classified objects.The built configuration is then applied to classify all pixels within each segment into glaciers and non-glaciers such as vegetation, bare rock, river and cloud using SVM and the linear kernel.
The classified images are manually corrected by visual comparison with images acquired in different years.Finally, the classified objects are merged and exported to vector polygons for further visually checking and manual edition in ArcMap, eliminating misclassified pro-glacial water and shadow areas by overlaying with DEM data and Google Earth images [20].The unidentified and cloud-blocked areas (~3% of total glacier areas) in any time period are excluded for the area change analysis for the four periods.

Glacier Combination from Different Tiles
The time of summer minimum snow cover in these areas was in August and early September based on our analysis using the improved daily cloud-free MODIS snow cover product [39].In order to obtain the optimal images for minimum cloud and snow cover, the selected four Landsat images (tiles) covering the study areas were acquired on different dates and even in different years (Table 1).The four images have about 10% overlap among their boundaries, and image 147/031 and image 146/031 have most glaciers (Figure 1).Therefore, each Landsat image is first subjected to glacier classification separately, and then the four classification images are combined into one image to cover the study area based on a priority scheme.Images of 147/031 (Path/Row), 146/031, 146/030 and 147/030 have the first, second, third and fourth priority in combination, respectively.In other words, the final glacier map is primarily (85%) based on the image 147/031, and the gap areas are complemented by the image 146/031, 146/030 and 147/030 subsequently.

Assessment of Glacier Outlines
The object-based classification approach is first validated by comparing glacier outlines derived from a high resolution GeoEye01 image acquired on 29 July 2013 and those derived from Landsat images on 31 July 2013 in the upper Koxqar Baqi Glacier Catchment (Figure 1), since there was no high spatial resolution image available on the date of the Landsat images, i.e., on 24 August 2007 or 15 August 2015 for the image 147/031.The selected validation images were constrained by the limited cloud-free high spatial resolution images and the Landsat images on the same summer date in the study area.The GeoEye01 image was also processed through the object-based image classification method.The derived glacier outlines are illustrated for detailed comparisons with CGI2 and GLIMS in the Tuomuer and Mushketoy Glacier Catchments.Glaciers are also compared with those of CGI2/GLIMS on the entire Central Tianshan.This study only maps clean ice, while both GLIMS and CGI2 include debris-covered glaciers.For better comparison, the debris-covered glaciers are cut off from GLIMS and CGI2 by manually delineating the lower edge of clean ice over Landsat images (bands combination of 5, 4, 3) [9].All glacier comparisons are carried out by spatial intersection and are summarzied in an error matrix table.
The uncertainty of glacier position caused by the comparison among different image cources is limited [40].The spatial resolutions are 28.5 m for Landsat images and 4 m for GeoEye01 images, and the mean error of images geolocation is 19.5 m.According to the uncertainty computation equation used in Ye et al. [41], the uncertainty for glacier comparison between Landsat images and GeoEye01 images is computed as 0.012 km 2 , and 0.016 km 2 for those between Landsat images in different time.These uncertainties due to different image sources are relatively small compared to the large glacier areas of about 4000 km 2 [41,42].Glacier areas are rounded to 0.1 km 2 for validation and 1 km 2 for change analysis on Central Tianshan in different periods.

Computation of 3D Surface Areas
The solid 3D earth surface is usually represented by a raster model like Digital Elevation Model (DEM) or a Triangulated Irregular Network (TIN) model.DEM is a continuous raster surface interpolated from discrete grid elevation by bilinear interpolation [43].The 3D surface area in each grid can be estimated from its 2D area and slope (called raster areas), or directly calculated from 8 3-dimensional triangles that are generated from the elevation of a 3 × 3 grid window by connecting each grid center point with the center points of the 8 adjacent grids [44].TINs are continuous vector surfaces, and the 3D surface area can be precisely measured from the facets of TINs (called polygon areas) and clipped to a certain polygon.The 3D surface areas calculated from a 3 × 3 grid window of DEM are consistently close (error < 0.1%) to the TIN-based area value at grid counts above 250 [44].The raster surface area is lower than that computed from a 3 × 3 grid window and the TIN-based polygon areas by −1.38% for a 2D area of 162 km 2 and a slope of 28.5 • , and the relative difference is even larger for smaller areas [45].
Although the raster surface areas are much faster in computing and are better for neighborhood analysis, they are not accurate in clipping to a specific polygon, especially for small polygons.Therefore, glacier 3D surface areas in this study are computed by the TIN surface for the entire study area, the 500-m elevation bands from 2500 to 6000 m, the 10-degree slope zones from 0-90 degree, the eight aspect fields, and the natural glacier catchments (Figure 1), i.e., Tekes (North), Kumeric/Aksu (West), Tailan, Muzart and Karasu Rivers (South), respectively.The TIN surface is established using the ASTER GDEM V2 data based on a Delaunay triangulation with an elevation tolerance of 5 m, ~1/1000 of the elevation range from 2500 to 7400 m [43,44].The glacier outlines are first clipped into different elevation bands, slope zones, aspect fields and river catchments.The elevation bands are directly grouped from GDEM data and are converted into polygons.Slope zones and aspect fields are reclassified from the slope raster and aspect raster computed from GDEM data and are converted into polygons.Then, the 3D polygon areas are computed by overlapping the clipped glacier polygons onto the TIN surface via the Polygon Volume tool in software ArcMap 10.3.
The elevation tolerance of 5 m in building TIN surface from GDEM data is relatively large compared to values used in the literature [43,44].It is constrained by the allowed maximum number of computing points (~20 million) in ArcMap 10.3 because of the small grid size (30 m) and large study area (23,000 km 2 ).Nevertheless, taking the glacier polygon in 2015 as an example, the glacier 3D areas computed using the TIN surface built by tolerances of 5, 10, 50 and 100 m are 5982, 5973, 5848 and 5774 km 2 , respectively.The 3D areas from the latter three tolerances are only −0.2%, −2.2% and −3.5% smaller than those by a tolerance of 5 m.This indicates that the TIN surface built by the tolerance of 5 m (and up to 10 m) can obtain a a relatively stable surface extent value, which is likely attributed to the large elevation range (4900 m) and a huge number of computing points/grids (>20 million).

Glacier 2D Area Validation and Comparisons
The object-based classification approach utilized in this study is first validated by comparing the glacier areas derived by this approach from a high spatial resolution GeoEye01 image (29 July 2013) and a Landsat image (31 July 2013) in the upper Koxqar Baqi glacier catchment (Figure 1).Their neighboring dates can minimize the impact of seasonal snow cover on glacier identification.Near half of the selected GeoEye01 image area is covered by perennial snow and glacier, and the glacier areas have an agreement of 90.8% by spatial intersection (Table 3).The commission error of the object-based classification approach is 16.7%, which is likely related to the transition pixels and the coarser resolution of Landsat images.The glacier outlines derived from the GeoEye01 image have a similar spatial pattern with those from the Landsat image, and can better capture the smaller glaciers and especially the glaciers' lower edges/terminus due to their fine spatial resolution (Figure 3), which is also reflected by the omission error of 9.2% from Landsat images.
Compared to CGI2, the glaciers classified in this study only include the clean ice/perennial snow, and a large part of the debris-covered glaciers (18.9 km 2 ) are not included in this study's glacier classification (Figure 3).Excluding the debris-covered glaciers in CGI2 and by spatial intersection, the 2D areas of clean ice in CGI2 have an agreement of 84.0% with those derived from Landsat images in this study (Table 3).CGI2 has a larger glacier area than that derived in this study, and the omission error is 16.0%.Part of the area difference may be attributed to the glacier melting during the 6 years between 2007 and 2013, while other factors such as classification methods, the existence of the seasonal snow cover on different dates and the identified debris-covered glaciers also contribute Water 2018, 10, 555 9 of 23 to the discrepancy.The commission error is 5.1%, which is much less than that (16.7%) against the GeoEye01 image.) and a Landsat image (31 July 2013) in the upper Koxqar Baqi glacier catchment (Figure 1).Their neighboring dates can minimize the impact of seasonal snow cover on glacier identification.Near half of the selected GeoEye01 image area is covered by perennial snow and glacier, and the glacier areas have an agreement of 90.8% by spatial intersection (Table 3).The commission error of the object-based classification approach is 16.7%, which is likely related to the transition pixels and the coarser resolution of Landsat images.The glacier outlines derived from the GeoEye01 image have a similar spatial pattern with those from the Landsat image, and can better capture the smaller glaciers and especially the glaciers' lower edges/terminus due to their fine spatial resolution (Figure 3), which is also reflected by the omission error of 9.2% from Landsat images.Compared to CGI2, the glaciers classified in this study only include the clean ice/perennial snow, and a large part of the debris-covered glaciers (18.9 km 2 ) are not included in this study's glacier classification (Figure 3).Excluding the debris-covered glaciers in CGI2 and by spatial intersection, the 2D areas of clean ice in CGI2 have an agreement of 84.0% with those derived from Landsat images in this study (Table 3).CGI2 has a larger glacier area than that derived in this study, and the omission error is 16.0%.Part of the area difference may be attributed to the glacier melting during the 6 years between 2007 and 2013, while other factors such as classification methods, the existence of the seasonal snow cover on different dates and the identified debris-covered glaciers also contribute to the discrepancy.The commission error is 5.1%, which is much less than that (16.7%) against the GeoEye01 image.The glacier areas derived from this study are further compared to CGI2 and GLIMS in the Mushketoy and Tuomuer Glacier Catchments (Figure 4 and Table 4).The Mushketoy Glacier Catchment mainly faces west in Kyrgyzstan.They are typically hanging glaciers.The total glacier areas are 236 km 2 , including 23 km 2 (9.8%) of debris-covered glaciers according to GLIMS.The glacier outlines derived in this study have a good spatial agreement with GLIMS (Figure 4a) except for the debris-covered glaciers at the lower edges and about 30 small isolated hanging glaciers labeled in GLIMS and identified in Google Earth images at the northern sides.There are 222 km 2 of clean ice derived in this study using the same Landsat image acquired on 18 August 2002 as GLIMS, and their agreement is 97.7% by spatial intersection (Table 4).The Tuomuer Glacier Catchment lies in China and mainly faces south and southwest (Figure 4b).They are valley glaciers and have 116 km 2 (18.8%) of debris-covered glaciers among the total 613 km 2 glaciers according to CGI2 (Table 4).The glacier outliness derived in this study have an agreement of 93.4% by spatial intersection with the clean ice of CGI2 after excluding the debris-covered glacier, while their total clean-ice areas have an agreement of 98.8% (Figure 4b).
In the entire Central Tianshan, the clean ice areas derived in this study in 2007 had an agreement of 88.6% by spatial intersection with those of CGI2 and GLIMS after excluding the debris-covered glaciers, which were 729 km 2 and 13.5% of the total glacier areas (Table 4).A similar ratio of debris-covered glaciers (14%) was reported in this region according to CGI2 (Guo et al., 2015).Their total clean-ice areas in 2007 were almost the same (Table 4, Figure 5).The complex topography, shadowing, seasonal snow and cloud cover together make it challenging to accurately delineate the glacier outlines at the regional scale.

Glaciers Changes
There are over three thousands of glaciers in Central Tianshan.The small glaciers dominate the number (i.e., count), while the large glaciers have the most areas (Figure 6).Those with areas < 1 km 2 contribute 75.5% in terms of the number of glaciers, but only 11.2% in terms of area.Those with areas > 10 km 2 contribute 2.4% in terms of the number of glaciers and cover 48.3% of the area.Both the glacier number and area have been changing dynamically.According to GLIMS+CGI1 and

Glaciers Changes
There are over three thousands of glaciers in Central Tianshan.The small glaciers dominate the number (i.e., count), while the large glaciers have the most areas (Figure 6).Those with areas < 1 km 2 contribute 75.5% in terms of the number of glaciers, but only 11.2% in terms of area.Those with areas > 10 km 2 contribute 2.4% in terms of the number of glaciers and cover 48.3% of the area.Both the glacier number and area have been changing dynamically.According to GLIMS+CGI1 and GLIMS+CGI2, the glacier areas have shrunk −25.1% from 7233 km 2 in the 1970s to 5415 km 2 in 2007 in Central Tianshan (Table 5).Since both glacier inventories are derived from different sources and methods, it is difficult to assess the exact changes of glaciers by their direct comparison.
GLIMS+CGI2, the glacier areas have shrunk −25.1% from 7233 km 2 in the 1970s to 5415 km 2 in 2007 in Central Tianshan (Table 5).Since both glacier inventories are derived from different sources and methods, it is difficult to assess the exact changes of glaciers by their direct comparison.According to the glacier outlines derived from Landsat images by using the same approach in this study, the glacier count (>0.01 km 2 , clean ice) dramatically declined from 3629 to 2920 (count), and glacier 2D area decreased by −8.1% from 4997 km 2 to 4593 km 2 between 2015 and 1989 (Figures 5  and 6, Table 5).The number of small glaciers with area < 0.05 km 2 have the largest reduction from 1989 to 2002, and the total areas of the small glaciers also shrunk at a much large rate than the large ones (Figure 6).During those four periods, the total glacier areas in Central Tianshan were relatively stable (−1.8%) during the first 13 years between 2002 and 1989, and decrased by −3.8% during the 5 years between 2007 and 2002, and −2.7% during 8 years between 2015 and 2007.The overall retreating trend and varying rates of glacier areas can be explained by the increase trend and fluctuations of air temperature and precipitation (Figure 7).It was a relatively wet and cool decade in the 1990s, especially the record low air temperature and high precipitation 1996, leading to stable glaciers at the regional scale between 2002 and 1989; in contrast, the record high air

Glacier 3D Area Computation
Glacier 3D areas in this study are computed by the most accurate TIN 3D surface model [44].Zhang et al. [45] reported that the raster 3D area is −1.4% smaller than the polygon 3D areas computed from a TIN surface model.Our investigations confirm that the raster 3D areas tend to

Glacier 3D Area Computation
Glacier 3D areas in this study are computed by the most accurate TIN 3D surface model [44].Zhang et al. [45] reported that the raster 3D area is −1.4% smaller than the polygon 3D areas computed from a TIN surface model.Our investigations confirm that the raster 3D areas tend to produce smaller 3D areas in two steps (Table 6).In the first step of converting the glacier polygon into glacier raster, e.g., in 2015 over Central Tianshan, the converted raster 2D areas (4341 km 2 ) were −5.5% (−252 km 2 ) smaller than the polygon 2D areas (4593 km 2 ).Consequently, in the second step, the raster 3D areas computed from the raster 2D areas (4341 km 2 ) were −5.9% (−352 km 2 ) smaller than the polygon 3D areas.This results in a much smaller relative difference (RD = 22.6%) of the raster 3D areas compared to their initial 2D polygon areas (4593 km 2 ), while having slightly smaller RD (29.7%) than the TIN polygon 3D areas if they are compared to the converted 2D raster areas (4341 km 2 ).A much larger area loss occurred during the polygon to raster conversions for steep terrains and for smaller polygons due to the fact that grids located cross the polygon boundary are ignored during the conversion [45].Therefore, it is evident that the polygon 3D area based on the TIN surface model in this study is more accurate than the raster 3D areas, especially for glacier analysis clipping into different elevation bands, slope zones and aspect zones.Note: * The second ratio is computed against the original glacier polygon 2D areas, while the first ratio is compared to the converted raster 2D areas.

Glacier Area Change in Slope Zones
The ratio of the glacier 3D surface extent to the 2D planar areas is determined by and increases with slopes (Figure 8c).Geometrically, the 3D area is 15% larger than the 2D area when the slope angle exceeds 30 • and even doubles as it exceeds 60 • .For example, there are 44.8% of glaciers' 2D areas distributed on slope angles larger than 30 • , resulting in total of 30.3% larger 3D areas than the 2D area in Central Tianshan in 2015 (Table 5, Figure 8a-c).The shrinkage rates of the 2D areas (−8.1%) were slightly larger than those of the 3D areas (−7.9%) between 2015 and 1989, although their absolute shrinkage areas (−404 km 2 ) were 112 km 2 smaller than the 3D area reduction (Figure 8d).The relative difference (RD) between 3D and 2D areas slightly increased from 30.1% in 1989 to 30.3% in 2015 since most of the ice melt occurred at the lower elevation areas with gentler slopes, which is supported by the much smaller RD (27.8%) for the melted ice 3D and 2D areas.The dominant reductions occurred in slope zones of 10-20 • and 40-50 • .

Glacier Area Change in Elevation Bands
Most (~93%) glaciers were distributed in the elevation bands ranging from 3000 m to 5000 m, especially 39.1% in the 3500-4000 m band (Figure 9a,b), where the 3D areas are 5%, 11%, 22% and 38% larger than their 2D areas within each 500-m elevation band (Figure 9c).From 1989 to 2015, the glaciers demonstrated an overall declining trend in those areas lower than 4500 m.As expected, glaciers in the lower elevation bands had larger shrinkage rates, with the largest retreating areas in the bands of 3000-3500 m.Within elevation zones higher than 4500 m, the glacier areas (18%) were quite stable during the study periods.5, Figure 8a-c).The shrinkage rates of the 2D areas (−8.1%) were slightly larger than those of the 3D areas (−7.9%) between 2015 and 1989, although their absolute shrinkage areas (−404 km 2 ) were 112 km 2 smaller than the 3D area reduction (Figure 8d).The relative difference (RD) between 3D and 2D areas slightly increased from 30.1% in 1989 to 30.3% in 2015 since most of the ice melt occurred at the lower elevation areas with gentler slopes, which is supported by the much smaller RD (27.8%) for the melted ice 3D and 2D areas.The dominant reductions occurred in slope zones of 10-20° and 40-50°.

Glacier Area Change in Elevation Bands
Most (~93%) glaciers were distributed in the elevation bands ranging from 3000 m to 5000 m, especially 39.1% in the 3500-4000 m band (Figure 9a,b), where the 3D areas are 5%, 11%, 22% and 38% larger than their 2D areas within each 500-m elevation band (Figure 9c).From 1989 to 2015, the glaciers demonstrated an overall declining trend in those areas lower than 4500 m.As expected, glaciers in the lower elevation bands had larger shrinkage rates, with the largest retreating areas in the bands of 3000-3500 m.Within elevation zones higher than 4500 m, the glacier areas (18%) were quite stable during the study periods.

Glacier Area Change in Aspect Fields
Aspect fields are closely related to the intensity of surface solar radiation, which determines the melting rate of glaciers to some degree.The northern aspect receives the least surface solar radiation, leading to dominant existing glaciers (29.1% in 2015) in the long run, while the southern aspect fields have fewer glaciers (Figure 10a).The northern field also has the largest ratio (3D/2D = 1.35), i.e., the steepest glacier slope, while the western field has the smallest ratio (1.26) or the

Glacier Area Change in Aspect Fields
Aspect fields are closely related to the intensity of surface solar radiation, which determines the melting rate of glaciers to some degree.The northern aspect receives the least surface solar radiation, leading to dominant existing glaciers (29.1% in 2015) in the long run, while the southern aspect fields have fewer glaciers (Figure 10a).The northern field also has the largest ratio (3D/2D = 1.35), i.e., the steepest glacier slope, while the western field has the smallest ratio (1.26) or the gentlest slope (Figure 10c).This is not distinguishable for glacier 2D and 3D areas during the four periods, which indicates they have a similar distribution pattern within each aspect field/orientation.The northern and eastern fields had larger area reduction than the western fields (Figure 10d).The glacier areas in the southern, southwestern and western fields in 2002 were even slightly larger than those in 1989, which is consistent with the identified glacier surges, for example, Qingbintan Glacier No. 74 and Qongterang Glaciers [24].

Glacier Area Change in River Catchments
The study area is divided into five river catchments, i.e., Tekes River (North), Kumeric/Aksu River (West), Tailan River, Muzart River and Karasu Rivers (South) (Figures 1 and 11).Topographic and glaciers characteristics of the five glacier catchments within the study area of Central Tianshan are summarized in Table 7.
The Tekes River catchment has the largest land area of 8142 km 2 , the Kumerik/Aksu River catchment had the highest glacier area of 1883 km 2 in 1989, and the Muzart River catchment has the steepest terrain, with a mean slope of 32.6° (Table 7, Figure 11).The mean slope for glacier surface in the five catchments is 28.8°, leading to a 30.3% larger glacier 3D area than the 2D area.The mean glacier coverage (2D area) was 20.6% in 2015, with the largest value of 37.4% in the Muzart River catchment and the lowest of 9.3% in the Karasu River catchment.The Tekes River catchment had the second lowest glacier coverage (2D area) by 12.4% in 2015, while it had the largest shrinking rate of −17.0% between 1989 and 2015.The second largest shrinking rate occurred in the Karasu River catchments, i.e., −14.2%, followed by the Muzart River catchment: −7.7%.The glacier area in

Glacier Area Change in River Catchments
The study area is divided into five river catchments, i.e., Tekes River (North), Kumeric/Aksu River (West), Tailan River, Muzart River and Karasu Rivers (South) (Figures 1 and 11).Topographic and glaciers characteristics of the five glacier catchments within the study area of Central Tianshan are summarized in Table 7.
The Tekes River catchment has the largest land area of 8142 km 2 , the Kumerik/Aksu River catchment had the highest glacier area of 1883 km 2 in 1989, and the Muzart River catchment has the steepest terrain, with a mean slope of 32.6 • (Table 7, Figure 11).The mean slope for glacier surface in five catchments is 28.8 • , leading to a 30.3% larger glacier 3D area than the 2D area.The mean glacier coverage (2D area) was 20.6% in 2015, with the largest value of 37.4% in the Muzart River catchment and the lowest of 9.3% in the Karasu River catchment.The Tekes River catchment had the second lowest glacier coverage (2D area) by 12.4% in 2015, while it had the largest shrinking rate of −17.0% between 1989 and 2015.The second largest shrinking rate occurred in the Karasu River catchments, i.e., −14.2%, followed by the Muzart River catchment: −7.7%.The glacier area in the Kumerik and Tailan River catchments near the Tuomuer region remained relatively stable (−2%) during this period.

Glacier Mapping
Manual correction is an important step to obtain relatively accurate glacier outlines in rugged terrains.The glacier outline accuracy is affected by many factors, such as cloud, shadow, seasonal snow and debris cover.At the catchment scale, manual delineation of glacier outlines is widely utilized with the help of high resolution images and topographic data [5,39,46,47].The assessment of the GLIMS framework confirmed that artificial interpretation remains the best tool for extracting glaciers from satellite images, especially when the mapping is conducted by the same person using a combination of different types of imagery [39,48].The glacier outlines in this study were corrected by comparing the images acquired in different years, seasons and time, on a 3D surface view for shade verification and on GoogleEearth.The clean ice accuracy (90.8%) against the high resolution GeoEye01 image in 2013 is similar to that reported in the literature, e.g., 89.3% in the Muzart Glacier Catchment [20], 93% for clean ice in the Mount Manaslu in Nepal using object-based classification from Landsat images [35], ~95% for clean ice by the band ratio method in European Alps and Chugach Mountains, Alaska [49].The relatively larger difference (~10%) in this study is mainly

Glacier Mapping
Manual correction is an important step to obtain relatively accurate glacier outlines in rugged terrains.The glacier outline accuracy is affected by many factors, such as cloud, shadow, seasonal snow and debris cover.At the catchment scale, manual delineation of glacier outlines is widely utilized with the help of high resolution images and topographic data [5,39,46,47].The assessment of the GLIMS framework confirmed that artificial interpretation remains the best tool for extracting glaciers from satellite images, especially when the mapping is conducted by the same person using a combination of different types of imagery [39,48].The glacier outlines in this study were corrected by comparing the images acquired in different years, seasons and time, on a 3D surface view for shade verification and on GoogleEearth.The clean ice accuracy (90.8%) against the high resolution GeoEye01 image in 2013 is similar to that reported in literature, e.g., 89.3% in the Muzart Glacier Catchment [20], 93% for clean ice in the Mount Manaslu in Nepal using object-based classification from Landsat images [35], ~95% for clean ice by the band ratio method in European Alps and Chugach Mountains, Alaska [49].The relatively larger difference (~10%) in this study is mainly attributed to the coarser Landsat image (30 m) against the higher spatial resolution of the GeoEye01 image (4 m), which can detect more small glaciers, glacier edges and terminus than the former.Compared to the clean ice outlines from CGI2 and GLIMS by spatial intersection (Table 4), their agreements range from 88.6% over the entire Central Tianshan to 97.7% in the Mushketoy glacier catchment.Their total glacier areas have a very similar value.Overall, the derived glacier outlines in this study achieve similar accuracy as CGI2 and GLIMS.They are suitable to be used to demonstrate the general status and change rates in different river catchments and over the entire Central Tianshan from 1989-2015.

Glacier 2D and 3D Area Difference
Glacier 2D areas are the virtual/uniform surface projected to the horizontal earth surface although they are close to the actual 3D surface extents in flat terrains, such as the large ice sheets in Greenland and Antarctica.Geometrically, the ratio between 3D area and 2D area increases with slope angles at the reciprocal of cosine (slope).The ratio enlarges dramatically when the slope angle is larger than 20 • .In the steep Muzart Glacier Catchment with a mean slope angle of 32.6 • , glacier 3D areas are 39% larger than the 2D areas (Table 7).In the entire Central Tianshan analyzed in this study, the mean value of the glacier slope angle is 28.8 • , resulting in 30.3% larger 3D areas than 2D areas.Such large area differences also exist for other alpine glaciers in other steep mountains, such as the eastern, northern and western Tianshan Mountains, Altai Mountains, Pamir Plateau, Qilian Mountains and Tibetan Plateau [16].
Glacier 2D areas had been adopted in the glacier community for a long period [25], and many documents, statistics and applications had been established based on them.For example, the empirical Volume-Area (V-A) power law equations have been applied to estimate glaciers total volumes at a catchment scale worldwide by using glacier 2D areas since there were few in situ measurements of ice volume using modern techniques, such as sounding echo, ground radar or gravity methods [50,51].The V-A functions often produce large errors that may exceed 50% in volume estimation of individual glaciers [52].The glacier 2D area change does not closely correspond to ice thickness changes, e.g., an increase in the accumulation zone and a decrease in the ablation zone, resulting in even larger errors, especially in estimating the ice volume changes by using glacier 2D areas in different years [1,20].Thus, a potential application of glacier 3D areas could be applied to improve the V-A functions for a better ice volume estimate since 3D areas have glacier surface altitude information [50,51,53,54].That is under our consideration but might divert the focus of this study.
In addition, glacier 3D areas in the four periods are computed using glacier 2D areas in the four periods and the TIN 3D model built from the same GDEM data that were derived from ASTER images in 2009.The bed rock/ice surface slope is the main contribution to the difference between glacier 3D and 2D areas.It is impossible to derive glacier 3D area changes caused by glacier thickness changes using one-time GDEM data, while the TIN 3D model can describe glacier 3D area changes caused by glacier terminal retreat.There were few ice surface elevation data with high accuracy and fine resolution (30 m of GDEM) at the global/regional scale before GDEM data.Our preliminary analysis also indicates that the glacier 3D area difference due to ice thickness change between 2000 and 2009 is much smaller than the uncertainties caused by the vertical accuracy of GDEM (2009) and SRTM data (2000).It is an ideal situation that the real ice surface extent is computed by the ice surface elevation data within each investigation period.That may be possible in the near future with the rapid development of satellite lidar and radar survey and mapping techniques.

Glacier 2D Areas Changes
There were 16,000 glaciers with a total area of 15,400 km 2 before the 1970s in the Tianshan Mountains according to glaciers of Union of Soviet Socialist Republics (USSR) and the first Chinese Glacier Inventory They were relatively stable from the late 1950s to the early 1970s [10].In recent years, with the help of airborne photos and satellite images, glaciers in this region have been intensively investigated with relatively high confidence at different spatial scales and between various temporal periods (Table 8).The glacier number and areas show an overall declining trend with large spatial heterogeneity and temporal variations during the past half century (Table 8, Figure 11).
Using the same object-based classification method and Landsat images, a reduction rate of −0.31% a −1 was found over the entire Central Tianshan between 1989 and 2015.The largest retreat rate (−0.76% a −1 ) occurred between 2002 and 2007, when the annual mean air temperature increased from 11.4 • C to 12.2 • C and annual precipitation decreased from 107 mm to 29 mm (lowest record since 1986) at the southern foot (Aksu and Paicheng) of Central Tianshan (Tables 5 and 8, Figure 7).The air temperature and precipitation changes were traced back to 1996, when it recorded the highest annual precipitation (186.2 mm) since 1955 and lowest temperature (9.7 • C) since 1986.Air temperature and precipitation are the two dominant factors controlling the advance or ablation rate of glaciers in Central Tianshan, and the area reductions always lag behind the air temperature increase for a few years [47].
The regional mean shrinkage rate (−0.31% a −1 ) considers all situations in Central Tianshan, while large spatial heterogeneity exists in each catchment (Table 7 and Figure 11).The northern Tekes River catchment had the largest shrinkage rate, followed by the southeastern Karasu and Muzart River catchments, and those areas had larger changes and risk in the future perennial snow/glacier fed runoff [2,3].In contrast, glaciers in the Kumerik and Tailan River catchments in the Tuomuer regions did not show evident change.

Summary
glacier outlines derived by the object-based classification approach in this study have good agreement with CGI2 and GLIMS.They are suitable to be used to investigate the general status and change rates during different periods over the entire Central Tianshan.The glacier number and areas had an overall declining trend with large spatial heterogeneity and temporal variations during the past half century.Glacier 2D areas were reduced by −8.1% between 1989 and 2015 in Central Tianshan.They were relatively stable during 1989-2002, and a large reduction occurred during 2002-2007, which could be explained by the increasing air temperature and decreasing precipitation.The northern Tekes River catchment had the largest shrinkage rate of −17.0% between 1989 and 2015, followed by the southeastern Karasu River −14.2% and Muzart River −7.7% catchments.In contrast, glaciers in the Kumerik and Tailan River catchments in the Tuomuer regions did not show evident change (−2%).
Glacier area distributions and changes were closely related to topographic factors such as slope, elevation and aspects.The northern aspect fields receive the least surface solar radiation, leading to dominant existing glaciers.The glaciers were near evenly distributed in slope zones of 0 • to 50 • , and only 4% occurred in areas with slopes larger than 60 • .The area-weighted mean slope angle is 28.8 • , resulting in 30.3% larger 3D areas than 2D areas.The slope zones of 10-20 • and 40-50 • had dominant area reduction.Two thirds of the glacier area occurred in elevation bands of 3500-4500 m, and 85% reductions occurred in the elevation bands of 3000-4000 m.Glaciers located above 4500 did not show evident change.The northern aspect field had the largest glaciers (29.1%), while the southeastern fields had the smallest (7.5%).The northern and eastern aspect fields had larger area reduction than the western fields.

acquired on 29 24 Figure 1 .
Figure 1.Study areas in the Central Tianshan Mountains.The analyzed areas are constrained by the 2500-m elevation contour (dark blue line).The crossing straight lines are the outlines of four Landsat images (background, composite bands 5, 4, 3) used in 2006-2007.Two glacier catchments (white polygons: Mushketoy and Tuomur) and the red polygon in the Tuomur Glacier Catchment are used to compare the glacier areas derived from Landsat images and those from GLIMS, CGI2 products and GeoEye images, respectively.The five major river catchments (red circles) are applied for glacier area change analysis.

Figure 1 .
Figure 1.Study areas in the Central Tianshan Mountains.The analyzed areas are constrained by the 2500-m elevation contour (dark blue line).The crossing straight lines are the outlines of four Landsat images (background, composite bands 5, 4, 3) used in 2006-2007.Two glacier catchments (white polygons: Mushketoy and Tuomur) and the red polygon in the Tuomur Glacier Catchment are used to compare the glacier areas derived from Landsat images and those from GLIMS, CGI2 products and GeoEye images, respectively.The five major river catchments (red circles) are applied for glacier area change analysis.

Water 2018 ,
10, x FOR PEER REVIEW 11 of 24

Figure 4 .
Figure 4. Comparison of glacier outlines of GLIMS and CGI2 (white polygon) against those derived in this study (yellow polygons and background Landsat images in composite bands 5, 4, 3) in the (a) Mushketoy and (b) Tuomuer Glacier catchments from Landsat images acquired on 18 Aug 2002 and 24 August 2007, respectively.The debris-covered glaciers were 23 km 2 (9.8%) and 115 km 2 (18.8%) in both glacier catchments according to GLIMS and CGI2.

Figure 4 .
Figure 4. Comparison of glacier outlines of GLIMS and CGI2 (white polygon) against those derived in this study (yellow polygons and background Landsat images in composite bands 5, 4, 3) in the (a) Mushketoy and (b) Tuomuer Glacier catchments from Landsat images acquired on 18 Aug 2002 and 24 August 2007, respectively.The debris-covered glaciers were 23 km 2 (9.8%) and 115 km 2 (18.8%) in both glacier catchments according to GLIMS and CGI2.

Figure 4 .
Figure 4. Comparison of glacier outlines of GLIMS and CGI2 (white polygon) against those derived in this study (yellow polygons and background Landsat images in composite bands 5, 4, 3) in the (a) Mushketoy and (b) Tuomuer Glacier catchments from Landsat images acquired on 18 Aug 2002 and 24 August 2007, respectively.The debris-covered glaciers were 23 km 2 (9.8%) and 115 km 2 (18.8%) in both glacier catchments according to GLIMS and CGI2.

Figure 5 .
Figure 5. Comparisons of glacier 2D area (clean ice) from inventories and those derived in this study from Landsat images over Central Tianshan from 1989 to 2015.The debris-covered glaciers in CGI2 and GLIMS are manually separated following the method of Haireti et al. [9].

Figure 5 .
Figure 5. Comparisons of glacier 2D area (clean ice) from inventories and those derived in this study from Landsat images over Central Tianshan from 1989 to 2015.The debris-covered glaciers in CGI2 and GLIMS are manually separated following the method of Haireti et al. [9].

Figure 6 .
Figure 6.The count (a) and 2D area (b) histograms of glaciers derived from Landsat images in this study for four periods (1989, 2002, 2007 and 2015) over Central Tianshan within different area intervals.The numbers above each column are the four-period mean frequency percentages of glaciers count (a) and total areas (b) within each area category, respectively.

Figure 6 .
Figure 6.The count (a) and 2D area (b) histograms of glaciers derived from Landsat images in this study for four periods (1989, 2002, 2007 and 2015) over Central Tianshan within different area intervals.The numbers above each column are the four-period mean frequency percentages of glaciers count (a) and total areas (b) within each area category, respectively.
temperature and low precipitation in 2007 resulted in the fast reduction of glaciers between 2007 and 2002.With increasing precipitation and declining air temperature as in 2010 and 2013, the shrinkage rates of glaciers slowed down during 2007 and 2015.This further confirms the reliability of the glacier 2D areas derived from Landsat images by the object-based classification in this study.Water 2018, x FOR PEER REVIEW 13 of 24 temperature and low precipitation in 2007 resulted in the fast reduction of glaciers between 2007 and 2002.With increasing precipitation and declining air temperature as in 2010 and 2013, the shrinkage rates of glaciers slowed down during 2007 and 2015.This further confirms the reliability of the glacier 2D areas derived from Landsat images by the object-based classification in this study.

Figure 7 .
Figure 7. Annual mean air temperature and total precipitation from 1955 to 2015 in Paicheng and Aksu at the southern foot of the Central Tianshan Mountain.

Figure 7 .
Figure 7. Annual mean air temperature and total precipitation from 1955 to 2015 in Paicheng and Aksu at the southern foot of the Central Tianshan Mountain.

Figure 8 .
Figure 8. Histograms of Glacier 2D area (a), 3D area (b), area differences between 3D and 2D (c) and area retreat between 1989 and 2015 (d) based on the glacier outlines derived from Landsat images in this study for the four periods (1989, 2002, 2007 and 2015) over Central Tianshan within different slope angle zones.The numbers above columns in (a,b) are the frequency percentages in 2015, area ratio of 3D to 2D in 2015 (c), and (d) shrinkage rates ((2015-1989)/1989) within each slope zone.

Figure 8 .
Figure 8. Histograms of Glacier 2D area (a), 3D area (b), area differences between 3D and 2D (c) and area retreat between 1989 and 2015 (d) based on the glacier outlines derived from Landsat images in this study for the four periods (1989, 2002, 2007 and 2015) over Central Tianshan within different slope angle zones.The numbers above columns in (a,b) are the frequency percentages in 2015, area ratio of 3D to 2D in 2015 (c), and (d) shrinkage rates ((2015-1989)/1989) within each slope zone.

Figure 9 .
Figure 9. Histograms of Glacier 2D area (a), 3D area (b), area differences between 3D and 2D (c), and area retreat between 1989 and 2015 (d) based on the glacier outlines derived from Landsat images in this study for four periods (1989, 2002, 2007 and 2015) over Central Tianshan within different elevation bands.The numbers above columns in (a,b) are the frequency percentages in 2015, area ratio of 3D to 2D in 2015 (c), and shrinking rates ((2015-1989)/1989) within each elevation band.

Figure 9 .
Figure 9. Histograms of Glacier 2D area (a), 3D area (b), area differences between 3D and 2D (c), and area retreat between 1989 and 2015 (d) based on the glacier outlines derived from Landsat images in this study for four periods (1989, 2002, 2007 and 2015) over Central Tianshan within different elevation bands.The numbers above columns in (a,b) are the frequency percentages in 2015, area ratio of 3D to 2D in 2015 (c), and shrinking rates ((2015-1989)/1989) within each elevation band.

Water 2018 , 24 Figure 10 .
Figure 10.Glacier 2D area (a), 3D area (b), area differences (km 2 ) between 3D and 2D (c), and area retreat between 1989 and 2015 (d) based on the glacier outlines derived from Landsat images in this study for four periods (1989, 2002, 2007 and 2015) over Central Tianshan within eight aspect fields.The numbers within each aspect field in (a,b) are their frequency percentages in 2015, area ratio of 3D to 2D (c), and shrinkage rates ((2015-1989)/1989).

Figure 10 .
Figure 10.Glacier 2D area (a), 3D area (b), area differences (km 2 ) between 3D and 2D (c), and area retreat between 1989 and 2015 (d) based on the glacier outlines derived from Landsat images in this study for four periods (1989, 2002, 2007 and 2015) over Central Tianshan within eight aspect fields.The numbers within each aspect field in (a,b) are their frequency percentages in 2015, area ratio of 3D to 2D (c), and shrinkage rates ((2015-1989)/1989).

Water 2018 , 24 Figure 11 .
Figure 11.Glacier 2D and 3D area (km 2 ) in four years of 1989, 2002, 2007 and 2015 and their shrinkage rates between 1989 and 2015 in the five river catchments.

Figure 11 .
Figure 11.Glacier 2D and 3D area (km 2 ) in four years of 1989, 2002, 2007 and 2015 and their shrinkage rates between 1989 and 2015 in the five river catchments.

Table 1 .
Landsat image tiles and their glacier 2D area/percentage used in this study.

Table 2 .
Description of glacier products investigated in this study.

Table 3 .
Error matrix of glacier 2D areas (clean ice) derived from Landsat images in this study and from the CGI2 and GeoEye01 image in the upper Koxqar Baqi Glacier catchment.The debris-covered glaciers in this catchment were 18.9 km 2 according to CGI2.

Table 4 .
Comparisons of glacier 2D areas derived from this study against GLIMS and CGI2 in Mushketoy and Tuomuer Glacier Catchments and Central Tianshan by spatial intersection.The source images of GLIMS and CGI2 were acquired on 18 August 2002 and 24 August 2007.This study uses the same Landsat images but different classification methods to determine glaciers.According to CGI2 and Landsat images, the debris-covered glaciers were 23.0 km 2 , 116.0 km 2 and 729 km 2 in the three respective regions.Intersection Agreement = Intersection/Ref.,Total area Agreement = This study/Ref., Ref. represents the clean ice area of GLIMS or CGI2.

Table 5 .
Statistics of glacier 2D and 3D areas (km 2 ) based on GLIMS/CGI2 and those extracted from Landsat images in Central Tianshan during 1989-2015.

Table 5 .
Statistics of glacier 2D and 3D areas (km 2 ) based on GLIMS/CGI2 and those extracted from Landsat images in Central Tianshan during 1989-2015.

Table 6 .
Comparisons of the glacier 3D area computed from DEM raster and TIN polygon surface models using glacier polygons in Central Tianshan in 2015.

Table 7 .
Topographic and glacier (Gla) characteristics of the five river catchments in Central Tianshan (TS).

Table 7 .
Topographic and glacier (Gla) characteristics of the five river catchments in Central Tianshan (TS).

Table 8 .
The glacier 2D area reported from different studies surrounding Central Tianshan (TS).* Area 1 and Area 2 are the glacier 2D areas in the beginning and end years of investigations, respectively.The capital letters after Central TS are their orientation relative to Central Tianshan.Bold fonts are those reported in this study for comparison.