The study area, Central Asia (35° -55° N, 45° -95° E; Fig. 1), includes Xinjiang Uygur Autonomous Region of China and five central Asian countries: Kazakhstan, Uzbekistan, Kyrgyzstan, Turkmenistan and Tajikistan. It covers an area of approximately 5.6× 10⁶km² and extends to the Caspian Sea to the west and the Gobi deserts of China to the east, Russia to the north, and Pakistan, Iran and Afghanistan to the south. The mountainous regions are mainly situated in the southeastern Central Asia and account for about 17.3% of the total area. The Kunlun Mountains, Pamir Plateau, Tianshan Mountains, and Altay Mountains are considered as the mountainous regions.

Fig. 1

	Figure Option ViewDownloadNew Window
	Fig. 1 Overview of Central Asia with MODIS tile grids (parallelograms shown as dashed lines) and major lakes (blue polygons). Data source for administrative information was cited from Hu et al. (2014), and the locations of lakes were obtained from ArcGIS 10.2.

Most part of Central Asia are dominated by typical temperate continental arid and semi-arid climate. Only southeastern Pamir Plateau belongs to alpine climate (Lioubimtseva and Henebry, 2009; Yin et al., 2016). The western and northwestern plain regions are open to cold northerly and northwesterly air inflows and also influenced by westerly moist influx from the North Atlantic. The southeastern mountainous regions e.g., Pamir Plateau, Tianshan Mountains, and Kunlun Mountains) are almost completely isolated from the moist air masses originating from the Indian Ocean. Therefore, the spatial distribution of precipitation shows strong gradients from north to south and also from mountains to lowlands (Lindzen and Nigam, 1987; Liu et al., 2010). For the entire Central Asia, the snowfall events mainly occur from the preceding November to succeeding April in one hydrological year (from 1 September to 31 August next year).

The 8-day snow products MOD10A2 (V005) and MYD10A2 (V005) (Hall et al., 2001) provided by the National Snow and Ice Data Center (NSIDC) were selected to generate the time series of cloud-free and high-accuracy snow-cover dataset to characterize the spatial and temporal variations of snow cover in Central Asia. For one hydrological year, there are 46 scenes of snow-cover data. NDSI values were used to identify snow pixels with a series of threshold values (Hall et al., 2002). Eleven categories of pixel cover types are listed in Table 1.

Table 1

Table 1

Table 1 Pixel cover types of the MOD10A2 and MYD10A2 products Pixel value Pixel cover types Pixel value Pixel cover types 0 Missing data 50 Cloud 1 No decision 100 Lake ice 11 Night 200 Snow 25 Land without snow cover 254 Detector saturated 37 Lake 255 Fill 39 Ocean

Table 1 Pixel cover types of the MOD10A2 and MYD10A2 products

Altitude information was required for the cloud removal algorithm and the analyses of snow-cover variations. We adopted the Digital Elevation Model data with the spatial resolution of 90 m (Version 4.0) derived from SRTM (Wang et al., 2008). The data are of global coverage and can be downloaded from the website, http://srtm.csi.cgiar.org.

Three Landsat ETM+ scenes with cloud coverage less than 5% and 30-m spatial resolution (Fig. 1) were selected to represent one plain (L1) and two mountainous (L2 and L3) regions for validating the accuracy of the cloud-free snow-cover dataset. More details are shown in Table 2. The preprocessing of Landsat ETM+ data, including radiometric calibration and atmospheric corrections, were implemented in ENVI software.

Table 2

Table 2

Table 2 Information of the Landsat ETM+ data for validation Scene Sensor Imaging time Path Row Cloud coverage (%) L1 ETM+ 19 Nov, 2014 167 25 2 L2 ETM+ 28 Dec, 2014 152 31 1 L3 ETM+ 15 Dec, 2014 141 33 5

Table 2 Information of the Landsat ETM+ data for validation

Additionally, global ecological regionalization data with 867 ecological regions (derived from National Earth System Science Data Sharing Infrastructure, National Science & Technology Infrastructure of China) were implemented in this study. In Central Asia, we selected the key ecological regions for plain regions (Badghyz and Karabil (BK) semi-desert, southern Central Asian (SCA) desert, northern Central Asian (NCA) desert, Kazakh semi-desert, Kazakh steppe, Kazakh upland, and Kazakh forest steppe) and for mountainous regions (Kunlun Mountains, Pamir Plateau, Tianshan Mountains, and Altay Mountains) to investigate the difference of SCR and variations among these regions.

