Snow Cover Variations and Controlling Factors at Upper Heihe River Basin, Northwestern China

Abstract

Snow is an important water resource and greatly influences water availability in the downstream areas. In this study, snow cover variations of the Upper Heihe River Basin (UHRB) during hydrological years (HY) 2003–2013 (September through August) is examined using the flexible multiday-combined MODIS snow cover products. Spatial distribution and pattern of snow cover from year to year for the basin is found to be relatively stable, with maximum snow cover area (SCA) and snow cover days occurring in HY2004, HY2008 and HY2012. A method, based on correlation coefficients between SCA and climate factors (mainly air temperature and precipitation), is presented to identify the threshold altitude that determines contributions of climate factors to SCA. A threshold altitude of 3650 ± 150 m is found for the UHRB, where below this altitude, both air temperature (T_air) and precipitation are negative factors on SCA, except in the winter season when both are positive factors. Above the threshold altitude, precipitation acts as a positive factor except in summer, while T_air is a negative factor except in autumn. Overall, T_air is the primary controlling factor on SCA below the threshold altitude, while precipitation is the primary controlling factor on SCA above the threshold altitude.

1. Introduction

As an important part of the hydrologic cycle, snow cover extent and variations through time are crucial to hydrologic modeling, water balance, and climate change [1]. Studying snow cover variation is of significant importance in monitoring and maintaining water management for ecosystem processes and irrigation practices [2]. For example, in Western China, snow cover is a very important water resource, since water from melted snow forms the headwaters for several of the largest rivers in Asia, including the Yangtze River and Yellow River. In semi-arid regions where precipitation is rare, snowmelt water is dominant among all water sources [3]. Consequently, snow can significantly influence agriculture and animal husbandry; for example, heavy snow could bring serious impact to animal husbandry [4,5]; spring drought resulting from snow scarcity can cause lower harvests in agriculture [6]; and fluctuation of albedo in response to change in the snow cover can cause significant change in the soil temperature, and is unfavorable for crop growth [7]. More than one sixth of the global population relies on water from mountainous snowmelt [8,9].

As a typical inland river basin in an arid area of Western China, the Heihe River Basin (HRB) is characterized by its distinct cold and arid landscapes: glaciers, frozen soils, alpine meadows, forests, irrigated crops, and desert, distributed among the upstream area, the middle stream area, and the downstream area [10,11,12,13,14]. The Upper HRB (UHRB) is rich in snow water resources, which is the majority of water resources for the entire HRB. A slight change in snow cover in the UHRB would greatly affect the water budget of the Heihe River, thus bringing a chain-reaction to agriculture and vegetation in the whole basin area [10]. Due to limited meteorological stations, in-situ snow cover observation data are insufficient to reflect the spatial and temporal distributions of snow cover in the UHRB [15]. The biggest advantage of remote sensing technologies is that they can provide information where ground stations are not available or ground station density is sparse.

The most recent and advanced snow cover product is produced by the MODerate resolution Imaging Spectroradiometer (MODIS) flown on the Earth Observing System (EOS) Terra and Aqua platforms, launched in 1999 and 2002, respectively [16]. Maurer et al. [17] reported that one of the main difficulties in estimating snow cover area (SCA) from visible remote sensing sensors (such as GOES and AVHRR) is cloud. In a validation study of MODIS snow cover images over Austria, Parajka and Blöschl [18], found that 63% of the areas were cloud-covered, although on cloud-free days, MODIS snow mapping accuracy was about 95%. A multi-day flexible combination of MODIS Terra/Aqua snow cover products was developed to effectively reduce cloud cover and was widely applied in many areas, such as in the USA and China [19,20,21,22].

It has been suggested that temporal and spatial distributions of snow cover are primarily influenced by air temperature (T_air) and precipitation [23]. Low T_air (below zero) and precipitation are the two major weather conditions required for snowfall. The SCA in mountainous regions is also highly correlated with altitude because temperature usually decreases with height (stratosphere is a special case). Topography determines the uplift of moist air that can trigger condensation and precipitation [24].

Comparing with areas at low altitude, high snow cover in mountain areas is not as vulnerable to atmospheric warming, if temperature remains sufficiently low to enable snowfall and keep snowpack from melting during winter and spring. Previous studies [25,26,27,28] reported different trends of SCA change at low and high altitude sites with climate warming. It is very important to know the threshold altitude that separate the different effects of climate factors on snow cover distribution in a mountainous range. Beniston (2012) [25] recently found a threshold altitude at approximately 1500–2000 m a.s.l. in the European Alps, without implementing an objective method for identifying the threshold value. Morán-Tejeda et al. [24] reported an altitude dependence between climate factors and snowpack variability in Switzerland by using a multivariate regression model. They found a threshold altitude of ~1400 m where temperature is the main explanatory variable for snow cover below 1400 m and precipitation is a better predictor of snowpack variability above 1400 m. They also indicated that as the climate warms, the threshold altitude would increase.

One purpose of this study is to examine if such a threshold altitude exists in the UHRB and what its value is, if it indeed exists. Our study finds that an altitude of 3650 ± 150 m in the UHRB defines the climate factor contributions and that such contributions fluctuate over four different seasons. In Section 2 and Section 3, the study area, data, and method used in this study are described. In Section 4, the temporal and spatial variations of snow cover, and their correlation with air temperature and precipitation are shown. In Section 5, we discuss different factors, including altitude, seasonality, and sublimation, in determining the correlations between temperature/precipitation and snow cover. Conclusions are given in Section 6.

