A Review on the Little Ice Age and Factors to Glacier Changes in the Tian Shan, Central Asia A Review on the Little Ice Age and Factors to Glacier Changes in the Tian Shan, Central Asia

Mountain glaciers are a reliable and unequivocal indicator of climate change due to their sensitive response to changes in temperature and precipitation. The importance of mountain glaciers is best reflected in regions with limited precipitation, such as arid and semi-arid central Asia. High concentration of glaciers and meltwater from the Tian Shan contribute considerably to the freshwater resource in Xinjiang (China), Kyrgyzstan and nearby countries. Documenting glacier distribution and research on glacier changes can provide insights and scientific support for water management in central Asia. As the most recent glacial event, the Little Ice Age (LIA, approximately AD 1300–1850) signifies the cold periods prior to the warming trend in the twentieth century. Here we present an overview of topics recently studied on the modern and LIA glaciers in the Tian Shan of the central Asia. With data sets of the Glacier Inventory of China and the presumed LIA glacial extents, we applied statistical models in a case study of the eastern Tian Shan to examine the impact of local topographic and geometric factors on glacier area changes. The findings of glacier size and elevation as key local factors are representative and con sistent with other studies.


Introduction
Extending the timescale of climate variations from the past century to the past millennium allows us to see the broader context of the unprecedented global warming today. Prior to the twentieth century warming, the cold extremes of the Little Ice Age (LIA) appeared to be remote in human's memory and often easily ignored. The term 'Little Ice Age' was first introduced by Matthes in 1939 in Report of Committee on Glaciers published on Eos [1]. Although his original term to describe 'an epoch of renewed but moderate glaciation that already has lasted about 4000 years' has been overtaken as 'neoglacial', in this report, it was noted that '… the glacier-oscillations of the last few centuries have been among the greatest that have occurred during the 4000-year period…' [1]. The term has now been more formally and widely adopted to describe the period approximately from AD 1300 to 1850, characterized by lower temperature over most of the globe and growth of glaciers to a more advanced extent than that of prior and post the period [2].
Following the warm period of the 'Medieval Climatic Optimum', a dramatic series of glacier advances and retreats symbolized the cooling events of the LIA. However, despite the evidence of similar events from many places around the world, the occurrence of the LIA cannot be assumed synchronous in time and uniform across space. A large amount of evidence of the LIA was assembled initially in the regions like Europe, Greenland and the Arctic, around the North Atlantic, e.g. [3][4][5][6][7][8]. The expanded stages of mountain glaciers, cycles of excessive cold, severe droughts, hot summers, or unusual rainfall were captured in abundant historical documents such as artistic paintings, ship logs and agricultural records or diary from Mediterranean, alpine Europe and further north. For examples, see Refs. [4,6,9]. An increased variability of the climate, as well as other LIA-type events induced significant impacts on human society, as evidenced in existing historical documents. As noted in Ref. [10], the beginning of the LIA is marked by the heavy rains, severe winters and harvest fails in 1315-1319 widespread in England. But, it needs to be noted that firstly, the timing of cold conditions occurred differently from region to region, and secondly, throughout the span of the LIA, the climate was never monotonically cold or always favourable to glacial expansion but instead with sometimes disastrous shifts between warmth and coldness at centennial, decadal, or even annual scales [11]. Although the observational record is less available outside of the North Atlantic region, it is well studied that mountain glaciers advanced far beyond their modern limits in highlands of Asia, the Andes of South America, New Zealand, western North America and other ranges, e.g. [12][13][14][15][16][17][18]. Grove's book [2] provides a summary of the LIA-type events from all major regions over the world.
In this chapter, we focus on the Tian Shan range in central Asia, a less-studied region compared to other key mountains on the Earth, to synthesize (1) the critical regional settings of the Tian Shan glaciers; (2) current documentation and identification of modern glaciers and the LIA glacial events; and (3) the influence of climate and local factors to glacier change since the LIA. With LIA glacial chronologies established at a few sites across the Tian Shan, it makes the examination of the spatiotemporal pattern of the LIA glacial advances possible. The spatial variation of glacier changes reflects its response to climatic shifts as well as the local topography and geometry. What modifications in climate systems occurred here during the LIA? What local factors be attributed to explain the glacier variability? This chapter provides a holistic review of such questions based on current literature. At the end, we took an example of a dataset of 865 glaciers in the eastern Tian Shan to examine the significant local factors to glacier changes using random forest model, and we compared the results with previous findings.

