Glacier Area and Snow Cover Changes in the Range System Surrounding Tarim from 2000 to 2020 Using Google Earth Engine

: Glacier and snow are sensitive indicators of regional climate variability. In the early 21st century, glaciers in the West Kunlun and Pamir regions showed stable or even slightly positive mass budgets, and this is anomalous in a worldwide context of glacier recession. We studied the evolution of snow cover to understand whether it could explain the evolution of glacier area. In this study, we used the thresholding of the NDSI (Normalized Difference Snow Index) retrieved with MODIS data to extract annual glacier area and snow cover. We evaluated how the glacier trends related to snow cover area in ﬁve subregions in the Tarim Basin. The uncertainty in our retrievals was assessed by comparing MODIS results with the Landsat-5 TM in 2000 and Landsat-8 OLI in 2020 glacier delineation in ﬁve subregions. The glacier area in the Tarim Basin decreased by 1.32%/a during 2000–2020. The fastest reductions were in the East Tien Shan region, while the slowest relative reduction rate was observed in the West Tien Shan and Pamir, i.e., 0.69%/a and 1.08%/a, respectively, during 2000–2020. The relative glacier stability in Pamir may be related to the westerlies weather system, which dominates climate in this region. We studied the temporal variability of snow cover on different temporal scales. The analysis of the monthly snow cover showed that permanent snow can be reliably delineated in the months from July to September. During the summer months, the sequence of multiple snowfall and snowmelt events leads to intermittent snow cover, which was the key feature applied to discriminate snow and glacier.


Introduction
The Tarim Basin is the largest inland drainage basin in China, located in the northwest of the Tibetan Plateau and far away from the ocean. The basin is surrounded by the mountain ranges of the Tien Shan in the north, Pamir in the west, West Kunlun in the south and East Kunlun in the southeast [1]. These high mountains are the ideal place for glacial development. Glaciers and snow meltwater account for about 75% of the runoff [2], thus being the main contributor to water security in the Tarim Basin, particularly for irrigated agriculture along the fringes of the Taklamakan Desert [3][4][5]. Therefore, changes in upstream glaciers have an effect on the whole Tarim Basin runoff [6,7]. In the context of global warming, the annual average temperature has increased at a rate of 0.2 • C per decade in the Tarim Basin [8]. Extreme hydrological events have intensified, such as glacier lake outburst floods (GLOFs) and glacier surging [9][10][11]. Thus, it is essential to investigate glacier variability and the potential impacts on the long-term sustainability of water supply.
Taking into account that debris-covered ice cannot be detected by NDSI, we only distinguished snow and clear ice in this study. The usability of Moderate-Resolution Imaging Spectroradiometer (MODIS) images is improved by the sub-daily temporal resolution and global coverage [44]. The Terra and Aqua satellites orbit the earth and cross the equator about 3 h apart, increasing the likelihood of acquiring cloud-free images [45]. In order to diminish data volumes and to eliminate cloud-affected observations, composite products are produced from observations collected over the course of 8 days [46]. At present, the MODIS snow cover products are generated and made available by the National Aeronautics and Space Administration (NASA). The snow products generated with the MODIS data acquired from Terra (MOD10A2) and Aqua (MYD10A2) are an eight-day combined dataset available at a spatial resolution of 500 m, which contains maximum snow cover information over the compositing period [47].
Snow and glacier ice are characterized by high reflectance in the visible bands and low reflectance in the near-infrared [48], which explains why thresholding NDSI is widely used for snow cover mapping [36,49]. Mishra et al. [50] evaluated the range of NDSI values of snow in the Himalayan region and found that NDSI can vary between 0.04 and 0.92 as the amount of snow cover increases. Zhang et al. [51] suggested that the NDSI threshold of 0.35 is the optimal threshold for use in the Tibetan Plateau based on daily snow-depth measurements during 2003-2013. For global analyses, an NDSI threshold of 0.4 is still recommended [44]. In this study, we also use the standard NDSI threshold of 0.4.
Google Earth Engine (GEE) is a platform that provides a large collection of satellite data and products and allows users to visualize and analyze data through a web-based system [52]. The platform is an application programming interface, which uses a combination of JavaScripts and functions to process and analyze the data that are provided to the user. GEE is a powerful geospatial processing and analysis tool not only because it is computationally time-efficient but also due to the expanding availability of many open-source geospatial datasets. It has been widely used to map rock glaciers [53], water areas [54], wildfire detection [54], among other applications. To explore the snow cover relationship with changes in glacier area, we carried out this study using Google Earth Engine.
We used MODIS reflectance products and spatial analysis techniques to investigate the evolution of glacier area from 2000 to 2020 in mountain ranges in the Tarim Basin. According to the drainage basin delineation from the Hydrological data and maps based on Shuttle Elevation Derivatives at multiple Scales (HydroSHEDS), we divided the whole Tarim Basin into five regions: East Tien Shan, West Tien Shan, Pamir, West Kunlun, and East Kunlun. We examined the evolution of snow cover in relation to its effects on the glacier area in different regions in the Tarim Basin, and further explored whether the Pamir-W Kunlun glacier anomaly was related to the regional evolution of snow cover. We evaluated the MODIS glacier area using higher spatial resolution glacier area retrieved from Landsat 5-Thematic Mapper (TM) in 2000 and L8-Operational Land Imager (OLI) in 2020. Finally, we combined our results with digital elevation models and previous studies of glaciers to evaluate local glacial changes and their effects on water resources.