2. Study Area

As the second largest inland river basin in China, the HRB has become an important research interest in the past decades due to its typical topography, climate, and environmental conditions [29]. The basin is located at the northeast edge of the Tibetan Plateau (TP). Owing to its unique natural environment, the HRB can serve as an ideal laboratory for studying hydrological processes in cold and arid regions [30].

Figure 1. Drainage basin of the upper Heihe River, Northeast edge of the Tibetan Plateau (TP) (shadow area of the inset map), with two weather stations: Qilian (2787 m) and Yeniugou (3320 m) denoted as black stars. The basin is segmented into four elevation zones. Zone a: 1708–3000 m, Zone b: 3000–3500 m, zone c: 3500–4000 m, and Zone d: 4000–4787 m. The 31 grid cells (25 km × 25 km) are extracted from APHRODITE precipitation data.

Figure 1. Drainage basin of the upper Heihe River, Northeast edge of the Tibetan Plateau (TP) (shadow area of the inset map), with two weather stations: Qilian (2787 m) and Yeniugou (3320 m) denoted as black stars. The basin is segmented into four elevation zones. Zone a: 1708–3000 m, Zone b: 3000–3500 m, zone c: 3500–4000 m, and Zone d: 4000–4787 m. The 31 grid cells (25 km × 25 km) are extracted from APHRODITE precipitation data.

The UHRB (Figure 1), ranging from 38.6°N to 39.1 °N and 98.7°E to 101.2°E, is about 10,009 km² area, with an altitude range from 1674 to 5108 m. Temporary snow usually exists under 2700 m; patchy snowpack is found in ranges from 2700 to 3400 m; and perennial snowpack exists above 3400 m. Alpine steppe and desert steppe exist under 2800 m; alpine meadow and bushes are in ranges from 2800 to 4000 m; bare ground exists above 4000 m [29]. As shown in Figure 1, the study area is segmented into four elevation zones: (a) 1708–3000 m; (b) 3000–3500 m; (c) 3500–4000 m; and (d) 4000–4787 m.

3. Data and Methodology

3.1. Data

3.1.1. SRTM DEM

The Shuttle Radar Topography Mission (SRTM) obtained elevation data in February 2000 on a near-global scale to generate the most complete high-resolution digital topographic database of Earth. The SRTM consisted of a specially modified radar system that flew onboard the Space Shuttle Endeavour during an 11-day mission. The surface elevations used in this study are from SRTM digital elevation model (DEM) data.

3.1.2. Snow Cover Data

Standard MODIS snow products from the Terra satellite have been available since February 2000. The products include daily and 8-day composite 500 m resolution swath and tile products (which include both binary and fractional snow cover (FSC) and snow albedo), and 0.05° resolution products on a climate-modeling grid (CMG) (which also include both binary and FSC). These snow products (from Collection 4 reprocessing) have been validated and are available to order through the National Snow and Ice Data Center [16]. The overall absolute accuracy of the well-studied 500 m resolution swath (MOD10 L2) and daily tile (MOD10A1) products (binary products) is ~93%, with variations by land-cover types and snow conditions [16]. The binary snow cover products from the MODIS instruments on both Terra and Aqua satellites of the NASA EOS are used in this study, including MODIS daily (MOD10A1, MYD10A1) and 8-day maximum composite (MOD10A2) from HY2003 to HY2013. Only one MODIS tile h25v05 is needed to cover the study area. Binary snow cover products (snow or no snow) are used in this study, which may, of course, cause some uncertainties. The MODIS factional snow cover products would be used for future study, although their accuracy is still a concern as they currently lack validation.

3.1.3. Temperature

Air temperature from two weather stations (Figure 1), Qilian station (2787 m) and Yeniugou station (3320 m), and Terra MODIS 8-day land surface temperature (LST) (MOD11A2) product are used to build empirical equations for producing air temperature for each grid cell for the study. For MODIS LST products, there is a published uncertainty of ±1 degree K [31]. Yu et al. (2011) validated the MODIS 1 km daily LST products using ground based longwave radiation observations from four automatic weather stations in the HRB and found that the mean bias was less than 1 K and the mean absolute error was about 2–3 K [32].

3.1.4. Precipitation Data

Asian Precipitation Highly-Resolved Observational Data Integration Toward Evaluation (APHRODITE) daily gridded precipitation (25 km cell size) [33] is the only long-term (1951 onward) continental-scale daily product for Asia including the Himalayas, South and Southeast Asia and mountainous areas in the Middle East [33,34]. It is based on gauge observations from the Global Telecommunication System, precompiled datasets, such as the Global Historical Climatology Network, the National Climate Data Center, the Food and Agriculture Organization of the United Nations, and the national hydro-meteorological services. The spatial interpolation scheme takes into account the effect of mountain ranges by giving a high weight to gauges on slopes inclined to the target location and a low weight to gauges on the leeward side behind a mountain ridge. Based on a comparison between co-located 17 gauges and APHRODITE precipitation in the HRB from 1961 to 2007 (not shown), the correlation coefficients of the two datasets range from 0.79 to 0.96 (p < 0.01), with the root mean squared error from 0.33 to 1.10 mm/day. This indicates the overall good quality of the APHRODITE precipitation product.

3.2. Methods

3.2.1. Deviation of Air Temperature from LST Data