Although the MODIS 8-day products have better quality than the daily products (Zhou et al., 2005; Ault et al., 2006; Wang et al., 2008), it is still not realistic to apply the data directly without removing the cloud-covered pixels. The cloud cover is the major obstacle influenceing the accuracy of the analytical results (Zhou et al., 2005; Hall and Riggs, 2007). Wang et al. (2008) reported a mean accuracy of 19% as the winter cloud mask used in MOD10A2 in northern Xinjiang of China. In this study, a five-step algorithm was designed to estimate the snow status of the cloud-covered pixels.

The first step combines Terra (MOD10A2) and Aqua (MYD10A2) satellite data observed at different times on the same day (Wang and Xie, 2009) to obtain the maximum SCR, and the overall accuracy of this step is 92% (Dietz et al., 2013). The second step is to spatially combine the four directly “ side-bordering” neighboring pixels of the cloud-covered pixels. A cloud-covered pixel is defined as a snow-covered pixel when three of its four directly neighboring pixels are snow. The overall accuracy of this method is 92.38% (Gafurov and Bá rdossy, 2009). The third step involves integrating the environmental lapse rate of temperature variations along the altitude gradients. If the eight neighboring pixels are cloud-free, the center pixel that was covered with clouds will be examined. The accuracy of this step can reach 94% (Gafurov and Bá rdossy, 2009). The fourth step, proposed as the Snowl algorithm by Parajka et al. (2010), is intended to eliminate cloud-covered pixels at extremely low or high altitude. The basic concept is to calculate the minimum and maximum altitude of snow-covered pixels. Based on the altitude of the upper snow line, the cloud-covered pixels are assigned as snow when they are above the upper snow line but are assigned as land when below the lower snowline. The accuracy of this step can reach 95% (Parajka et al., 2010). The fifth step is designed to further check the snow-covered pixels based on the daily status of each pixel in one hydrological year. It is important to note that the dataset of Eight Day Snow Cover Chronology (EDSC) contained in the MOD10A2 and MYD10A2 is also needed to obtain daily snow-cover status of each pixel. The overall accuracy of this step is higher than 78% (Gafurov and Bá rdossy, 2009; Dietz et al., 2013).

To obtain reliable assessment of snow-cover variations, accurate evaluation of the cloud-free snow-cover dataset is essential for the utilization of the datasets. However, it is rather difficult to conduct in-situ measurements in these regions. Therefore, we applied indirect methods to validate the cloud-free snow-cover dataset in this study. The indirect validation method was proposed by Gafurov and Bá rdossy (2009) and improved by Dietz et al. (2013). The overall accuracy depends on two factors: the overall accuracy under clear-sky conditions and the inherent accuracy of the algorithm. Previous studies suggested that the accuracy of cloud-free pixels in MOD10A2 can reach 93% under clear-sky conditions. The overall accuracy of algorithms can be estimated by the weighted average accuracy based on the accuracy of algorithm adopted in each step. The overall accuracy of the cloud-free snow-cover dataset thus can be calculated by Equation 1.

${{T}_{ai}}=[(Pixel{{s}_{clear-sky}}\times 93\%+Pixel{{s}_{step1}}\times 92\%+Pixel{{s}_{step2}}\times 92.38\% \\ +Pixel{{s}_{step3}}\times 94\%+Pixel{{s}_{step4}}\times 95\%+Pixel{{s}_{step5}}\times 78\%)/{{T}_{p}}]\times 100\%$

where T_ai denotes the overall accuracy of the cloud-free snow-cover dataset under the indirect validation method, and T_p is the total number of pixels in a given dataset. Pixels represents the number of cloud-covered pixels that were removed.

