Impact Analysis of Climate Change on Snow over a Complex Mountainous Region Using Weather Research and Forecast Model (WRF) Simulation and Moderate Resolution Imaging Spectroradiometer Data (MODIS)-Terra Fractional Snow Cover Products

Abstract

Climate change has a complex effect on snow at the regional scale. The change in snow patterns under climate change remains unknown for certain regions. Here, we used high spatiotemporal resolution snow-related variables simulated by a weather research and forecast model (WRF) including snowfall, snow water equivalent and snow depth along with fractional snow cover (FSC) data extracted from Moderate Resolution Imaging Spectroradiometer Data (MODIS)-Terra to evaluate the effects of climate change on snow over the Heihe River Basin (HRB), a typical inland river basin in arid northwestern China from 2000 to 2013. We utilized Empirical Orthogonal Function (EOF) analysis and Mann-Kendall/Theil-Sen trend analysis to evaluate the results. The results are as follows: (1) FSC, snow water equivalent, and snow depth across the entire HRB region decreased, especially at elevations over 4500 m; however, snowfall increased at mid-altitude ranges in the upstream area of the HRB. (2) Total snowfall also increased in the upstream area of the HRB; however, the number of snowfall days decreased. Therefore, the number of extreme snow events in the upstream area of the HRB may have increased. (3) Snowfall over the downstream area of the HRB decreased. Thus, ground stations, WRF simulations and remote sensing products can be used to effectively explore the effect of climate change on snow at the watershed scale.

1. Introduction

According to the Intergovernmental Panel on Climate Change (IPCC) Fifth Assessment Report, historical temperature records show that atmospheric temperatures increased up to 0.83 degrees Celsius from 1880 to 2012 [1]. Temperature has risen unequivocally over the last 100-plus years, and researchers have theorized that climate change has dramatically altered snow levels because water changes easily from solid to liquid in response to minor temperature changes [1,2].

Climate change is expected to accelerate the hydrologic cycle: as temperatures increase, a greater fraction of precipitation falls as rain rather than snow and more accumulated snowpack melts, resulting in greater runoff in the winter, less snowpack in the spring, and less runoff in the summer [2]. The enhanced hydrological cycle is also expected to increase snowfall amounts through increased moisture availability [3]. Both Burnett et al. [4] and Wright et al. [5] have indicated that lake-effect snowfall increases along with temperatures over the Great Lakes, although they used different methods to study the sensitivity of snowfall to climate change; the former used historical records of oxygen isotopes from the twentieth century, while the latter used the Advanced Research Weather Research and Forecasting Model (WRF-ARW). However, station observations of snow in spring at lower elevations in the Northern Hemisphere show decreasing snowfall trends. Snowfall records in the United States show that total snowfall over the 48 contiguous states has decreased by approximately 0.19 percent per year since 1930. Snowfall over southern China decreased obviously from 1979 to 2006 [6]. This decline has also been found across northeastern China since 1951 [7,8]. Snowfall records across the Tibetan Plateau of China indicate a more complicated pattern; snowfall in this region increased obviously from 1960 to 1998 [9] but has decreased since 2000 [10]. Across the Heihe River Basin (HRB), on the northeastern Tibetan Plateau, the change in snowfall is spatially heterogeneous [11]. Based on the above studies, the following research questions are relevant at the local scale: (1) Related to regional temperature increases, is there an elevation tipping point, below which snowfall decreases and above which snowfall increases? (2) For an entire mountainous region, how is the total volume of snowfall flux affected by rising temperature?

According to the National Oceanic and Atmospheric Administration, the increase in temperatures over the most recent 15-year period nearly equaled the increases that occurred over the entire last half of the 20th century. Therefore, analysis of the effects of climate change over the most recent years of the 21st century is commonplace. The regional hypotheses, questions and statements above provided the basis for an evaluation of long-term snowfall rates at high spatial and temporal resolutions using a WRF-ARW model of the mountainous region in the Heihe River Basin (HRB, shown in Figure 1) from 2000 to 2013. The spatial and temporal distribution of high-resolution snowfall data simulated by the weather research and forecast model (WRF) model is similar to observations across northeastern China [12,13]. Lee et al. [14] indicated that a heavy snowfall event generated by the WRF model matched the observations well.