Since air temperature (T_air) is a better variable than LST to snow cover [31,32,35], we correlate T_air measured at the two weather stations (Qilian and Yeniugou stations) and co-located LST derived from MODIS in two different cases, snow cover and snow free, determined by the co-located and concurrent MODIS snow cover product. It is assumed that, in a given elevation zone, the T_air and LST follow a linear relationship [31]. Below are the four steps used to derive air temperature from MODIS LST data.

(i)

MOD10A2 8-day snow products from HY2002 to HY2007 are used to separate snow and no snow periods for the two station pixels.

(ii)

MOD11A2 8-day LST data from HY2002 to HY2007 are extracted for the two stations and are then separated into snow and no snow periods based on the MOD10A2 snow product analysis.

(iii)

The regression equations between 8-day mean air temperatures measured at a station and 8-day MODIS LST at the station pixel for the same period HY2002 to HY2007, for snow and no snow cases, are derived, as shown in Figure 2. These equations (unit K) are: (1) y = 148.97 + 0.4209x; (2) y = 73.882 + 0.6955x; (3) y = 107.6 + 0.559x; and (4) y = 53.78 + 0.757x. There are many more data in the no-snow case than the snow case for both stations, with R² values of the no-snow case larger than those of the snow case. The R² values at Yeniugou station are larger than R² values at Qilian station. This might indicate that the T_air is correlated better with LST in higher altitude (Yeniugou station) than at lower altitude (Qilian station). The difference of land-surface and air temperatures dependence is related to the aerodynamic resistance to heat transfer [36,37]. At higher altitude, the air temperature is lower. Therefore the heat transfer is also lower, resulting in smaller differences between air temperature and land-surface temperature. This might be the reason that the correlation coefficient between air temperature and LST is higher at higher altitude than that at lower altitude.

(iv)

Simulate T_air of each cell based on the MODIS LST, using the above equations in Figure 2. The Equations (1) and (2) (snow and no snow cases) derived at the Qilian station (2787 m) are used to simulate T_air in the elevation zone (a) 1708–3000 m. Equations (3) and (4) (snow and no snow cases) derived at the Yeniugou station (3320 m) are used to simulate T_air for elevation zones (b) 3000–3500 m; (c) 3500–4000 m; and (d) 4000–4787 m. Using the equations derived from the Yeniugou station (in Zone b) for the zones (c) and (d) might cause some errors; however, the difference should be minor and should not impact the overall patterns of T_air change.

Figure 2. Relationships (Equations (1)–(4)) between T_air and LST for Qilian station for (a) snow and (b) no-snow observations and Yeniugou station for (c) snow and (d) no-snow observations from hydrological years 2002 to 2007, with linear regression lines and R² values also included.

Figure 2. Relationships (Equations (1)–(4)) between T_air and LST for Qilian station for (a) snow and (b) no-snow observations and Yeniugou station for (c) snow and (d) no-snow observations from hydrological years 2002 to 2007, with linear regression lines and R² values also included.

3.2.2. Gridding Partition

Since the resolution for APHRODITE (25 km × 25 km) is much coarser than that of snow product (500 m × 500 m) and LST (1 km × 1 km), 31 APHRODITE grid cells covering the whole UHRB are used to represent the altitude differences (Figure 1), in examining the interrelationships between altitude, precipitation, T_air and SCA. These cells are marked as P1 to P31 in Figure 1. The central altitude of each cell is treated as the altitude for the cell in this study.

3.2.3. Snow Cover Mapping Methods

One of the most effective ways to maximize the snow cover retrieval is by minimizing cloud cover [18,19,38,39]. In this study, the multi-day flexible combination of MODIS Terra/Aqua snow cover products is used to remove cloud cover, with two thresholds, 10% in maximum cloud cover and 8 days of maximum composition, to control the multiday combinations [19,20,21,22].

Snow covered days (SCD) maps are then derived from the flexible multiday combinations of MODIS images (less than 10% cloud cover). Typically, snow covered days (SCD) is defined as the total number of days that a pixel is covered with snow in a hydrological year [40]. For each pixel, at a given day, a value is given depending on the snow cover condition, 1 as snow and 0 as no snow. The sum of values for a given pixel of an entire year is defined as the SCD of the pixel cell. The calculation equation is shown as:

$$SC{D}_j={\displaystyle \sum }_{i=1}^{365}Cel{l}_i\left(Cel{l}_i=0 or 1\right)$$

(5)

where SCD_j is the snow cover days in a pixel j and i is the ith day. Since there are still pixels with cloud cover, the maximum SCD of a pixel could be less than 365 days.

3.2.4. Correlation between Climate Factors and Snow Cover Area (SCA)

Four seasons, autumn (Septempber–November), winter (December–February), spring (March–May), and summer (June–August), are used to differentiate the seasonal effects. SCA and air temperature are spatially averaged in each grid cell (25 km × 25 km). Accumulative precipitation of each cell from daily APHRODITE product is calculated with an 8-day time interval to match with the SCA and T_air (8-day).

Pearson’s correlation coefficient between two variables is defined as the covariance of the two variables x and y divided by the product of their standard deviations as shown in the following equation.

$$R=\frac{\text{n}({\displaystyle \sum }\text{xy})-({\displaystyle \sum }\text{x})({\displaystyle \sum }\text{y})}{\sqrt{[\text{n}{\displaystyle \sum }{\text{x}}^2-{({\displaystyle \sum }\text{x})}^2][\text{n}{\displaystyle \sum }{\text{y}}^2-{(\text{y})}^2]}}$$

(6)