A direct validation method was employed to further validate the reliability of the cloud-free snow-cover dataset, . Firstly, we obtained snow-cover maps through visual interpretation of Landsat ETM+ data in the selected cloudless regions. The snow-cover maps derived from Landsat ETM+ and cloud-free 8-day snow-cover dataset were referred to as SCM-ETM+ and SCM-MODIS, respectively. Through the intersection tool in ArcGIS, we obtained the common snow area (C_snow) and common land area (C_land) between SCM-ETM+ and SCM-MODIS. Additionally, the snow area and land (i.e., snow-free) area in SCM-ETM+ were referred to as T_snow and T_land, respectively, which were considered as the truth value to validate the overall accuracy as Equation 2:

${{T}_{ad}}=\left( \frac{{{C}_{snow}}+{{C}_{land}}}{{{T}_{snow}}+{{T}_{land}}} \right)\times 100%$, (2)

where T_ad denotes the overall accuracy of the cloud-free snow-cover dataset under the direct validation method.

The cloud-free snow-cover dataset, which was processed by the five-step approach to remove the effects of the cloud and validated through the direct and indirect validation methods, was adopted to analyze the spatial and temporal variability over the period 2000-2015. The snow-cover frequency (SCF) is defined by Equation 3.

$SCF=\sum\limits_{{}}^{{}}{{{S}_{(x, y)}}}/N$, (3)

where N is the total number of images for cloud-free snow-cover dataset in one hydrological year, and S_{(x, y)}is the status of pixel cover types. When this index is 1, it indicates that the pixel is covered by snow, and 0 indicates non-snow-covered pixels.

The SCR is defined as the percentage of snow-covered area in a single 8-day image data S_area to the total area of Central Asia T_Area in Equation 4.

$SCR=({{S}_{area}}/{{T}_{area}})\times 100%$, (4)

whereS_area(km²) stands for the area of snow cover in a single 8-day image, andT_Area(km²) is the total area of Central Asia. A 5-year running mean and the corresponding confidence interval were calculated to further examine the trend of the SCR during the 15 hydrological years, .

Empirical orthogonal function (EOF) was employed to analyze the anomaly field of the snow-cover frequencies (SCFs) in the spatial and temporal dimensions. A small set of uncorrelated variables, instead of the original information produced by linear transformation for data dimensional reduction, was applied in this method (Xu et al., 2009). This method is widely used for multivariable analysis in meteorology, geology, and geography (Ryu and Jenkins, 2005) with clear advantages of no need for specific functional form, fast convergence rate and powerful information enrichment capability, and decomposing stationary data with an irregular distribution over finite regions (Singh, 2004). It is necessary to analyze the characteristics of the eigenvectors, especially when the number of actual points in space is larger than samples. The actual physical meanings of these eigenvectors can be determined using the method proposed by North et al. (1982), and the errors δiof eigenvalue ican be estimated by Equation 5.

$\delta {{\lambda }_{i}}\cong {{\lambda }_{i}}\sqrt{\frac{2}{N}}$, (5)

$\Delta {{\lambda }_{i}}=\left| {{\lambda }_{i}}-{{\lambda }_{i-1}} \right|$, (6)

where N stands for the number of samples, and λ _i denotes the i^th eigenvalues. λ _i and λ _i-1 represents the latter and the former eigenvalues, respectively. If δ λ _i/Δ λ _i is less than 1, the difference between these two eigenvalues is greater than the errors, implying that λ _i and λ _i-1 are fairly well separated.

Additionally, in a specific EOF mode, the positive (negative) anomaly denotes above (below) the average value of SCR during 2000-2015 in this study.

Based on the ecological regionalization data, we investigated the fluctuating characteristics of the mean SCR every 8 d during 2000-2015 over the different ecological regions (plain and mountainous regions). Furthermore, the seven ecological regions of plain regions were divided into three latitude levels (i.e., low-latitude regions, mid-latitude regions, and high-latitude regions) to clarify the characteristics of snow-cover variations with latitude. Meanwhile, the characteristics of snow-cover variations over four mountainous regions were also analyzed. Eventually, the inter-annual variations of SCR over plain regions with different elevation bends and mountainous regions (Pamir Plateau, Tianshan Mountains, Kunlun Mountains, and Altay Mountains) were investigated using linear least square method during 2000-2015.

