Landsat 5/7 images from 2001 and Landsat 5 images from 2011 were obtained from the GEE platform to extract the glaciers in the Tianshan in 2001 and 2011. In 2021, Landsat 8 was utilized to supplement the validation process by filling in missing pixels for certain Sentinel 2 pixels where cloud interference could not be eliminated. These orthorectified images represent the surface reflectance values calculated using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) algorithm for Landsat 5 and 7 (Masek et al., 2006), and they have a spatial resolution of 30 m.

Sentinel-1 contains all of the ground range detected (GRD) scenes. Each scene has four band combinations (corresponding to scene polarization), and the possible combinations are single-band VV or HH and dual-band VV+VH or HH+HV. Sentinel-2 is a wide-swath, high-resolution, multispectral imaging mission supporting Copernicus Land Monitoring studies. The Sentinel-2 L2 data on GEE were downloaded from scihub and were computed by running sen2cor. Sentinel series data that had a spatial resolution of 10 m were used for glacier classification in 2020, and some data from 2021 were used to replace missing 2020 data.

The global 25 m PALSAR/PALSAR-2 (The Phased Array type L-band Synthetic Aperture Radar) mosaic is a seamless global SAR image created by mosaicking strips of SAR imagery from PALSAR/PALSAR-2. For each year and location, the strip data were selected through visual inspection of the browse mosaics available over the period, with those showing a minimum response to surface moisture were preferentially used. Its polarization backscattering coefficient contains HH and HV, and this dataset was used as the SAR data source for 2011.

Two sources of digital elevation and surface data were used: NASADEM and ALOS World 3D-30 m (AW3D30). NASADEM reprocesses SRTM data with improved accuracy by incorporating auxiliary data. AW3D30 is a global DSM dataset with a horizontal resolution of approximately 30 m (1 arcsec mesh). The DEM data used to extract ridge lines and calculate glacier length were AW3D30, and the DEMs used to calculate the surface elevations, slopes, and aspects of glaciers were NASADEM and AW3D30.

The Randolph Glacier Inventory (RGI 2.0 and 6.0, representing the glacier boundaries before 2000 and in approximately 2010, respectively) was applied as a reference to evaluate the glacier delineation results. RGI glacier outlines and a distance-based buffer of 100 m were used as the baseline representation of glacier area in this study. Some glaciers were missed in version 6.0, so they were supplemented with version 2.0. The glacier inventory was uploaded to the GEE Assets.

To analyse regional climate forcing, the ERA5-Land Monthly Averaged by Hour of Day reanalysis data was analysed to investigate changes in climate throughout the study. The monthly 2-m air temperature and total monthly precipitation from 1991 to 2021 were considered, and this Climate Reanalysis dataset had a horizontal resolution of 0.1°× 0.1° (the resolution on GEE platform is marked as 11,132 m).

The names of the datasets used on the GEE and their respective time scales, quantities, and resolutions can be found in Table 1.

A pixel-based machine learning classification method was used to identify and classify ice, snow, debris-covered ice and other features (water, land, vegetation, etc.) on GEE, and the whole process consisted of several significant processing steps (Figure 2). First, satellite images of selected years were filtered, and annual image composites were generated after masking clouds and cloud shadows. Second, manual visual interpretation was performed based on the annual image composites from the previous step to select feature sample points and generate training and validation datasets. Third, multiple indices were generated from the multi-sensor data and fused into the original image in the form of multiple bands together with elevation and slope data. PCA compressed many bands in the original image into fewer uncorrelated bands. The first three significant bands which could explains 95% of the overall variance were selected as the input data for the subsequent supervised classification. Fourth, four machine learning algorithms (RF, SVM, GTB, CART) were applied and evalu- ated and the one with the highest accuracy was selected. The study area was divided into 0.5°×0.5° geographic grid cells, and the classifier was trained to extract glaciers and other features based on the classification features of the training dataset. Fifth, the classification results were post-processed to eliminate the ‘salt-and-pepper noise’, and the results were compared with the validation dataset to evaluate the classification accuracy. Sixth, the changes in glacier areas were analysed and the influence of climatic and topographic factors were explored.