where x and y are two variables, R is the correlation coefficient. Since correlation between SCA and T_air/P are determined in the same time period from HY2003 to HY2013, n equals to 11 here. By plotting the correlation coefficients against the altitudes of individual grid cells in different seasons, we are able to assess the correlations between climate factors and SCA, in terms of seasons and altitudes.

4. Results

4.1. Characteristics of Temporal Variations of Snow Cover Areas (SCA)

Snow cover parameters within a hydrological year include the timing of snow cover onset/disappearance in a snow cycle, SCA, and maximum SCA. Figure 3 presents the snow cover time series for the study area. In most of the hydrological years, there are two major snow periods, autumn and spring seasons, while in between (winter months), there is a period of low relative SCA. For example, in HY 2005 and HY 2011, the SCA even completely dropped down below 40% during the winter months. The minimum SCA values in summer would be the glaciers and perennial snow cover of the basin.

Figure 3. Time series of snow cover area (SCA) percentage (%) over the basin area from flexible multiday combined MODIS images, in the Upper Heihe River Basin during HY2003–HY2013.

Figure 3. Time series of snow cover area (SCA) percentage (%) over the basin area from flexible multiday combined MODIS images, in the Upper Heihe River Basin during HY2003–HY2013.

We define a SCA peak as the percentage of SCA that is greater than 40% of the basin area, as seen in Figure 3 and summarized in Table 1. The HY2011 shows the smallest number (3) of peaks (larger than 40%), while in the HY2004 there are a total of 13 peaks. The highest peak of 85% appeared in HY 2008, the greatest SCA during the 11-year period. HY 2004, HY 2008, and HY 2012 are the three hydrological years with relatively higher SCA in the past decade. In HY2004, the largest SCA peak appeared in autumn, while in HY 2008, the largest SCA peak appeared in winter. In HY2012, there are two peaks in November and February, respectively.

Table 1. Number of Peaks with Snow Cover Percentages of More Than 40% in Autumn, Winter Spring, and Summer, During Hydrological Years 2003–2013.

Hydrological Year	Number of Peaks
	Autumn	Winter	Spring	Summer
2003	2 (55%, 60%)	3 (45%, 45%, 40%)	5 (40%, 40%, 40%, 65%, 60%)	0
2004	6 (40%, 50%, 60%, 80%, 55%, 80%)	2 (55%, 50%)	5 (75%, 70%, 60%, 60%, 40%)	0
2005	5 (40%, 40%, 40%, 40%, 50%)	1 (50%)	4 (40%, 50%, 60%, 40%)	0
2006	3 (40%, 40%, 40%)	1 (40%)	2 (50%, 60%)	0
2007	4 (40%, 70%, 40%, 40%)	2 (55%, 40%)	3 (50%, 80%, 60%)	0
2008	4 (50%, 60%, 60%, 40%)	2 (75%, 85%)	5 (40%, 40%, 70%, 40%, 80%)	0
2009	4 (40%, 85%, 40%)	2 (55%, 65%)	3 (45%, 50%, 40%)	1 (40%)
2010	3 (50%, 40%, 70%)	0	3 (60%, 70%, 70%)	0
2011	2 (45%, 85%)	0	1 (80%)	0
2012	5 (50%, 60%, 55%, 45%, 40%)	1 (75%)	3 (55%, 50%, 40%)	0
2013	2 (40%, 40%)	1 (40%)	2 (60%, 60%)	0

Table 1. Number of Peaks with Snow Cover Percentages of More Than 40% in Autumn, Winter Spring, and Summer, During Hydrological Years 2003–2013.

Hydrological Year	Number of Peaks
	Autumn	Winter	Spring	Summer
2003	2 (55%, 60%)	3 (45%, 45%, 40%)	5 (40%, 40%, 40%, 65%, 60%)	0
2004	6 (40%, 50%, 60%, 80%, 55%, 80%)	2 (55%, 50%)	5 (75%, 70%, 60%, 60%, 40%)	0
2005	5 (40%, 40%, 40%, 40%, 50%)	1 (50%)	4 (40%, 50%, 60%, 40%)	0
2006	3 (40%, 40%, 40%)	1 (40%)	2 (50%, 60%)	0
2007	4 (40%, 70%, 40%, 40%)	2 (55%, 40%)	3 (50%, 80%, 60%)	0
2008	4 (50%, 60%, 60%, 40%)	2 (75%, 85%)	5 (40%, 40%, 70%, 40%, 80%)	0
2009	4 (40%, 85%, 40%)	2 (55%, 65%)	3 (45%, 50%, 40%)	1 (40%)
2010	3 (50%, 40%, 70%)	0	3 (60%, 70%, 70%)	0
2011	2 (45%, 85%)	0	1 (80%)	0
2012	5 (50%, 60%, 55%, 45%, 40%)	1 (75%)	3 (55%, 50%, 40%)	0
2013	2 (40%, 40%)	1 (40%)	2 (60%, 60%)	0

Figure 4a shows seasonal mean variations of SCA in a hydrological year. Autumn and spring are the two primary seasons with greater SCA, while summer has the lowest SCA. SCA in winter season varies a lot, although in most of the years, it is less than that of either autumn or spring or both seasons, with only the one exception that it is the highest in HY2008. There is a suggested 3–4-year cycle of annual SCA change. During a typical cycle, the annual SCA rises, drops and rises again. Taking the cycle during 2003 to 2006 as an example, there is an overall annual SCA increasing from 2003 to 2004 and then decreasing from 2004 to 2006. The same pattern appears from 2006 to 2010, in which the SCA increases from 2006 to 2008, and decreases from 2008 to 2010. However, with only 11 years of MODIS data, it is hard to reach any definitive conclusion on a linear time trend of the SCA for this area. Consequently, a longer time series of data are needed to reach more definitive conclusions.