Three MODIS tiles were selected to represent the plain (h21v03) and the mountainous regions (h23v04 and h24v05) to validate the accuracy of the cloud-free snow-cover data. These scenes in Figure 2 were randomly selected from two hydrological years (2003/2004 and 2014/2015) after cloud removal. Bars represent the percentage of pixels processed by each step in cloud removal algorithm, and the curves are the monthly average overall accuracy of the cloud-free dataset. The accuracies of the second and third steps in the cloud removal algorithm were combined because a notably smaller percentage of the pixels was processed compared with the other steps.

Fig. 2

	Figure Option ViewDownloadNew Window
	Fig. 2 Diagrams of the monthly average overall accuracy vs. the percentage of the processed pixels in each step of the five-step algorithm for the tiles of the plain (h21v03) and the mountainous areas (h23v04 and h24v05) in 2003/2004 and 2014/2015

It can be seen that the cloud-covered pixels are mainly removed by the first step that combines the Terra and Aqua data and the fourth step that employs the Snowl algorithm. The two steps accounted for at least 30% of the cloud-covered pixels, and the maximum fraction of cloud-covered pixels (92.38% in h21v03) was removed in May of 2003/2004. Generally, the overall accuracies of the three MODIS reference tiles (h21v03, h23v04 and h24v05) fall in the range of 90%-94%.

For the three reference tiles in each month, the overall accuracy of the cloud-free dataset ranges from 92.79% to 93.34%, and the accuracy in winter is slightly lower than those in other seasons. The annual mean accuracy of the cloud-free dataset in the plain areas (93.01%) is slightly higher than that in mountainous areas (92.97%). The overall accuracy of the cloud-free dataset is consistently greater than 92% throughout the selected year, even in December 2003 when cloud coverage reached the maximum (approximately 40.8%).

As complementary of indirect validation method, the direct validation method is also performed. Figure 3 shows the comparison between SCM-ETM+ and SCM-MODIS, including the false color composite image derived from Landsat ETM+ using bands 5, 4, and 3 (Figs. 3a-1, a-2, and a-3), the spatial distribution of snow cover through visual interpretation based on Landsat ETM+ (Figs. 3b-1, b-2, and b-3), and the cloud-free snow-cover dataset based on MODIS (Figs. 3c-1, c-2, and c-3). A relatively good consistency is observed between the snow-cover maps. Through the application of Equation 2, we obtained the overall accuracy of the cloud-free snow-cover dataset (Table 3). The overall accuracies of SCM-MODIS were greater than 90% in both the plain and the mountainous regions. But, the SCM-MODIS in the plain region (L1) showed a slight higher accuracy than that in the mountainous regions (L2 and L3), being consistent with the conclusion obtained from the indirect validation method.

Fig. 3

Figure Option ViewDownloadNew Window

Fig. 3 A comparison of snow-cover maps obtained through visual interpretation of Landsat ETM+ (2nd column) and through cloud removal of MODIS 8-day snow-cover products (3rd column). Left to right: false color composite image derived from Landsat ETM+ using bands 5, 4, and 3 (a-1, a-2, a-3), snow-cover maps obtained through visual interpretation of Landsat ETM+ (b-1, b-2, b-3), and snow-cover maps obtained through the removal of clouds from the MODIS 8-day snow-cover products (c-1, c-2, c-3).

Table 3

Table 3

Table 3 Results comparison between SCM-ETM+ and SCM-MODIS Scene Csnow(km2) Cland(km2) Tsnow(km2) Tland(km2) Tad(%) L1 2495.11 4530.26 2506.78 5058.77 92.86 L2 6941.32 1239.12 7026.56 1989.68 90.73 L3 2473.50 10, 059.53 2502.60 11, 124.70 91.97 Note: Csnow, common snow area; Cland, common land area; Tsnow, snow area and land (i.e., snow-free) area in SCM-ETM+; Tad, overall accuracy of the cloud-free snow-cover dataset under the direct validation method.

Table 3 Results comparison between SCM-ETM+ and SCM-MODIS

We generated a mean SCF map for Central Asia during the study period (Fig. 4a) based on the cloud-free snow-cover dataset of Central Asia. In general, the SCF in the plain regions was lower than that in the mountainous regions; the perennial snow cover accounted for only approximately 0.34% of the study area and mainly distributed in the Pamir Plateau, the Alay Range (Fig. 1), and the middle/eastern Tianshan Mountains. In general, The SCF is dependent on altitude.