Study Area
The range system surrounding Tarim Basin (34)(35)(36)(37)(38)(39)(40)(41)(42)(43)(44) • N, 75-92 • E) covers an area of about 1,020,000 km 2 , with elevations ranging between 773 and 7500 m. The mountain streams originate from snow and glacier meltwater at high elevations. The basin is mainly fed by four major streams: Aksu River from West Tien Shan to the north, the Karakax River and Yurungkax originating in the north slope of Kunlun Shan to the south, and the Yarkant River from the Pamir Mountains [55] (Figure 1). The Taklamakan Desert is situated in the middle of the Tarim Basin, encircled by mountains and glaciers. It is characterized by an extremely arid desert climate, with little precipitation and high potential evaporation [56]. Due to the extremely arid desert climate in the lower reach of these rivers, most of the rivers disappear in the Taklamakan Desert. The climate of the Tarim Basin is continental warm-temperate, with a mean annual precipitation of approximately 50 mm and potential evapotranspiration of up to 2200 mm annually [57]. Under the dominance of westerlies, the precipitation at high elevation is above 300 mm and almost always in the form of snow. Temperature ranges between −3 °C in the winter and 40 °C in summer [58].

Data
MODIS 8-day 500 m Surface Reflectance product (MOD09A1) was used in this study to extract glacier and seasonal snow cover. To map annual glacier area, six tiles of MODIS, h23v4, h23v5, h24v4, h24v5, h25v4, and h25v5, were employed to construct mosaic images. The MODIS product version V006 is based on improved snow, cloud and cloud shadow detection algorithms [59]. These MODIS products were imported in the Google Earth Engine Application Programming Interface (GEE API). Data acquired between July and September, when snow cover was minimum, were used to delineate glacier area. Temporal composite image was created by applying the GEE reducer function to obtain the median radiance over all valid samples from July to September for each pixel. This composite image was then used to calculate the NDSI and then to extract the annual glacier outline by applying the prescribed threshold. Table 1 describes the properties of the MODIS spectral bands [59].  The climate of the Tarim Basin is continental warm-temperate, with a mean annual precipitation of approximately 50 mm and potential evapotranspiration of up to 2200 mm annually [57]. Under the dominance of westerlies, the precipitation at high elevation is above 300 mm and almost always in the form of snow. Temperature ranges between −3 • C in the winter and 40 • C in summer [58].

Data and Method
MODIS 8-day 500 m Surface Reflectance product (MOD09A1) was used in this study to extract glacier and seasonal snow cover. To map annual glacier area, six tiles of MODIS, h23v4, h23v5, h24v4, h24v5, h25v4, and h25v5, were employed to construct mosaic images. The MODIS product version V006 is based on improved snow, cloud and cloud shadow detection algorithms [59]. These MODIS products were imported in the Google Earth Engine Application Programming Interface (GEE API). Data acquired between July and September, when snow cover was minimum, were used to delineate glacier area. Temporal composite image was created by applying the GEE reducer function to obtain the median radiance over all valid samples from July to September for each pixel. This composite image was then used to calculate the NDSI and then to extract the annual glacier outline by applying the prescribed threshold. Table 1 describes the properties of the MODIS spectral bands [59]. Surface Reflectance retrievals from L5-Thematic Mapper (TM) and L8-Operational Land Imager (OLI) data are provided by the United States Geological Survey (USGS) and are available in GEE (see Table 1). Here, we used the USGS L5-TM and L8-OLI Surface Reflectance Tier 1 dataset to extract glacier areas to evaluate our results of the glacier delineation from MODIS data in 2000 and 2020 accordingly. Cloudless L5-TM in 2000 and L8-OLI in 2020 were selected for validating MODIS glacier delineation. As the reflectance of snow and ice is high in the visible band and near zero in the short-wave infrared band, bands 2 and 5 of L5-TM and bands 3 and 6 of L8-OLI were used to delineate snow and glaciers. The method to delineate glaciers was similar to the one applied to MODIS data.

Boundary of the Tarim Basin
The boundary of the Tarim Basin was defined using the HydroSHEDS data to determine the boundaries of the sub-basins and of the entire Tarim Basin. HydroSHEDS data provide hydrographic information on global and local scales [60] available at https: //hydrosheds.org (accessed on 1 October 2021). The HydroSHEDS database provides a suite of geo-referenced datasets, including watershed boundaries, drainage directions, channel networks, flow accumulations, and river topology details [61]. The data are a good source to depict the catchment areas or the watershed boundaries on different scales.

Glacier Inventory
The Randolph Glacier Inventory (RGI 6.0) were applied as a reference to evaluate the glacier delineation based on MODIS data. The RGI 6.0 provided debris-covered information in our study area ( Table 2). It was estimated that 3.91% of glaciers are covered by debris in Tarim Basin [42]. The RGI 6.0 dataset was motivated by the Fifth Assessment Report of the Intergovernmental Panel on Climate Change (IPCC AR5) [62], which utilized most of the Landsat Thematic Mapper (TM) and Enhanced TM Plus (ETM+) images between 2003 and 2010 [63] to delineate glaciers in China.

SRTM DEM
The Shuttle Radar Topography Mission (SRTM) DEM is a 30 m resolution void-filled digital elevation model, produced by NASA based on the STRM Interferometric Synthetic Aperture Radar (InSAR) data. This dataset was created and made available by the USGS and the German Aerospace Center (DLR). The accuracy of SRTM DEM is regarded as better than 20 m in the horizontal and 16 m in the vertical direction expressed as linear error at a 90% confidence level [64]. The SRTM DEM dataset was generated using data in February 2000, while our glacier area was determined for each year between February 2000 and February 2020. In this study, we employed the SRTM DEM to acquire the elevation attributes of glaciers, as many earlier studies where the SRTM DEM was applied to support analyses on the distribution of glaciers with elevation [65,66].