The seasonal and annual accumulative precipitations of HY2003–HY2007 are shown in Figure 4b (APHRODITE precipitation data is not available after 2007). The seasonal precipitation is following the same trend: summer > autumn > spring > winter. Obviously, precipitation in winter is the lowest, while SCA in winter is not always the lowest. For example, it is even the highest in 2008. This is because the SCA in autumn can sometimes survive through the winter, but, of course, not all the SCA in autumn can survive through the winter. Otherwise, the SCA in winter would always be higher than that in autumn. That this is not the case, therefore, suggests, either sublimation or melting or both would happen in winter seasons in the study area. These possibilities are discussed further in the later sections.

Figure 4. Mean seasonal SCA in HY2003–HY2013 (a) and cumulative precipitation in HY2003–HY2007 (b) in the upper Heihe River Basin.

Figure 4. Mean seasonal SCA in HY2003–HY2013 (a) and cumulative precipitation in HY2003–HY2007 (b) in the upper Heihe River Basin.

4.2. Characteristics of Spatial and Temporal Variations of Snow Cover Days (SCD)

Figure 5 shows the spatial distribution of SCD in the UHRB during HY2003–HY2013. The SCD values range from 0 (blue) to 365 (red) for each pixel. The red areas indicate areas with longer days of snow cover. The variations in SCD are coincident with the variation in SCA. HY2004, HY2008 and HY2012 are three years with both the largest SCA and SCD. Specifically, in HY2004, the basin area in the northwest had an average SCD larger than 200 days, much higher than that in any other years. The areas with SCD greater than 60 are usually considered as stable snow-cover areas in China, they are the major containers of snow water resources, and where many large rivers originate. Comparing with the DEM in Figure 1, the areas with higher SCD (larger than 250 days) are consistent with the elevation zone of 4000 m and higher, corresponding well with the huge mountains such as Qilian Mountains. The most persistent SCAs (with SCD greater than 300 days) are concentrated in or near these huge mountaintops and ridges, mostly central-west area, northeast area, and southwest area. These three areas are glaciers and perennial snow cover areas of the basin.

Figure 5. Spatial distribution of snow cover days (SCD) in the upper HRB during HY2003–HY2013. The SCD values range from 0 (blue) to 365 (red) days.

Figure 5. Spatial distribution of snow cover days (SCD) in the upper HRB during HY2003–HY2013. The SCD values range from 0 (blue) to 365 (red) days.

4.3. Correlation between Air Temperature and Snow Cover

Figure 6 shows the air temperature from MODIS LST, based on the four equations shown in Figure 2, for the four elevation zones: (a) 1708–3000 m; (b) 3000–3500 m; (c) 3500–4000 m; and (d) 4000–4787 m. In a typical cycle of T_air in a year, from September to January, the T_air decreases to the lowest level and then increases to the highest level in July. Thereafter, the T_air can be relatively stable for several months before another cycle. Generally, the annual variation pattern at different elevations remains the same, while the absolute T_air is decreasing with the increasing elevation from zones (a) to (d). This suggests that the equations derived from the Yeniugou station (in Zone b) are reasonably good for T_air estimations for the zones (c) and (d).

Because Qilian (2787 m) and Yeniugou (3320 m) are at elevation zones (a) 1708–3000 m and (b) 3000–3500 m, respectively, the average SCA percentage from HY2003 to HY2013 in elevation zones (a) and (b) are compared with the corresponding mean air temperatures measured at the two stations (Figure 7), respectively. Similar to in Figure 3 for the entire area in each year, there are two SCA peaks in both zones (a) and (b), one before January and one after January. These two peaks split a hydrological year into four separated snow cover periods, with slight time differences due to the elevation and temperature differences in two elevation zones, while the separations between periods II and III of the two elevation zones are basically in the similar time, i.e., early January, usually the minimum T_air, or close to the minimum T_air. Overall, as temperature decreases from September to January, the SCA increases (Period I) and then decreases (Period II), with the peak SCA in early November (Zone a), or end of October (Zone b). As temperature increases from January to July, the SCA shows increase (Period III) and then decrease (Period IV), with the maximum SCA in late March (zone a) or early April (zone b).

Figure 6. The simulated annual cycle of T_air (k) for four different elevation zones for HY 2003 to HY 2013.

Figure 6. The simulated annual cycle of T_air (k) for four different elevation zones for HY 2003 to HY 2013.

The overall SCA in Zone b is larger than in Zone a due to the overall colder and higher elevation in Zone b than in Zone a. The resultant time period (Period II and Period III) is longer in Zone b than in the Zone a because of the overall higher altitude and lower temperature in Zone b. It is expected that the resultant time (Period II and Period III) in the Zone d is the longest of the four elevation zones.