Fig. 4

	Figure Option ViewDownloadNew Window
	Fig. 4 Spatial distribution of the annual mean snow-cover frequency (SCF) derived from the snow-cover frequencies from 2000/2001 to 2014/2015 (a), annual mean SCF changes with latitude in the plain regions (Profile line 1, b) and mountainous regions (Profile line 2, (c)), and the linear relation between annual mean SCF and altitude (d)

Profile lines (Fig. 4), which were inspired by Dietz et al. (2013), can be used to analyse the SCF variations along a geographical gradient. Profile line 1 (Fig. 4b) exhibits a typical distribution of the SCF with a north-south gradient, revealing that the SCF decreased with decreasing latitude at a rate of approximately 0.03 per degree in latitude (R²=0.97, P< 0.001) between 39° N and 54° N. However, a slight increasing trend was observed with a rate of 0.015 per degree in latitude (R²=0.85, P< 0.001) between 35° N and 39° N.

Profile line 2 (Fig. 4c) is drawn across mountainous regions (Tianshan Mountains and Alay Range) and the intermontane Fergana Valley. Along this line (i.e., line 2), the SCF did not vary with latitude. Rather, it increased with increasing altitude (Yatagai et al., 2012).

For the entire Central Asia, a strong correlation between the mean SCF and altitude was found with a rate of approximately 0.12 per kilometer (R²=0.81, P< 0.001; Fig. 4d) and the strongest SCF’ s dependency on altitude occurs between 4000 and 6000 m in the Kunlun Mountains and the eastern Pamir Plateau.

Figure 5 presents the temporal variations of the mean SCR and SCF. The SCR exhibits a slight and insignificant increasing trend (R²=0.07, P=0.36; Fig. 5a), suggesting that the mean snow-cover extent in Central Asia has not changed significantly during the study period. The SCF also showed a slight insignificant increasing trend (R²=0.12, P=0.21; Fig. 5b), also suggesting that the SCF has not changed significantly during the study period.

Fig. 5

	Figure Option ViewDownloadNew Window
	Fig. 5 Changes in annual mean snow-cover rate (SCR; a) and annual mean SCF (b) over Central Asia from 2000/2001 to 2014/2015

Figure 6 illustrates the spatial distribution of the slope of SCF (Fig. 6a) and the corresponding P values (Fig. 6b). It can be seen that the majority of Central Asia (approximately 65% of the total pixels) experienced an increase of SCF. But, only 4.3% of these SCF-increasing pixels is statistically significance and they were distributed in western and central-southern Kazakhstan and northwestern and southern Turkmenistan. By contrast, about 35% of Central Asia underwent a decrease in SCF, of which 5.2% of the area had a significant decrease. The significant decrease mainly occurred in western Tarim Basin and southern and eastern Xinjiang of China. The SCF tended to increase in most of the plain areas (e.g., Kazakhstan, Uzbekistan and Turkmenistan), although the increasing trend was not statistically significant.

Fig. 6

	Figure Option ViewDownloadNew Window
	Fig. 6 Spatial change trends in the SCF (a) and the statistical Pvalues associated with the spatial change trends in the SCF (b) in Central Asia during the period 2000/2001-2014/2015

EOF analysis was employed to detect the spatial and temporal variations of snow cover using a small number of variables (Table 4). The first two EOF modes (EOF-1 and EOF-2) represented the predominant change patterns and contributed 39.95% of the spatial anomaly variability. These first two modes were thus selected to analyze the anomaly variations of snow cover to reveal the main physical processes behind the snow-cover variations in Central Asia.

The first mode EOF-1 (Fig. 7a) exhibits the SCF variation with positive or negative anomalies in 27.65% of the total variance. Positive anomalies account for 77.5% of the total pixels, mostly in western and southern Kazakhstan. The remaining 22.5% reflects the negative anomalies predominately distributed in the Pamir Plateau, Tianshan Mountains, Kunlun Mountains, Altay Mountains, and northeastern Xinjiang in China. The positive anomalies in SCF variability exceed 2.0 for some pixels, whereas the negative anomalies are mostly lower than 1.0. The corresponding orthogonal time coefficients (Fig. 7b) tend to increase although insignificantly (R²=0.14, P=0.16). However, it is noteworthy that the SCF anomaly variability is significant in the years of 2001/2002 and 2011/2012 at a significance level of 0.05, being correspondent to the results suggested in Figure 5b. All these analyses confirm that the majority of Central Asia experienced an increase in SCF during the study period. However, in western and southern Kazakhstan, about 35% of the area underwent a decrease in SCF. It mainly occurred in the mountainous areas.