As previously described, 30% of the 15,000 sample points for each year were randomly selected as the validation dataset to evaluate the accuracy of the RF model. After noise removal, the OA in 2001, 2011, and 2021 were improved compared to those before filtering, with values of 93.6%-99.7%, 96.0%-99.6%, and 94.7%-99.4%, respectively, and the kappa coefficients ranged from 0.889-0.996, 0.932-0.987, and 0.903-0.987 (Figure 6). Among the sub-regions, the lowest overall accuracies were found in the Central Tianshan, with values of only 93.6%, 96.0%, and 94.7%, respectively, and the kappa coefficients were 0.889, 0.938 and 0.912. In contrast, the Eastern Tianshan had the highest overall accuracy and the highest kappa coefficients. Among the local feature classifications, the overall accuracy of ice and snow identification was the highest, at almost 100%, although there was some confusion between ice and snow. Debris were identified with the lowest accuracy, and although some debris were identified accurately, there were cases where they were identified as ice or other features; however, overall, debris could be better distinguished from other features (Table 2).

Based on the vectorized glacier profile, the uncertainty of glacier area was obtained by buffering, and we calculated the statistical error of glacier area from 2001 to 2021 for the three periods as ±353.87 km², ±310.80 km², and ±200.90 km². Specifically, 2021 had a finer pixel resolution (10 m), and therefore, the error was minor compared to that in 2001 and 2011 (pixel resolution of 30 m).

The glacier areas in 2001, 2011, and 2021 for the four sub-regions (Central Tianshan, Dzhungarsky Alatau, Eastern Tianshan, and Northern/Western Tianshan) of the Tianshan were calculated on GEE using multisource remote sensing data and analysed their changes (Figure 7). The largest regional glacier area was in Central Tianshan, which accounted for 60% of the total glacier area in the Tianshan. The glacier area of Dzhungarsky Alatau was relatively small, only approximately 7% of the size of the glacier area in the Central Tianshan. The average area of a single glacier in the Central Tianshan was also larger than that in other regions. Across the Tianshan, glacier area decreased from 13,610.33 ± 353.87 km² to 12,431.83 ± 310.80 km² from 2001 to 2011, a loss of 1178.50 km² at a rate of 0.87%/a, while glacier area decreased to 10,573.93 ± 200.90 km² from 2011, which was a rate of 1.49%/a, showing an accelerated retreat overall compared to the previous decade.

Although glacier area has decreased in all sub-regions since 2001, they have had different trend patterns. The decrease in glacier area in Dzhungarsky Alatau from 2011 to 2021 was significantly lower than that from 2001 to 2011, i.e., a ‘slower retreat.’ Although glacier area retreat has accelerated in the other three sub-regions in the last decade, the inter-annual rate of change in the central and Northern/Western Tianshan was more stable and increased less, while the rate of glacier retreat in the Eastern Tianshan increased significantly in 2011-2021 compared to the previous decade. Specifically, the glacier areas in Dzhungarsky Alatau, Central Tianshan, Northern/Western Tianshan, and Eastern Tianshan decreased by 57.29 km², 566.44 km², 442.44 km², and 112.32 km² during 2001-2011, at rates of 0.98%/a, 0.70%/a, 1.86%/a, and 0.43%/a, respectively, while the loss rates for 2011-2021 were 0.27%/a, 0.96%/a, 2.14%/a, and 2.84%/a, indicating glacier area reductions of 14.28 km², 724.24 km², 412.54 km², and 706.83 km², respectively. It is important to note in particular that the Northern/Western Tianshan has consistently exhibited a high rate of retreat since 2001, with a value of almost 2.00%/a (Table 3).

The evolution of glacier areas within each sub-region also showed significant differences. The Tianshan was divided into 0.5°×0.5° geographic grid cells and analysed the trend of glacier area change in different grid cells (Figure 8). Overall, the majority of glacier areas within the geographic grid cells show shrinkage at diverse rates of magnitude. Glacial retreat was more severe in mountain peripheries, river valleys, and gentler areas than in higher elevation mountains. The areas with the most dramatic changes in area values were located at the junction of the western part of the Central Tianshan with the eastern part of the Northern/Western Tianshan, the western part of Eastern Tianshan and the western part of Northern/Western Tianshan. Although the values of area reduction were substantial, the relative rates of change were not necessarily more significant than those in other regions, particularly in the Central Tianshan (approximately 80°E). Specifically, in the Tomur-Khan Tengri region, a positive glacier area growth of 0.2%/a-0.4%/a was observed between 2001 and 2011, while a shrinkage rate of -0.02 to -0.4%/a between 2011 and 2021, indicating an anomalous steady state. The glacier area loss in the interior of the Dzhungarsky Alatau exceeded that in its surrounding areas, yet the rate of area retreat in the surrounding areas surpassed that in the interior. A region of significant positive glacier area growth is present in the Eastern Tianshan near 88°E, and the rate of retreat in other regions has accelerated over the past decade. This acceleration in the rate of retreat is also widespread in various geographic grid cells of the northwestern Tianshan, exhibiting a more extensive distribution.