Table 2 shows the Pearson correlation coefficients between the in-situ T_air and the corresponding SCA of the station pixel for the two stations in the four periods. The high negative correlations (at 99.9% confidence level) in Period I and Period IV indicate that T_air is a negative factor in snow formation and snow melting periods. This is easily understandable, since decreasing T_air in autumn opens the opportunity for snowfall and increasing T_air in spring causes snowmelt. However, in the snowpack duration period (Period II and Period III), significantly positive correlations (at 99.9% confidence level) between T_air and SCA may indicate that decreasing temperature (below 0 degrees) would cause no snowfall or decreasing snowfall, while the continuously strong sublimation would cause SCA decrease (Period II), since below 0 degrees T_air would not have snowmelt. Increasing temperature (still below 0 degrees), however, would potentially increase snowfall, resulting in SCA increase (Period III). The effects of T_air and precipitation on SCA against the elevation are examined in greater detail in the section below.

Figure 7. SCA (red) at elevation zones (a) 1708–3000 m and (b) 3000–3500 m and corresponding mean air temperature (blue) (T_air) at Qilian station (a) and Yeniugou (b) from HY2003 to HY2013. A year is divided into four periods I, II, III, and IV (purple dashed lines).

Figure 7. SCA (red) at elevation zones (a) 1708–3000 m and (b) 3000–3500 m and corresponding mean air temperature (blue) (T_air) at Qilian station (a) and Yeniugou (b) from HY2003 to HY2013. A year is divided into four periods I, II, III, and IV (purple dashed lines).

Table 2. Pearson correlation coefficients (p < 0.001) between average snow cover area and mean air temperature during the four different snow cover periods (indicated in Figure 7) at Qilian and Yeniugou stations from HY2003 to HY2013.

Snow Period	I	II	III	IV
Qilian	−0.72	0.77	0.54	−0.74
Yeniugou	−0.70	0.62	0.38	−0.52

Table 2. Pearson correlation coefficients (p < 0.001) between average snow cover area and mean air temperature during the four different snow cover periods (indicated in Figure 7) at Qilian and Yeniugou stations from HY2003 to HY2013.

Snow Period	I	II	III	IV
Qilian	−0.72	0.77	0.54	−0.74
Yeniugou	−0.70	0.62	0.38	−0.52

4.4. Seasonal and Altitude Effects on Correlations between Snow Cover and Climate Factors

Figure 8 shows the Pearson correlation coefficients between different climate factors and SCA from the 31 grid cells over the four seasons, with the zero correlation line and linear regression lines included to differentiate the effects of T_air and precipitation on snow cover as the altitude rises. The vertical line in the figure is the threshold altitude determined by a trial and error method. First, we select an initial threshold altitude around 3700 m, by plotting the regression lines of correlation coefficients between SCA and T_air and precipitation, in areas below and above the threshold altitude, respectively; two crossing points (one from regression lines between T_air and SCA and one from regression lines between precipitation and SCA) obtained could be exactly on, below, and/or above the initial threshold altitude. If both points are on the line, the initial threshold altitude is defined as the threshold altitude. Otherwise, a new threshold altitude is set, and the same process is repeated. The process only stops when both crossing points are on the same vertical line, which is then defined as the altitude threshold line for the season. As seen in Figure 8, the effects of T_air, precipitation, and elevation on SCA are quite different in the four seasons and threshold altitude also differs slightly from season to season, but within the range of 3650 ± 150 m.

Figure 8. Altitude effects on correlation coefficients between climate factors (air temperature (T_air) and precipitation (P)). T_air and P at low altitude (black and red) and high altitude (blue and pink) are labeled with different colors, respectively. The purple solid line defines the threshold altitude in each season. The cyan dashed line represents the correlation coefficient of zero.

Figure 8. Altitude effects on correlation coefficients between climate factors (air temperature (T_air) and precipitation (P)). T_air and P at low altitude (black and red) and high altitude (blue and pink) are labeled with different colors, respectively. The purple solid line defines the threshold altitude in each season. The cyan dashed line represents the correlation coefficient of zero.

In autumn, both T_air and precipitation are negative factors on SCA below the threshold altitude and their impacts decrease as altitude increasing and become to near zero at or near the threshold altitude. Both T_air and precipitation turn out to be mostly positive factors on SCA above the threshold altitude and their effects increase as altitude increases. T_air and precipitation are almost comparable in terms of correlation coefficient values with SCA at all altitude levels, although T_air is more negative below the threshold altitude and precipitation is more positive above the threshold altitude, i.e., T_air has a stronger effect on SCA in areas below the threshold altitude, while precipitation has a stronger effect on SCA in areas above the threshold altitude.

In winter, precipitation is a more important positive factor than T_air in areas below the threshold altitude, although the overall impact of precipitation on SCA shows a slightly decreasing tendency as altitude increases, while the T_air shows a fast increasing tendency as altitude increases. In areas above the threshold altitude, the effect of precipitation on SCA is attenuating to zero correlation as altitude increases, although it is still a positive factor on SCA. T_air, however, becomes mostly a negative factor in areas above the threshold altitude. It is noticeable that the blue line for the winter panel is not applicable since during winter the correlations between air temperature and SCA at high altitude vary a lot and there is no clear change trend.

In spring, both T_air and precipitation are negative factors on SCA below the threshold altitude, although T_air is the primary factor on SCA and its impact becomes more important at and near the threshold altitude. In areas above the threshold altitude, precipitation, however, becomes a positive factor and the primary factor on SCA, and its impact increases as altitude increases. In contrast, the impact of T_air on SCA above the threshold altitude is reducing as altitude increases.

In summer, both precipitation and T_air are negative factors on SCA at all altitudes, except that T_air seems to have no correlation with SCA below the threshold altitude.