Glaciers in the Tian Shan
As located in the arid and semi-arid environment, glaciers in high mountains in central Asia are sensitive indicators of climate change, as well as important freshwater storages in the region. Major glacier-covered mountains in central Asia include the Altay, the Urals, the Tian Shan, the Kunlun, the Karakoram and the Himalayas (Figure 1), and they represent a wide range of climatic and hydrological conditions. The Tian Shan is approximately located at 40.5°N to 43.5°N and 75.5°E to 94.8°E.

Geographical and climatic settings
'Tian Shan' means 'Heavenly Mountains' (Shan = Mountains) in Chinese. Formed by the collision of the Indian and Eurasian continental plates about 40 to 50 million years ago, the Tian Shan is a ~2500 km long mountain series stretching from the western boundary of Kyrgyzstan across most of the Xinjiang Uyghur Autonomous Region, China (Figure 1). It is the largest mountain chain in the world's midlatitude arid region. Its highest peak is Tumur (other names: Victory Peak, Jengish Chokusu, or Pobeda), 7439 m above sea level (a.s.l.), on the border between easternmost Kyrgyzstan, Kazakhstan and China. Bounded by the Taklimakan Desert to the south and the Gurbantonggut Desert to the north, the Tian Shan indispensably serves The Tian Shan with glacier data downloaded from the GLIMS Glacier Database (http://www.glims.org/download/); black dots represent the study sites with the LIA moraine identified by numerical dating methods [56][57][58][59]. (c) Study area of the eastern Chinese Tian Shan, with modern and LIA glacier data in three sub-regions.  (Figure 1).
The geographical location at the centre of the Earth's largest continent determines its typical continental, temperate climate in the Tian Shan. The confluence of major climate systems makes central Asia a transitional region, particularly sensitive to changes in the spatial pattern of climate systems. The overall continentality is characterized by sharp local differences. The elevation is the main factor influencing air temperature distribution rather than longitude or latitude. The contrasts in seasonal and diurnal temperature are high. The annual amplitude in monthly temperature of up to 38°C is common in the Kyrgyz Tian Shan, but regional differences exist [19]. Precipitation varies both longitudinally from west to east and from the northern slope to the southern slope. Such features are closely associated with two large-scale atmospheric circulations dominated in this area: midlatitude westerlies and the Siberian High. The moisture from the Atlantic Ocean and closed drainage basins such as the Aral, Black and Caspian seas is carried by the westerlies which delivers abundant precipitation to the western end of the mountains [20][21][22], while the orographic effect of the mountain barrier reduces it to the minimum amount at the eastern end. Sorg et al. [23] summarized that on the windward northwestern slopes, the annual precipitation can be up to 1500 to 2000 mm in high elevations, whereas to the east in the interior regions, it can be as low as 100 mm. In general, the Siberian anticyclonic circulation (high pressure) dominates the range in winter season with low temperatures helpful to maintain glaciers; after March, the westerlies associated with mid-latitude cyclones convey precipitation in form of snow to favour the growth of glaciers.

