Located in the arid and semi-arid Central Asia, the Tian Shan glaciers are critical as “water tower”, providing water resource to millions of populations     . Among many near-parallel, west-easting trending mountain ranges of its 2500 km span, the central Tian Shan (CTS) encompasses the highest massif, the Tomur Massif (a.k.a. Jengish Chokusu Massif, its Peak at 7439 m asl.), marking the center of the contemporary glaciation with some of the largest glaciers developed in the sculpted valleys   . Local people adjacent to the massif is much dependent on snow and ice meltwater for civil, agricultural, and hydroelectrical uses. Due to its transboundary nature, freshwater scarcity in CTS could be a serious issue to trigger political dispute among neighboring countries (e.g. China, Kyrgyzstan, and Kazakhstan)  . Therefore, knowledge about glacier status and dynamics in this region is of fundamental scientific value to water management and related decision making   .
Observations from remote sensing data serve a great opportunity for glacier monitoring and change detection in recent decades. Several studies have shown the undergoing significant shrinkage of glaciers in the TS (e.g.  -  ). However, this time span of the glacial history is limited by the availability of satellite or aerial imagery data sources. In order to better understand the response of glaciers to the changing climate, efforts need to be made to reconstruct a longer history of former glacial periods. The well-preserved geomorphological features in the formerly glaciated landscapes can be the valuable proxy for this need. During the past glacial and interglacial cycles in Quaternary, the high elevations of the TS were repeatedly and extensively modified by erosional and depositional work of glacier expansions and recessions. Numerical dating of such glacial landforms has been essential to establish chronologies of Quaternary glacial advances across ranges of the TS (e.g.  -  ). In areas around the Tomur Massif in CTS, paleoglacial reconstruction has been primarily conducted on the southern flank of the mountains     . Stroeven et al.  provides an extensive glacial geomorphology map for the TS to show a comprehensive picture of the former glaciation. However, less attention has been paid to investigate the glacial history in the late Holocene (e.g. the Little Ice Age), which connects the satellite based modern era and the ancient Pleistocene glacial history. The cooling event of the Little Ice Age (LIA) around 1450 - 1850 AD marks the most recent glacier advancing and follows the switch to retreating under the warming climate   . The relatively fresh-looking depositional and erosional land features in the foreland and along the margins of the modern glaciers are natural archives of glacial dynamics at this centennial time scale. In this study, we mapped glacial landforms of the newly deglaciated landscapes at two headwater glacier basins of Mt. Harajoriha in the CTS, and the mapping results will lay a detailed geomorphological framework for the forthcoming age constraints of the late Holocene glacial activities.
The linkage between the late Holocene glacial land features and modern glaciers is often directly illustrated with the highly dynamic glacio-fluvial and geomorphic processes in a short distance (a few kilometers) extending from the glacier snout. Knowing the surface features is important to evaluate glacier status and to predict future glacier change. Glacier surface represents a complex of different types of glacier facies and materials, such as ice, crevasse, snow, slush, supraglacial debris cover and others, and reflects the atmospheric and topographic environment in forming them. In general, the accumulation zone of a glacier is mainly dominated with snow facies, while the ablation zone is characterized by exposed glacial ice. In the CTS area, the debris-cover is often present in the lower reach of glaciers. It is estimated 20% - 60% coverage of debris in ablation zones for debris-covered glaciers in the TS  . The presence and the thickness of the debris cover greatly influence the energy exchange between the atmosphere and the ice, and have been the focus of many studies in High Asia (e.g.  -  ). The classification of satellite imagery acquired at the end of the melting season has long been utilized to examine the areal variations on a glacier and to estimate the equilibrium line    . Here, we adopted an automated mapping method from Bhardwaj et al.  to test the capability of extracting supraglacial debris and other glacier facies using Landsat 8 images on the studied glaciers.
With a fusion of field investigation, satellite imagery and digital elevation models, this study aims to present the geomorphological evidence for the late Holocene glaciation and to characterize the glacier surfaces in the glaciated areas of the Mt. Harijoriha in CTS. Due to the scarcity of field data and the remoteness of the location, a lack of research exists in the study area and this work will fill the gap by presenting the first-hand observations and mapping results. The detailed glaciological and geomorphological information will enrich our knowledge of glaciers in the less-studied region of the Tian Shan and provide insights for reconstruction of young glacial events such as the LIA and the response to the climatic forcing.
2. Study Area
The Tian Shan is subject to active tectonics evidenced by several intracontinental deformations in the Asia’s interior   . The most recent uplift following the India-Asia collision resulted the Tian Shan with some of the highest mountains outside of the Himalayas. Under compressive tectonic forces, several W-E trending ranges formed across the region and are separated by intermountane basins  . The central part of the TS mainly comprises the Tomur, Khan-Tengri, and Haditor-Harajoriha Ranges (Figure 1), with average altitudes of 6600 m, 6000 m, and 5400 m above sea level (a.s.l.), respectively  . Along with the high elevations, the interplay of major atmospheric circulations, namely the westerlies, the Siberian high, and Asian monsoon, prescribes the perfect conditions for glacier formation and evolution in the CTS. The massif around the Peak Tomur (7439 m a.s.l.), the highest peak of the TS, is recognized as the
Figure 1. Regional setting of the mountain ranges in central Tian Shan and the location of the study area in the Harajoriha Range (red box). The terrain background on the main map was derived from GTOPO30 (~1 km resolution) developed by US Geological Survey; the dataset of glaciers was extracted from the Randolph Glacier Inventory Version 6.0; the upper-left inset map was created with ArcGIS basemap of the World Topo Map (Sources: Esri, DeLorme, HERE, TomTom, Intermap, increment P Corp., GEBCO, USGS, FAO, NPS, NRCAN, GeoBase, IGN, Kadaster NL, Ordnance Survey, Esri Japan, METI, Esri China (Hong Kong), swisstopo, MapmyIndia, and the GIS User Community); and country boundaries were downloaded from Natural Earth at www.naturalearthdata.com.
center of modern glaciation in the TS   . The orographic blocking effect on the westerlies causes the spatial variation of precipitation on the windward side and leeward side of the mountains  . The annual precipitation at the equilibrium line altitude (~4500 m) of modern glaciers in the CTS can be as much as 1000 mm  , while it decreases at a lapse rate of 30 mm per 100 m descending to piedmont at valley mouth  . Snow or rain fall primarily between May and September; therefore, glaciers experience both ablation and accumulation in summer   
The Harajoriha Range is located to the northeast of the Peak Tomur and is within the Xinjiang Uyghur Autonomous Regions of China. The mapped area is about 489 km2 and is focused on the upper watershed of two rivers in this range: the Xiata river (a.k.a. the North Muzart River) to the northern slope of the mountain ridge and the Muzart River to the south (Figure 2). The Muzart Pass at 3582 m a.s.l. saddles across the two drainage basins and it has been a historical route for travelers to abridge north Xinjiang and south Xinjiang  . On the northern slope, the glacier catchment at the head of the Xiata valley (referred as the Xiata basin hereafter) supplies snow/ice meltwater to the Xiata River which
Figure 2. Studied glaciers and glacier basins in a Sentinel-2 scene (acquisition date: 9/17/2017, bands 7-4-3). Glacier outlines are obtained from the Randolph Glacier Inventory Version 6.0 (www.glims.org/RGI/).
flows northward to join the Tekes River and then the Ili River; on the southern slope, the Muzart glacier catchment (the Muzart basin hereafter) is the source area of the Muzart River which flows southward and then eastward to ultimately join the Tarim River. Both river systems deliver vital water resource to the socioeconomic and agricultural development in the downstream villages.
A total of 70 glaciers exist in these two headwaters. They are distributed at elevations ranging from 2700 m a.s.l. to more than 6000 m a.s.l.. The highest peak reaches up to 6553 m a.s.l. (Figure 2). The modern equilibrium line altitude (ELA) is estimated between 4200 m and 4500 m for glaciers at the Muzart basin   . The largest glaciers in the Xiata and Muzart basins are the Aerqialeteer Glacier and the Muzart Glacier, respectively. Both are compound valley glaciers with several tributaries extending from cirques along the mountain ridgelines. Individual cirque glaciers, hanging glaciers and small valley glaciers are also present. Unlike the anciently glaciated downstreams, upper valleys are characterized with newly sculpted steep mountain sides, showing fresh signs of erosion and associated mass wasting processes and subsequent sediment deposition in the valley bottom. As a consequence of the abundance of the erodible materials and the topography, several glaciers, especially the large valley glaciers, are covered with various amount of supraglacial debris in the ablation zone. One study  has focused on mapping and dating the past glaciers at the Muzart valley, proposing a LIA age for the moraine at the glacier front, while no studies have been conducted in the Xiata basin on the northern slope.
3. Data and Methods
Both the late Holocene glacial geomorphology and the modern glacier surface classification were mapped based on a combination of various remote sensing data sources and field validation. In the end of July in 2017, a field campaign was set out to the Aerqialeteer Glacier, the largest glacier in the Xiata basin. Owing to the challenging logistics and terrains, we were only able to reach the foreland of the glacier snout on foot. Field observations are taken at the marginal LIA moraines and within the proglacial area. From there, a DJI® Mavic Pro drone was deployed to capture images and videos from an aerial perspective. Information collected from the field reconnaissance provides opportunities of ground truthing to validate and refine the mapping derived from remote satellite analyses of glacial landforms and surface types in the study.
3.1. Geomorphological Mapping
The primary imagery used for glacial landform mapping was an ESA Sentinel-2 image acquired on 9/17/2017 and freely downloaded from USGS EarthExplorer. This scene was selected because it covers the entire study area with no clouds and minimum seasonal snows and the date is not far from the field campaign. Compared to other freely accessible, commonly-used satellite imagery (e.g. Landsat TM, ETM+) in geomorphological mapping (e.g.      ), the Sentinel-2 imagery with a spatial resolution of 10 m in the RGB and near-infrared bands provides a better choice for years after 2015. The false-color composite (7-4-3 bands) image was imported to the ArcGIS environment where each geomorphological feature was then digitized/vectorized on screen. High resolution Google Earth imagery was frequently consulted for cross-checking and better identification of the land features with its flexible 3D visualization. Some small-scale landforms like moraine ridges were directly mapped in Google Earth and then exported and converted from the kml format to shapefiles in ArcGIS.
The NASA Shuttle Radar Topography Mission (SRTM) Global 1-arc second (~30 m resolution) Digital Elevation Model (DEM) was applied with hillshade effect to enhance the visualization of the terrain. It was used to help identify and correct the boundary of some glacial landforms with relief changes. It was also the input data for the derivation of basin boundary, contours, summits, and river channels in ArcGIS. Ridgelines were extracted automatically from the inverse DEM with the terrain analysis tool in SAGA. All of these topographic features were checked visually with the overlay of the Sentinel-2 image, and some were smoothened, truncated, and removed to create a clear representation of the topographic context. The mapped glacial geomorphological features include glacial valleys, moraine complex, moraine ridges, and trimlines, and are overlain on the grey-scale hillshaded DEM. All datasets and vectorized features were projected to the Universal Transverse Mercator (UTM) coordinate system Zone 44N. In the Xiata basin, GPS waypoints and panoramic photos from the field serve as validations for some on-screen interpretations of land features, especially when features have a dimension smaller than 10 m. In areas without field data, a combination of all available materials, namely the Sentinel-2 image, DEM, Google Earth, and maps from previous publications, was carefully examined to help interpret and delineate the landforms.
3.2. Classification of Surface Features
Relying on the differences in spectral and thermal reflectance characteristics, many methods have been developed to detect glaciers facies (e.g.     ) and supraglacial debris cover, where the latter has been challenging due to its similarity with the surrounding terrain (e.g.      ). In this study, the classification of glacier surfaces, including clean ice, dirty ice, snow, slush zone, and debris cover, was attempted by applying an automated mapping algorithm proposed by  . The method utilizes 3 optical bands and 1 thermal band of Landsat 8 data and then calculates a band ratio for thresholding different types of the surface. The remote sensing analysis steps are outlined in a flowchart (Figure 3). Prior to the band math, 16-bit unsigned DN (digital number) values of the Landsat 8 optical and thermal bands need to be converted to the ρλ, “TOA (top of atmosphere) reflectance” and T, “at-satellite brightness temperature” with parameters from the metadata file of the Landsat scene. The conversion formulae are provided and updated by the USGS Landsat Missions, and the most recent version of the formulae was applied to obtain ρλ for optical bands of 1, 5, 6, 9 and T for the thermal Band 10 in this case (Figure 3). A late summer Landsat 8 OLI/TIRS scene acquired on 9/4/2017 was selected because of its minimal cloud cover and no snow cover to ensure only capturing the glacier features. Glacier boundaries were obtained from the Randolph Glacier Inventory Version 6.0 (www.glims.org/RGI/) and were used for subsetting the ratio image. Threshold values for dividing various surface types were slightly modified from the original divisions in  depending on this particular scene and are provided in Table 1.
4. Results and Discussion
In this section, a case study was first produced for the proglacial landscapes in the Aerqialeteer Glacier in the Xiata basin. With field observations, this helps to exemplify the glacial and morphological characteristics of the study area. Then, a composite map of the study area was presented to show the late-Holocene glacial landforms including glacial valleys, moraine complexes and ridges, and trimlines, as well as various glacier facies and debris cover in modern glacier surfaces.
Figure 3. The method flowchart of Landsat 8 image processing algorithm for mapping glacier surface types. (Equations to derive ρλ and T are provided by USGS at https://www.usgs.gov/land-resources/nli/landsat/using-usgs-landsat-level-1-data-product; modified from  ). where: = Top of atmosphere (TOA) planetary reflectance; = TOA planetary reflectance, without correction for solar angle; = Local sun elevation angle; = Pixel values; = Band-specific multiplicative rescaling factor; = Band-specific additive rescaling factor; = TOA brightness temperature (Kelvin); = TOA spectral radiance; = Band-specific thermal conversion constant.
Table 1. Classification of various surface types based on threshold values modified from  .
(* and **: the original values are 0.001 and 5000 respectively in  ).
4.1. Proglacial Geomorphology in Xiata Basin
The geomorphological configuration was investigated in detail at the Xiata site where field work was conducted. GPS tracks show the spread of the footprint at the top of the sparsely vegetated marginal moraines on both sides of the valley and across the outwash plain in front of the debris covered glacier snout (Figure 4(a)). The terminal moraine is located approximately 2 km down the clean ice
Figure 4. Geomorphology and field photos at frontal area of the Aerqialeteer glacier in the Xiata Basin. (a) Geomorphological map of the proglacial zone and glacier fronts based on field experience and satellite imagery (Sentinel-2 and Google Earth); (b) the foreland viewed in Google Earth with an image date of 9/7/2007; (c) the corresponding view from Sentinel-2 image (bands 7-4-3) of 9/17/2017; (d) a panorama view from the southern bank of the LIA moraine, showing both the clean ice glacier front and the debris-covered tongue, and other features at the proglacial area; (e) subtle moraine ridges on the northern bank of the valley; (f) several meters high ice cliff under thick debris cover, with meltwater undercuts; (g) large boulders siting on the lateral moraine at the northern bank.
edge of the glacier. Based on geomorphological analysis and numerical dating of similar moraines in the TS, previous studies have suggested the timing of such formation during the LIA (e.g.        ). These LIA moraines not only mark the extent of the most recent glacier advances, but are also used to define the boundary of a proglacial environment in which intense and variable glacio-fluvial dynamics occur (Figure 4(b), Figure 4(c);  ). The proglacial area of the Aerqialeteer glacier covers about 0.8 km2, and the outer wall of the LIA terminal moraine reaches down to 2720 m a.s.l., below the local treeline (~2900 m a.s.l.). Dense forests, dominantly composed of Picea schrenkiana, are distributed beyond the edge of the LIA moraines, whereas on top of the morainic terrains, a few pioneer species have sparsely established colonization after the vegetation elimination and creation of the moraine from the last glacier advance (Figure 4(d)). Such a contrast marks a distinct boundary on the satellite imagery. Other prominent features such as small proglacial lakes and subtle lateral moraine ridges can only be distinguished with help of higher resolution image, i.e. Google Earth and/or field panoramic photos (Figure 4(d), Figure 4(e)). The stream channel pattern is short-lived, indicating the frequent redistribution of glaciofluvial sediments and a high instability in the proglacial system. It is well reflected by comparing the 9/7/2007 Google Earth imagery with the 9/17/2017 Sentinel-2 image, which demonstrates a disappearance of a large proglacial lake (about 26000 m2) within the period (Figure 4(b), Figure 4(c)). Such dramatic changes could impose a potential glacial lake outburst flood (GLOF) hazard to the downstream. The Aerqialeteer Glacier has two tributaries contributing meltwater to the runoff in the proglacial area, but the mechanism and processes are remarkably different due to the difference in surface type (discussed in Section 4.3). The southern branch is a debris-covered glacier with stagnant ice bodies buried underneath the heavy amount of debris at the lower end, whereas the northern branch is a clean-ice active glacier. Several exposed ice cliffs are witnessed in the field at various parts underneath the surface debris cover which extends a few hundred meters further than the tongue of the clean ice glacier branch (Figure 4(d), Figure 4(f)). For ice cliffs next to the main drainage channel, the disintegration of the ice bodies is prominent at the base, showing active melting in the late summer (Figure 4(f)). The thickness of the debris-cover is visually estimated to be about 10 m - 20 m across various parts. The land surface patterns and processes within the proglacial environment are rapidly changing and represent a revenue that displays a continuity between the late Holocene glacial heritages and the present-day glacier characteristics.
4.2. Descriptions of Glacial Landforms
Glacial valleys are delineated to confine the extent of the U-shaped cross-profile valleys worked by the LIA glaciers (Figure 5). The large body of viscous ice is powerful to abrade the bedrock to broaden and deepen the valley, leaving the eroded materials deposited as till or moraines at the valley floor. The resulting steep, smooth valley walls and gentle valley bottom characterize glacial valleys at the highest mountain areas, whereas transitioning to the lower terrain, fluvial valleys with V-shaped cross-sections often signify the non-glacial condition. In
Figure 5. The composite map of the late-Holocene glacial landforms and surface classification in the study area.
our cases, the extent of glacial valleys is constricted to the recently deglaciated area not far (<5 km) from the margins of the modern glaciers. The maximum extent of the late Holocene portion of the glacial valleys was delineated based on the location of the fresh morainic deposition and the sharp change in slope transitioning to the relatively older deposition. Small valley channels outward the cirques are included as glacial valleys as well, and they often connect to the mainstream valley as truncated spurs (Figure 6). These young glacial valleys are primarily identified and mapped based on the DEM and Sentinel-2 imagery.
Moraine complexes and moraine ridges are positive relief formed by deposition of unstratified till materials from past glacial advances. For the purposes of describing the late Holocene glacial geomorphology, moraine complexes and ridges mapped in the two glacier basins are the first set of such landforms from the head of the valley, and the sequence of older moraines down in the valley is not discussed here. The key component in the moraine complexes category is the LIA moraines which clearly delimit the margins of the LIA glaciers. These moraines are composed of a large number of erratic boulders, gravels, and other debris, and are sparsely vegetated with limited soil development. Large boulders
Figure 6. Examples of small tributary glaciers to the mainstream Muzart Glacier as represented in (a) the SRTM DEM (with semi-transparent hillshade effect), (b) Sentinel-2 of 9/17/2017 in false color composite of band 7, 4, 3, and (c) mapping results of confined glacial valleys, moraines, and glacier facies of snow, clean ice, and dirty ice, and supraglacial debris (the same legend as in Figure 5). Location of the area is marked in Figure 5.
with height more than 2 m can be commonly found at the surface of latero-frontal moraines (Figure 4(g)). The LIA terminal moraines can stand up to 30 - 40 m above the adjacent lower terrain and are often in arc-shape encompassing a moraine-dammed proglacial lake or incised by meltwater streams. Lateral moraines are sharp-crested and linearly located along the valley sides. Moraine ridges are outlined for the most prominent crests exhibited on the well-formed marginal moraines. From the field, several parallel, shallow ridges are found on morainic deposit at the steep slope of the northern valley side of the Aerqialeteer Glacier (Figure 4(e)). Moraine ridges are good approximation of the former stillstands of the glacier during the retreat. Hummocky moraines are also included in the moraine complexes category. Across the Muzart Pass which lies between the two basins, thick till deposits form an undulating surface of ground moraine. Such complexes are likely to suggest a stagnating status of former extensive valley glaciers at this location.
Trimline is a near-horizontal linear feature along the sidewall of a glacial valley that marks the higher stand to which the former glacier was able to erode. It is an indicator of the previous thickness of the glacier. During the LIA maximum, a thicker glacier advanced down the valley and impacted the surface that was in contact with the moving ice. In the studied basins, erosional actions removed surface vegetation, and as the LIA glacier retreated and its surface was lowered, a clear boundary is exposed to separate the vegetated surface and the eroded, vegetation-free eroded region. The presence of a sharp contrast in color and texture on valley slopes allows easy identification of trimlines in the late-Holocene glacial valleys. However, LIA trimlines do not always form on valley slopes or are not discernible especially when in conjunction with moraine formation. At the field site, a section of distinct trimline was found extending about 100 m on one side of the valley, clearly illustrating the abrasive marks where the former glacier encountered the resistant bedrock (Figure 4(a), Figure 4(d)). Other trimlines are mapped from inspection of satellite imagery and Google Earth.
4.3. Glacier Surfaces
Because the dynamics of a glacier are largely determined by the mass and energy exchange between the overlying atmosphere and the glacier surface, it is an important task to discriminate zones of different surface materials. Due to the high terrain relief in the Harijoriha glacier basins, direct field-based observations or measurements from the surface are extremely difficult. Fortunately, previous studies have long explored the feasibility to distinguish surface zones using remote-sensing data and have shown applicability of such methods in good agreement with field validation (e.g.       ). The employment of a recently developed mapping algorithm on the studied glaciers using a Landsat 8 satellite image acquired near the end of the ablation season was able to generate a pixel-based classification of different glacier surface types, including clean ice, snow, dirty ice, slush zone, and debris covers  .
The output pattern of each distinct type of glacier surfaces reflects the general zonation of accumulation and ablation in glacier mass balance. As shown on the classified map (Figure 5), glaciers in both the Xiata basin and Muzart basin are predominantly composed of clean ice at the surface, indicating the ablation-induced mass loss exceeding the mass gain through the accumulation of snow. This can be explained from the acquisition date (9/7/2017) of the selected Landsat 8 image. This time of the year corresponds to the end of the ablation season, when the accumulation of snow from the previous winter has been mostly melted and the solid bare ice gets exposed to the surface before the new accumulation season starts. The snow zone where the winter snow does not melt throughout the summer is only present at limited extents, particularly at higher altitudes (above 5000 m a.s.l.) and north-facing slopes in the Muzart basin (Figure 5). In contrast, the highest elevations of the Xiata glaciers are below 5000 m a.s.l. and snow facies are almost completely absent at their surfaces, suggesting the lack of permanent accumulation zone at the end of ablation season. On the one hand, this observation implies that local topographic factors, such as elevation and aspect orientation, can strongly influence the evolution of glacier mass gain or loss over time. On the other hand, it may suggest other potential sources for accumulation other than the persisting solid precipitation that account for the mass input to the glaciers. More investigation is needed to examine the role of wind-blown snow or avalanching in contributing the mass gain to these shrinking glaciers. During the ablation season, the slush zone may appear as an area completely saturated with water and graded up glacier into the wet snow facies, as defined in  . The discrimination of this zone in the Harijoriha shows a small area along the margin of snow zones (Figure 5). This agrees with its concept as an ephemeral transition zone from the ice facies and the snow facies  . Crevasses were extracted when the mapping algorithm applied to glaciers in the Indian Himalaya  , however, they cannot be derived in our study area. Another prominent feature distinguished from the automated mapping is the supraglacial debris. Using the pre-defined glacier boundaries from the glacier inventory, extensive debris-covered regions were detected in the ablation areas of some valley glaciers (Figure 5, Figure 6). Dirty ice pixels, which delineate the ice mixed with debris, are often found in between the debris areas and the clean-ice. Large patches of dirty ice are clearly seen at the lowest portion of the two main glaciers, the Aerqialeteer and the Muzart Glaciers (Figure 6). From satellite images and field observations, this part of the glacier demonstrates a thick and rugged debris surface with dead ice buried underneath (Figure 4(d), Figure 4(f)), whereas the debris surface with active ice up glacier is smoother and shows more streamlined flow features (Figure 5).
4.4. Quality and Uncertainty
The results are considered to provide a detailed record of the late-Holocene glacial geomorphology and glacier surface features in the studied basins of the central Tian Shan. Remote sensing data were carefully examined when applied in the manual delineation of glacial landforms and in the automated classification of surface types. The optimal selection of Landsat 8 and Sentinel-2 images with no disturbance of snow and clouds were acquired at the end of ablation season, representing a complete and accurate picture for inspection and interpretation of the glaciers and the surrounding landscape. The acquisition dates (9/17/2017 for Sentinel-2 and 9/4/2017 for Landsat 8) are close to the days of field reconnaissance of last days in July 2017, which enables a plausible comparison between field checks and image extractions. Investigation in the field is strictly limited by the difficulties of the terrain, and using panoramic photos derived from drone images lends a broader view than the sights obtained from the ground. Uncertainties from different sources may influence the accuracy of the mapped glacial imprints. First, the resolution of the DEM and satellite imagery limits the size of land features that can be identified for mapping, and some small glacial landforms, such as trimlines and narrow moraine ridges, may not be distinguishable if present. Higher resolution images, if available by Google Earth, helped to detect these features as much as possible. The mapped glacial valleys and moraine complexes are larger features (>100 m) and less affected by the resolution of data sources. Second, the automated classification based on optical and thermal bands of Landsat 8 imagery could be an effective method, but the choice of threshold values to differentiate glacier surface types is somewhat empirical and requires multiple rounds of modifications where uncertainties might be introduced. The third source for ambiguities lies in the glacier boundaries which were used for subsetting the ratio image. Uncertainties or inaccuracies in the RGI glacier outlines may cause misrepresentation of the actual glacier extents, and when processing for surface classifications, the results may contain some errors. For examples, the rock outcrops at the northern slopes of the Muzart basin were not excluded from the glacier extent and were consequently misclassified as debris cover due to the similarity of the two in the spectral analysis. Increasing the accuracy level of the prerequisite glacier data offers a way to improve the accuracy of the mapping results.
The contemporary glaciers and their recently deglaciated proglacial landscapes demonstrate dynamic and responsive environments in the backdrop of climate change and continuous retreat since the LIA maximum glacier advance. This work presents a complete picture of the late Holocene glacial geomorphology and modern surface characteristics in such environments for two glacier basins in the remote Harijoriha Mountains in central Tian Shan, where no previous detailed mapping existed. The Muzart basin on the southern slope and the smaller Xiata basin to the north preserve similar glacial landforms, including glacial valleys, terminal and lateral moraines, and trimlines. The mapped extent down in the valley is bounded by the young moraine ridges likely deposited during the LIA when the glacier was more extensive. The distinction of the vegetation cover between outer walls of the LIA moraines and the proglacial areas facilitates the delineation of the landforms from the remote sensing imagery. The timing of the moraine formation is yet to be determined, and this mapping work provides a framework for our forthcoming chronological study to improve the understanding of the late-Holocene glaciers in the central Asia. By automated mapping using Landsat 8 imagery, various glacier facies and supraglacial debris were generated for each glacier. The spatial distribution and comparison of snow and ice facies at two basins are indicative to an overall negative mass balance (ablation/ice zone overwhelms accumulation/snow zone) and the role of non-climatic factors such as elevation and/or aspect in surface characterization. To refine the classification result, the accuracy of the prerequisite glacier outlines needs to be improved. Field check was particularly helpful when the lower reach of the glacier is obscured by dead ice with a thick layer of debris cover. Despite the limitation and uncertainties, the dataset produced by this study presents the glaciological and geomorphological characteristics of the intensely glaciated environment and showcased the applicability of remote sensing mapping in understudied places like central Tian Shan.
This work was supported by the funds from the Natural Science Foundation of China (No. 41501067) and Scholarly Excellence grant from South Dakota State University. We are grateful to Xiulin Shen, Changshan Li, Shengchun Lu, Changwei Li and Xinli Li for their devoted assistance during the fieldwork. We also thank TianXi Forestry Administration and the Xiata Geologic Park Administration for granting us permission to access and survey the field sites. Last but not the least, we appreciate the two anonymous reviewers for their comments and suggestion to improve the manuscript.
 Yao, T., Wang, Y., Liu, S., et al. (2004) Recent Glacial Retreat in High Asia in China and Its Impact on Water Resource in Northwest China. Science in China Series D, 47, 1065.
 Zhao, J., Song, Y., King, J.W., et al. (2010) Glacial Geomorphology and Glacial History of the Muzart River Valley, Tianshan Range, China. Quaternary Science Reviews, 29, 1453-1463.
 Khromova, T.E., Dyurgerov, M.B. and Barry, R.G. (2003) Late-Twentieth Century Changes in Glacier Extent in the Ak-Shirak Range, Central Asia, Determined from Historical Data and ASTER Imagery. Geophysical Research Letters, 30.
 Bolch, T. (2007) Climate Change and Glacier Retreat in Northern Tien Shan (Kazakhstan/Kyrgyzstan) Using Remote Sensing Data. Global and Planetary Change, 56, 1-12.
 Kutuzov, S. and Shahgedanova, M. (2009) Glacier Retreat and Climatic Variability in the Eastern Terskey-Alatoo, Inner Tien Shan between the Middle of the 19th Century and Beginning of the 21st Century. Global and Planetary Change, 69, 59-70.
 Li, K., Li, Z., Gao, W. and Wang, L. (2011) Recent Glacial Retreat and Its Effect on Water Resources in Eastern Xinjiang. Chinese Science Bulletin, 56, 3596-3604.
 Petrakov, D., Shpuntova, A., Aleinikov, A., et al. (2016) Accelerated Glacier Shrinkage in the Ak-Shyirak Massif, Inner Tien Shan, during 2003-2013. Science of the Total Environment, 562, 364-378.
 Narama, C., Kondo, R., Tsukamoto, S., et al. (2007) OSL Dating of Glacial Deposits during the Last Glacial in the Terskey-Alatoo Range, Kyrgyz Republic. Quaternary Geochronology, 2, 249-254.
 Koppes, M., Gillespie, A.R., Burke, R.M., et al. (2008) Late Quaternary Glaciation in the Kyrgyz Tien Shan. Quaternary Science Reviews, 27, 846-866.
 Kong, P., Fink, D., Na, C. and Huang, F. (2009) Late Quaternary Glaciation of the Tianshan, Central Asia, Using Cosmogenic 10Be Surface Exposure Dating. Quaternary Research, 72, 229-233.
 Li, Y.K., Liu, G.N., Kong, P., et al. (2011) Cosmogenic Nuclide Constraints on Glacial Chronology in the Source Area of the Urumqi River, Tian Shan, China. Journal of Quaternary Science, 26, 297-304.
 Lifton, N., Beel, C., Hättestrand, C., et al. (2014) Constraints on the Late Quaternary Glacial History of the Inylchek and Sary-Dzaz valleys from in Situ Cosmogenic 10Be and 26Al, Eastern Kyrgyz Tian Shan. Quaternary Science Reviews, 101, 77-90.
 Chen, Y., Li, Y., Wang, Y., et al. (2015) Late Quaternary Glacial History of the Karlik Range, Easternmost Tian Shan, Derived from 10Be Surface Exposure and Optically Stimulated Luminescence Datings. Quaternary Science Reviews, 115, 17-27.
 Blomdin, R., Stroeven, A.P., Harbor, J.M., et al. (2016) Evaluating the Timing of Former Glacier Expansions in the Tian Shan, A Key Step towards Robust Spatial Correlations. Quaternary Science Reviews, 153, 78-96.
 Zhao, J., Liu, S., He, Y. and Song, Y. (2009) Quaternary Glacial Chronology of the Ateaoyinake River Valley, Tianshan Mountains, China. Geomorphology, 103, 276-284.
 Zhao, J., Wang, J., Harbor, J.M., et al. (2015) Quaternary Glaciations and Glacial Landform Evolution in the Tailan River Valley, Tianshan Range, China. Quaternary International, 358, 2-11.
 Haeberli, W., Müller, P., Alean, P. and Bösch, H. (1989) Glacier Changes Following the Little Ice Age—A Survey of the International Data Basis and Its Perspectives. In: Oerlemans, J., Ed., Glacier Fluctuations and Climatic Change, Springer, Dordrecht, 77-101.
 Huang, L., Li, Z., Han, H., et al. (2018) Analysis of Thickness Changes and the Associated Driving Factors on a Debris-Covered Glacier in the Tienshan Mountain. Remote Sensing of Environment, 206, 63-71.
 Scherler, D., Bookhagen, B. and Strecker, M.R. (2011) Spatially Variable Response of Himalayan Glaciers to Climate Change Affected by Debris Cover. Nature Geoscience, 4, 156-159.
 Benn, D.I., Bolch, T., Hands, K., et al. (2012) Response of Debris-Covered Glaciers in the Mount Everest Region to Recent Warming, and Implications for Outburst Flood Hazards. Earth-Science Reviews, 114, 156-174.
 Juen, M., Mayer, C., Lambrecht, A., et al. (2014) Impact of Varying Debris Cover Thickness on Ablation, a Case Study for Koxkar Glacier in the Tien Shan. The Cryosphere, 8, 377-386.
 Collier, E., Maussion, F., Nicholson, L.I., et al. (2015) Impact of Debris Cover on Glacier Ablation and Atmosphere-Glacier Feedbacks in the Karakoram. The Cryosphere, 9, 1617-1632.
 Xiang, Y., Yao, T., Gao, Y., et al. (2018) Retreat Rates of Debris-Covered and Debris-Free Glaciers in the Koshi River Basin, Central Himalayas, from 1975 to 2010. Environmental Earth Sciences, 77, 285.
 Rabatel, A., Dedieu, J.-P. and Vincent, C. (2005) Using Remote-Sensing Data to Determine Equilibrium-Line Altitude and Mass-Balance Time Series, Validation on Three French Glaciers, 1994-2002. Journal of Glaciology, 51, 539-546.
 Bhardwaj, A., Joshi, P., Snehmani, et al. (2015) Applicability of Landsat 8 Data for Characterizing Glacier Facies and Supraglacial Debris. International Journal of Applied Earth Observation and Geoinformation, 38, 51-64.
 Yang, S., Li, J. and Wang, Q. (2008) The Deformation Pattern and Fault Rate in the Tianshan Mountains Inferred from GPS Observations. Science in China. Series D, Earth Sciences, 51, 1064-1080.
 Jolivet, M., Dominguez, S., Charreau, J., et al. (2010) Mesozoic and Cenozoic Tectonic History of the Central Chinese Tian Shan, Reactivated Tectonic Structures and Active Deformation. Tectonics, 29, TC6019.
 Han, H., Liu, S., Wang, J., et al. (2010) Glacial Runoff Characteristics of the Koxkar Glacier, Tuomuer-Khan Tengri Mountain Ranges, China. Environmental Earth Sciences, 61, 665-674.
 Xu, X., Kleidon, A., Lee, M., et al. (2010) Late Quaternary Glaciation in the Tianshan and Implications for Palaeoclimatic Change: A Review. Boreas, 39, 215-232.
 Sidjak, R.W. (1999) Glacier Mapping of the Illecillewaet Icefield, British Columbia, Canada, Using Landsat TM and Digital Elevation Data. International Journal of Remote Sensing, 20, 273-284.
 Fu, P., Heyman, J., Hättestrand, C., et al. (2012) Glacial Geomorphology of the Shaluli Shan Area, Southeastern Tibetan Plateau. Journal of Maps, 8, 48-55.
 Loibl, D.M. and Lehmkuhl, F. (2015) Glaciers and Equilibrium Line Altitudes of the Eastern Nyainqêntanglha Range, SE Tibet. Journal of Maps, 11, 575-588.
 De Jong, T., Copland, L. and Burgess, D. (2018) Changes in Glacier Facies Zonation on Devon Ice Cap, Nunavut, Detected from SAR Imagery and Field Observations. The Cryosphere, 1-28.
 Paul, F., Huggel, C. and Kääb, A. (2004) Combining Satellite Multispectral Image Data and a Digital Elevation Model for Mapping Debris-Covered Glaciers. Remote Sensing of Environment, 89, 510-518.
 Bolch, T., Buchroithner, M., Kunert, A. and Kamp, U. (2007) Automated Delineation of Debris-Covered Glaciers Based on ASTER Data. In: Gomarasca, M.A., Ed., GeoInformation in Europe, Mill Press, Rotterdam, 403-410.
 Shukla, A., Gupta, R.P. and Arora, M.K. (2010) Delineation of Debris-Covered Glacier Boundaries Using Optical and Thermal Remote Sensing Data. Remote Sensing Letters, 1, 11-17.
 Bhambri, R., Bolch, T. and Chaujar, R.K. (2011) Mapping of Debris-Covered Glaciers in the Garhwal Himalayas Using ASTER DEMs and Thermal Data. International Journal of Remote Sensing, 32, 8095-8119.
 Racoviteanu, A., Williams, M.W., Racoviteanu, A. and Williams, M.W. (2012) Decision Tree and Texture Analysis for Mapping Debris-Covered Glaciers in the Kangchenjunga Area, Eastern Himalaya. Remote Sensing, 4, 3078-3109.
 Solomina, O., Barry, R. and Bodnya, M. (2004) The Retreat of Tien Shan Glaciers (Kyrgyzstan) since the Little Ice Age Estimated from Aerial Photographs, Lichenometric and Historical Data. Geografiska Annaler Series A Physical Geography, 86, 205-215.
 Yi, C., Liu, K., Cui, Z., et al. (2004) AMS Radiocarbon Dating of Late Quaternary Glacial Landforms, Source of the Urumqi River, Tien Shan—A Pilot Study of 14C Dating on Inorganic Carbon. Quaternary International, 121, 99-107.
 Xu, X. and Yi, C. (2014) Little Ice Age on the Tibetan Plateau and Its Bordering Mountains, Evidence from Moraine Chronologies. Global and Planetary Change, 116, 41-53.
 Zeng, Q., Cao, M., Feng, X., Liang, F., Chen, X. and Sheng, W. (1984) A Study of Spectral Reflection Characteristics for Snow, Ice and Water in the North of China. In: Hydrological Applications of Remote Sensing and Remote Data Transmission, Proceedings of the Hamburg Symposium, IAHS Publication, London, 451-462.
 Engeset, R.V., Kohler, J., Melvold, K. and Lundén, B. (2002) Change Detection and Monitoring of Glacier Mass Balance and Facies Using ERS SAR Winter Images over Svalbard. International Journal of Remote Sensing, 23, 2023-2050.