The glaciers in the Tianshan are distributed at 2200-7500 m a.s.l., with more than 95% of the glacier area mainly distributed at 3400-5400 m a.s.l., i.e., near the glacier equilibrium line. Dzhungarsky Alatau is at a lower altitude, so its glacier distribution altitude band is also lower than those of the other three regions, and its glacier terminus can reach a minimum altitude of approximately 2200 m a.s.l. (Figure 9).

However, the glaciers distributed at approximately 2200-2400 m a.s.l. in 2001 completely disappeared after 2011. The Central Tianshan has the highest average elevation and, therefore, the most elevation bands spanned by glacier distribution in its region. The glacier equilibrium line in Dzhungarsky Alatau was much lower than that in the Central Tianshan. Our results show that, in terms of altitudinal band distribution, the most drastic reduction in glacier area has been in the range of 3400-3600 m a.s.l., with 671.38 km², 569.21 km², and 417.99 km² of glacier area in the three phases, respectively, and a loss of 253.39 km² at a rate of 1.89%/a. The four sub-regions with the most severe loss of glacier area were 3800-4000 m a.s.l. in the Central Tianshan, 3000-3200 m a.s.l. in the Eastern Tianshan, 3200-3400 m a.s.l. in the Northern/Western Tianshan, and 3000-3200 m a.s.l. in Dzhungarsky Alatau. Part of the glacier elevation was more significant than 5800 m a.s.l. Changes occurred relatively steadily, where the loss was slight, and the shrinkage rate was only 0.3%/a or less, and the proportion of glacier area to the total area in this region even slightly increased. Only approximately 1.6%-1.7% of the glacier area was distributed at elevations below 3400 m a.s.l.

The median elevation is commonly used to estimate the equilibrium line height, so it describes the height position in glacier topographic analysis. Our results showed that the median elevation of glaciers of different area classes was above 4000 m a.s.l., and it is increased with the glacier area. Among the different glacier size categories, glaciers in the smallest area class (0-1 km²) had a more significant median elevation distribution but a narrower elevation range, while glaciers in the larger area class had a more comprehensive elevation range but a closer median elevation distribution (Figure 10a).

In terms of slope distribution, northwest, north, and northeast were the main directions of glacier distribution in the Tianshan, where the area of glaciers with north as the main direction accounted for approximately one-third of the glacier area, while those with the northwest and northeast directions accounted for 17.6% and 15.7% of the area, respectively (Figure 10b). It should be noted that in the area of the prominent peak of the Tianshan, the main directions of some massive valley glaciers, such as the Tomur Glacier, the North-South Inilchek Glacier, and the Tugai Berezi Glacier, were southwest, west, and east, respectively. However, most of the glaciers at the edge of the mountain range have developed towards the north, and most of the glaciers with more area loss are also in this direction. The average slope of glaciers in the Tianshan is 25°, and the larger the glacier area is, the gentler its average slope is. In terms of slope gradation, more than 5500 km² of glaciers have an average slope between 10° and 20°, accounting for approximately 40% of the glaciers, and more than 5000 km² of glaciers have an average slope gradient between 20° and 30° (Figure 10c). A tiny percentage of glaciers have average slopes less than 10° and greater than 40°, and the remaining < 1000 km² of glaciers have average slopes steeper than 30°. A significant reduction in glacier area occurred in the less than 10° and greater than 40° slope bands (Figure 10c).

Discussing the causes of uncertainty generation is essential to improve the accuracy of glacier identification and accurately verify the glacier area variation. The uncertainty in our results had two main components: the classification error when the machine learning algorithm was run on GEE and the error generated when the glacier profile was depicted.