Status of current glaciers
A high concentration of glaciers in the Tian Shan is known as the "Water Tower of Central Asia" [23]. They play an important role in water cycle in this arid environment. The meltwater runoff into transboundary rivers, such as the Syr Darya River, the Ili River, the Kaidu River and the Urumqi River, contributes substantially to the freshwater resource in downstream ecosystems and for over 50 million populations in the central Asia region [23][24][25]. While melting glaciers in a first phase release an increasing amount of water, the undergoing reduction of glacier volume will eventually reduce water availability and induce potential political disputes on water allocation issue [26]. Documenting the distribution of present-day glaciers and studying the history of glacier changes in the Tian Shan are fundamental and critical for current and future development in this region.
Existing records of glaciers in the Tian Shan have been derived from aerial photographs, topographic maps and more common nowadays, satellite images, e.g. [27,28]. The key role of remote sensing in glacier monitoring has been widely recognized, especially in areas like remote mountain ranges in the Tian Shan where the field-based glaciological survey is difficult to conduct. Two main parameters of glaciers obtained from remote-sensing data are the surface area and elevation, which have been commonly used to derive changes of glacial extents and thickness through time. Two international projects, World Glacier Monitoring Service (WGMS) [29] and Global Land-Ice Measurements from Space (GLIMS) [30], have been long devoted to create worldwide glacier inventories and to provide easy access to standard data for the public.
In the Tian Shan, glacier data have been collected both from the Kyrgyz side and the Chinese side. Aizen et al. [31] mapped glaciers in the Kyrgyz Tian Shan based on remote-sensing data and identified 7590 glaciers with a total area of 13,271 km 2 and an estimated 1840 km 3 volume of ice. The Chinese Academy of Sciences compiled the First Glacier Inventory of China (GIC) using topographic and aerial photographs acquired during the 1950s-1980s [32] and updated to the Second GIC in 2014 based on mostly Landsat images acquired between 2006 and 2010 [33]. This dataset is accessible via the Cold and Arid Regions Science Data Center at Lanzhou (http://card.westgis.ac.cn/), as well as through the GLIMS Glacier Database (http:// glims.colorado.edu/glacierdata/). The GIC data followed the GLIMS guidelines and contain attributes such as glacier name, glacier identifier, drainage code, glacier area, absolute and relative accuracies, debris-covered area and many other parameters. According to the Second GIC, the total number of glaciers in the Chinese Tian Shan is 7927, and they cover an area of 7256 km 2 and an ice volume of 781 km 3 . Such inventory data sets provide great convenience to know the status of present glaciers in the Tian Shan. Due to the fact that glaciers are changing constantly, successive inventories are anticipated to keep tracking the status of glaciers, under the potential of better available data and methods in future.

Evidence of LIA glaciers
Two distinct time periods, the 'Medieval Warm Period' and the 'LIA', represent unique but opposite climate conditions during the last millennium, and they set contemporary glacier changes into a long-term context [6]. Arguments and inconsistency exist on the accurate definition of the two climate episodes. Relatively abundant records in Europe and North America help identify such millennial or centennial glacier activities, but meanwhile, reflect a bias in the spatial representation when limited data are available in other key regions like central Asia. Until recent decades, studies on reconstructing past climate and glacier variations in the Tian Shan started to sprout, thanks to the development of techniques and better accessibility to the place.
In the Tian Shan, the terminal moraine characterized with no or little vegetation cover, massive piling of loose tills, sharply-crested ridge and a location a few hundred meters apart from glacier front is often believed as the mark of the LIA maximum extent [17,[34][35][36] ( Figure 2). Numerical dating using proxy records is the approach for assigning the formation ages of such moraines and thus the timing of glacier fluctuations. These proxy materials include erratic boulders, buried wood, trees and lichens. Lichenometry is the dating method that derives moraine age from measuring the diameter of the largest/oldest lichen and uses the lichen's growth rate to infer the time of the moraine's formation. Although this method has been adopted since the 1950s, the difficulties such as species identification, establishing accurate growth curve and lack of knowledge of its colonization time on boulders creates many uncertainties and limits its use, e.g. [37][38][39][40]. Radiocarbon dating of organic matters, such as buried tree trunks that have been incorporated in moraine sediments or transported tree logs that are deposited on till plain or outwash plain, can help identify the maximum age of moraine formation. The drawback of this method is the large uncertainties associated with the age [41][42][43]. In the case of the Tian Shan, lack of organic materials within or nearby moraines often prevents the application of radiocarbon dating. Dendrochronology, which uses tree-ring crossdating technique, is another method to date glacial moraines and can often provide precise age control [44][45][46][47]. The year of the innermost ring of living trees that grow on the moraine can indicate the minimum age of the moraine. Tree-ring dated outermost ring of glacially killed trees can be used to cross-validate radiocarbon results. But, the challenge is also the availability of the dating materials: trees are rarely present at margins of glaciers in the Tian Shan. Since the 1990s, cosmogenic exposure dating technique has allowed great improvement in dating glacial landforms. It measures the concentration of the targeted isotope that is only produced by the interaction of cosmic rays with minerals in rocks and then calculates the surface exposure time of the rock [48][49][50][51]. Nuclide isotope 10 Be is the most widely used one in the application of reconstructing glacial chronology [52]. The advancement in techniques has made this method dating young event on timescales of 100 years possible, for examples, in New Zealand [16], Greenland [53], Switzerland [54] and other places including the Tian Shan [55][56][57].
Several studies focused on the direct dating of the LIA moraines in the Tian Shan are marked on Figure 1. As the climatic conditions shifted with cycles of coldness and warmth during the entire LIA, glaciers in the Tian Shan have been found with two or three moraine ridges that could be identified as LIA moraines [34], indicating different periods of glacial stagnations as glaciers respond to climate changes. In the Chinese Tian Shan, the Chinese Academy of Sciences established the Tian Shan Glaciological Station near the Urumqi Glacier No. 1 in 1959, and thereafter, the Glacier No.1 has become one of a few benchmark glaciers in central Asia (Figures 1 and 2a). The dating of the LIA moraines at this site has been conducted by different groups of researchers with various dating methods. Using lichenometry, Chen [58] dated three moraine ridges with minimum formation ages of 1538 ± 20 AD, 1777 ± 20 AD and 1871 ± 20 AD, respectively. Using AMS radiocarbon dating of inorganic carbonate coating of glacial boulders on the outermost moraine, Yi et al. [59] obtained two ages of 450 ± 120 yr and 480 ± 120 yr (calibrated in Ref. [17]), which are averaged and converted to 1535 ± 120 AD, consistent with Chen's results. Our group (authors and collaborators as mentioned in Acknowledgments) conducted cosmogenic exposure dating of the boulders from the fresh moraines at the same site and another site to the west, the Haxilegen Pass (Figures 1b and  2d). The ages of seven samples from the Glacier No. 1 site clustered around 430 ± 110 yr, and the outer moraine at the Haxilegen site was dated to 430 ± 40 yr [57]. Similar ages from different methods provide reliable evidence of the LIA glacial extent, suggesting that the glacier retreated from its LIA maximum extent, 700-900 m from glacier front, approximately 430 years ago at sites in the Chinese Tian Shan. In the Kyrgyz Tian Shan, extensive work of dating LIA moraines has been done using lichenometry. Solomina et al. [36] studied the retreat of 293 glaciers across the Kyrgyz Tian Shan, and lichenometric dates of moraines indicated nearly identical maximum glacier advances occurred in the seventeenth to the mid-nineteenth centuries. Using cosmogenic exposure dating, Koppes et al. [56] measured two 10 Be ages for one boulder from the outer moraine and one from the inner moraine in the Ala-Archa Valley, the Kyrgyz Front Range. They are recalculated in [17] to 612 ± 111 yr and 284 ± 75 yr which both belong to the LIA period.
The chronology evidence of the LIA maximum extents at the selected study sites better allows us to use morphologic appearance to interpret potential LIA extent in less-studied areas in the Tian Shan. Indeed, more moraine chronological data across different mountain ranges are needed to examine the timing and the extent of the LIA glaciers and to depict the spatial heterogeneity of LIA glacial events in the region.

Climatic driving
The retreat or advance of a glacier is attributed to the amount of snow accumulation and ablation and is a result of an integrated response to climate. Temperature, precipitation, or a composite of climate inputs acts as the major driving force on glacial fluctuations: increased snowfalls or decreased temperature could favour a positive mass balance, making a glacier grow; otherwise, glaciers may recede. Although it is well known that glaciers are a sensitive indicator of climate change, the complexity of the glacier-climate relationship is not an easy puzzle to solve because glaciers are a component linked with many other components in the natural systems and the forcings that influence glacier distribution and changes operate at different scales. It can take several decades for a glacier to respond to some change in climate, and this time lag varies non-linearly and due to many other non-climatic factors such as glacier types and topography [60,61].
Like in many mountain ranges around the world, the Tian Shan glaciers are retreating or disappearing in response to the increasing atmospheric temperature in the past a few decades. Almost all meteorological stations have observed a warming trend since the 1970s in central Asia [23]. The IPCC Fifth Assessment Report AR5 [62] noted that the warming was particularly strong in winter (November to March), with a rate of 2.4°C per 50 years, in the semi-arid area of Asia. The estimate of the increase in mean annual temperature is about 0.1-0.2°C per decade [62]. Such changes in temperature cause less persistent snow accumulation and prolonged melting season for glaciers. This is likely associated with a weakening or spatial shift of the Siberian high pressure further to the east over the continental Asia in the past century [20,23].
Precipitation changes driven by the zonal and meridional atmospheric circulation patterns do not show spatially coherent trends in central Asia [28,63]. Instrumental data of recent decades revealed that the mean annual precipitation generally shows an increasing trend in both northern ranges of the Kyrgyz Tian Shan and the Chinese Tian Shan [64][65][66], but a decreasing trend in the interior range, the central Tian Shan [23]. Aizen et al. [20] argued that a strengthened westerly flow associated with the warm season precipitation (accounts for most of the annual total) could be the explanation, but the spatial distribution of precipitation decreasing from northwest to southeast reflects the effect of mountain blocking. The regions with increased precipitation undergo glacier retreating too, indicating the impact of current warming on glaciers is not compensated by more precipitation or the sensitivity of glaciers is higher to temperature change in central Asia [28,67].
The variability and long-time change of the climate system in central Asia are closely related to the interconnections with the large-scale oceanic-atmospheric circulations, such as the El Nino/Southern Oscillation (ENSO), the Atlantic Multidecadal Oscillation (AMO) and the North Atlantic Oscillation (NAO). Several studies have reconstructed millennium-long records of these oscillations and found that the La Niña condition, warm phase AMO and positive NAO are correlated with the above average winter or annual precipitation in arid central Asia [6,68,69]. Although a lack of quantity and quality of direct observations hampers the confidence in the assessment of past changes, inferred long-term trends of the indices provided possible explanations of the climatic forcing/mechanisms during the LIA. Evidence derived from climate proxy data is another common approach to extend information on the past climate conditions. Through reconstructions using tree rings, e.g. [70][71][72], ice core, e.g. [18,73], eolian sediment, e.g. [74] and lake sediment, e.g. [75], the regional climate during the LIA was depicted as cold and wet in arid central Asia. For example, the Dunde ice core record from the Qilian Mountain, ~500 km southeast of the Tian Shan, identified three long cold periods, around mid-fourteenth century, mid-sixteenth to seventeenth century and the nineteenth century during the LIA [73], which correspond well with the series of glacial advances dated in the moraine chronology [57][58][59]. IPCC AR5 [62] reported prevailed wetter conditions in central Asia during the LIA compared to the Medieval Climate Anomaly and the twentieth century. On a regional scale, the wet and cold conditions are the primary climatic driving for the LIA glacier expansions in the Tian Shan, but on a more local scale, the non-uniform response to similar climate changes implies other influential factors that need to be taken into account.

Local factors
Many non-climatic factors influence the development of each glacier and determine how each glacier responds to climate change individually and differently in future. For instance, glacier geometry and topography play an important role in the spatial pattern of glacier change variability. Such factors include glacier size, elevation range, hypsometry (areal distribution by elevation), surface orientation (aspect), slope, shape and surface characteristics (e.g. debriscover). They vary from region to region and from a glacier to another glacier. Understanding the role of the factors that lead to disparate glacier responses at a local scale will improve our knowledge of glacier-climate interaction and help predict future glacier changes.
Previous studies have explored the impact of a single topographic factor or a combination of factors on the glacier behaviour. For example, Pratt-Sitaula et al. [76] applied cosmogenic 10 Be dating on moraines at ice extent maxima in five valleys at Annapurna, Nepal, and found that the glacial asynchrony at neighbouring valleys under no spatial differences in climate is attributed to the effect of hypsometric characteristic. They argued that mountain glaciers with source areas at higher altitude are more likely to lead to glacial advance, while glaciers with lower maximum altitudes tend to retreat when facing a uniform climate change from a coolerdrier late glacial to a warmer-wetter early Holocene [76]. Numerous studies utilized satellite remote sensing data to assess the sequential changes of glaciers during past decades and revealed that in all mountain regions where glaciers exist today, glaciers shrunk considerably over the past 150 years, after the end of the LIA, and many small glaciers have disappeared [32]. Glaciers of different sizes are proven to have strong correlations with the rate of glacier changes, and small glaciers are undergoing stronger retreat compared to large-sized glaciers, e.g. [77][78][79][80]. The presence of many small glaciers today are formed due to the rapid disintegration of medium-sized glaciers and are likely to disappear if the warming trend of climate continues [80]. Topographic conditions also matter to the sensitivity and the response time of glaciers to an instantaneous change in climate. Large glaciers that have multiple tributaries and low slope gradient take a longer time to adjust their extent, whereas smaller glaciers respond faster with a higher changing rate; thus, small glaciers are more sensitive to climate change. The influence of the aspect factor can be interpreted in terms of solar radiation receipt, as well as the wind effects to areas with moderate relief [77]. According to the analysis using large data sets from the World Glacier Inventory, Evans [81] measured slope aspect of a total of 66,084 glaciers in 51 regions over the world and found a broadly consistent poleward aspect in middle and high latitudes.
A limited number of studies have been conducted in the Tian Shan to examine the relationship between local factors and glacier changes. In northern Kyrgyz Tian Shan, Bolch [78] compared glacier retreats from 1955 to 1999 among several valleys. He mentioned that the heterogeneity of glacier changes is dependent on the size, as well as the climatic regime divided among northern slopes and southern slopes. Our previous study [82] in the central Chinese Tian Shan was the first to examine the importance of topographic factors to the spatial variability of changes in glacier area and equilibrium line altitude (ELA) since the LIA. The results showed that glacier size and mean elevation range are the two key factors and could explain up to 64% of the relative area change, but the ELA change cannot be well explained by the regression of the local factors [82]. Debris cover is another factor that can have either positive or negative impact on surface ablation because it modifies ice melt rates and the heat conduction in the ablation zone [83]. The thickness of the supraglacial debris was considered to be related to the ablation rate [84]. However, how debris cover influences glaciological processes is quite complex and varies case by case [85,86]. Most glaciers in the central Tian Shan are the debriscovered type, and glaciers in the eastern Tian Shan are mostly clean-ice type.
A glacier advances or retreats as a result of both local conditions, i.e. topography and its geometry and climate conditions. Understanding these internal and external factors on glacier changes is helpful to discuss the sensitivity of glacier response and is important to interpret paleoclimatic conditions and future glacier evolution. With nearly uniform changes in climate at local scales, more research needs to be focusing on comprehensively investigating the role of local factors played in resulting non-uniform glacier responses.

A case study: using statistical models to evaluate the importance of local factors on glacier changes in Eastern Tian Shan
In this case study, we aim to investigate the relationship between the local factors and glacier changes since the LIA in the eastern Tian Shan. The study area (42°40′N-44°00′N, 85°40′E-94°50′E) is entirely bounded within Xinjiang, China, and contains several mountain ranges trending west-east (Figure 1). Glaciers exist beyond 3000 m a.s.l., where monthly mean temperatures are −16 to −13°C in winter and 3-5°C in summer [87]. We assumed that at the regional scale, the group of glaciers was impacted by climate change at a similar magnitude since the LIA, and the varied behaviour of any individual glacier reflects the influence from non-climatic factors, including the local topography and glacial geometry.

Data
The glacier change is defined as the ratio of the areal difference between LIA and the modern glaciers to the area of LIA extents. The extent of modern glaciers was obtained from the Second Glacier Inventory of China [33]. The vector shape file of targeted glaciers in the study area contains attributes of glacier area, elevations and many others following the GLIMS standards. The LIA extent of glaciers is based on the manual delineation of the fresh moraines showing similar features with the sites already dated as LIA moraines in the Tian Shan (see Section 2.3 above). Further details about the delineation can be found in Ref. [88]. In total, 640 LIA glacial extents with corresponding 865 modern glaciers are included in the study. Three sub-regions were defined based on their geographical ranges: the Boro-Eren (BE) range, the Bogeda (BG) range and the Karlik (KL) range. The relative area change for each individual glacier is the dependent variable and seven local factors (glacier area, median elevation, slope, aspect, solar radiation, longitude and latitude) are taken into account as independent variables ( Table 1). The 1-arc second (~30 m) Shuttle Radar Topography Mission (SRTM) DEM was used to calculate topographic factors that are derived based on elevation. The summary statistics of both the dependent and independent variables used in our model are shown in Table 2

Methods
To minimize the possible bias due to the non-Gaussian distribution of the data, a logarithmic transformation was applied to glacier area. To incorporate circular data of aspect into the calculation, the first harmonic of Fourier transformation was used to convert the aspect degrees to the cosine and sine components. Such a transformation on aspect data was suggested in the past statistical analysis of glaciers: the sine term represents the difference along a westeast direction, whereas the cosine term measures north-south differences [81]. The sub-region codes, 'BE', 'BG' and 'KL', were included as a dummy variable in the statistical regression.
Pearson's pairwise correlation was used to examine if each of the independent variables is significantly positively or negatively correlated with the relative glacier area change and with the others. Considering the multicollinearity between some of the independent variables as well as the possible non-linear relationship between our independent variable and the dependent variable, we used the random forest regression model to identify the factors that have the greatest impact on glacier change. Random forest is similar to regression trees in that a hierarchical classification is used to recursively split data into smaller subsets based on dependent variables. The original data set will be partitioned based on the relationships between the dependent and independent variables, forming different splitting conditions (nodes), and the process of partitioning will not stop until a perfect tree is created, or until pre-defined conditions are met for terminating the expansion of the tree. The random forest does not make an assumption regarding the distribution of the data, and as the bootstrapping procedures are used, the predictability of the random forest model is improved meanwhile the possibility of over-fitting is minimized [90]. The random forest model used the difference in mean square error (MSE) for a test sample with and without a certain variable randomly permuted to determine the importance of each variable. The difference is averaged among all the trees generated and is considered a measure of how the predictive ability is reduced when a certain variable is excluded from the model. Traditional theories concerning how the aforementioned local factors interact with glacier change are often based on scientific abduction, with limited field/lab experiment supporting them, as these types of experiments are labour-intensive, expensive, and time-consuming. We used the randomForest package [90] along with other packages in R Language (http://www.r-project.org) to do all statistical analyses in this study.

Results and discussion
Among nine variables ( Table 2), the glacier area, the longitude ('x'), the latitude ('y'), the sine aspect ('sina') and cosine aspect ('cosa') are not normally distributed (Figure 3). The longitude variable reflects the distribution of glaciers in the study area, and three clusters can be identified. The cosine component of glacier aspect shows a highly skewed distribution, as a lot of glaciers are north-facing in our study area (cosine ~ 1), whereas the sine component shows a rather uniform distribution, suggesting that glaciers in our data set do not show any clustering pattern concerning east-west facings.
Pearson's pairwise correlation shows that the glacier change is highly correlated with glacier area, elevation, longitude, solar radiation, slope and aspect at 0.01 significance level ( Table 3).
For example, the correlation coefficient between the cosine aspect and the solar radiation is as low as −0.75 (p < 0.0001), indicating a high collinearity between the variables. For traditional regression models, the highly correlated variables might result in a potential redundancy using the two variables. The receipt of solar radiation already reflects the aspect factor of a glacier: in the northern hemisphere, a north-facing glacier, which has a higher value of cosine aspect, receives less solar energy compared to a south-facing glacier. If both solar radiation and aspect factors are taken into account, the collinearity between the variables might result in biased estimation of the coefficients in traditional regression models. The random forest regression model calculates an R 2 of 0.541, which suggests that in our study area, 54.1% of the variance in glacier changes can be explained by these local factors. We used 500 trees, the default value for the random forest regression, and the model stabilized when the number of trees exceeds 200 (Figure 4). We used a linear regression between the predicted glacier changes calculated from the random forest model and the observed glacier changes to assess the model's predicative ability, and the linear regression had an R 2 of 0.54 (p < 0.0001; Figure 4). There was no apparent sign of the model over-or under-predicting the glacier change in our study area, and the standard error for our model stayed stationary for different observations.
The relative importance of each independent variable in determining the glacier change is shown in Figure 4. The importance of each independent variable in our model is determined by testing how the accuracy of the model will be affected if any individual variable is randomly permuted. The increase in mean squared error (MSE) is calculated for any individual variable when the variable is randomly permuted in the model. The importance of each variable in each tree is calculated, averaged and then normalized using the standard error.  [82] and partial least squares regression [89]. It is commonly recognized that glaciers situated at higher elevations are exposed to colder temperature, thus they tend to recede less; small-sized glaciers tend to be more sensitive to the similar level of climate change compared to larger ones, especially as they expose more proportion of its area under the ELA. Our model ranked the slope as the third most important factor controlling glacier change, followed by latitude, longitude, regional variation, and the solar radiation.  The random forest model also produces the representation of the hierarchical classifications which result in high or low glacier changes in our study area. The final regression tree created from the random forest predictions is shown in Figure 5. Based on the result, the indication of the mechanism driving glacier changes in our study area can be inferred from the split conditions at each node. We only used split conditions that are significant at the 0.05 level in this figure, since the tree can be morphologically complex when every single node is extended. Also, as the number of nodes becomes larger, our capability to interpret each node becomes very limited. Figure 5 illustrates the split conditions (nodes) in the random forest model. For example, the first split (node 1) occurs at the elevation of 3795 m; for glaciers with elevation either lower or higher than 3795 m, the next split (node 2 and 17) is based on the size of the glacier.
Low values of glacier retreat occur when the elevation is greater than 3795 m and the area of the glacier is larger than 4.039 km 2 (node 35, glacier change = 20.284%). Larger glaciers that are situated at high altitudes tend to be less sensitive to the impact of climate change, as the majority of the mass of such glaciers are experiencing accumulation, and it is likely that only limited area of such glaciers is experiencing more ablation than accumulation. For glaciers situated at an elevation higher than 3795 m with an area between 2.356 km 2 and 4.039 km 2 , the glacier change is also relatively small (28.926%). For glaciers situated higher than 3795 m but with a smaller area (≤2.356 km 2 ), a smaller slope (≤24.962°) helps limit glacier retreat. However, if the elevation of the glacier is between 3795 m and 3971 m, the latitude of glacier can make a difference in glacier area variation (41.604% for the south of 43.7° and 29.22% for the north of 43.7°.) Extreme high values of glacier retreat are observed in node 4 (79.112% retreat). This node represents glaciers in the KL region with an elevation lower than 3795 m and smaller than 0.884 km 2 . It should also be noted that even for glaciers located higher than 3795 m, high glacier retreat (67.816%, node 27) is observed. This node represents the glaciers in KL region with an area smaller than 2.356 km 2 , slope greater than 24.96° and elevation lower than 3884 m. The random forest model is also able to pick up at three split conditions (node 3,13,26) where the branches with KL region all show greater glacier changes than the other group, although the glaciers in this region show an overall smaller retreat.

Summary
Located in central Asia, one of the most extreme continental regions on Earth, the Tian Shan, is particularly important because it possesses a high concentration of mountain glaciers which contribute a significant amount of freshwater to populated areas in the lowlands across the bordered countries. These glaciers are also sensitive to climate change and respond to temporal variations in the dominance of major climate systems, such as the mid-latitude westerlies and the Siberian high pressure. The Tian Shan has received little attention compared to other key mountains in the world, partly because of the accessing challenges.
As a WSW-ENE trending ~2500 km arc of mountain system, the Tian Shan includes many ranges across the Kyrgyz Tian Shan in the west and the Chinese Tian Shan in the east. The climatic settings show a high contrast of local differences. On the whole mountain scale, the temperature distribution is quite uniform regardless the change along the elevation, while the precipitation pattern shows strong variations from west to east and from northern slopes to southern slopes. Through some international and national efforts of glacier monitoring, glacier data set/inventories have been collected in the Kyrgyz and Chinese parts using recent remote-sensing data. Overall, there are more than 15,000 glaciers in the whole range, and large glaciers are mostly found in the west, whereas a high number of smaller ones with much less area and volume of ice are found in the east.
Our focus of the Little Ice Age (LIA) allows us to place the contemporary glacier changes under the warming climate into a broader context of glacier history. Although most evidence of the LIA climate and glaciers were assembled in Europe and North America, more and more studies have been conducted in an area like the Tian Shan to reconstruct the past climate conditions and decipher the timing and magnitude of the late Pleistocene glaciations. To date, a few studies have successfully constrained the chronology of the maximum or the successional LIA glacial advances in the Tian Shan, and the dating methods include lichenometry, radiocarbon dating and most recently, cosmogenic nuclide exposure dating. LIA glacial advances occurred around 200 years ago, 450 years ago and 650 years ago based on the moraine chronologies at some sites in the Chinese Tian Shan, but the sequence of the moraine/sub-moraines are not necessarily well-preserved at the front of every glacier. Climate reconstructions from the climate-proxy data, such as tree rings, lake sediments and ice core, revealed a cold and wet climate during the LIA in arid central Asia, and the timing generally agreed with the advances of glaciers. It is well recognized that glaciers respond to changes in climatic conditions, i.e. temperature and precipitation, and during the LIA, the cold, wet climate served as the major driving force to glacier expansions. Admittedly, the mechanism of glacier-climate relationship is complex and difficult to untangle, however, even if assuming similar climate changes at local scales, the variability of glacier changes still exists and implies the influence of other non-climate driving factors, for example, the topographic and geometric settings.
In a case study, we found that the selected local factors (elevation, area, slope, aspect, solar radiation and location) can explain more than 50% of glacier changes since the LIA in the eastern Chinese Tian Shan. This high level of explained variance indicates that local geomorphometric setting is important in determining glacier behaviour and its response to climate change. Among considered factors, glacier size and elevation are the two dominant factors, and this result is in agreement with findings in previous studies. The importance of the other factors is much less than that of these two. At a sub-regional scale, a decreasing trend is observed along a west-east gradient in the correlation between longitude and area change. Such a spatial variation of glacier changes reflected by the three sub-regions might not only be attributed to the difference in their local settings, but also to the gradient in climate systems (i.e. mid-latitude westerlies and the Siberian high pressure). Random forest model is helpful for us to better understand the mechanism between glacier retreat and local factors, and based on the split conditions in random forest model, scenarios resulting in greater or lower glacier retreat are identified.
Future studies in the Tian Shan should consider focusing on several aspects: (1) utilizing high-resolution satellite images to document glacier status and recent glacier changes; (2) developing new sites and improvements in dating techniques to create more well-constrained glacial chronologies; (3) filling the spatial gaps in past climate information by adding more site-specific climate reconstructions across the range and (4) modelling the response of glaciers to climate changes with considerations of local geomorphometric factors.