Snow cover not only affects land-surface albedo and, consequently, the global water and energy balance, but is also considered an important indicator of climate change. The inverse relationship between hemisphere-averaged monthly mean snow cover and temperature fluctuations has been fine-tuned by many researchers since the 1990s [15,16,17]. Fractional snow cover (FSC) is the fraction of a grid cell covered by snow and is one of the key variables related to snow cover. Furthermore, FSC estimates have a higher spatial resolution than do snow water equivalent (SWE) products because the former are retrieved from data acquired by visible to near-IR sensors, while the latter are acquired using passive-microwave sensors [18]. Due to their higher spatial resolution, FSC measurements are more suitable for analyzing the regional effects of climate change on snow. The accuracy of Moderate Resolution Imaging Spectroradiometer Data (MODIS)-Terra seasonal FSC compared with in situ Chinese snow observations is approximately 9 0% across the Tibetan Plateau [19]. Wang et al. [20] also found that the MODIS/Terra FSC data is exceptionally accurate in areas with high snow depth in Northern Xinjiang, China.

Thus, MODIS/Terra FSC data can provide a reliability guarantee when diverse snow variables from WRF simulations are used to evaluate the impact of climate change on snow over a complex mountainous region, because precipitation (including rainfall and snowfall) values simulated by the WRF model are not as reliable as values generated for other near-surface atmospheric elements. Given the enormous potential to supplement ground observations and their ability to serve as a source of reference data, the inclusion of remote sensing data in our analysis was necessary to supplement the WRF simulation results and evaluate the effects of climate change on snow levels.

This study used both snowfall at an hourly 0.05° resolution simulated by the WRF-ARW model and FSC measurements at a monthly 0.05° resolution calculated from MODIS-Terra data to evaluate the performance of the climate model on snow across the HRB. We use the satellite-retrieved snow cover product supplemented by temperatures simulated by the WRF model to determine whether climate change has occurred across the HRB and to evaluate the relationship between snowfall and climate change across the HRB with the goal of answering the research questions listed above. Section 2 describes the research region, observational data, the WRF model’s configuration, the MODIS-Terra FSC product, the Empirical Orthogonal Function analysis and the Mann-Kendall/Theil-Sen Trend analysis. The results are presented in Section 3, followed by a discussion and conclusions in Section 4 and Section 5, respectively.

2. Materials and Methods

2.1. Research Region

The HRB (shown in Figure 1) is the second largest inland river area in the arid and semi-arid northwestern region of China. It is located at 37°45’–42°45’ N and 97°–102° E and covers an area of approximately 140,000 km² [21], encompassing the Qinghai, Gansu and Inner Mongolia provinces. Large topographical differences from upstream to downstream have resulted in elevation changes from approximately 5500 m to 1000 m, respectively; and distributed landscapes within the region are diverse, including: glaciers, frozen soils, alpine meadows, forest, irrigated crops, riparian ecosystems, deserts, and gobis. The HRB is comprised of three distinct regions: the mountainous zone, including the headwaters, tributaries, and up-stream reaches of the Heihe River; the alluvial-diluvial fan zone along the central, or midstream reaches, of the Heihe River; and the northern Alxa high plains region in the downstream area. The mountainous zone runs through the northern Qilian Mountains region, where the elevation ranges from 2000 to 5500 m and precipitation ranges from 250 mm in the low mountain zone to 500 mm in the high mountain zone [22]. The mountainous zone includes the headwaters of the river and shapes the strong drainage and water retention zones in the midstream zone. The midstream zone is located in the central part of the Hexi Corridor and is fairly flat, with elevations ranging from 1000 to 2000 m. This area comprises the alluvial-diluvial fan zone, where 91% of the agricultural areas in the entire basin are located; in addition, runoff from the mountainous zone is captured and used here. Mean annual precipitation in the midstream zone can be as low as 100 mm. The downstream portion of the HRB incorporates the northern Alxa high plains and has a mean elevation of less than 1000 m. The area is comprised of diluvial, alluvial and lacustrine plains and is dominated by the Gobi, desert, and bare land. Overall, the mountainous zone of the HRB is the cradle that supports the cultivated areas and civilization of the midstream zone, as well as uses in the downstream areas. The HRB has suffered from water and ecological stress since 1949, when rapid population growth in the region began [23]. Large areas of woodland and natural oases throughout the downstream zone have been threatened by reduced water flows caused by increasing water withdrawals for agriculture irrigation and municipal water supplies across the midstream zone [24].

2.2. Precipitation from the China Meteorological Administration Station Data