Since both the Qilian and Yeniugou stations are located below the threshold altitude, 3650 ± 150 m, it is clear that the case shown in Figure 7 and Table 2 is actually a subset of the bigger picture discussed here in Figure 8. It is also clear that Period I corresponds to the autumn season, periods II and III to the winter season, and period IV to the spring and summer seasons. Although the T_air is a positive factor on SCA in the periods II and III, precipitation is a much important factor (also positive) on SCA in these altitudes (below the threshold altitude), as indicated in Figure 8.

In summary, below the threshold altitude, both T_air and precipitation are negative factors on SCA, except in the winter season when both are positive factors. Above the threshold altitude, precipitation acts as a positive factor except in summer when it is a negative factor, while T_air is a negative factor except in autumn when it is a positive factor. In the text below, we try to provide some physical explanations of these phenomena.

In areas below the threshold altitude, it is easily understandable for both T_air and precipitation as negative factors for SCA with T_air as the primary factor of the two for the autumn and spring seasons, since decreasing T_air enables snowfall in autumn and increasing T_air causes snowmelt in spring. There is almost no relation between T_air and SCA in summer, since there is not much snow cover left in summer after the spring melt. In winter season, precipitation as a positive factor is also understandable since more precipitation results in higher SCA, and vice versa. However, T_air as a positive factor on SCA in winter season is more complicated. As discussed earlier (for Figure 4 and Figure 7), it is suggested that the sublimation could be the direct and primary reason to reduce SCA as T_air decreases, since decreasing temperature (at T_air below 0 degrees) could not cause snowmelt, while the continuously strong sublimation could cause decreasing SCA. As indicated from Liston and Sturm (1998) [41] and Pomeroy et al. (1998) [42], sublimation is stronger at higher altitudes with extremely low temperature, low relative humidity, more dry winds, and less air pressure. Increasing T_air in winter season would potentially increase precipitation (snowfall), resulting in higher SCA. This effect is therefore understandable.

In areas above the threshold altitude, precipitation is the primary factor and positive factor for SCA in all seasons except in summer. Considering that T_air would be mostly below 0 degrees at high altitudes (i.e., above 3650 ± 150 m) in spring, autumn, and winter seasons, more precipitation generally means more snowfall then. In summer time, however, precipitation is a negative factor, as warmer rain would wash away any remaining snow cover. It is also understandable that T_air is a negative factor for SCA in winter and spring seasons, since T_air is already far below 0 degree in such high altitude area and increasing T_air would decrease snowfall, or vice versa. It is also understandable that T_air is a negative factor for SCA in summer since increasing T_air would speed up the remaining snow cover melt in summer in higher altitudes. The positive role of T_air playing on SCA in autumn season would be explained from two aspects. First, as seen in Figure 8, although both T_air and precipitation have a positive role on SCA, the precipitation is the primary factor; second, decreasing T_air might correlate with decreasing snowfall, while continuously strong sublimation in such high altitudes would reduce the SCA as mentioned before.

5. Discussion

In this study, we use correlation analysis to explain the different roles of temperature and precipitation on snow cover variability in a large scale in UHRB. Using gridded precipitation and remotely sensed snow cover and LST (converted to T_air) data, we are able to perform multivariate statistical analysis on a wide range of altitudes. Our results show that altitude significantly determines the ability of these two climate factors to influence the snow cover formation, melting, and duration. On the other hand, in different seasons, the role of climate factors varies significantly and the threshold altitude could fluctuate in a range of 300 m. As a whole, the interplay of altitude and season affects the precise placement of the threshold altitude at UHRB.

5.1. Altitude Factor

Based on our findings, in areas below the threshold altitude, the relationship between T_air and snow cover is always negative, except in the winter season; in areas above the threshold altitude, the relationship between precipitation and snow cover is always positive, except in the summer season This result is consistent with that found by Morán-Tejeda et al. [24] that temperature is a controlling factor for snow cover below the threshold altitude, while precipitation is the controlling factor for snow cover above the threshold altitude. However, our results provide much more detailed information. For example, temperature can become a positive factor in winter season (below the threshold altitude) and autumn season (above the threshold altitude) when sublimation in this low temperature, high wind, and low pressure environment could be a major factor for SCA decrease. Unlike T_air, precipitation can function either as a water source in the snow formation period or as a washing away force in the snow melting period. Below the threshold altitude, when temperature is above the freezing point, precipitation acts a negative factor to melt snow. Although the dependence of snowpack on temperature/precipitation has been suggested by indirect observations of Scherrer [28], Beniston [43] and Marty [44], and subsequently quantified by Morán-Tejeda et al. [24], our results indicate different trends of linear relationships between snow cover and climate factors below and above the threshold altitude. In the paper by Morán-Tejeda et al. [24], the correlation coefficients between climate factors and snow depth at different altitude levels follow a linear trend along the whole altitude level (same slope). However, in our results, the slopes of the trend lines of correlation coefficients along altitude are different below and above the threshold altitude. The use of different linear trends allows us to estimate the threshold altitudes where either precipitation or temperature becomes a better predictor of snow cover variability than the other.

5.2. Seasonal Factor

The threshold altitude to distinguish the different roles of climate factors on snow cover is found to seasonally fluctuate in the study area, i.e., 3700 m, 3670 m, 3500 m, and 3800 m, for autumn, winter, spring, and summer, respectively. This is the first time that such detailed study is conducted and that such different threshold altitudes in different seasons are found. The reason for such a difference is not clear. One of the reasons could be related to the coarse resolution of data used in the study. Since the only precipitation data available for the study is in 25 km resolution, other parameters, SCA and T_air are also averaged to the same cell size; the central altitude of each 25 km cell is used to represent the altitude of each cell. Nevertheless, this study provides the basic and significant information for further studies on this topic. It is also expected that such threshold altitude will be changed as global warming continues.