Some studies have similarly reported that the most intense glacier retreat in the Tianshan is generally distributed in the outer mountain ranges and lower elevations near densely populated plains or river valleys, e.g., Narama et al., (2010) examined four regions and found that the most significant changes occurred in the periphery, with a 19% change (~0.63%/a) between 1970 and 2000. Bolch (2007) used Soviet glacier inventories and Landsat data to examine the Zariiski and Kungui Alato valleys in the northern Tianshan and found a >32% decrease (~0.72%/a) between 1955 and 1999, while in the more immense glaciers and higher elevation interior alpine areas, studies found that they have shown relative stability in recent years (Milner et al., 2017). These results corroborate the validity and accuracy of our work. In terms of sub-region scale, it was found that glaciers in the Eastern Tianshan and Dzhungarsky Alatau have most significant negative mass balance rate (Barandun et al., 2021). In contrast, the annual mass balance variability was lower in the Central Tianshan (Barandun et al., 2021), and the glacier area shrinkage rate in the Northern/Western Tianshan had spatial variability ranging from 0.11%/a to 3.7%/a (Zhang et al., 2021), which was more consistent with our results. Huang et al. (2021) suggested that the glacier area shrinkage rate in the Eastern Tianshan was 0.8%/a between 1990 and 2018, but our results show that the glacier loss rate in the Eastern Tianshan consistently exceeded 2%/a. Our results indicate a stable and minimal increase in glacier area in this region between 2001 and 2011, followed by a subtle decreasing trend between 2011 and 2021, which aligns with the findings of other researchers. Zhou et al. (2021) found that this region encompasses nearly all the surge glaciers and advancing glaciers in the Tianshan, with advances ranging from several hundred meters. These glaciers are larger in size, and their advancing portions contribute more significantly to glacier area change compared to smaller glaciers. Li et al. (2018) similarly concluded that the glacier mass balance in this region from 2000 to 2011 was much lower than in other regions of the Tianshan. Table 5 offers comparisons of our results in various regions with the glacier area shrinkage rate from previous studies. An accelerating trend of glacier retreat has been observed in recent years (Wan et al., 2020), while earlier studies have reported less glacier change in the past decade. Based on this comparison, the glacier area shrinkage rate derived from our study for the first and second decades of the 21st century is consistent with those found in other studies. For some smaller glaciers or those located at the periphery of mountain systems, they are more sensitive to climate change. Although their absolute area change values may not be larger than those of larger glaciers, their relative glacier area change rates will be more pronounced. The above comparisons show that our work on identifying and classifying glaciers at large regional scales with high spatial resolution on the GEE is reliable. Such a process is helpful for annual glacier mapping to obtain a quick and accurate picture of their changes.

2 m temperature and total precipitation from ERA5-Land Monthly Averaged by Hour of Day reanalysis data were chosen to analyse the distance-level changes in summer and winter mean temperature and precipitation in the Tianshan for 2001-2011 and 2011-2021 compared to the previous decade (Figure 12), where June-August were defined as summer and December-February were defined as winter. The temperature changes in the Tianshan since 2001 showed inconsistent spatial variability. Summer temperatures showed varying degrees of upwards trends, while winter temperatures showed only weak negative pitch levels from 2001 to 2011 in parts of Dzhungarsky Alatau and the Eastern Tianshan, but from 2011 to 2021, the negative pitch levels increased significantly, particularly in mountainous regions. The variation in summer precipitation exhibited considerable spatial variability. High mountainous areas of the Central Tianshan and its adjacent regions, including the Northern/ Western Tianshan and Eastern Tianshan, as well as the mountainous regions of the Dzhungarsky Alatau, demonstrated an increase in precipitation between 2001 and 2011. Meanwhile, the remaining areas displayed various degrees of precipitation decline. Precipitation in the Northern/Western Tianshan continued to exhibit a negative pitch level during the summer of 2011-2021. However, most of the mountainous regions of the Eastern Tianshan experienced increased precipitation during the same period. The trend of winter precipitation remains unclear, except for the Northern/Western Tianshan, where it transitioned from progressively wet to progressively dry. In general, summer precipitation increases slightly in mountainous areas and decreases gradually in the surrounding valleys and plains.

Summer temperatures mainly control glacier ablation and winter precipitation mainly regulates glacier accumulation (Bonekamp et al., 2019), and precipitation in the Tianshan region is mainly concentrated in summer (Guan et al., 2021a, 2021b). Although summer precipitation has increased, its relatively small proportion of precipitation has had a much smaller impact on glaciers than summer temperatures. The continued increase in summer temperatures can be considered an important driving force contributing to the continued shrinkage of glaciers in the Tianshan. Specifically, augmented summer precipitation and negligible alterations in summer temperature in the regions near 80°E in the Central Tianshan and 88°E in the Eastern Tianshan may provide a plausible explanation for the observed expansion of glacier areas in these regions.