Table 4

Table 4

Table 4 Eigenvalues λ , variance contributions R, and the results of the North test δ λ /Δ λ for the EOF analyses EOF analysis EOF-1 EOF-2 EOF-3 EOF-4 EOF-5 EOF-6 λ (10-4) 9.43 4.20 3.17 2.29 2.07 2.04 R (%) 27.65 12.30 9.31 6.71 6.08 5.99 δ λ /Δ λ - 0.29 1.14 0.15 0.17 - Note: EOF, empirical orthogonal function; numbers with the acronym EOF represent EOF modes 1-6.

Table 4 Eigenvalues λ , variance contributions R, and the results of the North test δ λ /Δ λ for the EOF analyses

EOF-2, which contributes 12.30% of the total variance, displays a strong contrast between high latitude/mountainous regions and low-latitude/basin regions in SCF. That is, positive anomalies of SCF occurred in low-latitude regions including southern Kazakhstan, Uzbekistan, Turkmenistan, and Tarim Basin (Fig. 7c), whereas negative anomalies were found in northern Kazakhstan and most mountainous areas with high latitude or altitude. In addition, the fluctuations of the time coefficients corresponding to EOF-2 were close to the central value of 0 (Fig. 7d), suggesting that the SCF anomalies revealed by EOF-2 were representative.

Fig. 7

	Figure Option ViewDownloadNew Window
	Fig. 7 Empirical orthogonal function model 1 (EOF-1) of the annual SCF anomalies (a) and the corresponding time series of EOF-1 (b) and EOF-2 of the annual SCF anomalies (c) and the corresponding time series of EOF-2 (d)

Since the spatial distribution of snow cover also correlates with the ecological regions, the ecological regions in Central Asia are further divided into plain and mountainous regions (Fig. 8). The snow-cover variations in the different plain regions (Fig. 9a) indicated that they were robustly and positively correlated with latitude. According to the characteristics of snow-cover variations in the plain regions, the seven ecological regions located in the plain regions were divided into the low-latitude regions (consisting of the BK semi-desert and the SCA desert), the mid-latitude regions (consisting of the NCA desert and the Kazakh semi-desert), and the high-latitude regions (consisting of the Kazakh steppe, Kazakh upland, and Kazakh forest steppe) to further investigate the characteristics of snow cover. In winter (Julian day 337 to 57), the SCR (snow-cover rate) in the low-latitude, mid-latitude, and high-latitude regions ranges from 3.41% to 26.84%, from 38.29% to 73.54%, and from 83.15% to 99.60%, respectively. At the same time, the snow-cover duration prolonged with the increasing latitude in the plain areas. The SCR variation in the mountainous regions (Fig. 9b) also shows a similar latitudinal dependency.

Fig. 8

	Figure Option ViewDownloadNew Window
	Fig. 8 Spatial distribution of the key ecological regions in Central Asia

Fig. 9

	Figure Option ViewDownloadNew Window
	Fig. 9 Dynamic variations of the mean SCR in the different ecological regions during the period 2000-2015 for the plain areas (a) and the mountainous areas (b)

To further investigate the characteristics of snow cover in the different ecological regions, we computed the annual mean SCR (Fig. 10). The mean SCR variations in the low-latitude, mid-latitude, and high-latitude regions of the plain areas during the past 15 years indicated that SCR increases with the increasing latitude. The inter-annual variation of SCR in the mountainous areas suggests that Altay Mountains had the maximum SCR during the period 2000-2015, followed by Pamir Plateau, Tianshan Mountains, and Kunlun Mountains, suggesting that SCR increases with increasing latitude in the mountainous areas.

Fig. 10

	Figure Option ViewDownloadNew Window
	Fig. 10 Inter-annual variations of SCR in the different ecological regions of the plain areas (a) and the mountainous areas (b)