5.3. Sublimation Factor

Sublimation, involving complex mass and energy exchanges, has been recognized as an important hydrological process in high altitudes regions [45]. The UHRB is a typical inland area in China with arid and cold climate where sublimation occurs more readily [41,42,45,46,47]. There is a general agreement that higher sublimation and melt rates are driven by greater exchanges of latent and sensible heat, and higher net radiation [48]. During the process of sublimation, a lot of energy is needed to vaporize the snow and/or ice (720 calories per one gram of ice), thus the T_air also decreases during the large sublimation event, which takes away energy from the environment including radiation and heat [49].

As mentioned earlier, we suggest that sublimation would be the major reason for SCA decreasing in winter season in the study area, since snowmelt is not possible due to below zero degrees T_air during the winter period. Li et al. (2009 and 2014) [50,51] studied the sublimation behavior in UHRB. Based on their ground-based measurement data in Dadongshu Station, the sublimation was 1.65–1.90 times the precipitation in January and February, and became less than the precipitation thereafter (Table 3) [51]. Their studies clearly support our interpretation that sublimation is the reason of decreasing SCA in winter of the UHRB since no surface runoff before end of March was observed [50,51,52]. Previous studies also indicated that large amount of snow could be removed by sublimation given a high wind speed [53,54]. Wind speed in the Tibetan Plateau could be up to 16.9 m per second at the altitude above 4100 m [55]. As an important factor in sublimation, high-speed wind can also blow the snow away and change the distribution of snow cover in the local area.

Table 3. The monthly precipitation and sublimation measured in Dadongshu Station in 2008 (modified from [51]).

Month	Sublimation (mm)	Precipitation (mm)	Ratio (S/P)
January	59.9	31.6	1.90
February	52.7	32.0	1.65
March	24.5	34.1	0.72
1 April–12 May	3.7	106.9	0.03

Table 3. The monthly precipitation and sublimation measured in Dadongshu Station in 2008 (modified from [51]).

Month	Sublimation (mm)	Precipitation (mm)	Ratio (S/P)
January	59.9	31.6	1.90
February	52.7	32.0	1.65
March	24.5	34.1	0.72
1 April–12 May	3.7	106.9	0.03

5.4. Regional Factor

The average threshold altitude 3650 m in UHRB is much higher than that found in Switzerland (1400 m) [24] for several reasons. First, latitude and elevation levels in UHRB area of China (latitude: 38.9 °N; altitude: 1674 m to 5108 m) are quite different from those in Switzerland (latitude: 47.0 °N; altitude: 316 m to 2690 m). The higher elevation and more arid and cold climate in UHRB area also facilitate the sublimation phenomenon, which is one major factor in determining the snow cover variations in UHRB. However, no sublimation effect is observed or considered in the Switzerland snow area [24]. Second, the accuracy of remote sensing data in this study is obvious lower than that of the ground station data used in the Switzerland study. Especially, the resolution of precipitation in this study is as low as 25 km × 25 km, thus only 31 grid sites are studied in our area. This resolution could bring non-negligible influence to the results. In order to confirm our findings, more accurate climate data and a larger amount of test sites are needed to perform more careful statistical analysis. More study areas are needed to reveal the threshold altitudes at different locations.

6. Conclusions

The flexible multiday combination with two thresholds of maximum cloud percentage (10%) and maximum composite days (8) removes cloud coverage significantly, allowing efficient and accurate mapping of snow cover. The resulting mean cloud percentage of 4.06% in UHRB is much less than the daily products of 32.8%. The resulted SCA data are analyzed in detail and are used to examine the snow covered days (SCD) change for the periods of 2003 to 2013. Although there are general differences in snow cover onset and melt from year to year and from season to season, the spatial distribution and pattern of SCD from year to year for the basin seem relatively stable, with the HY2004, HY2008 and HY2012 having the maximum SCA and SCD.

The resulting SCA data, together with temperature, precipitation, and altitude, are examined in detail for their interrelations. It is found that an altitude of 3650 ± (150) m defines the contributions of climate factors, with seasonal fluctuation. During autumn, the snowfall is controlled by both T_air and precipitation negatively below the threshold altitude and positively above the threshold altitude. Winter is the snowpack duration season; the snow cover is mainly controlled by precipitation positively and is also controlled by sublimation effects that decrease the SCA. Spring is the snowmelt season, thus snowfall is primarily controlled by T_air negatively below the threshold altitude and positively controlled by precipitation above the threshold altitude. Seasonal snow cannot survive in summer season, when both T_air and precipitation are negative factors for snow cover.

The method presented in the study to locate the threshold altitude based on correlations between snow cover and climate factors would be very useful for long-term monitoring of threshold altitude level on a large scale. This method can overcome the ubiquitous drawback in snow cover research for lack of sufficient meteorological stations in high mountain areas. Moreover, it might provide alternative understanding of influence of global climate change to snow cover variations. Threshold altitude can be used as an input for climate models as a basis for future predictions. The threshold altitude is believed to change over time as the climate continues to warm. In this study, the only period studied is HY2003 to HY2013. Using the HY2003 to HY2013 as a baseline, we should be able to observe the trend of threshold altitude change over time and potentially be able to predict the threshold altitude in future decades.