Distinguishing between solid and liquid precipitation, whether in liquid or solid form, is important to the sustainable development of the HRB. The launching of meteorological satellites has resulted in rich and varied liquid precipitation products extracted from remote sensing data over the years. This refined-grid liquid precipitation data can be used at regional scales by assimilating the remote sensing-based precipitation data into regional climate models (RCMs). However, corresponding snowfall-related products are scarce; thus, corresponding evaluations and analyses are also relatively poor. Nevertheless, the evaluation and analysis of solid precipitation cannot be ignored because solid precipitation is an essential supplement to the runoff of the entire arid and semi-arid region. Furthermore, snowfall has a complex relationship with climate change.

Long-term snow observation data are scarce, both in the HRB and throughout China. Daily precipitation data with type information are available only from 1951 to 1979 and must be acquired from the China Meteorological Administration (CMA). However, after 1979, no precipitation type information was recorded as a supplement to precipitation flux data. Ding et al. [25] explored the relationship between precipitation type and surface elevation and meteorological variables, including daily mean air temperature, daily mean relative humidity, daily mean surface pressure, daily total precipitation, and precipitation type based on observational data from 1951 to 1979 using 824 CMA stations distributed throughout China to rebuild long-term precipitation type data for all stations in China from 1979 to 2010.

As a supplement to the CMA station data, the Dabenying automatic meteorological station (99°53′ E; 38°16′ N; 2980.2 m) located in the Hulu watershed in the upstream region of the HRB was established in August 2010. A Chinese standard precipitation gauge with a Double-Fence International Reference has been set up in this station. The Dabenying station has recorded precipitation every half hour since September 2010. Precipitation phases (snow, rain and mixed) are distinguished using the CMA standard [26], and the volume of snowfall flux has been corrected through a comparison experiment in the Qilian Mountains [27]. Therefore, in addition to CMA station data, the Dabenying station data [28] are used to test snowfall simulation results from the WRF model. Figure 1 shows the locations of the CMA stations and the Dabenying site used in this study.

2.3. Snow-Related Variables Simulated by WRF-ARW and Model Configuration

In this study, snowfall was simulated by the WRF model. The WRF model is a new generation mesoscale numerical weather forecasting system [29,30,31,32] that serves both the operational and research communities. The WRF-ARW model uses a terrain-following hydrostatic pressure coordinate system that permits vertical grid stretching [33]. Arakawa C grid staggering [34] is used for horizontal discretization. As described in Pan et al. [35], after comparing 9 different microphysical parameters sensitive to analysis and 2 initial and boundary conditions with other fixed parameters, one set of physical configurations was selected as the most suitable for simulating rainfall in the HRB. Because both rainfall and snowfall belong to the same moisture variable, the WRF-ARW model physical configuration in this study is the same as that used by Pan et al. [35]: the Kessler scheme [36] was used as the Microphysics parameterization, the Kain-Fritsch Scheme [37] as the Cumulus parameterization, the Yonsei University scheme [38] as Planetary Boundary parameterization, the Noah LSM [39] as the land-surface parameterization, the Dudhia scheme [40] as the shortwave radiative parameterization, and the Rapid Radiative Transfer Model (RRTM) scheme [41] as the longwave radiative parameterization. In our research, WRF-ARW was initialized by real boundary conditions using the National Centers for Environmental Prediction’s (NCFP) Final (FNL) Operational Global Analysis data from the Global Forecast System, which has a resolution of 1° × 1°. Two-way nested computational domains of 44 × 56 × 27 and 120 × 130 × 27 grid points were set with horizontal resolutions of 0.25° and 0.05°, respectively. The center longitude is 100 E and the center latitude is 40.28 N. For more detailed information about the WRF-ARW configuration in this study, we refer readers to Pan et al. [35].

Snowfall simulated by the WRF model was evaluated using monthly observational data recorded at the CMA stations from 2000 to 2010 and daily observational data recorded in situ at the Dabenying station from 15 September 2010 to 14 September 2012. To match station positions using the WRF data, bilinear interpolation of the gridded snowfall data simulated by the WRF model was performed. Three precipitation types are recorded by the Dabenying and CMA stations: snowfall, rainfall and a mixture of both. However, no ratio of snowfall to rain exists in the mixed records, making it difficult to extract the actual volume of snowfall from the mix. Therefore, the ratios of mixed snow/rainfall yielded by the WRF model were used to extract the snowfall from the mixed records in these stations’ observational data. This ratio has not been tested or verified by current observational data; however, fortunately, the amount of mixed precipitation recorded at Dabenying is only a small fraction of the total precipitation (approximately 7%). A comparison between the observed and simulated data can be found in Figure 2.

