1. Introduction
Snow cover is an important component of the cryosphere, as well as an indispensable variable in the study of earth system science and climate change. Snow is highly dynamic and has a great influence on land surface energy budgets and atmospheric circulation patterns due to high albedo and good thermal insulation [1,2,3,4,5]. Under the background of global warming, changes in snow cover will lead to ecological and environmental problems such as a reduction of water availability, extreme weather events, and frequent disasters, profoundly affecting the ecosystem and sustainable socio-economic development of the countries and regions concerned [6,7], which has received widespread attention from these countries. The snowline is identified as the boundary separating snow-covered areas from snow-free areas [8,9,10]. Due to the patchiness of the snow cover edge, it is also defined as a narrow belt that represents a zone of ∼50% snow coverage [9,11,12]. As an important indicator of climate change, snowline can reflect the increase and reduction of snow, and it can also comprehensively reflect the basic state of the climate and environment of mountains, plateaus, and polar regions that lack meteorological stations [13].
The snowline altitude (SLA) is a valuable metric that integrates the competing effects of snow accumulation and melt, and it can be used to assess future changes in snow cover [1,14]. Regional SLA and its inter- and intra-annual variability are key characteristics that indicate temporal variation in snow cover and duration of snow melt, and help in assessing the hydrologic cycle balance [15]. For the estimation of snow-covered area and its temporal evolution, SLA can be used as an input for hydrological modeling or validation of snow model simulations [16,17]. Besides, SLA estimates can be applied to remove clouds from satellite snow cover products [18,19]. The SLA at the end of the melting season can also serve as a good proxy for the equilibrium line altitude (ELA) on glaciers, where much of the remaining end of summer snow cover is located, and therefore for glacier mass balance [20,21,22,23,24,25].
Central Asia is the core area of the Silk Road economic belt. As the places most distanced from the oceans in the Eurasian continent, the Tienshan Mountains (TS) are called the water tower of Central Asia. They consist of a series of high mountains, basins, and valleys, boasting one of the most developed glacier mountains in the world [26,27]. Major rivers recharged by glacier/snow melt water originating from the TS (e.g., the Ili River, Syr Darya River, Amu Darya River, Tarim River, and Chu River) feed the lowland areas of Kazakhstan, Kyrgyzstan, Uzbekistan, and China’s north-western Xinjiang Uyghur Autonomous Region, which together form one of the largest irrigated areas in the world [26,28,29]. Spatiotemporal variations of the snow cover in the TS exert a significant impact on the changes of water source for its surrounding arid regions and the hydrological and biological processes [30,31,32]. Therefore, the detection and exhaustive analysis of spatiotemporal variation of the SLA over the entire TS are quite essential for the protection and utilization of local water resources.
Satellite remote sensing offers the opportunity to monitor seasonal snow dynamics in the inaccessible areas with rugged terrain and hostile climate. With the advantages of high spectral resolution, high temporal resolution, wide spatial coverage, and being freely available, the remote sensing snow cover products generated from the Moderate Resolution Imaging Spectroradiometer (MODIS) provide an excellent opportunity to study the snow cover for global or large-scale areas. Hence, it is also suitable to assess snowline changes in a continuous time series and space for quite a large-scale area, based on the MODIS snow cover products. The MODIS snow cover products [33] have been widely used to depict the spatiotemporal patterns of the seasonal or transient SLA in mountainous areas [15,19,34,35,36]. Besides, some evaluation studies have suggested a high accuracy of MODIS snow cover products under clear skies when compared with in situ observations and other higher resolution satellite data at both regional and global scales [37,38,39].
Seasonal snow cover in high latitude and altitude regions in the northern and southern hemispheres is very sensitive to climate change [40,41]. Previous studies indicated that nearly all regions have experienced a warming trend, and that the average global temperature has increased 0.85 °C (from 0.65 to 1.06 °C) over the period from 1880 to 2012 [42]. As one of the most sensitive and prominent areas responding to global climate changes, the TS have especially experienced this obvious warming trend over the past few decades [43,44]. The temperature experienced a sharp increase in 1997, and since then has been in a state of high variability [45]. Dramatic changes in temperatures have a strong impact on mountainous hydrological processes and water resources, and increase the runoff in the TS due to the accelerated glacier/snowmelt [46,47]. Therefore, under the background of such climate change, examining the correlation between SLA and the meteorological factors can help to understand the dynamic behavior of snow cover.
This study systematically investigated the SLA dynamics by taking advantage of long-term, continuous observations from MODIS snow cover products over the TS. The objectives of this study are to (1) determine the SLA and generate the daily SLA dataset over the TS from 2001–2019; (2) analyze the spatiotemporal variations of SLA; (3) investigate the potential influence of topographic factors (slope gradient and aspect) on SLA, and the correlation between SLA and meteorological factors (precipitation, temperature, and solar radiation).
2. Data and Methodology
2.1. Study Area
As the largest mountain system in Central Asia, the TS stretches approximately 2500 km long and 250–350 km wide within 38°–47° N and 67°–95° E, spanning regions from Uzbekistan to Kyrgyzstan, and from southeastern Kazakhstan to Xinjiang of China (Figure 1) [48]. The TS is bordered by the Junggar Basin to the north and the Tarim Basin to the south. The major peaks in the TS stand from about 4000 to 6000 m above sea level—with the highest peak, Thomuer, at 7439 m, and the second highest peak, Khan Tengri, at 6995 m—and make up an important water tower in this in semiarid region [49]. The TS has abundant precipitation due to westerly circulation and unique topography, exhibiting heavily glaciated and snow-covered regions; annual precipitation in areas with altitude greater than 2000 m mostly exceeds 500 mm and can even reach more than 2000 mm; there are 15,953 glaciers with a total area of 15,416 km2 [26,50,51]. Based on individual mountain ranges, drainage, and climate features, the TS are geographically divided into four regions [1,26,28,52]: the Western Tienshan Mountains (WTS), the Northern Tienshan Mountains (NTS), the Central Tienshan Mountains (CTS), and the Eastern Tienshan Mountains (ETS) (Figure 1). The climate characteristics and cryosphere change show obvious differences in these four subregions [26,48,49,52]. The total area of the study area is 135.5 × 104 km2; and the areas of these subregions account for 13.9% (CTS), 34.7% (WTS), 20.3% (NTS), and 31.1% (ETS), respectively (Figure 1).
2.2. Data
2.2.1. MODIS Fractional Snow Cover (FSC) Data
In this study, the MOD10A1 data for 2001–2019 distributed by the National Snow and Ice Data Center (NSIDC) [53] are utilized. There are MOD10A1 data of two versions (Collection 5 and Collection 6) through the NSIDC platform. The MOD10A1 C5 for 2001–2016 and MOD10A1 C6 for 2017–2019 are used in our work. MOD10A1 C5 provides both binary snow cover (the pixels are either snow or no-snow) and fractional snow cover (FSC) [33]. The MODIS FSC mapping algorithm, developed by Salomonson and Appel [54], is based on a statistical-linear relationship developed between the normalized-difference snow index (NDSI) from MODIS and the true subpixel fraction of snow cover as determined using Landsat scenes. Evaluation studies have demonstrated that MODIS FSC data are highly accurate (mean absolute error less than 0.1) [38,39,54,55].
FSC data are no longer provided in the MOD10A1 C6, therefore, for example, the NDSI_Snow_Cover data is used instead. In order to keep the data consistent in this study, the FSC data was obtained based the C6 NDSI_Snow_Cover product using the linear regression method [56,57]. The FSC calculation formula for MOD10A1 C6 is as follows:
FCSMOD= (−0.01 + (1.45× NDSI)) × 100(1)
The original MOD10A1 FSC data has a sinusoidal projection, which is mosaiced and resampled from the original 463.3 m pixel size to 500 m, and converted to a WGS84-UTM geographical projection using the MODIS reprojection tool (MRT). The final mosaiced images are converted to GeoTIFF file format.
2.2.2. Meteorological Data
The ERA reanalysis, from the European Center for Medium-Range Weather Forecasts (ECMWF), is produced with a sequential data assimilation scheme, advancing forward in time by using 12-hourly analysis cycles [58,59]. The meteorological data in the TS from 2001 to 2019 are derived from the most recently released fifth generation of ECMWF reanalysis, ERA5. ERA5 reanalysis extends back to 1979 and provides hourly data at a high resolution (0.25 degree ~ 31 km). To facilitate many climate applications, the monthly-mean averages have been calculated as well in the datasets of “ERA5 monthly averaged data on single levels from 1979 to present” [60]. This dataset has been regridded to a regular lat-lon grid of 0.25 degrees for the reanalysis. To survey the linkages between SLA changes and climate variations, the monthly temperature, precipitation, and surface net solar radiation data in the TS from 2001 to 2019 are extracted from this dataset.
2.2.3. Other Data
The Shuttle Radar Topography Mission (SRTM) is a set of freely available global digital elevation data covering over 80% of the globe. In this study, the data are available at 1 arc-second (∼30 m) resolution, with vertical accuracy reported as less than 16 m from the U.S. The grid size of the DEM is resampled to 500 m in order to remain consistent with the MODIS snowline pixels, and 500 m DEM is used to derive the SLA and calculate the slope gradient and aspect within the study area.
In order to facilitate the analysis of SLA dynamics and the extraction of ERA5 meteorological grid data (about 31 km spatial resolution), the entire study area is divided into 1715 grids (i.e., the numbers of grids for the NTS, CTS, ETS, and WTS are 350, 240, 526, and 599, respectively) with a cell size of 30 km. Eighty-two Landsat Operational Land Imager (OLI) images (Table 1) covering the four grids for validation are used as reference data for the SLA evaluation. The selection of these Landsat images primarily depends on cloud cover.
2.3. Methodology
Taking the 130th day of 2016 as an example, the process of SLA determination from MODIS FSC products is presented in Figure 2, including the cloud removal from MODIS FSC data, extracting the largest lake mask, and determination of snowline and SLA.
2.3.1. Cloud Removal from MODIS FSC Data
The extensive cloud cover over snow-covered region is a critical issue that greatly limits the applications of the MODIS snow cover product to monitor the snow cover and snowline [33,38]. A validation study of MODIS snow cover images over the Tibetan Plateau found that about 47.3% of the areas were cloud covered [39]. Previous studies suggest that the FSC maps more accurately represent the gradual changes of snow cover in each pixel than the binary snow cover maps [33,54], thus the use of the FSC data could be better for the removal of the cloud cover by temporal filtering [39]. Following the cloud removal method for MODIS FSC products developed by Tang et al. [39], the daily cloud-removed FSC data in the study area are produced from 2001 to 2019. The cloud removal algorithm is based on the cubic spline interpolation algorithm (temporal filtering). Details on the cubic spline interpolation cloud removal method and the relevant accuracy evaluation strategies can be found in the works of Tang et al. [39], and the method is also similar in principle to that proposed by Dozier et al. [61]. From the applications of the cloud removal method in the Tibetan Plateau and the TS [30,39,62,63], the cloud removal method is efficient in retrieving the FSC information of the cloud-covered pixels in these areas, with the overall mean absolute error less than 0.1. There is a high consistency between MODIS-derived snow-covered days (SCD) and the in-situ observed SCD, with a mean consistency over 85%, and a mean absolute error of less than 4.2 days [30,39]. The higher consistency between MODIS-derived SCD and in situ SCD indicates that the cloud-removed MODIS FSC data have a high accuracy to monitor the snowline in the TS.
2.3.2. Extracting the Largest Lake Area Mask
There are around 1667 lakes in the TS, with total area 96.5 ± 14.2 km2 of the TS. The number of lakes along with the total lake area have been in a state of continuous increase in the past few decades [52]. In general, the inland water bodies have been automatically identified (the code is 237) in MODIS snow products. However, due to the existence of lakes, it is easy to misclassify snow in the cold season. Seasonally frozen lakes (lake ice) are often erroneously identified as snow-covered areas due to the similar spectral characteristics of snow and ice. In such cases, the low-elevation lake boundaries are extracted as snowline, which causes incorrect results.
To solve this problem, the largest lake area mask (Figure 2b) is produced from a composite of 19 years (2001–2019) of MOD10A1 FSC data (Figure 2a) to eliminate the influence of lake ice. The inland water that is identified in each pixel in the 19 years of MOD10A1 FSC data for more than one day is taken as the lake pixel. Simultaneously, some snowline pixels located in the largest lake area mask were eliminated in the process of snowline extracting. Figure 3 presents the spatial distribution of the largest extracted lake area mask. Meanwhile the mean SCD is used as the base map, which is derived from the cloud-removed daily MODIS FSC datasets during 2001–2019.
2.3.3. Determination of Snowline and SLA
The determination of snow cover through the availability of the cloud-removed MODIS FSC data is an important step. The snowline is theoretically equivalent to the snow cover outline, which is at the edge of the snow-covered areas. A binary map of snow is tuned to classify a pixel as snow if its coverage is greater than 50% (i.e., FSC ≥ 50%) and snow-free otherwise (Figure 2c). Thus, the pixels at the edge of the snow-covered areas (snowline pixels) are extracted; for each snow-covered pixel, if there is one or more of its surrounding pixels that are no-snow and the pixel is not located in the largest lake area mask, we marked it as a snowline pixel (Figure 2d).
Once the snowline pixels have been marked, the SLA (Figure 2d) can be determined by overlaying the snowline image on the 500 m DEM. Then the determination of SLA for every 30 km grid was carried out by means of the SLA of 500 m snowline pixels.
2.3.4. SLA Dynamics Analysis
In this study, the application of the Mann–Kendall (M–K) test is to assess the temporal significance of trends in meteorological factors and SLA from 2001 to 2019 on a grid-by-grid basis. The M–K test is a nonparametric method of monotonic trends that has successfully been applied to detect trends in a time series [64,65,66,67]. The statistical significance of temporal trends is assumed by Z value, and the significance of the trends is evaluated with a confidence level of p < 0.05 [68]. Sen’s slope is used to estimate the magnitude of trends in meteorological factors and SLA from 2001 to 2019 in this study. The advantage of Sen’s slope is to deal with data of non-normal distribution, and it is used instead of a linear regression because it limits the influence that the outliers have on the slope [69,70]. In addition, correlation analysis was performed to analyze the relationship between SLA, temperature, precipitation, and solar radiation on a grid-by-grid basis.
3. Results
3.1. Comparison of SLA Derived from MODIS FSC and Landsat OLI Images
To evaluate the method, we compared the SLA of four grids extracted from Landsat images (Table 1) with the extraction results from MODIS. The process of SLA determination for a grid from Landsat images can be sketched as follows. First, Landsat images are classified as snow or non-snow using the SNOWMAP approach [71]. Then, the boundaries (snowline pixels) between snow cover and land are extracted. Lastly, the SLA of the grid is determined as the average elevation of the multiple snowline pixels by combination with the 30 m SRTM DEM. The mean absolute error and root mean square error are employed to evaluate the reliability of the MODIS-derived SLA (Table 2). In the four validation grids, the mean absolute error of MODIS-derived SLA compared with the Landsat-derived values is between 7.8 and 14.6 m, and the root mean square error is between 9.7 and 16.7 m. We believe that the MODIS-derived SLA with such accuracy in the 30 km grids can be applied to investigating the spatiotemporal patterns of SLA in the TS.
3.2. Spatiotemporal Patterns of SLA
Figure 4 presents the spatial pattern of multiyear (2001–2019) mean SLA for different months over the TS. The spatial distribution of the SLA over the entire TS exhibits a large spatial heterogeneity and corresponds well with the distribution of snow duration (SCD in Figure 3). Apparently, in the peripheral mountainous area, SLA is relatively low. By contrast, SLA is high in the vast interior area (especially in the CTS). Furthermore, Figure 5 and Figure 6 show the time series of SLA for the four subregions from 2001 to 2019, and the annual cycles of SLA over the entire TS and the four subregions, respectively. Strong seasonal variations in SLA are found over the whole area of the TS (Figure 4, Figure 5 and Figure 6). The seasonal variations range of SLA is from about 2000 to 4100 m. The SLA from the end of November to early March of the following year is lower than 2400 m, with relatively large standard deviations reflecting the high interannual variability, whereas in the June–September period, the average SLA is greater than 3900 m and with a relatively small interannual variability (Figure 6b). Apparently, the SLA over the TS progressively increases, beginning in March due to the snowpack ablation, and reaches the highest in August. A steady decrease of SLA begins in September along with the accumulation of snow, and reaches relatively low (average SLA of 2209 m) in winter (Figure 6b).
However, for different subregions, the differences in the seasonal fluctuations of SLA exhibit distinct climatological characteristics (Figure 6a). The NTS is situated under the strong influence of the Siberian anticyclonic circulation, which brings precipitation mainly in the form of snowfall in the cold season; simultaneously, it is also influenced by frontal cyclonic circulation and northern jet stream, which bring considerable precipitation even in the cold season [49]. As a result, the lowest SLAs occur in the NTS (Figure 6a). In winter, the WTS is weakly influenced by the Siberian anticyclonic circulation, and the southwest cyclonic circulation is moderate in winter, bringing in warm moist air and maximum precipitation; however, winter precipitation is relatively rare for the CTS and ETS [30,49]. This may explain why the SLA in the WTS is relatively lower than that in the CTS and ETS (Figure 6a). The CTS has high mountains, so the transport of water vapor from the Atlantic is blocked. The CTS is also adjacent to the Taklamakan Desert, and has dry air, unlike the other subregions of the TS. In addition, the topographic-altitude-controlled spatial pattern of the SLA could be ascribed to the mass elevation effect, which is essentially the result of the thermodynamic effect of mountain masses and has been recognized as a significant contributor to the vertical distribution of mountain snowline and timberline [72,73]. Thus, SLA in the CTS is dramatically higher than the other regions.
To explore the variations of SLA in detail, the slope (Sen’s slope) of the trend in SLA for the whole study area and the counts of SLA grids at the 5% significance level were calculated on a grid-by-grid and month-by-month basis during the period 2001–2019 (Table 3 and Table 4). Besides, the change trend (slope) of SLA and its significance level in August for the 19 years are depicted in Figure 7 as an example. On the whole (excluding the null value grids), the SLA over the entire TS shows a rising trend in August (the average slope of the grids is 1.95 m/a), although a large number of grids (53.5%) are characterized by weak trends in SLA (−1.5 < slope < 1.5 m/a) (Figure 7). Overall, 58.1% of the grids are increased (0 < slope < 9.7 m/a), but only 28% of the grids are characterized by decreasing trends (−10.1 < slope < 0 m/a). For the different subregions, the average slopes of NTS, ETS, WTS, and CTS are 1.36 m/a, 2.79 m/a, 2.08 m/a, and 1.29 m/a, respectively. As shown in Table 3 and Table 4, the SLA over the entire TS from 2001 to 2019 generally displays an increasing trend. In particular, the increasing trends are notable both at the start of the melt season SLA (March and April) and the end of melt season SLA (July and August) for the four subregions. Despite the number of significantly decreased SLA grids (93) exceeding the number of the significantly increased grids (10), the SLA in November over the entire TS is rising. A slight decrease of the SLA in June is found for the entire TS (−1.17 m/a) during 2001–2019, with most of the decrease stemming from the ETS.
3.3. The Influences of Topographic Factors on SLA
In high mountain regions, like the TS, SLAs in areas with similar climatic conditions might respond differently given their different and complex topographies. In order to investigate the variation characteristic of SLA in different topographic conditions, the aspect and slope gradient calculated by DEM are used in this study. The aspect of the TS is set as the eight ranges: north (0°–22.5°, 337.5°–360°), northeast (22.5°–67.5°), east (67.5°–112.5°), southeast (112.5°–157.5°), south (157.5°–202.5°), southwest (202.5°–247.5°), west (247.5°–292.5°), and northwest (292.5°–337.5°). Likewise, the slope gradient is divided into four ranges (0°–10°, 10°–20°, 20°–30°, and >30°).
From the SLA over the TS relative to different slope gradients located in four aspect ranges (Figure 8), it is evident that the SLA increases with the steeper slope gradient. As the slope gradient gradually increases, the SLA changes greatly. SLA within the slope gradient of 0°–10° is considerably lower than other slope gradients in winter. In addition, areas having a southerly aspect receive more solar radiation than the north, resulting in more heat absorption and higher SLA (Figure 8). In contrast, SLAs in the Himalayan region are lower on the southerly aspect because of the influence of two atmospheric circulations: the Indian Ocean and the summer southwest monsoon [36]. The two circulation systems, combined with the huge topographic landform, exert climatic controls on the distribution of SLA in the Himalayan region. Due to significant shielding of water vapor transport by high topography as atmospheric travels north across the region, the southern aspects are enriched and the northerly aspect depleted in their abundance of precipitation.
3.4. Spatiotemporal Characteristics of Meteorological Factors Over the TS
Monthly distribution of precipitation, temperature, and solar radiation in the TS during 2001–2019 shows obvious seasonal variations (Figure 9). The annual mean precipitation across the TS is 469.4 mm, which is concentrated between April and July; the annual mean temperature is 7.5 ℃; and the annual mean radiation is 143.1 W m−2 (Figure 9a). The temperature and solar radiation tend to be remarkably consistent in their changes during the year. Spatial distributions of the precipitation, temperature, and solar radiation are shown in Figure 9b–d. The WTS and NTS have the greater annual precipitation, and are characterized by a relatively humid climate, while the CTS and ETS have a typical continental climate. Besides, areas peripheral to the TS are characterized by less precipitation. In contrast, the temperature and solar radiation are low in the central region and high in the marginal region.
3.5. The Influences of Meteorological Factors on SLA
To explore possible mechanisms for SLA changes of the TS, we examined the correlation between SLA and the variations of three important meteorological factors (temperature, precipitation, and solar radiation) by correlation analysis. In this study, the monthly averaged ERA5 climate reanalysis data were selected to investigate the relationship between SLA, temperature, precipitation, and solar radiation at monthly scales. Taking October as an example, the Pearson correlation coefficients between SLA, temperature, precipitation, and solar radiation during 2001–2019 are calculated on a grid-by-grid basis (Figure 10). As can be seen from Figure 10, the correlation between SLA, temperature, and solar radiation shows obvious consistency. The SLA increases with the increased temperature and radiation in most areas of the TS (particularly in the WTS and CTS). It is evident that the SLA decreases with increased precipitation, except for a small part of the ETS (Figure 10). Additionally, correlation coefficients between SLA, temperature, and solar radiation vary in different regions and months (Table 5). In March–June, it shows a significant correlation between SLA, temperature, and solar radiation in the entire TS, while the correlation between SLA and precipitation is insignificant. From September to November, these meteorological factors have a great influence on SLA changes, with solar radiation being the dominant one. In the winter period (i.e., from December to February of the following year), the SLA shows no obvious correlation with temperature, precipitation, or solar radiation. Because of the temperature being well below freezing in winter, there is little temperature-related snow melting. The increased sublimation of the snow may cause the reduction of winter snow at high altitudes where it is frequently accompanied by high winds. More than half of the snow mass in the Tibetan plateau is lost by sublimation in winter [74].
4. Discussion
The TS is the main water source and ecological barrier of Central Asia, with unique climatic and hydrological conditions. Based on MODIS snow cover products, this study extracts spatiotemporal continuous SLA of the entire TS for 19 years. The uncertainty of the MODIS-derived SLA may come from the following: (1) the errors that occurred due to the pixel size of the remote sensing images, slope, and aspect of the terrain, affecting the accuracy of the georeferencing and the quality of the DEM [75,76]; or (2) the errors that occurred in the MODIS FSC mapping algorithm [38,77] and cloud removal method [39], although the MODIS SCD threshold is calibrated in this method. For MODIS FSC mapping, Masson et al. has shown that NDSI linear regression produced better results when used alongside atmospheric and topographic correction [78]; a new linear regression model could be defined if a large series of validated and very accurate ground truth data were to be established. In addition, the linear spectral unmixing approach is a more reliable method for estimating snow cover fraction, and tends to produce more stable results than the NDSI linear regression method [78]. Therefore, future research should examine the SLA method presented here with the MODSCAG and MODImLAB approaches [79,80]. Given the high resolution, the Landsat and Sentinel-2 images are suitable for high-resolution SLA extraction. However, the cloud cover and longer revisit period greatly limit their ability. Google Earth Engine (GEE) is a cloud-based platform that makes it easy to access high-performance computing resources for processing very large geospatial datasets [81]. The GEE platform integrates the Landsat and Sentinel-2 images from over the last few decades, which can be accessed and processed by the Earth Engine application programming interface (API). Therefore, the GEE platform provides a possible method for SLA extraction in a continuous time and space for a large-scale area.
The MODIS-derived SLA method and gridded approach in this study can be efficiently applied to SLA extraction for other large-scale areas. For a catchment scale, however, the regional snowline elevation can be better obtained by the method developed by Krajčí et al. in the case of cloud cover (<70%) [15]. However, this method must estimate an elevation (i.e., the regional snowline elevation) where the number of snow-covered pixels below and the number of snow-free pixels above are minimized, and only obtains one SLA value for the entire basin [15]. If we want to estimate SLA dynamics and their long-term trends at the catchment level, the recent study gives an available method of random forest regression [82] combining high resolution satellite imagery and climate reanalysis data.
Under the background of global warming, and the state of high variability in temperature of the TS since 1997 [45], the cryosphere of the TS has been changing rapidly [28,83,84,85]. Like many glaciers worldwide, the glaciers in the TS have generally been losing mass over the past decades [28,86,87]. In this work, we find the SLA in the TS generally shows a rising trend in the 19 years (Table 3 and Table 4). This increasing trend is especially significant at the end of melting season (July and August), and especially in August, with average slopes of NTS (1.36 m/a), ETS (2.79 m/a), WTS (2.08 m/a), and CTS (1.29 m/a). These rising trends of SLA at the end of melting season indicate a decrease in the annual mass balance of glaciers across the 19 years. Glaciers in the TS are generally losing mass (i.e. average annual mass balance is already negative) [28,85]. These results are consistent with the conclusions drawn from previous studies on the relationship between the glacier mass balance and SLA at the end of melting season [23,25,88]. It would be more beneficial to reconstruct the annual mass balance time series if higher resolution data could be used for estimating spatial and temporal continuous SLA over a large-scale area.
Due to the specific geographical location and climatic conditions, the snow cover and mountain glaciers in the TS are very vulnerable to climate change. Our results show a strong correlation between SLA, temperature, precipitation, and solar radiation (Figure 10), although the absolute values of the correlation coefficients vary due to the different subregions. That is, the SLA increases with increased temperature/radiation and decreased precipitation. Solar radiation is the dominant climatic factor affecting the changes of SLA in the TS (Table 5). Due to the reduction in surface albedo caused by the decrease in glacier and snow coverage, more solar radiation will be absorbed, and will then amplify cryosphere warming, further intensifying ablation and snowmelt [89,90]. Nevertheless, the effect of climate change on SLA is very complex. The 19 years of available MODIS information is insufficient for definitive conclusions about climate change. Longer time series of data need to be examined in further studies to obtain some more definitive conclusions about temporal trends of SLA and their relationship with climate change.
5. Conclusions
In this study, a large-scale SLA investigation methodology is developed using the cloud-removed MODIS FSC products (MOD10A1). The 500 m resolution of a daily SLA dataset is generated over the TS for the period 2001–2019. The spatiotemporal patterns of SLA over the entire TS during 2001–2019 are investigated using the gridded approach, with specific attention to the four subregions. We also explore the potential influence of topographic factors (slope gradient and aspect) on SLA and the correlation between SLA and meteorological factors (temperature, precipitation, and solar radiation) on a grid-by-grid and month-by-month basis. The main findings are summarized as follows:
(1) The large-scale SLA monitoring method is efficient in monitoring the spatiotemporal patterns of the SLA in the TS. Compared with the Landsat-derived SLA of the four validation grids (30 km), the mean absolute error of MODIS-derived SLA is between 7.8 and 14.6 m, and the root mean square error is between 9.7 and 16.7 m.
(2) Our results show strong seasonal fluctuations of SLA over the TS from 2001 to 2019, ranging from about 2000 (in late December) to 4100 m (in early August). The distribution of SLA over the TS shows a large spatial heterogeneity due to complex topography and geomorphology. The SLA increases with the steeper slope gradient. The SLA of the northerly aspect is generally less than that of the south, due to receiving more solar radiation.
(3) The SLA over the entire TS from 2001 to 2019 shows a rising trend. Except for a slight decrease in June, the SLA increased in all other months, especially at the start of the melting season (March and April) and the end of the melting season (July and August). The SLA increases with increased temperature/radiation in the TS (particularly in the WTS and CTS) and decreases with increased precipitation. Solar radiation is the dominant climatic factor leading the SLA dynamics, and temperature has a greater influence on SLA than precipitation. This study could enrich the understanding of the response of snow cover dynamics to climate, and provide scientific support for eco-environment sustainable management in the high mountain region. If global warming continues, the melting of snow and peak runoff in the TS will increase. The early snowmelt may lead to a decrease in summer flow, which in turn changes the flow regimes and water availability, thereby affecting the ecosystem, agriculture, and water resources in the densely populated downstream areas.
Author Contributions
All authors contributed to the design and development of this manuscript. Conceptualization, Z.T.; methodology, G.D. and Z.T.; validation, G.D., G.H., and J.W.; formal analysis, G.D. and Z.T.; investigation, G.D. and G.H.; data curation, G.S., J.L., and J.W.; writing—original draft preparation, G.D.; writing—review and editing, G.D., Z.T., and G.H.; supervision, Z.T.; project administration, Z.T.; funding acquisition, Z.T. All authors have read and agreed to the published version of the manuscript.
Funding
This study was supported financially by the National Natural Science Foundation of China (No. 41871058), the State Key Laboratory of Cryospheric Science, Northwest Institute of Eco-Environment and Resources, Chinese Academy Sciences (No. SKLCS-OP-2020-08); the Foundation for Innovative Research Groups of the Natural Science Foundation of Hunan Province, China (No. 2020JJ1003); and the Scientific Research Foundation of Hunan Education Department, China (No. 20B227).
Institutional Review Board Statement
Not applicable.
Informed Consent Statement
Not applicable.
Data Availability Statement
Data sharing not applicable.
Acknowledgments
We would like to thank the editors and reviewers (especially Keith A. Brugger) for their helpful comments and constructive suggestions. MODIS fractional snow cover data were obtained from the National Snow and Ice Data Center (NSIDC) (
Conflicts of Interest
The authors declare no conflict of interest.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figures and Tables
Enlarge this image.Figure 1. Location and the extent of the TS, and boundaries of the four subregions.
Enlarge this image.Figure 2. Process of SLA determination, including figures of (a) a subset of MODIS FSC image (130th day of 2016), (b) the cloud-removed MODIS FSC and the largest lake area mask image, (c) the binary snow image, and (d) final snowline image with the SLA values.
Enlarge this image.Figure 3. Distribution of mean MODIS-derived SCD, and the largest lake area over the TS.
Enlarge this image.Figure 6. The annual cycle of SLA for the four subregions (a), and the whole region of the TS (b). The SLA values in (a, b) are averages of 19 years from 2001 to 2019. The error bars in (b) show the standard deviation, indicating the interannual variations of SLA from 2001 to 2019.
Enlarge this image.Figure 7. Change trend (slope) of SLA and its significance level in August during 2001–2019 over the TS.
Enlarge this image.Figure 8. SLA over the TS relative to different slope gradients (0°–10°, 10°–20°, 20°–30°, and >30°) of different aspects (north, northeast, east, southeast, south, southwest, west, and northwest) in spring (a), summer (b), autumn (c) and winter (d). The SLA values are averages of 19 years, from 2001 to 2019.
Enlarge this image.Figure 9. (a) Monthly distribution of precipitation, temperature, and solar radiation in the TS during 2001–2019. Annual spatial distributions of the precipitation (b), temperature (c) and solar radiation (d) over the TS in the period 2001–2019.
Enlarge this image.Figure 10. Spatial distribution of the Pearson correlation coefficients between SLA, temperature, precipitation, and solar radiation in October during 2001–2019.
Information about Landsat Operational Land Imager (OLI) images used in validation of SLA—snowline altitude.
| Year | Grid 1 |
Grid 2 |
Grid 3 |
Grid 4 |
|---|---|---|---|---|
| 2013 | 6 August, 7 September, 23 September | 27 July, 12 August, 28 August, Sep 29 | 19 July | 9 September, 25 September |
| 2014 | 10 September, 26 September | 14 July, 15 August, 31 August | 6 July, 6 July, 7 August, 23 August | 10 July, 26 July, 11 August, 12 September |
| 2015 | 11 July, 12 August, 29 September | 17 July, 18 August, 3 September, 19 September | 9 July, 11 September, 27 September | 13 July |
| 2016 | 14 August, 15 September | 4 August, 5 September, 21 September | 11 July, 27 July, 28 August, 29 September | 1 September, 17 September |
| 2017 | 2 September, 18 September | 6 July, 22 July, 7 August | 30 July, 31 August, 16 September | 2 July, 19 August, 4 September, 20 September |
| 2018 | 19 July, 4 August, 20 August | 10 August, 26 August, 27 September | 1 July, 3 September, 19 September | 21 July, 6 August |
| 2019 | 7 August, 23 August, 8 September, 24 September | 12 July, 28 July, 13 August, 29 August | 4 July, 5 August, 21 August, 22 September | 24 July, 25 August |
| Total | 19 | 24 | 22 | 17 |
Comparison of SLA extracted from Landsat OLI images and MODIS.
| Grids | Mean Absolute Error (m) | Root Mean Square Error (m) |
|---|---|---|
| 1 | 14.6 | 16.7 |
| 2 | 10.2 | 11.5 |
| 3 | 7.8 | 9.7 |
| 4 | 9.5 | 10.7 |
The change trend and the counts of SLA grids with significant trends for the entire TS and the four subregions during the period 2001–2019.
| Months | Entire TS | NTS | ETS | WTS | CTS | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| P | N | Trend | P | N | Trend | P | N | Trend | P | N | Trend | P | N | Trend | |
| Jan | 40 | 35 |
|
3 | 6 |
|
13 | 15 |
|
13 | 9 |
|
11 | 5 |
|
| Feb | 36 | 18 |
|
7 | 6 |
|
14 | 7 |
|
8 | 3 |
|
7 | 2 |
|
| Mar | 15 | 5 |
|
3 | 1 |
|
8 | 3 |
|
3 | 1 |
|
1 | 0 |
|
| Apr | 15 | 3 |
|
6 | 2 |
|
4 | 1 |
|
1 | 0 |
|
4 | 0 |
|
| May | 2 | 1 |
|
2 | 0 |
|
0 | 0 | 0 | 1 |
|
0 | 0 | ||
| Jun | 1 | 14 |
|
0 | 0 | 0 | 13 |
|
1 | 1 | 0 | 0 | |||
| Jul | 9 | 0 |
|
0 | 0 | 6 | 0 |
|
1 | 0 |
|
2 | 0 |
|
|
| Aug | 9 | 3 |
|
1 | 0 |
|
4 | 0 |
|
3 | 1 |
|
1 | 2 |
|
| Sep | 10 | 6 |
|
1 | 1 | 4 | 3 |
|
0 | 2 |
|
5 | 0 |
|
|
| Oct | 11 | 9 |
|
0 | 0 | 0 | 5 |
|
3 | 3 | 8 | 1 |
|
||
| Nov | 10 | 93 |
|
1 | 33 |
|
2 | 48 |
|
4 | 9 |
|
3 | 3 | |
| Dec | 67 | 57 |
|
10 | 22 |
|
16 | 18 |
|
32 | 15 |
|
9 | 2 |
|
Notes: The NTS, ETS, WTS and CTS indicate the Northern Tienshan Mountains, the Eastern Tienshan Mountains, the Western Tienshan Mountains and the Central Tienshan Mountains, respectively. P and N indicate the counts of SLA grids with statistically significant positive and negative trends at the 5% significance level, respectively; when P is higher than N, a red up arrow indicates an increasing trend, and when P is lower than N, a blue down arrow indicates a decreasing trend.
Table 4The slope of the change trend for the entire TS and the four subregions during the period 2001–2019.
| Months | Entire TS | NTS | ETS | WTS | CTS | ||||||
|---|---|---|---|---|---|---|---|---|---|---|---|
| SlopeA | Slopep | Slopen | Slopep | Slopen | Slopep | Slopen | Slopep | Slopen | Slopep | Slopen | |
| Jan | 4.74 | 22.49 | −15.54 | 10.51 | −11.58 | 27.63 | −16.84 | 17.93 | −15.96 | 25.09 | −15.63 |
| Feb | 12.21 | 24.83 | −13.03 | 32.52 | −13.59 | 28.06 | −16.31 | 23.24 | −4.22 | 12.51 | −13.07 |
| Mar | 13.53 | 20.70 | −7.99 | 19.47 | −2.90 | 28.17 | −8.14 | 6.07 | −12.64 | 8.47 | 0 |
| Apr | 11.98 | 17.06 | −13.45 | 18.00 | −16.19 | 15.88 | −7.95 | 18.29 | 0 | 16.53 | 0 |
| May | 9.94 | 18.75 | −7.69 | 18.75 | 0 | 0 | 0 | 0 | −7.69 | 0 | 0 |
| Jun | −1.17 | 14.43 | −8.93 | 0 | 0 | 0 | −8.89 | 14.43 | −9.46 | 0 | 0 |
| Jul | 1.31 | 9.03 | 0 | 0 | 0 | 9.95 | 0 | 10.10 | 0 | 5.77 | 0 |
| Aug | 1.95 | 7.85 | −8.42 | 7.81 | 0 | 7.18 | 0 | 8.14 | −10.10 | 9.67 | −7.58 |
| Sep | 4.05 | 14.86 | −13.97 | 26.75 | −13.73 | 7.14 | −13.68 | 0 | −14.52 | 18.67 | 0 |
| Oct | 1.12 | 17.04 | −18.34 | 0 | 0 | 0 | −18.55 | 14.65 | −18.17 | 17.93 | −17.79 |
| Nov | 1.40 | 15.47 | −22.52 | 15.61 | −26.41 | 13.04 | −22.43 | 14.29 | −11.99 | 18.61 | −12.81 |
| Dec | 2.94 | 19.82 | −17.25 | 14.43 | −16.40 | 15.06 | −16.68 | 22.20 | −19.04 | 25.80 | −17.15 |
Notes: Slope indicates Sen’s slope estimator; SlopeA indicates the average slope of the SLA grids; Slopep and Slopen indicate the average slope of the SLA grids with statistically significant positive and negative trends at the 5% significance level, respectively.
Table 5Pearson correlation coefficients between SLA, temperature, precipitation, and solar radiation in different subregions at monthly scales for the period of 2001–2019.
| Months | Temperature | Precipitation | Radiation | |||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| NTS | ETS | WTS | CTS | NTS | ETS | WTS | CTS | NTS | ETS | WTS | CTS | |
| Jan | −0.008 | −0.004 | −0.025 | −0.024 | 0.014 | 0.016 | −0.037 | 0.032 | 0.032 | 0.043 | −0.021 | 0.023 |
| Feb | −0.022 | −0.041 | −0.141 | 0.010 | −0.025 | 0.027 | −0.080 | −0.162 | 0.015 | −0.025 | −0.038 | 0.148 |
| Mar | 0.413 | 0.300 | 0.641 ** | 0.450 | −0.133 | −0.073 | −0.174 | −0.307 | 0.430 * | 0.347 | 0.627 ** | 0.545 ** |
| Apr | 0.678 ** | 0.559 ** | 0.586 ** | 0.624 ** | −0.205 | −0.180 | −0.408 | −0.412 | 0.681 ** | 0.616 ** | 0.580 ** | 0.715 ** |
| May | 0.619 ** | 0.523 * | 0.710 ** | 0.678 ** | −0.155 | −0.130 | 0.004 | −0.036 | 0.627 ** | 0.624 ** | 0.610 ** | 0.687 ** |
| Jun | 0.421 | 0.336 | 0.717 ** | 0.554 ** | −0.220 | −0.245 | −0.251 | −0.155 | 0.640 ** | 0.492 * | 0.671 ** | 0.511 * |
| Jul | 0.279 | 0.439 | 0.135 | 0.402 | −0.137 | −0.151 | −0.022 | −0.147 | 0.220 | 0.250 | −0.045 | 0.377 |
| Aug | 0.244 | 0.375 | 0.292 | 0.296 | −0.318 | −0.227 | −0.010 | −0.383 | 0.321 | 0.300 | 0.120 | 0.406 |
| Sep | 0.569 ** | 0.557 ** | 0.520 * | 0.641 ** | −0.599 ** | −0.468 * | −0.385 | −0.555 ** | 0.773 ** | 0.652 ** | 0.598 ** | 0.708 ** |
| Oct | 0.564 ** | 0.568 ** | 0.490 * | 0.498 * | −0.604 ** | −0.382 | −0.650 ** | −0.432 * | 0.769 ** | 0.615 ** | 0.768 ** | 0.592 ** |
| Nov | 0.501 * | 0.390 | 0.611 ** | 0.399 | −0.496 * | −0.315 | −0.623 ** | −0.391 | 0.642 ** | 0.413 | 0.753 ** | 0.485 * |
| Dec | −0.086 | −0.038 | 0.194 | 0.055 | −0.084 | −0.043 | −0.375 | −0.212 | 0.121 | 0.093 | 0.484 * | 0.268 |
Note: * and ** indicate statistical significance at the 0.05 and 0.01 level, respectively.
© 2021 by the authors.