Some topographic factors, such as elevation, slope, and slope direction, influence glacier variability (Liu et al., 2013). It has been confirmed that increases in temperature and precipitation lead to differences in elevation-related snow cover area expression (Kunkel et al., 2016), with a decrease in snow cover area at lower elevations and an increase at higher elevations, similar to the decrease in glacier area at lower elevations and the stabilization in glacier area at higher elevations shown in our results. Glaciers in the Tianshan are primarily located in the lower elevation zone (3600-4600 m a.s.l.), a region with higher temperatures and less precipitation than those observed in higher elevations, where glaciers continue to retreat, a result that is close to the results provided by Shangguan et al., (2009), who found that the most significant reductions in the glacier area occurred at elevations of 4100-4500 m a.s.l. Our results found that glaciers in smaller area classes (0-5 km²) had much larger relative area loss rates, although they were smaller than those in larger area classes (5-10 km² and >10 km²). Many small glaciers have small accumulation areas, receive limited accumulation and are more sensitive to the temperature rise, so their retreat is more severe. The convergence of several branch glaciers forms some large composite valley glaciers, and their accumulation areas are higher in elevation and widely distributed. Some accumulation areas are gentler than small glaciers, and the higher elevation areas are thinner but stable in area change when precipitation increases in higher mountain areas. However, the terminal elevations of these large valley glaciers are usually lower, and thus their glacier terminals are subject to significant ablation by increasing temperatures. The orientation of the glacier is closely related to the intensity of surface solar radiation it receives, so the average slope orientation determines the accumulation and melting rates of the glacier to some extent (Olson and Rupper 2019). The surface solar radiation received on the south side of the mountain range is much higher than that on the south side of the mountain range, and thus, more glaciers are developed on the north side of the mountain range than on the south side (i.e., glaciers with N, NW, and NE as the average slope orientation are more predominant). In addition, glaciers on the south side of the watershed in some areas have declined more than those on the north side, reflecting the difference in solar radiation received by the north and south sides of glaciers. The mean slope of glaciers varies significantly from region to region. It is mainly related to the mean slope of the catchment, with some larger glaciers having a more comprehensive altitudinal range and, therefore, a more significant mean slope; moreover, smaller overhanging glaciers may have a more significant mean slope and may be distributed over a wide range of altitudes. Some studies have also indicated that the debris coverage of some glaciers in Central Asia has increased as the glaciers continue to shrink (Scherler et al., 2011; Kääb et al., 2012; Milner et al., 2017). Debris, which usually covers the end of the ice tongue and may be continuous from higher elevations to the glacier below, also has a significant influence on the glacier ablation. Thin debris accelerates glacier melting by reducing the albedo of the glacier surface, but debris thicker than a few centimetres may inhibit glacier ablation by isolating heat exchange between ice and air. Because the statistics of glacier area were not continuous in time, the effect of surface debris is not discussed further here, although it was clear from the area change that some surface debris-covered glaciers have retreated at a faster rate than glaciers not covered by surface debris.

The analysis of the causes of glacier changes in the Tianshan is similar to the conclusions of others in that changes in the degree of glacier retreat depend on climatic and topographic conditions to some extent (Sorg et al., 2012; Farinotti et al., 2015; Kraaijenbrink et al., 2017; Pritchard, 2019). Studies have shown that since the 1980s, the Central Asian region where the Tianshan is located has experienced a transition from warm-dry to warm-wet, i.e., an increasing trend in temperature and precipitation (Zhang et al., 2019a; Yao et al., 2021). The high mountain ranges of the Tianshan intercept the moisture carried by the air currents from the mid-latitude westerly wind belt to generate precipitation, and the convergence between the westerly wind belt disturbances and the atmospheric circulation from the Siberian high pressure has a triggering effect on the precipitation in the mountains (Yao et al., 2016). It has been suggested that the Scandinavia (SCAND) teleconnection pattern represents important circulation variability affecting Tianshan summer precipitation, and the vigorous high pressure over the Ural Mountains and the low pressure over Central Asia during the SCAND negative phase in summer jointly lead to enhanced moisture transport from the Arctic Ocean to the Tianshan (Guan et al., 2021b). It has also been suggested that climate forcing is the primary driver of sensitivity to uneven changes in material balance in alpine Asia, explaining up to 60% of the spatial variation in glacier variability (Sakai and Fujita 2017), while glacier morphology has been found to explain up to 36% of the spatial mass balance variability in the Tianshan (Brun et al., 2019).