Meteorological Data
In order to analyze regional climate forcing, the China Meteorological Forcing Dataset (CMFD) was downloaded from the National Tibetan Plateau Data Center (http://data.tpdc. ac.cn/zh-hans/data/8028b944-daaa-4511-8769-965612652c49, accessed on 3 January 2021). The CMFD is a high spatial-temporal resolution gridded near-surface meteorological dataset developed specifically for studies on land surface processes [67]. The dataset was generated through the fusion of remote sensing products, reanalysis data and in situ weather station data. The dataset includes 2 m air temperature and precipitation rate with a temporal resolution of three hours and a spatial resolution of 0.1 • [68]. The annual mean temperature and precipitation from 1979 to 2018 were calculated for each sub-basin (see Section 5.2) in our study.

Method
Google Earth Engine is suitable to develop and implement long-term and large area time series remote sensing applications; its high-performance computation capabilities offer a state-of-the-art platform to explore and process earth observation data without actually downloading such large datasets [52].
To explore the relationship of snow cover with glacier area in different regions of the Tarim Basin, we evaluated snow cover based on the MODIS 8-day composite images between 18 February 2000 and 25 January 2021. Summer is the best season to separate glaciers from snow cover through satellite images because snow not contributing to glacier mass balance melts in summer due to very few snow events occurring in summer in High Mountain Asia [49]. However, glacier mapping is challenging because clouds frequently appear during summer in the Tarim Basin. To mitigate this problem, we developed a summer composite image to delineate glaciers, as done in [69]. We use summer observations to separate the glacier and snow cover area based on the minimum snow cover during summer, i.e., from July to September, to delineate glaciers [49].
The method primarily includes two parallel workflows: (1) glacier extraction and (2) snow cover extraction. The glacier extraction workflow includes: selection and retrieval of datasets, generating a summertime composite image, annual glacier extraction, water mask generation to exclude water pixels from glacier, and glacier area delineation. Most of the water pixels were located close to the glacier terminal or to the glacier boundary and can be excluded using NDSI, but there were still a few misidentified pixels. To ensure the accuracy of estimated glacier changes, we masked off the water pixels. The snow cover extraction workflow includes: retrieval of 8-day MODIS composite datasets, identifying pixels with end of summer snow cover and delineation of intermittent snow cover by thresholding the median NDSI (see Section 3.2.1 for details), water mask generation excluding water pixels from snow cover, and snow cover area calculation. The glacier area thus delineated includes end of summer snow cover, which contributes to ice accumulation.

Extracting Glacier Area
The boundaries of the Tarim Basin and its sub-basins were defined based on Hy-droSHEDS data (see Section 3.1.3). We uploaded the Boundary shapefile to the GEE driver and defined it as Region of Interest. The Date Range was used to select datasets within a specific time period.
The Datasets applied to annual delineation of glaciers were MODIS (called 'MODIS/006/MOD09A1′ in GEE). The coarse spatial resolution of MODIS can only capture glaciers larger than 0.25 km 2 and it is likely to overestimate the rate of decrease in small glaciers. Glaciers were extracted using summer composites of MODIS surface reflectance data. For example, to extract the glacier area in 2000, using the GEE Image Collection function, we selected, accessed and processed the MOD09A1 images from 1 July 2000 to 30 September 2000, i.e., the summer period to minimize the impact of snow cover. We called this step "median surface in-band reflectance", which was applied to generate a summer composite image for each year. All reflectance retrievals for each pixel contributed to determine the median reflectance during summer. This compositing procedure image yielded combinations of in-band reflectance such that NDSI was lower when observing rapidly varying surface features, such as clouds and seasonal snow cover, and higher for stable snow and glaciers. Next, we calculated the NDSI based on the composite image.
NDSI can be applied to discriminate glacier/snow from many other land cover features and decreases the influence of atmospheric effects and viewing angles. This is achieved by normalizing the spectral difference between MODIS band 4 and band 6 to their sum [70]:

Extracting Glacier Area
The boundaries of the Tarim Basin and its sub-basins were defined based on Hy-droSHEDS data (see Section 3.1.3). We uploaded the Boundary shapefile to the GEE driver and defined it as Region of Interest. The Date Range was used to select datasets within a specific time period.
The Datasets applied to annual delineation of glaciers were MODIS (called 'MODIS/ 006/MOD09A1' in GEE). The coarse spatial resolution of MODIS can only capture glaciers larger than 0.25 km 2 and it is likely to overestimate the rate of decrease in small glaciers. Glaciers were extracted using summer composites of MODIS surface reflectance data. For example, to extract the glacier area in 2000, using the GEE Image Collection function, we selected, accessed and processed the MOD09A1 images from 1 July 2000 to 30 September 2000, i.e., the summer period to minimize the impact of snow cover. We called this step "median surface in-band reflectance", which was applied to generate a summer composite image for each year. All reflectance retrievals for each pixel contributed to determine the median reflectance during summer. This compositing procedure image yielded combinations of in-band reflectance such that NDSI was lower when observing rapidly varying surface features, such as clouds and seasonal snow cover, and higher for stable snow and glaciers. Next, we calculated the NDSI based on the composite image.
NDSI can be applied to discriminate glacier/snow from many other land cover features and decreases the influence of atmospheric effects and viewing angles. This is Remote Sens. 2021, 13, 5117 8 of 24 achieved by normalizing the spectral difference between MODIS band 4 and band 6 to their sum [70]: Glaciers were extracted by applying a threshold to the NDSI calculated with median summertime values of the at-surface reflectance in MODIS bands 4 and 6 (see Equation (1)). Pixels with NDSI ≥ 0.4 are identified as being glacier or stable snow [42,44]. A threshold of 0.11 on an NIR band is widely used to prevent water pixels from being classified as glaciers [37,48,71]. Furthermore, the Primary Glacier Boundary ( Figure 2) is generated by applying these threshold settings. This operation is repeated for each year.
The snow products are affected, however, by multiple limitations, such as the misclassification of water pixels [44]. For this reason, we removed water pixels when re-generating glacier and snow cover areas. The water pixels were identified by applying a threshold to the Normalized Difference Water Index (NDWI). In addition, since melting occurs at the glacier surface during summer, liquid water may also exist on glaciers and will give reflectance spectra similar to liquid water. In order to further diminish the impact of water pixels, we generated a water mask based on the summertime composite image. The water pixels were identified by applying the Normalized Difference Water Index [72], due to the strong water absorption in the near-infrared band and relatively higher reflectance in the green band. The NDWI [72] is defined as: Pixels are identified as water when NDWI > 0.2 and the near-infrared reflectance in band 2 is smaller than 0.2 [37,48]. This step is to generate a water mask according to the summer composite images and further remove the water pixels. The annual glacier area was calculated after removing the water mask. The annual glacier boundary was generated from 2000 to 2020 using the composite images during summer.
As we mentioned in Section 3.1.2, the L5-TM and L8-OLI Surface Reflectance product (namely, 'LANDSAT/LT05/C01/T1_SR' and 'LANDSAT/LC08/C01/T1_SR' in GEE) data were used to evaluate the delineation of glaciers based on the lower resolution MODIS data. Due to the limited number of cloudless L5-TM scenes, we used L5-TM scenes from 1 January 2000 to 30 December in 2000 to ensure the coverage of study area. Cloudless L8-OLI scenes from 1 July 2020 to 30 September in 2020 were mosaiced in GEE.
Next, the glacier area is obtained after removing the water mask in 2000 and 2020, respectively. We compared the glacier area in 2000/2020 generated by MODIS with the glacier area in 2000/2020 generated by L5-TM/L8-OLI data to evaluate our glacier delineation based on MODIS data ( Figure 3).

Extracting Snow Cover
In this study, we separated the glacier and intermittent snow cover mainly based on the stability of summer snow cover using frequent observations. We considered summer as the season with the smallest snow area, and thus considered the area delineated by thresholding the median NDSI during summer as the glacier and stable snow cover area. Intermittent snow cover was calculated using data acquired at 8-day intervals during winter, i.e., from 18 February 2000 to 25 January 2021, and using the MODIS data accessed via GEE. Snow cover was delineated by applying the same threshold, as done with glaciers, to the 8-day NDSI images calculated with each 8-day MODIS in-band at-surface reflectance in bands 4 and 6. The snow cover generation method is similar to the process of glacier generation. The primary snow cover was generated when NDSI ≥ 0.4 and band 2 > 0.11, then the intermittent snow cover was the total number of pixels identified as snow in all available images during winter. The water mask was also generated when NDWI > 0.2 and the near-infrared reflectance in band 2 < 0.2. Next, the water mask was excluded from the first estimate of intermittent snow cover. The calculated intermittent snow cover area was output in the GEE table format. Finally, the glacier area, determined with the summertime data, was excluded from the first estimate of the intermittent snow cover area within the same year.

Accuracy Assessment
Glaciated areas may be misclassified as intermittent snow due to similar spectral characteristics [73]. In our proposed method, intermittent snow is identified by evaluating the temporal variability of snow cover for each pixel and relating it to the median surface spectral reflectance. Estimating the glacier area accuracy is challenging, as reference data are limited. Here, glacier area delineated with the MODIS images by applying the NDSI threshold (see Section 3) was evaluated against the area of five subregions with the boundaries of glaciers in 2000 delineated with L5-TM and as with L8-OLI image data in 2020 ( Table 2). The delineation with MODIS data of glaciers in different subregions was also compared with previously published results (see Section 5.1).
Glacier delineation by a spectral band combination is supposed to be one of the most efficient methods applicable to debris-free glaciers, but it is not appropriate to debriscovered glaciers [74]. As mentioned in the Introduction, the debris-covered glacier portions were not taken into account in our study due to the limited incidence of changes in the extension of debris facies on the total change in glacier area of the study region [42,43,75,76]. According to the published RGI 6.0, debris-covered glaciers take up 3.91% of the total glacier area in the Tarim Basin, i.e., errors of estimate have limited effect on the glacier area result.
The uncertainty in the MODIS 2000 glacier area was higher in East Tien Shan and East Kunlun, when compared with the estimates based on the L5-TM image data ( Table 2). The results for the Pamir and West Kunlun regions were the most accurate, with −6.76% and −7.94% relative differences when comparing the MODIS 2000 estimates with the L5-TM 2000 data. The comparison of L8-OLI delineated glacier area and MODIS glacier area in 2020 indicated a higher uncertainty in East Tien Shan and Pamir. The results for the West Tien Shan and East Kunlun regions were the most accurate, with 2.55% and 6.19%.
To estimate the rate of change of glacier area, we fitted a linear relationship to the annual values in Figures 4 and 5, then we took the slope of this relationship as estimated rate of change. To evaluate the accuracy of these estimates, we calculated the relative Root Mean Square Deviation as: where x i is the actual value, x i is the linear regression value, and N is the number of non-missing data points. The Relative Root Mean Square Deviation was calculated is then: The RRMSD has been added in Figures 4 and 5. The error of estimate on glacier area in a region was −6.48% (weighted average of all MODSI cases in Table 2).

Variability of Glacier Area in the Tarim Basin
The annual glacier area in the five subregions of the Tarim Basin was retrieved from MODIS data from 2000 to 2020 by applying the method described in Section 3 ( Figure 4). In the whole Tarim Basin, the glacier area decreased by 7975.71 km 2 from 2000 to 2020, i.e., 25.11% of the initial glacier area, at a rate of 0.94%/a, and it decreased fastest during 2016-2020 at a rate of approximately 1.82%/a. In other words, the decreasing rate accelerated after 2016.

Variability of Glacier Area in the Tarim Basin
The annual glacier area in the five subregions of the Tarim Basin was retrieved from MODIS data from 2000 to 2020 by applying the method described in Section 3 ( Figure 4). In the whole Tarim Basin, the glacier area decreased by 7975.71 km 2 from 2000 to 2020, i.e., 25.11% of the initial glacier area, at a rate of 0.94%/a, and it decreased fastest during 2016-2020 at a rate of approximately 1.82%/a. In other words, the decreasing rate accelerated after 2016.
The evolution of glacier area in the five subregions of the Tarim Basin was further analyzed separately ( Figure 5). It was found that glacier area decreased in all the five subregions during the 2000-2020 period ( Figure 5), but there were clear differences in the evolution across the five subregions.
The rate of change in clean glacier area in E Tien Shan, W Tien Shan, Pamir, W Kunlun, and E Kunlun was 2.98%/a, 1.07%/a, 0.5%/a, 1.04%/a, and 0.81%/a, respectively, during 2000-2020. These rates imply that the total reduction in 2000-2020 was largest, i.e., 50.86%, i.e., 154.13 km 2 , in the E Tien Shan area. The area of glaciers in W Tien Shan decreased by 973.76 km 2 (13.02%). In the Pamir region, the area of glacier decreased by 20.52%. The reduction in the W Kunlun and E Kunlun regions were 30.98% and 36.49%, respectively.

Changes in Glacier Area at Different Elevations
Most glaciers in the Tarim Basin are located at elevations in the range 3000-7500 m (Figures 6 and 7). Many small glaciers are scattered at these elevations in the Tien Shan region. In the Pamir and West Kunlun, glacier elevation ranges from 5000 to 7500 m. West Kunlun is the largest glacial region in the Tarim Basin. In East Kunlun, glaciers are distributed above 6000 m. The terminus position of most glaciers is above 3000 m, while the snow line is between 3600 and 5700 m [17]. In general terms, temperature decreases and precipitation increases with elevation, but spatial variability is significant and glacier response to climate variability needs to be evaluated at different elevations.
Glacier area in the 5500-6000 range is the largest fraction of total glacier area in the Tarim Basin, and the area was 9150.77 km 2 in 2000, and decreased to 7434.62 km 2 in 2020. Our analysis shows that the decrease in glacier area was largest in the 3500-5500 m range: the glacier area was 18,358.29 km 2 in 2000, and decreased to 12,568.57 km 2 in 2020. Glacier area in the 6000-6500 range was 3926.68 km 2 in 2000, and decreased slightly to 3511.61 km 2 in 2020. However, the area of glaciers above 6500 remained stable between 2000 and 2020, and even increased slightly, although the glacier area above 6500 m accounts for 1.08% of total glacier area. The stability is most likely due to the low temperature at higher elevations, while air temperature increased at lower elevations.
Our analysis shows that the decrease in glacier area was largest in the 3500-5500 m range: the glacier area was 18,358.29 km 2 in 2000, and decreased to 12,568.57 km 2 in 2020. Glacier area in the 6000-6500 range was 3926.68 km 2 in 2000, and decreased slightly to 3511.61 km 2 in 2020. However, the area of glaciers above 6500 remained stable between 2000 and 2020, and even increased slightly, although the glacier area above 6500 m accounts for 1.08% of total glacier area. The stability is most likely due to the low temperature at higher elevations, while air temperature increased at lower elevations.

Snow Cover Area in Sub-Basins
Glacier area is determined by mass accumulation and ablation, with the dominant contributor to ice accumulation being snowfall [77]. Snowfall variability can be observed through intermittent snow cover. Snow cover can be detected by optical remote sensing because of its distinct spectral signature [48], but frequent cloud cover and mixed water pixels in high mountain regions need to be dealt with. As we mentioned in Section 3.2, we applied a water mask to minimize the water effect and used 8-day composite surface reflectance products to minimize the impact of cloud cover in the retrieval of intermittent snow cover. The temporal evolution of snow cover area showed a typical monthly cycle with lowest area observed in July and August (Figure 8). The fastest expansion of snow cover happened in October and the fastest reduction occurred in April (Figure 8). Overall, no significant trend could be observed in the time series at full temporal resolution in Tarim Basin from 2000 to 2020.
The seasonality of snow cover in each subregion was better captured by the mean monthly snow cover area (Figure 9). The snow cover area started increasing in August and grew until February in the following year. At a lower temperature, the water vapor carried by the westerlies precipitated in solid form in the high mountain regions, thus causing an increase in snow cover area. The snow cover area decreased from February to August, while snow melting already started at lower elevations in February. Snow cover area in the E Kunlun, however, already started decreasing in January, which might be related to the strengthening westerlies [31]. In the W Tien Shan region, the annual evolution of snow cover was slightly different, with the largest snow cover area being reached

Snow Cover Area in Sub-Basins
Glacier area is determined by mass accumulation and ablation, with the dominant contributor to ice accumulation being snowfall [77]. Snowfall variability can be observed through intermittent snow cover. Snow cover can be detected by optical remote sensing because of its distinct spectral signature [48], but frequent cloud cover and mixed water pixels in high mountain regions need to be dealt with. As we mentioned in Section 3.2, we applied a water mask to minimize the water effect and used 8-day composite surface reflectance products to minimize the impact of cloud cover in the retrieval of intermittent snow cover. The temporal evolution of snow cover area showed a typical monthly cycle with lowest area observed in July and August (Figure 8). The fastest expansion of snow cover happened in October and the fastest reduction occurred in April (Figure 8). Overall, no significant trend could be observed in the time series at full temporal resolution in Tarim Basin from 2000 to 2020.
Snow Index in W Kunlun was the largest among the five regions, which means the W Kunlun accumulated more snow. RSI does not indicate the area of remaining snow; it reflects the snow cover change within one year. For a given minimum snow cover area, the RSI is larger when the snow area is smaller in winter. Accordingly, an increasing trend of RSI does not imply that the area of remaining snow increases, only that the fraction of minimum snow area in summer over maximum snow area in winter is increasing. Further detailed discussion can be found in Section 5.2.  The seasonality of snow cover in each subregion was better captured by the mean monthly snow cover area (Figure 9). The snow cover area started increasing in August and grew until February in the following year. At a lower temperature, the water vapor carried by the westerlies precipitated in solid form in the high mountain regions, thus causing an increase in snow cover area. The snow cover area decreased from February to August, while snow melting already started at lower elevations in February. Snow cover area in the E Kunlun, however, already started decreasing in January, which might be related to the strengthening westerlies [31]. In the W Tien Shan region, the annual evolution of snow cover was slightly different, with the largest snow cover area being reached in February. Also worth noting was the lower rate of decrease, especially from April to July, in the W Kunlun when compared with other subregions. In addition, the snow cover rate of increase from August to December was higher than in other subregions.

Comparison with Previous Studies
Other authors have studied changes in glacier area in the Tarim Basin, and we evaluated our findings against this literature evidence. To date, regional investigations on changes in glacier area were limited in spatial coverage, so it is impossible to verify our results in detail. The comparisons are performed using reported results on the annual rate of change in glacier area within slightly different regions and period of time. In these studies, the glacier area shows a diversified variability in the W Tien Shan, E Tien Shan, Pamir, W Kunlun, and E Kunlun (Table 3).
A decrease in glacier area was recorded in the W Tien Shan during 1970-2018. The rate of shrinkage in glacier area varied from 0.11%/a to 3.7%/a within different sub-basins, with the Sary-Jaz River Basin having the largest reduction rate during 1990-2010 [78]. The reduction rate of the Sary-Jaz River Basin was slightly larger than our study result of 1.07%/a. Our study area included the Central Tien Shan, while other studies paid more attention to the China part of the Tien Shan [6,19,79]. Some disagreement in the rate of decrease estimated by different studies may be due to differences in the area of interest [12]. In addition, the large number of small glaciers (<1 km 2 ) [19] may be another reason for the discrepancy in estimated annual rate. Differences in the observation period may clearly lead to significant differences, especially between earlier and recent studies, since the glacier rate of decrease may have been accelerating due to increasing temperature and reduced precipitation in recent years [49]. In addition, most studies were focused on one To detect and evaluate a possible multi-annual trend in snow cover, we defined the Remaining Snow Index (RSI) as: where 'A max ' is the yearly maximum value of the monthly snow cover area, with the maximum value usually being reached in January or February (Figure 9). The 'A min ' is the yearly non-zero minimum value of the monthly snow cover area, usually being reached in July or August (Figure 9). The Remaining Snow Index is an indicator of the snow remaining annually at the end of the yearly cycle; thus, it is related to mass accumulation in glaciers, but it is a relative, not an absolute, indicator of snow area. Here, we used the Remaining Snow Index to analyze the snow variation in 2000-2020 in each subregion (Figure 10), which could not be captured by the data plotted in Figure 8. The Remaining Snow Index in W Kunlun was the largest among the five regions, which means the W Kunlun accumulated more snow. RSI does not indicate the area of remaining snow; it reflects the snow cover change within one year. For a given minimum snow cover area, the RSI is larger when the snow area is smaller in winter. Accordingly, an increasing trend of RSI does not imply that the area of remaining snow increases, only that the fraction of minimum snow area in summer over maximum snow area in winter is increasing. Further detailed discussion can be found in Section 5.2.

Comparison with Previous Studies
Other authors have studied changes in glacier area in the Tarim Basin, and we evaluated our findings against this literature evidence. To date, regional investigations on changes in glacier area were limited in spatial coverage, so it is impossible to verify our results in detail. The comparisons are performed using reported results on the annual rate of change in glacier area within slightly different regions and period of time. In these studies, the glacier area shows a diversified variability in the W Tien Shan, E Tien Shan, Pamir,

Comparison with Previous Studies
Other authors have studied changes in glacier area in the Tarim Basin, and we evaluated our findings against this literature evidence. To date, regional investigations on changes in glacier area were limited in spatial coverage, so it is impossible to verify our results in detail. The comparisons are performed using reported results on the annual rate of change in glacier area within slightly different regions and period of time. In these studies, the glacier area shows a diversified variability in the W Tien Shan, E Tien Shan, Pamir, W Kunlun, and E Kunlun (Table 3).
A decrease in glacier area was recorded in the W Tien Shan during 1970-2018. The rate of shrinkage in glacier area varied from 0.11%/a to 3.7%/a within different sub-basins, with the Sary-Jaz River Basin having the largest reduction rate during 1990-2010 [78]. The reduction rate of the Sary-Jaz River Basin was slightly larger than our study result of 1.07%/a. Our study area included the Central Tien Shan, while other studies paid more attention to the China part of the Tien Shan [6,19,79]. Some disagreement in the rate of decrease estimated by different studies may be due to differences in the area of interest [12]. In addition, the large number of small glaciers (<1 km 2 ) [19] may be another reason for the discrepancy in estimated annual rate. Differences in the observation period may clearly lead to significant differences, especially between earlier and recent studies, since the glacier rate of decrease may have been accelerating due to increasing temperature and reduced precipitation in recent years [49]. In addition, most studies were focused on one or two river catchments of the Tarim Basin, while our study covers the entire Tarim Basin. The glacier area shrinkage rate was 0.8%/a between 1990 and 2018 for the E Tien Shan in the study of Huang et al. (2021) [49], while the glaciers shrank 2.98%/a for the same area during 2000-2020 according to our results. In the Pamir and East Kunlun region, the rate of glacier decrease was smaller than in the East Tien Shan region. We only consider the part of the Pamir region draining into Tarim Basin. The sub-basin within Pamir shrank by 0.05% to 0.59%/a during 1990-2018 [49,[80][81][82], while our result of 1.01%/a is larger than previous studies. In Pamir and W Kunlun, the glacier snouts have been retreating, stable or advancing [49,80]. The smallest rate of glacier contraction was observed in the Pamir regions [31]. The deviation is possibly caused by the different observation period, i.e., the rate of decrease was 0.46 during 2000-2010, while it was 0.59 during 1990-2018 according to Huang et al. (2021). More comparisons in the five subregions and periods are shown in Table 3. Table 3. Comparison of glacier area from various studies surrounding the Tarim Basin. The subregions correspond to Figure 6.

Glacier Area Variation in Response to Climate Forcing and Snow Cover
Since the end of the 1980s, a transition from warm-dry to warm-wet conditions with both increasing temperatures and precipitation was observed in northwestern China [87]. Many factors have caused a reduction in glacier area, but precipitation and temperature were the most significant, where the summer temperature mainly regulated glacier melting and winter precipitation mainly regulated the accumulation of ice. Here, annual precipitation and temperature between 1979 and 2018 in each sub-basin are evaluated.
In E Tien Shan, both annual precipitation and temperature increased during 1979-2004, while they decreased during 2004-2018 (Figure 11a).   [19] showed that the temperature and precipitation in the Aksu River Basin (W Tien Shan) increased in 1975-2016, and the warming temperatures were the main reason for the glacier recession in Tien Shan [19]. A positive trend in annual precipitation and temperature was noted in W Tien Shan (Figure 11b). Increases in temperature and precipitation have led to an elevationdependent snow cover response, with decreased snow cover area at lower elevations and increased snow cover area at higher elevations [88]. The increase in precipitation has been induced by increasing water vapor transport by the westerly airflow [89]. Westerly-affected regions such as the Pamir and W Kunlun have extensive winter snow cover [90,91]. The smallest rate of glacier contraction was observed in the Pamir, with a reduction rate of 0.50%/a [31]. Positive mass balances observed from geodetic studies (using DEMs and ICESat data) in W Kunlun in 2006-2010 were reported by [12].
According to [55], the mean annual precipitation in the southern part of the Tarim Basin is lower than in the northern and western regions. Higher temperature and higher water vapor content in the atmosphere may lead to larger mass accumulation than in earlier years. In the E Kunlun, both temperature and precipitation increased from 1979 to 2004, and decreased from 2004 to 2018. We found that the glacier area decreased from 2000 to 2004 (Figure 5e), which is in consistent with the noted trends in climatic condition.
No trend was detected in the annual snow-covered area, with significant interannual variability highlighted by calculating the RSI (see Figure 10) and generally a small fraction of winter snow-covered area remaining after summer.
We also analyzed the normalized cross-correlation functions of glacier area and temperature, glacier area and precipitation, and glacier area and snow cover in the five Tarim subregions. The detailed results are presented in the Supplementary Materials ( Figure S1). We calculated the normalized cross-correlation function by shifting step by step the response time series with respect to the forcing time series ten years forwards and ten years backwards, using the matlab function (xcorr(x,y, 'coeff')). The results show that the cross-correlation is symmetric around a zero time-lag where the correlation is highest. This indicates that the response of glacier area to temperature, precipitation and snow area is within less than a year, which is the temporal resolution of the data applied in this analysis. Precipitation and snow cover show positive correlation with glacier area, as expected, while the correlation of glacier area with temperature is positive in W Tien Shan and negative in other subregions.  [67].
The mean elevation of the glacier area in the Tien Shan region is lower than in the W Kunlun region, which increases the vulnerability to climate variability of the glaciers in the former region. Fan et al. (2014) [92] concluded that the spatial variability of temperature and precipitation in the high elevation mountain zone is significant, while the lower elevation zone (3500 m a.s.l.) is uniformly affected by less precipitation, higher temperature and large irradiance [89]. Most glaciers in the Tien Shan region are distributed at elevation between 3000 and 5500 m and experienced a continuous glacier recession. In our study, the absolute loss is highest in W Kunlun, and the highest relative rate of decrease is found in the E Tien Shan. Shangguan et al. (2009) [17] found that the altitude of the maximum glacier reduction was 4900-5400 m a.s.l., while another peak in the decrease of glacier area was found at an altitude of 4100-4500 m a.s.l. Our results also show that the glacier area varies greatly between 4000 and 5500 m a.s.l. Shangguan et al. (2009) [17] also found little change in the area of glaciers at an altitude of 6100 m, which is not consistent The mean elevation of the glacier area in the Tien Shan region is lower than in the W Kunlun region, which increases the vulnerability to climate variability of the glaciers in the former region. Fan et al. (2014) [92] concluded that the spatial variability of temperature and precipitation in the high elevation mountain zone is significant, while the lower elevation zone (3500 m a.s.l.) is uniformly affected by less precipitation, higher temperature and large irradiance [89]. Most glaciers in the Tien Shan region are distributed at elevation between 3000 and 5500 m and experienced a continuous glacier recession. In our study, the absolute loss is highest in W Kunlun, and the highest relative rate of decrease is found in the E Tien Shan. Shangguan et al. (2009) [17] found that the altitude of the maximum glacier reduction was 4900-5400 m a.s.l., while another peak in the decrease of glacier area was found at an altitude of 4100-4500 m a.s.l. Our results also show that the glacier area varies greatly between 4000 and 5500 m a.s.l. Shangguan et al. (2009) [17] also found little change in the area of glaciers at an altitude of 6100 m, which is not consistent with our result.
Some studies used hydrologic modelling to describe the response of glaciers to temperature and precipitation [93][94][95]. However, the temperature and precipitation (liquid water) data used for hydrologic modelling mainly depend on weather stations scattered in the high elevation mountain area [93,94]. Additionally, elevation, monsoon, and terrain vary significantly in the Tarim region [96]. The widely noted changes in glacier area raised the awareness of the urgency of studying the mountain glaciers in the Tarim Basin [97,98]. The estimated response of glacier area by hydrologic modelling, however, is affected by the coarse spatial resolution and by uncertainty in the parameterization of terrain and elevation effects [96].

Limitations and Advantages of the Method
The MODIS data products take advantage of the temporal frequency. In this study, we used this data product to evaluate the annual change rate in glacier area and in the seasonal snow cover area. This is a systematic way to explore the relationship between glacier and snow cover areas. The MODIS coarse spatial resolution, however, limits the accuracy in the estimation of glacier area. The estimated total glacier area is affected by the minimal glacier size in relation to the spatial resolution of applied image data. The minimal size of detectable glaciers with MODIS data is 0.25 km 2 (500 m × 500 m). Our estimates of the area of small glaciers and of their distribution may differ from results based on high-or medium-resolution data. In addition, the examples of Muztag and Kungur Tagh were used to show the limitations and advantages of the proposed NDSI based extraction method. We used the standard NDSI threshold of 0.4 to extract glacier area and intermittent snow cover in our study. In Figure 3, the extracted glacier outline does not correctly capture the debris-covered region, which may cause an underestimation of glacier area [74]. In W and E Kunlun, the fractional debris-covered area is less than 5% of the glacier area, so the effect of debris is limited [42]. The limitation of our study is that we did not map the debris-covered region. A glacier is regarded as "perennial mass of ice, and possibly firn and snow, originating on the land surface by the recrystallization of snow or other forms of solid precipitation and showing evidence of past or present flow" [52]. "Firn is snow that has survived at least one ablation season but has not been transformed to glacier ice" [52]. However, glacier mapping is challenging because clouds frequently appear during summer in the Tarim Basin. To mitigate this problem, we developed a summer composite image to delineate glaciers, as done in [53]. We choose the data acquired in July and September to delineate the glacier area, but in high altitude regions, such as W Kunlun (more than 5500 m a.s.l.), the perennial snow ("snow that has survived at least one ablation season but has not been transformed yet to glacier ice" [99]) can remain for years at high elevation. We have evaluated the impact of perennial snow on our results as follows.
The perennial snow appears to be adjacent to the glaciers, as delineated in the RGI 6.0 ( Figure 12). The perennial snow was mapped annually in the entire Tarim Basin. The perennial snow that occurred in more than 5 years out of 20 was analyzed to determine its frequency distribution (Figure 13). The snow persistent for the entire observation period, i.e., 2000-2020, was already counted in our analysis of the evolution of glacier area and has no impact on our results. The snow with shorter persistence was not present continuously and cannot be counted as glacier area. Remote Sens. 2021, 13, x FOR PEER REVIEW 20 of 25

Conclusions
Glacier area extraction from satellite images is an important aspect of glacier mapping and studies on glacier changes. The selection of cloud-free images is difficult when using low revisit images such as the ones acquired by L8-OLI. Our method relies on the frequency of usable image pixels to minimize the impact of clouds and to detect intermittent snow. The major innovation in our method is the use of the median surface reflectance to composite several images pixel by pixel. This does not require that cloud-free observations of all pixels are available on the same day, just that a sufficient number of cloud-free samples is acquired for each pixel. The composite image improves the reliability of estimated change in glacier area, because the median is a robust statistic. The NDSI method combined with threshold setting was effective, leading to systematic estimates of the glacier area and monthly snow cover.
Changes in glacier and intermittent snow area between 2000 and 2020 were first evaluated for the entire Tarim Basin. The results show that the glacier area decreased by 7975.71 km 2 (25.11%), with the annual rate of decrease being 0.94%/a. Most glaciers are distributed between 5000 and 6000 m a.s.l., where the glacier area fluctuates most. Above 6000 m a.s.l., the glacier area was rather stable. The consistent spatial distribution of perennial snow suggests that the RGI 6.0 boundaries may need to be revisited.
The Tarim Basin glaciers show considerable heterogeneity in the five different subregions due to the different circulation of the atmosphere and environmental conditions. Overall, the five subregions show a decreasing trend during the 2000-2020 period. In the Tien Shan region, the glacier area decreased rapidly. The intra-annual variability of snow

Conclusions
Glacier area extraction from satellite images is an important aspect of glacier mapping and studies on glacier changes. The selection of cloud-free images is difficult when using low revisit images such as the ones acquired by L8-OLI. Our method relies on the frequency of usable image pixels to minimize the impact of clouds and to detect intermittent snow. The major innovation in our method is the use of the median surface reflectance to composite several images pixel by pixel. This does not require that cloud-free observations of all pixels are available on the same day, just that a sufficient number of cloud-free samples is acquired for each pixel. The composite image improves the reliability of estimated change in glacier area, because the median is a robust statistic. The NDSI method combined with threshold setting was effective, leading to systematic estimates of the glacier area and monthly snow cover.
Changes in glacier and intermittent snow area between 2000 and 2020 were first evaluated for the entire Tarim Basin. The results show that the glacier area decreased by 7975.71 km 2 (25.11%), with the annual rate of decrease being 0.94%/a. Most glaciers are distributed between 5000 and 6000 m a.s.l., where the glacier area fluctuates most. Above 6000 m a.s.l., the glacier area was rather stable. The consistent spatial distribution of perennial snow suggests that the RGI 6.0 boundaries may need to be revisited.
The Tarim Basin glaciers show considerable heterogeneity in the five different subregions due to the different circulation of the atmosphere and environmental conditions. Overall, the five subregions show a decreasing trend during the 2000-2020 period. In the Tien Shan region, the glacier area decreased rapidly. The intra-annual variability of snow

Conclusions
Glacier area extraction from satellite images is an important aspect of glacier mapping and studies on glacier changes. The selection of cloud-free images is difficult when using low revisit images such as the ones acquired by L8-OLI. Our method relies on the frequency of usable image pixels to minimize the impact of clouds and to detect intermittent snow. The major innovation in our method is the use of the median surface reflectance to composite several images pixel by pixel. This does not require that cloud-free observations of all pixels are available on the same day, just that a sufficient number of cloud-free samples is acquired for each pixel. The composite image improves the reliability of estimated change in glacier area, because the median is a robust statistic. The NDSI method combined with threshold setting was effective, leading to systematic estimates of the glacier area and monthly snow cover.
Changes in glacier and intermittent snow area between 2000 and 2020 were first evaluated for the entire Tarim Basin. The results show that the glacier area decreased by 7975.71 km 2 (25.11%), with the annual rate of decrease being 0.94%/a. Most glaciers are distributed between 5000 and 6000 m a.s.l., where the glacier area fluctuates most. Above 6000 m a.s.l., the glacier area was rather stable. The consistent spatial distribution of perennial snow suggests that the RGI 6.0 boundaries may need to be revisited.
The Tarim Basin glaciers show considerable heterogeneity in the five different subregions due to the different circulation of the atmosphere and environmental conditions. Overall, the five subregions show a decreasing trend during the 2000-2020 period. In the Tien Shan region, the glacier area decreased rapidly. The intra-annual variability of snow cover was largest for W Kunlun and smallest in the E Tian Shan. The analysis of the cross-