In addition to snowfall, WRF-ARW-simulated values for SWE, snow depth, and rainfall were extracted as a supplement to the moisture trend analysis, and WRF-ARW-simulated temperature values were extracted as a supplement to climate change events across the HRB.

2.4. MODIS-Terra FSC Remote Sensing Product

This product is based on a snow-mapping algorithm named SNOMAP [42,43,44] that uses MODIS from the Terra satellite’s instruments and employs a Normalized Difference Snow Index (NDSI), along with a criteria test to provide daily and global snow cover. The MODIS/Terra Snow Cover Monthly L3 Global 0.05Deg CMG (MOD10CM) dataset [45], which is calculated from the daily global snow cover products and is new in Version 5 (V005), was selected for this study. This dataset contains FSC and Quality Assessment (QA) data in Hierarchical Data Format-Earth Observing System (HDF-EOS) format, along with corresponding metadata from 2000 to present, organized in a 0.05-degree Climate Modeling Grid (CMG).

2.5. Empirical Orthogonal Function Analysis

Empirical Orthogonal Function (EOF) analysis is used by climate scientists as a powerful tool for compressing data and reducing its dimensions with the goal of expressing possible spatial variability modes (i.e., patterns) and to explore how these patterns change over time [46]. Gridded climate data can be described as follows:

$$F_{ijk}=F\left(t_i,{\theta }_j,{\varnothing }_k\right)$$

(1)

where $F_{ijk}$ is a three-dimensional function with discretized time $t_i,i=1,\dots ,n$, latitude ${\theta }_j,j=1,\dots ,p_1$, and longitude ${\varnothing }_k,k=1,\dots ,p_2$. We can reshape $F_{ijk}$ into a gridded space-time field $X\left(t,s\right)$ as follows:

$$X=\left(\begin{array}{ccc}{x}_{11}& \cdots & x_{1p}\\ ⋮& \ddots & ⋮\\ x_{n1}& \cdots & x_{np}\end{array}\right)$$

(2)

$$\overline{x}=\left({\overline{x}}_{·1},\dots ,{\overline{x}}_{·p}\right)$$

(3)

$$x_{ts}^{\prime }=x_{ts}-{\overline{x}}_s$$

(4)

$$X^{\prime }=X-{\left(1,\dots ,1\right)}^T\overline{X}=(I-\frac{1}{n}\left(1,\dots ,1\right){\left(1,\dots ,1\right)}^T)X$$

(5)

where $p=p_1\ast p_2$. $\overline{x}$ is the time average of field $X$ for spatial grid points, $x_{ts}^{\prime }$ is the anomaly field, and $I$ is the $n\times n$ identity matrix. The covariance matrix is defined by

$${\displaystyle \sum }=\frac{1}{n-1}{X}^{\prime T}{X}^{\prime }$$

(6)

Given a direction $\mathrm{a}={\left(a_1,\dots ,a_p\right)}^T$, $X^{\prime }a$ has maximum variability. The variance of the time series $X^{\prime }a$ is defined by

$$\mathrm{var}\left(X^{\prime }a\right)=\frac{1}{n-1}{\left|\right|X^{\prime }a\left|\right|}^2=\frac{1}{n-1}{\left(X^{\prime }a\right)}^T\left(X^{\prime }a\right)=a^T{\displaystyle \sum }\mathrm{a}$$

(7)

The goal of the following is to find the top most-weighted components for the time series gridded data:

$$\underset{a}{\max }\left(a^T{\displaystyle \sum }\mathrm{a}\right),s.t.a^{T}a=1$$

(8)

$${\displaystyle \sum }\mathrm{a}=\lambda a$$

(9)

$$\frac{100{\lambda }_k}{{{\displaystyle \sum }}_{k=1}^p{\lambda }_k}\%$$

(10)

where $\lambda a$ is a simple eigenvalue problem, the k-th eigenvector $a_k$ of the covariance matrix is the k-th EOF, and $c_k=X^{\prime }{a}_k$ is the k-th principal component.

$$c_k\left(t\right)={\displaystyle \sum }_{s=1}^p{x}^{\prime }\left(t,s\right)a_k\left(s\right)$$

(11)

In our study, six sets of $F_{ijk}$ are composed of gridded precipitation, temperature, snowfall, SWE, snow physical depth simulated by the WRF model and monthly satellite-retrieved fractional snow cover from 2000 to 2013. The returned eigenvalues are normalized such that the sum of the squares for each EOF pattern equals one. The principal component (PC) time series (eigenvectors) are normalized by the weights to obtain the time series of the mean areal amplitudes.

2.6. Mann-Kendall/Theil-Sen Trend Analysis

The Mann-Kendall trend test and the Theil-Sen estimator are used to derive the trend in the variables over time. The Mann-Kendall test was formulated by Mann [47] to determine non-linear trends and turning points.

The Mann-Kendall test statistic, S, is calculated by the following equations:

$$\mathrm{S}={\displaystyle \sum }_{i=1}^{n-1}{\displaystyle \sum }_{j=i+1}^{n}sgn\left(x_j-x_i\right)$$

(12)

$$sgn\left(x_j-x_i\right)=\{\begin{array}{c}+1, if x_j-x_i>0\\ 0, if x_j-x_i=0\\ -1, if x_j-x_i<0\end{array}$$

(13)

$$Var\left(S\right)=\frac{n\left(n-1\right)\left(2n+5\right)-{{\displaystyle \sum }}_{i=1}^m{t}_i\left(t_i-1\right)\left(2t_i+5\right)}{18}$$

(14)

$$Z_S=\{\begin{array}{c}\frac{S-1}{\sqrt{Var\left(S\right)}}, if S>0\\ 0, if S=0\\ \frac{S+1}{\sqrt{Var\left(S\right)}}, if S<0\end{array}$$

(15)

where $n$ is the total number of data points, $x_i$ and $x_j$ are the data values at time series $i$ and $j$ ($j>i$), respectively, $sgn\left(x_j-x_i\right)$ is the sign function, $m$ is the number of tied groups, $t_i$ denotes the number of ties of extent $i$, $Var\left(S\right)$ is the variance, and $Z_S$ is the standard normal test statistic. Positive values of $Z_S$ indicate an increasing trend, and vice versa.

The Theil-Sen estimator is a robust regression mechanism and an unbiased estimator of the true slope for simple linear values [48,49].

The Theil-Sen estimator for a simple linear regression is shown below:

$$Y_t=\alpha +\beta t+{\epsilon }_t$$

(16)

where $Y_t$ is the variable value at year $t$, $\alpha$ is the linear intercept, $\beta$ is the linear slope, and ${\epsilon }_t$ is the deviation from the linear straight line.

3. Results

3.1. Spatial Patterns

Figure 3 shows the spatial distribution of the average annual precipitation, 2m temperature, snowfall, SWE, snow depth and FSC; of the data simulated by the WRF model, FSC was retrieved from MODIS/Terra. This figure shows that the amounts of precipitation flux and snow-related variables decrease from upstream to downstream. Across the upstream region, precipitation flux in the eastern portion is larger than that in the western portion; in other words, precipitation flux peaks in the eastern portion of the upstream area. However, the snow-related variables show a reverse trend; these variables peak near the southwestern HRB, where the Qiyi Glacier is located and where annual 2m air temperatures are lowest. The spatial distribution and the peak of FSC extracted from the remote sensing data agree well with the variables simulated by the WRF model, which indicates that the satellite-retrieved data used in this study and that simulated by the model are both reasonable and effective.

Figure 4 shows the altitudinal profiles of mean annual precipitation, 2m temperature, snowfall, SWE, snow depth and FSC at various elevations from 2000 to 2013. The 2m temperature data follow the temperature lapse rate across the HRB near the land surface; in other words, the higher the altitude, the lower the temperature. The precipitation flux is high at elevations between 2500 and 4000 m, especially in the 2500–3000 m and 3500–4000 m ranges. Snowfall is greatest between 2500 and 3000 m, well above the snowfall rates at other elevations. The SWE, snow depth and FSC are greatest at the highest elevations and are larger in the 2500–3000 m range than in the 3000–3500 m range, which may have resulted in maximum snowfall being measured in the 2500–3000 m range.

3.2. Temporal Pattern

Figure 5 shows that the trends in accumulated annual precipitation across the upstream and midstream zones increased noticeably from 2000 to 2013. Precipitation increased at a rate of 1.82 mm/year across the entire HRB basin, and 8.57 mm/year, 2.01 mm/year and 0.36 mm/year for the upstream, midstream and downstream zones, respectively. The pace of these increases is consistent among the sub-basins. Temperature increased at a rate of 0.05 K/year across the HRB basin and 0.06 K/year, 0.05 K/year and 0.05 K/year for the upstream, midstream and downstream zones, respectively.

The accumulated annual snowfall across the upstream and midstream zones rose slightly between 2000 and 2013 but did not rise in the downstream region. The rates (coefficients) of increase in snowfall are 0.55 mm/year, 0.02 mm/year and −0.25 mm/year for the upstream, midstream and downstream zones, respectively. Figure 6 also shows that snowfall in the upstream zone is far higher than in the midstream or downstream zones. The accumulated annual snowfall in the upstream zone is nearly 4 times as high as in the midstream zone and more than 10 times as high as in the downstream zone.

Both the SWE and snow depth decreased across the upstream and downstream zones. FSC increased slightly across the upstream zone. Peaks were evident in the downstream curves in 2003 for water snow equivalent, snow depth and FSC, which again indicates that the data from satellite remote sensing and from the model simulation are both reasonable and effective. Moreover, these anomalies are not synchronized in 2009. FSC decreased at a rate of −0.09%/year across the entire basin and −0.14%/year, −0.10%/year and −0.08%/year for the upstream, midstream and downstream zones, respectively. Temperatures increased significantly across the entire basin.

Figure 6 shows that annual precipitation, 2m temperature and snow-related variables varied with altitude from 2000 to 2013. The precipitation flux at elevations over 4500 m decreased sharply from 2000 to 2013, while it increased markedly from 2500 to 3000 m, reaching its largest values in 2005. Snowfall increased at elevations above 1500 m. At elevations between 2500 and 3500 m and above 3500 m, SWE and snow depth showed prominent decreases, especially at elevations over 4500 m, where FSC also decreased noticeably; however, it increased slightly between 3500 and 4500 m.

Figure 7 shows the mean annual number of snow days across the HRB and its sub-basins from 2000 to 2013 based on the gridded snowfall data from the WRF model. The mean annual number of snow days decreased notably from 2000 to 2013 for the entire region from upstream to downstream. The maximum mean annual number of snow days across the upstream zone was 76 days in 2003, and that zone’s minimum mean annual snow days was 48 days in 2013. FSC decreased at a rate of −0.76 days/year, −0.44 days/year, and −1.05 days/year for the upstream, midstream and downstream zones, respectively, while it decreased at a rate of −0.85 days/year for the entire HRB region.

3.3. Long-Term Observational Snowfall

Figure 8 shows annual precipitation, snowfall and 2m temperature from the CMA observation data across the HRB from 1960 to 2010. From 1979 to 2010, no precipitation type was recorded; therefore, precipitation was divided into liquid, solid and mixed precipitation according to Ding [21], as described in Section 2.2. Clearly, 2m temperature values increased at all 16 CMA observation stations by approximately 2 K degrees over the last 50 years. Precipitation flux increased sharply at one station (ID: 52737) located at an altitude of 2981.5 m in the upstream western portion of the HRB. The increasing precipitation trend could also be seen at two high-altitude stations (IDs: 52633 and 52645). A slight decreasing precipitation trend was found at one station (ID: 52378). Snowfall increased at three stations (IDs: 52737, 52546 and 52645) and decreased at two stations (IDs: 52754 and 52633).

The Mann-Kendall/Theil-Sen trend analysis of the CMA observation data for the past 50 years yielded the following results: precipitation increased by 0.58 mm/year with a probability of 0.98, snow increased by 0.08 mm/year with a probability of 0.90, and temperature increased by 0.036 K/year with a probability of near 1.00.

3.4. EOF Analysis

The first spatial pattern of normalized EOFs (Figure 9, left two columns) explains 86.6%, 99.6%, 50.8%, 47.8%, 59.2% and 28.8% of the total variance for annual precipitation, 2m temperature, snowfall, SWE, snow depth and FSC, respectively. Except for SWE and FSC, the first EOFs for snowfall and snow depth explain more than 50.0% of the total variance. The first four EOFs of FSC explain 28.8%, 12.8%, 11.5% and 10.1% of the total variance, respectively. Therefore, the first four EOFs of FSC are effective. The EOF 1 of the precipitation is high across the upstream zone and decreases slowly from upstream to downstream; the highest value is in the southeastern upstream area. The EOF 1 of the 2m temperature shows a high center across the downstream zone and decreases slowly from downstream to upstream; the highest value is in the northeastern downstream area. The EOF 1 for snowfall shows alternating stripes of high values from upstream to midstream; mostly high values occur at mid-altitudes, while a negative value occurs in the eastern downstream area. The EOF 1 s for SWE and snow depth show low values across the entire HRB; the lowest value is in the western upstream area. The EOF 1 for FSC is low across the entire HRB, with minimums centered in the midstream zone. Some no-value spots appear at the EOFs for FSC, which indicate that water bodies were misclassified as snow cover across the HRB.

The first normalized PC time series (the right column in Figure 9) for precipitation, 2m temperature and snow-related variables were standardized to have a zero mean and unit variance. Precipitation increased during the rainy season, temperature increased during the cold season, and snow-related variables had large variations during the cold season over the past 14 years. Snowfall tended to decrease in December; however, it increased in October, November, January and February, which is why the peak in the third bar of Figure 9 (right column) is deflated and total snow flux is elevated.

3.5. Mann-Kendall/Theil-Sen Trend Analysis

The results of the Mann-Kendall/Theil-Sen trend analysis for precipitation and snow-related variables at different altitude ranges are shown in

Table 1. Probabilities and trends in precipitation, snowfall, SWE, snow depth and FSC at different altitude ranges.

Elevation (km)		<1.0	(1.0, 1.5)	(1.5, 2.0)	(2.0, 2.5)	(2.5, 3.0)	(3.0, 3.5)	(3.5, 4.0)	(4.0, 4.5)	≥4.5
Precipitation	Probability	0.00	0.26	0.94	0.96	0.96	0.56	0.77	0.73	0.92
	Trend (mm/year)	−0.10	0.14	3.142	11.02	14.01	−5.23	9.00	7.13	−15.12
Snowfall	Probability	0.77	0.81	0.09	0.42	0.26	0.34	0.34	0.17	0.09
	Trend (mm/year)	−0.07	−0.05	0.00	0.53	1.12	0.24	0.58	0.21	−0.16
Snow water equivalent	Probability	0.81	0.26	0.09	0.90	0.81	0.42	0.09	0.49	0.77
	Trend m-2/year)	−0.01	0.00	0.00	0.06	0.07	−0.01	0.00	−0.03	−0.10
Snow depth	Probability	0.87	0.81	0.00	068	0.68	0.09	0.56	0.87	0.99
	Trend (m/year)	−0.15	−0.07	0.00	0.28	0.45	0.02	−0.20	−0.57	−1.39
Fractional snow cover	Probability	0.62	0.56	0.00	0.77	0.87	0.09	0.42	0.00	0.73
	Trend (%/year)	−0.06	−0.13	−0.01	0.25	0.25	0.05	0.11	0.03	−0.68

. The largest increases in precipitation and snow-related variables occurred at elevations from 2500 to 3000 m with high probabilities. The negative trends for SWE, snow depth and FSC at elevations below 1.0 km and higher than 4.5 km had high probabilities. The largest decreasing precipitation trend occurred at elevations above 4.5 km with a high probability, while the largest decreasing snowfall trend also occurred at this altitude but had a low probability. At mid-elevations, the precipitation and snow-related variable trends were generally positive, while at the lower and upper elevations, they were generally negative.

4. Discussion

4.1. Response of Snow-Related Variables to Climate Change

Overall, the results unequivocally indicate a correlation between increasing temperature and precipitation as simulated by the WRF model across the HRB region: the positive trends for both variables were improved based on data reported in numerous studies [50,51,52,53,54]. However, the response of snow-related variables to increasing temperatures is complex, as shown in the EOF analysis and Mann-Kendall/Theil-Sen trend analysis results.

The effect of climate change (increasing temperatures) on FSC, SWE and snow depth across the upstream zone of the HRB is significant. In this area, these snow-related variables have decreased notably since 2000. Wang and Li [55] indicated that climate change had caused snowmelt runoff peaks to advance and had enlarged snowmelt runoff significantly in the upstream zone of the HRB based on observation data from the second half of the last century.

The effect of climate change on snowfall is complex. Snowfall tends to increase as temperature increases [56,57,58] across the upstream zone. This increased snowfall occurs mainly because (1) increasing temperatures contribute to an increase in the vapor ratio by increasing evaporation and (2) climate change impacts global water, energy balance, and redistribution values. Decreases in the number of snowfall days is primarily linked to temperature increases, which cause a delay in snow formation in autumn and cause snowfall to melt earlier in spring.

The reduced number of snowfall days coupled with a larger snowfall volume will increase the occurrence of extreme snow events across the upstream zone of the HRB. Peak spring runoff will be both weaker and occur earlier due to decreasing FSC, forcing agriculture to depend more heavily on irrigation in the midstream zone. In the downstream zone, increased evaporation resulting from increasing temperatures and reduced snowfall (this region experienced almost no increase in precipitation) will exacerbate impacts on an already weak ecological system.

4.2. The Importance of High Spatial Resolution Data in the Evaluation of the Impact of Climate Change on Snow

The elevation effect, especially in mountainous regions, is an additional and locally productive source of snow. The finer the spatial resolution of snow-related variables, the more effective such data are for evaluating the effects of climate change on these variables using detailed altitude information.

RCMs yield high-resolution meteorological data that can be used to evaluate the impact of climate change on regional water fluxes, especially those with high spatial resolution, which provides more detailed information concerning land-surface heterogeneity and finer forcing such as complex topography and land-surface variations [59,60,61].

4.3. Mutual Proofs for Remote Sensing Products and WRF Simulation

Remote sensing data has enormous potential to supplement ground observations. The mutual proof between remote sensing products and regional climate model simulations supports the reliability of regional scale climate change evaluations based on these data, especially across complex mountainous regions with sparse snow observation sites. In this study, a number of mutual proofs between the MODIS-Terra Fractional Snow Cover variable data (including snowfall, SWE and snow depth) and the WRF-simulated data were generated. For spatial distribution (Figure 3), the peak of mean annual FSC derived from the remote sensing data occurred near Qiyi Glacier in the western upstream zone of the HRB; this is the same location where peak mean annual snowfall, SWE and snow depth were simulated by the WRF model. For altitudinal profile, mean annual FSC at 2500–3000 m range was larger than that at the 3000–3500 m range; the same phenomenon was found for mean annual SWE and snow depth as well.

4.4. MODIS-Terra FSC Remote Sensing Product

This study used the monthly high spatial resolution MODIS-Terra FSC remote sensing product. The MODIS-Terra FSC data agrees well with the spatial patterns of SWE and snow depth over the last 14 years, and the peak location identified using the WRF simulation matches well with those two variables’ peaks. These data agreements validate one another: the satellite-derived data validates the simulated regional climate model data and vice versa. One of the MODIS data’s flaws across the HRB, as mentioned earlier, is that water bodies are sometimes misclassified as snow cover along the Heihe River, as was further demonstrated by the EOF analysis results. In the downstream zone of the HRB at the first eigenvalue of EOF analysis, no data exist for the water body grids because these have been misclassified as snow cover continuously over the past 14 years; therefore, there is no variation in these grids. Consequently, it is necessary to inspect and calibrate such misclassifications in future research.

5. Conclusions

This study evaluated the effects of climate change on snow based on both the long-term high spatial-temporal resolution snowfall simulated by the WRF model and the long-term high spatial resolution monthly FSC attained by the MODIS-Terra product over the past 14 years, supplemented by temperature and precipitation data. The results indicate that temperature has increased, that snowfall has increased, and that FSC, snow depth, SWE and the number of snowfall days has decreased. The effects of climate change on the HRB can be concluded as follows. (1) Due to increases in air temperature, the FSC, SWE, and snow depth across the HRB region have decreased over the past 14 years, especially at elevations above 4500 m; however, snowfall has increased in the mid-altitude ranges across the upstream zone of the HRB. (2) Across the upstream zone of the HRB, total snow flux has increased; however, the number of snowfall days has decreased over the past 14 years. Therefore, extreme snow events may increase across the upstream zone of the HRB. (3) Air temperature in the HRB across the downstream zone increased the most; however, snowfall across this zone has also decreased over the past 14 years. Thus, it is effective and reasonable to explore the effect of climate change on snow at the watershed scale using ground station data, WRF simulations and remote sensing data.

Given the effects of climate change on snow across the HRB as revealed in this study, our future work will assimilate snow-relative remote sensing products to improve the precision of snowfall forecasting through precipitation forecasting and, thus, be able to provide early warnings of extreme snowfall events in the upstream zone, predict irrigation recharge for the midstream zone and identify the need for ecological renewal in the downstream zone.