Attribution Analysis of Runoff Change in the Upper Reaches of the Kaidu River Basin Based on a Modified Budyko Framework

Abstract

The Kaidu River plays an important role in the water development and utilization in the Tarim River basin in northwestern China. In this study, we used a modified Budyko framework, which considered the snowmelt to analyze and attribute the runoff change in the upper Kaidu River basin based on the observations during the period of 1960–2010. The time series was divided into two periods: 1960–1995 and 1996–2010. The contribution rate of runoff change between these two periods and the elasticity coefficient of runoff were estimated to quantify the effect of climatic variables and landscape changes on runoff alteration. The results show that the increase in precipitation was the major cause of increase in runoff, whose contribution accounted for 81.42%. The contribution rate of the landscape change was lower than that of the precipitation change, accounting for 9.07%. The elasticity coefficient of runoff to precipitation was 1.24, and the elasticity coefficient of runoff to the landscape was −0.74. Compared with the original Budyko framework, without considering the snowmelt, the contribution rates of precipitation and potential evaporation to runoff change would decrease after considering the snowmelt in the modified Budyko framework, while the contribution rate of landscape would increase. The increased snow ratio would cause more fluctuations in the runoff. This study provides a valuable reference for the water resources management in the upper Kaidu River basin and deepens our understanding of the response of runoff to climate change in snowmelt-affected regions.

1. Introduction

Cold regions affected by snowmelt provide water resources for about one-sixth of the global population and influence most of the industrialized countries. However, they are expected to be particularly vulnerable to climate change [1]. As reported by the IPCC, the rise in global average air temperature was significant in the past 40 years, and it will reach 1.5 °C over the next two decades even under the low greenhouse gas emissions scenarios [2]. As a result, streamflow in snowmelt-affected catchments is undergoing significant changes due to climate change [3]. However, at the catchment scale, the contributions of snowmelt change to runoff alteration may vary significantly in a changing climate [4], which deserves great attention for water resources management in cold regions.

Snowfall is an important component of the water cycle in arid or semi-arid areas, and it has a great impact on freshwater supply and ecosystem [5]. Previous studies have reported a significant effect of snowfall change on runoff. Adam et al. [6] used a physical hydrological model to assess the impact of climate change on runoff, indicating that warm season runoff has decreased in most snow-affected areas, and the distribution of runoff at an annual scale has been altered. Nohara et al. [7] found that a warming climate would lead to earlier snowmelt and an increase in runoff, especially in middle- or high-latitude rivers. Barnhart et al. [8] suggested that future increases in regional air temperature may lead to earlier but slower snowmelt in mountainous areas, which might reduce the below-root zone percolation and streamflow, resulting in a corresponding reduction in runoff. In general, although the response of runoff to climate change in snowmelt-affected catchments has been extensively studied, attribution analysis of runoff changes remains a great challenge. This is because there are complex interactions between climate, vegetation, topography, and soils in snowmelt-affected areas [9].

Some methods are proposed for the attribution analysis of runoff changes in snowmelt-affected catchments. For example, the statistical method is widely used because of the low requirement for computation and great adaptability [10]. However, the statistical method simplifies or ignores the physical processes, which results in great uncertainties [11]. The physically based hydrological model, such as SRM and GE-SWAT, can simulate the runoff in detail, but it has high requirements on data and computation resources, with a large number of parameters to calibrate [12,13]. In addition, the water balance method combined with the Budyko framework has been considered an effective approach for estimating the contribution of climate variability and landscape changes to runoff alteration [14]. The Budyko framework is an analytical method that includes evapotranspiration, available energy, and available water supply, using a parameter that represents the basin characteristics [15,16,17]. The Budyko framework has been widely applied to quantitively study the effect of climate changes on runoff in different catchments [18,19]. However, the effect of snowmelt was rarely considered in the original Budyko framework, and the application of the Budyko framework in cold regions is not adequate.

The Kaidu River is located on the southern slope of the Tianshan Mountains in northwest China [20]. It provides 83.4% of inflow water to Bosten Lake [21]. Therefore, the Kaidu River plays an important role in the water resources development and utilization in the Yanqi Basin. In recent years, the climate in the Tianshan Mountains has exhibited a warmer and wetter trend [22], and the runoff of the Kaidu River has shown an obvious increasing trend in the past 40 years [20]. However, there are few studies on quantitatively evaluating the effect of different factors on runoff change, and the impact of snowmelt remains unclear [23]. Therefore, the response of runoff to climate change in the upper Kaidu River basin needs to be further examined.

Due to the complex relationship between changes in temperature and precipitation and their combined impact on snowmelt-affected runoff, trend detection and attribution analysis of runoff change are still difficult in cold regions [24]. In previous studies, as a common method, hydrologic models, such as SWAT, GR6J, and VHM, have been applied in different sub-regions in the Kaidu River basin to quantify the response of the hydrological cycle to climate changes [25,26,27]. Chen et al. [28] used Surfer 9.0 to obtain the quantitative runoff depth–precipitation–potential evapotranspiration relations of the upper Kaidu River basin, indicating that the graphic method was effective in quantifying the potential impacts of climate fluctuation on runoff. Zhuang et al. [29] developed an integrated multi-GCM-based stochastic weather generator and stepwise cluster analysis method to explore the nonlinear relationship between the input variables and outputs in basin systems. However, the contributions of snowmelt to the runoff are still poorly understood [30].

The main purpose of this study is to quantitatively evaluate the influence of different factors on the runoff change in the upper reaches of the Kaidu River basin using a modified Budyko framework incorporating snowmelt. The first-order hydro-climatic relationship (i.e., the relationship between runoff and snowmelt) is investigated [31].

2. Study Area and Data

Originating from the central Tianshan Mountains, the Kaidu River is one of the sources of the Tarim River. The upper reaches drain the headwater region of the river upstream of the Dashankou hydrological station, located within 41°30′ N to 43°21′ N latitudes and 82°58′ E to 86°55′ E longitudes (Figure 1a). The length of the upper reaches is 2222 km, and the drainage area is 18,827 km², with an average annual runoff amount of 3.49 × 10⁹ m³ [32]. The Kaidu River basin is dominated by a continental climate, and the air moisture mainly comes from the humid westerlies and the Arctic Ocean air in the northwest. The mean annual precipitation is about 332 mm, and the rainy season is from May to August. The upper reaches of the Kaidu River basin are located in the mountain area, and the mean annual pan evapotranspiration is more than 1100 mm. The evapotranspiration increases from north to south and from high altitude to low altitude [33]. The average annual air temperature in the basin is about −4.6 °C. The snow season generally lasts from November to March in the following year, lasting 139 days, and the maximum basin-averaged snow depth is 12 cm [34]. The basin is partially covered by glaciers with an area of 490.68 km², and the reserves are about 25.25 km³ [35]. Due to the limited rainfall, the basin is mostly covered by grassland, and the percentage of forest coverage is only 0.41%.

The measured daily climate data used in this study were obtained from 11 meteorological stations in and around the upper Kaidu River basin, and the observed river discharge data at the Dashankou hydrological station located at the outlet of the upper Kaidu River basin were collected from the Hydrological Yearbook in China. The locations of the stations are shown in Figure 1b, and detailed information about the stations is shown in

Table 1. The information on the hydrological and meteorological stations in the upper Kaidu River basin.

Station Number	Station Name	Latitude and Longitude	Elevation (m)	Station Type
51,241	Tuoli	(45.56 N, 83.36 E)	1077.8	Meteorological station
51,243	Kelamayi	(45.37 N, 84.51 E)	449.5	Meteorological station
51,334	Jinghe	(44.34 N, 82.49 E)	329.2	Meteorological station
51,367	Hutubi	(44.1 N, 86.51 E)	575.1	Meteorological station
51,469	Urumqi Mushizhan	(43.27 N, 87.11 E)	1930	Meteorological station
51,470	Tianchi	(43.53 N, 88.07 E)	1942.5	Meteorological station
51,477	Dabancheng	(43.21 N, 88.19 E)	1103.5	Meteorological station
51,542	Bayinbuluke	(43.02 N, 84.09 E)	2458	Meteorological station
51,567	Yanqi	(42.05 N, 86.34 E)	1055.3	Meteorological station
51,639	Shaya	(41.14 N, 82.47 E)	980.4	Meteorological station
51,656	Kuerle	(41.45 N, 86.08 E)	931.5	Meteorological station
705,700	Dashankou	(42.25 N, 85.73 E)	1270	Hydrological station

.

We collected daily meteorological data in a long series of 51 years from 1960 to 2010, including air temperature, relative humidity, wind speed, and sunshine hours. The 0.25° resolution gridded data set of precipitation was provided by the National Meteorological Information Center of China (http://cdc.nmic.cn, accessed on 10 April 2021). The hydrological station provided the annual runoff data from 1960 to 2010. The DEM data with a 90 m spatial resolution were obtained by the Shuttle Radar Topography Mission (SRTM) dataset (http://srtm.csi.cgiar.org, accessed on 20 October 2020). The land use data were derived from the remote sensing monitoring database of land use status in China provided by the Resource and Environmental Science and Data Center (http://www.resdc.cn, accessed on 26 March 2021).

Since the meteorological station did not record the types of precipitation, it is necessary to distinguish the precipitation into snowfall and rainfall. We used the parameterization scheme developed by Ding et al. [36] to determine the type of precipitation, and the snow and sleet quantity were counted in the snowfall amount:

$$type\left\{\begin{array}{l}snow, T_w\le T_{min}\\ sleet, T_{min}\le T_w\le T_{max}\\ rain, T_w\ge T_{max}\end{array}\right.$$

(1)

where T_w is the daily mean wet-bulb temperature, including air temperature, humidity, and atmospheric pressure information, which is more suitable for indicating the precipitation type than the air temperature. T_w is deduced as:

$$T_w=T_a-\frac{e_{sat}\left(T_a\right)\left(1-RH\right)}{0.000643p_s+\frac{\partial e_{sat}}{\partial T_a}}$$

(2)

where RH is the relative humidity; p_s is the atmospheric pressure (hPa); T_a is the air temperature (°C); e_sat (T_a) is the saturated vapor pressure at T_a, which is estimated as:

$$e_{sat}\left(T_a\right)=6.1078exp (\frac{17.27T_a}{T_a+237.3})$$

(3)

T_min and T_max are characteristic temperatures, representing the boundary between snow, sleet, and rain. They can be calculated using RH [36].

The study basin was divided into a grid system of 152 rows and 258 columns. Daily meteorological data averaged of the whole study basin were estimated by spatial interpolation based on the observations at meteorological stations, and the air temperature was interpolated considering its gradient of elevation. Potential evapotranspiration was calculated by the Penman-Monteith equation [37]. Then, the daily data were summarized into annual data.

3. Method

3.1. Mann-Kendall Nonparametric Test

The Mann-Kendall nonparametric test is used to detect trends in the hydrologic series and to determine the abrupt year of changes [38]. For one time series, x₁, x₂, …, x_n, the following parameters were defined as:

$$s_k={\displaystyle \sum }_{i=1}^k{r}_i\left(k=2,3,4,\cdots ,n\right)$$

(4)

$$r_i=\left\{\begin{array}{l}1,x_i>x_j\\ 0,\hspace{1em}\hspace{1em}or\end{array}\hspace{1em}\left(j=1,2,\cdots ,i\right)\right.$$

(5)

Assuming that the samples are random and independent, the statistic UF_k (UF_k is a normally distributed statistic) can be calculated as:

$$U{F}_k=\frac{s_k-E\left(s_k\right)}{\sqrt{Var\left(s_k\right)}}\left(k=2,3,4,\cdots ,n\right)$$

(6)

where UF₁ = 0; E(s_k) and Var(s_k) can be determined by the following formulas:

$$E\left(s_k\right)=\frac{n\left(n+1\right)}{4}$$

(7)

$$Var\left(s_k\right)=\frac{n\left(n-1\right)\left(2n+5\right)}{72}n$$

(8)

If UF_n is larger than the critical value of 1.96 (or lower than −1.96), the trend of the variable in the whole study period is considered to be significant at a 5% significance level. We can put the time series in reverse order and use Equations (4)–(8) to calculate a new statistic UB_k. At a 5% significance level, if UF_k is outside the range of critical values (α = ±1.96), the increasing trend (UF_k > 0) or decreasing trend (UF_k < 0) will be significant. The intersection of UB_k and UF_k is the abrupt point in the time series. Accumulative anomalies were also used to validate the result of the Mann-Kendall nonparametric test.

3.2. The Budyko Framework

In this study, we used the Yang-Choudhury equation [39] as the original Budyko framework:

$$E=\frac{E_{P}P}{{\left(P^n+E_P^n\right)}^{1/n}}$$

(9)

where P is precipitation (mm); E is actual evapotranspiration (mm); E_P is potential evapotranspiration (mm); n is a parameter representing the landscape of the basin. According to the water balance equation, neglecting the water storage change in the basin at the interannual scale, runoff (R) can be expressed as:

$$R=P-E$$

(10)

Equation (9) can be rewritten as follows:

$$\frac{P-R}{P}=\frac{E_P}{{\left(P^n+E_P^n\right)}^{1/n}}$$

(11)

Many previous studies have extended this framework. Du et al. [40] used equivalent precipitation to replace P to represent different available water resources for evapotranspiration. Here, we adopted a new modified framework proposed by Zhang et al. [41], considering the snowmelt.

We introduced an assumption that the snow cover sublimates barely and the melting snow is too little to produce runoff in the snow season; the melting snow water flows into rivers in the form of surface runoff without evaporation or infiltration loss in the melting season [41,42]. In other words, the entire process of snow, from accumulation to melting, does not participate in the evapotranspiration processes. Therefore, the equivalent precipitation (P_e) in this framework can be expressed as:

$$P_e=(1-r_s)·P$$

(12)

where r_s is the snow ratio, which is the ratio of snowfall to precipitation. The energy term also needs to be adjusted accordingly. Because part of the energy for evapotranspiration is used to melt the snow, the equivalent available energy is equal to the difference between E_P and melting heat equivalence R_m/L, where L is the latent heat of evaporation (2500 kJ/kg). R_m is the energy required for all snow to melt, and it is about 0.14r_s·P according to the estimations by Zhang et al. [41]. The equivalent available energy ($E_P^e$) can be calculated as:

$$E_P^e=E_P-R_m/L=E_P-0.14r_s·P$$

(13)

Then, we replaced P and E_P in the original framework by P_e and $E_P^e$. The original framework can be changed to the following form:

$$\frac{P-R}{\left(1-r_s\right)·P}=[1+(\frac{E_P}{\left(1-r_s\right)·P}-\frac{0.14r_s}{1-r_s}{)}^{-n}{]}^{-1/n}$$

(14)

The r_s in the study basin was estimated to be about 0.28. The energy correction term 0.14r_s/(1 − r_s) was much smaller than the original term E_P/((1 − r_s)·P). For example, using the long-time average values, the energy correction term and original term were estimated to be 0.05 and 4.44, respectively. The difference between them are two orders of magnitude away. Thus, the energy correction term can be neglected, and Equation (14) can be simplified as:

$$1-\frac{R}{P}=[{\left(1-r_s\right)}^{-n}+(\frac{E_P}{P}{)}^{-n}{]}^{-1/n}$$

(15)

3.3. Contribution Analysis of Runoff Change

Here, we give an additional assumption that the runoff only changes between steady states without any transient changes [41]. Therefore, an analytical framework can be introduced to quantify the contribution of different factors to runoff change between two periods. Equation (15) can be rewritten in the following form:

$$R=P-[P^{-n}{(1-r_s)}^{-n}+E_P{}^{-n}{]}^{-1/n}$$

(16)

We assumed that the factors are independent of each other. Then, the total differential equation of R can be expressed as:

$$dR=\frac{\partial R}{\partial P}dP+\frac{\partial R}{\partial E_P}d{E}_P+\frac{\partial R}{\partial r_s}d{r}_s+\frac{\partial R}{\partial n}dn$$

(17)

where

$$\left\{\begin{array}{l}\hspace{1em}\hspace{1em} \frac{\partial R}{\partial P}=1-\frac{P-R}{P}\frac{E_P^n}{{\left[P\left(1-r_s\right)\right]}^n+E_P^n}\\ \hspace{1em}\hspace{1em} \frac{\partial R}{\partial E_P}=-\frac{P-R}{E_P}\frac{{\left[P\left(1-r_s\right)\right]}^n}{{\left[P\left(1-r_s\right)\right]}^n+E_P^n}\\ \hspace{1em}\hspace{1em} \frac{\partial R}{\partial r_s}=\frac{P-R}{1-r_s}\frac{E_P^n}{{\left[P\left(1-r_s\right)\right]}^n+E_P^n}\\ \frac{\partial R}{\partial n}=-\frac{P-R}{n}(\frac{ln\left[P^n{\left(1-r_s\right)}^n+E_P^n\right]}{n}-\frac{{[P\left(1-r_s\right)]}^{n}ln\left[P\left(1-r_s\right)\right]+E_p^n ln E_p}{{[P\left(1-r_s\right)]}^n+E_P^n})\end{array}\right.$$

(18)

The contribution of different factors to runoff change can be expressed as follows:

$${\eta }_x=\frac{\partial R}{\partial x}·\frac{dx}{dR}=\frac{\partial R}{\partial x}·\frac{\Delta x}{\Delta R}$$

(19)

where x is a factor that can influence the runoff. Δx is the difference of this factor between the two periods. ΔR is the difference in the runoff between two periods.

3.4. Climate Elasticity Method

We used the climate elasticity method to analyze the sensitivity of runoff to different factors. The concept of elasticity was used to explore the reaction of runoff under varying climatic conditions [43]. The elasticity coefficient of runoff (ε) can be calculated as:

$${\epsilon }_x=\frac{\partial R}{\partial x}·\frac{x}{R}$$

(20)

where ε_x represents the ratio of proportional change in runoff and proportional change in an influencing factor, indicating the sensitivity of runoff to this factor [44].

4. Results

4.1. Trends and Variations of Precipitation, Snowfall, E_P, and Runoff

We firstly analyzed the trends in different climate factors, which helped us better understand their impacts on runoff. Figure 2 displays the long-term trends in annual precipitation, snowfall, E_P, and runoff from 1960 to 2010. The annual precipitation, runoff, and average air temperature all showed a significant increasing trend, and the trend in annual precipitation was more obvious during the period of 1990–2010. The average values of annual precipitation, runoff, and air temperature were 358.1 mm, 185.4 mm, and −2.43 °C, with an increase of 17.8 mm, 10.3 mm, and 0.03 °C per decade, respectively. E_P had a significant decreasing trend, and it declined by 16.6 mm every decade, and the average E_P from 1960 to 2010 was 1144.2 mm. The average annual snowfall was 99 mm during the period of 1960–2010, with a variation of 8.6 mm per decade. The snow ratio in the period of 1960–2010 also increased.

We used the Mann-Kendall nonparametric test to analyze the trend of runoff. Figure 3a shows that the abrupt change point was 1995, and this result passed the test of significance at a 95% confidence level. Figure 3b shows that the accumulative value of runoff was at the lowest point in 1995. Therefore, the time series can be divided into two periods: 1960–1995 and 1996–2010 for analyzing the change in the runoff. The period from 1960 to 1996 is defined as the base period, and the period from 1996 to 2010 is defined as the change period. The parameter n in the two periods can be calculated by Equation (15).

Table 2. The value of annual precipitation, snowfall, average air temperature, snow ratio, E_P, runoff, wind speed at 10 m, sunshine time, and relative humidity during the base period (1960–1995) and change period (1996–2010).

Runoff Period	Precipitation (mm)	Snowfall (mm)	Temperature (°C)	EP (mm)	Snow Ratio	Runoff (mm)	Basin Parameter	Wind Speed (m/s)	Sunshine Time (h)	Relative Humidity
1960–1995	339.9	92.7	−2.7	1151.3	0.275	171.4	0.735	3.4	7.9	0.53
1996–2010	401.8	114.0	−1.7	1127.0	0.284	219.0	0.710	2.9	7.7	0.55

shows the average value of annual precipitation, snowfall, air temperature, snow ratio, E_P, runoff, parameter n, wind speed at 10 m, sunshine time, and relative humidity in the upper Kaidu River basin during the two periods.

Table 3. The trend of precipitation, snowfall, E_P, runoff, air temperature, and snow ratio in the period from 1960 to 2010.

Variable	Trend	Linear Change Rate
Precipitation	increase	1.78 (mm/year)
Snowfall	increase	0.86 (mm/year)
EP	decrease	−1.43 (mm/year)
Runoff	increase	1.03 (mm/year)
Temperature	increase	0.03 (°C/year)
Snow ratio	increase	0.001

Note: The items marked in red indicated that the trend was significant at 95% confidence level.

summarizes the trends of these variables detected by the Mann-Kendall nonparametric test. Yin et al. found that the wind speed, net radiation, and relative humidity played a dominant role in the change of E_P from 1964 to 1993 in northwestern China, and the wind speed and relative humidity played a dominant role in the change of E_P from 1994 to 2018 [45]. A decrease in sunshine time and wind speed and an increase in relative humidity will lead to a decrease in E_P [46]. Therefore, in spite of the increase in air temperature, the decrease in wind speed and sunshine hours, together with the increase in relative humidity in the change period compared with that in the base period (Table 2), caused a reduction in E_P.

We also display the data in the two periods in the Budyko space. As shown in Figure 4, the data points in the base period and change period were in different Budyko curves. Compared with the base period, the dryness index (E_P/P) in the change period decreased significantly, indicating a tendency for a warmer and wetter climate in the upper Kaidu River basin.

4.2. Attribution Analysis of Runoff Change

The climate factors (i.e., precipitation and E_P) and landscape factor n in the two periods varied to different degrees (Table 2). Compared with the data in the base period, runoff, precipitation, and snow ratio in the change period increased by 28.76% (47.6 mm), 18.23% (62.0 mm), and 3.35% (0.009), respectively; E_P and parameter n decreased by 2.12% (24.4 mm), and 3.39% (0.025), respectively.

As shown in Figure 5a,b, the increase in precipitation caused the runoff to increase by 38.7 mm, among which the contribution of snow ratio is 1.3 mm, and the corresponding contribution rates were 81.42% and 2.64%, respectively. The decrease in E_P and n caused an increase in the runoff by 0.9 mm and 4.3 mm, respectively. Their contribution rates were estimated at 1.83% and 9.07%, respectively. The contribution rate of other factors to the runoff change was 7.68%. These results indicate that the change in precipitation dominated the increase in runoff in the upper reaches of the Kaidu River basin. Change in landscape represented by the alteration in parameter n was the secondary cause of runoff change, and the contribution rate of the landscape change was 9.07%. Change in E_P had a quite limited impact on runoff change, which corresponded to the characteristics of arid or semi-arid regions, which was the water-limited region.

Figure 5c displays the effect of proportional changes in precipitation, snow ratio, E_P, and parameter n on runoff change. During the base period, an increase of 1% in precipitation and snow ratio increased the runoff by 1.24% and 0.22%, respectively. In contrast, an increase of 1% in E_P and parameter n reduced the runoff by 0.24% and 0.74%, respectively. During the change period, an increase of 1% in precipitation and snow ratio increased the runoff by 1.23% and 0.19%, respectively, and an increase of 1% in E_P and parameter n, in contrast, reduced the runoff by 0.23% and 0.69%, respectively. The elasticity coefficient of runoff to precipitation was the largest among all the elasticity coefficients, and the parameter n had the second largest elasticity coefficient. Compared with the elasticity coefficients in the base period, the elasticity coefficients of runoff to precipitation, E_P, r_s, and parameter n in the change period were reduced by 0.81%, 4.17%, 13.64%, and 6.76%, respectively. This result indicates that the elasticity coefficient of runoff to precipitation did not change significantly in the two periods.

The results of the contribution rates and elasticity coefficients suggested there was a significant impact of precipitation change on runoff alteration, but the effects of other factors should not be ignored. Although the contribution rate of precipitation change was evidently larger than that of landscape change (parameter n), which seems to suggest that the precipitation had a much larger effect on runoff, we should note that the absolute value of their elasticity coefficients was not quite different. The relatively small difference between the elasticity coefficients of runoff to precipitation and landscape (parameter n) indicates that the sensitivity of runoff to landscape (parameter n) was close to that of precipitation. Meanwhile, the change rate of precipitation (18.23%) was larger than that of landscape (3.39%), which led to a further expansion in the difference in contribution rates of change in precipitation and landscape (Equation (19)). That was why the contribution rates and elasticity coefficients exhibited different conclusions when the sensitivity of runoff to different factors was similar, and this phenomenon also occurred between E_P and snow ratio. In addition, the landscape is easily affected by human activities. Some studies proposed that the value of parameter n, which represents the landscape, was determined by vegetation coverage, soil property, slope, and land use [47,48,49]. The land use and vegetation coverage could be changed by humans in the short term. Several studies even argued that the parameter n could reflect the influence of human activity [50,51]. Therefore, we suggest that the impact of landscape change on runoff alteration is also important for water resources management in the upper Kaidu River basin.

4.3. The Influence of the Snowmelt on Runoff Change

We compared the results of this study with that of the original Budyko framework without considering the snowmelt. The parameter n in the original Budyko framework during the base period and change period was 0.574 and 0.564, respectively. Figure 6 shows the contribution rates of different factors to runoff change in the original Budyko framework and the results of the modified Budyko framework in this study.

The contribution rates of precipitation, E_P, and their total contribution rate in the original Budyko framework were larger than that in the modified Budyko framework; however, the contribution rate of parameter n in the original Budyko framework was less than that in the modified Budyko framework. For Equation (19), if factor x is precipitation, we can get the following equation at the interannual scale:

$${\eta }_P^{\prime }=\frac{\partial E}{\partial P}·\frac{\Delta P}{\Delta E}=\frac{\partial \left(P-R\right)}{\partial P}·\frac{\Delta P}{\Delta P-\Delta R}=\left(1-\frac{\partial R}{\partial P}\right)·\frac{\frac{\Delta P}{\Delta R}}{\frac{\Delta P}{\Delta R}-1}=\frac{\frac{\Delta P}{\Delta R}-{\eta }_P}{\frac{\Delta P}{\Delta R}-1}$$

(21)

$${\eta }_P=\frac{\Delta P}{\Delta R}-(\frac{\Delta P}{\Delta R}-1){\eta }_P^{\prime }$$

(22)

If the factor x is not precipitation, we can get another equation:

$${\eta }_x^{\prime }=\frac{\partial E}{\partial x}·\frac{\Delta x}{\Delta E}=\frac{\partial \left(P-R\right)}{\partial x}·\frac{\Delta x}{\Delta P-\Delta R}=-\frac{\partial R}{\partial x}·\frac{\frac{\Delta x}{\Delta R}}{\frac{\Delta P}{\Delta R}-1}=-\frac{{\eta }_x}{\frac{\Delta P}{\Delta R}-1}$$

(23)

$${\eta }_x=-(\frac{\Delta P}{\Delta R}-1){\eta }_x^{\prime }$$

(24)

where η_x is the contribution rate of one factor to runoff change; ${\eta }_x^{\prime }$ is the contribution rate of one factor to evapotranspiration. When ΔP is greater than ΔR, η_P is negatively correlated with ${\eta }_P^{\prime }$, and the absolute value of η_x is positively correlated with that of ${\eta }_x^{\prime }$.

Considering that snowmelt did not participate in the evapotranspiration process, the water available for evapotranspiration was less than precipitation in the modified Budyko framework. At the same time, part of the energy was consumed by snowmelt, which reduced the available energy for evapotranspiration. The available energy is generally more than the available water in arid and semi-arid areas. The reduction of available water or energy increased the sensitivity of evapotranspiration to water and energy. If the available water was reduced nearly to zero and the available energy remained constant, evapotranspiration would be limited by the available water, and any increase in available water would cause an equivalent increase in evapotranspiration [17]. Similarly, if the available energy was reduced nearly to zero and the available water remained constant, the evapotranspiration would be limited by the available energy, and any increase in available energy would cause an equivalent increase in evapotranspiration (the water and energy limitation in Figure 4). However, if the available water and energy were reduced at the same time and the reduction of the available energy was much smaller than that of the available water, the change rate of evapotranspiration with available water ($\partial E/\partial P_e$) would increase, and that to available energy ($\partial E/\partial E_P^e$) would decrease in the original Budyko framework. When the snowmelt was considered, the available water and energy in the modified Budyko framework were less than that in the original Budyko framework, and the difference between available water was more than that between available energy. The value of $\partial E/\partial P_e$ and $\partial E/\partial E_P^e$ in the original Budyko framework were 0.331 and 0.049, respectively, while the values of $\partial E/\partial P_e$ and $\partial E/\partial E_P^e$ in the modified Budyko framework were 0.375 and 0.036, respectively. Therefore, although the corresponding factor changed by the same amount in the two frameworks, the contribution rate of precipitation change to evapotranspiration alteration in the modified Budyko framework was larger than that in the original Budyko framework, and the absolute value of the contribution rate of E_P change to evapotranspiration alteration in the modified Budyko framework was smaller than that in the original Budyko framework. According to Equation (22), if the contribution of precipitation to evapotranspiration (${\eta }_P^{\prime }$) increased, the contribution of precipitation to runoff (η_P) would decrease. Similarly, in Equation (24), if the contribution rate of E_P to evapotranspiration (η_Ep′) increased (the value of η_Ep′ was negative), the contribution rate of E_P to runoff (η_Ep) would decrease. That was why the contribution rate of precipitation and E_P in the original Budyko framework were larger than that in the modified Budyko framework.

The appearance of snow also influences the landscape characteristics (such as soil, vegetation coverage, etc.) [52,53], which was reflected in the differences between the parameter n in the original and modified Budyko framework. There was a complex relationship between snow and landscape characteristics [41]. The contribution rate of parameter n had a larger value in the modified Budyko framework, which might be the result of these multiple mechanisms.

4.4. The First-Order Hydro-Climatic Relationship

Figure 7 shows that the snow ratio elasticity coefficient of runoff has a positive correlation with the snow ratio. A larger snow ratio corresponded to a larger runoff elasticity coefficient. In the wet years, with a high snow ratio, the basin would get more runoff than the basin with a low snow ratio when the snow ratio increased. In dry years, the runoff loss of the basin with a low snow ratio would be less than that with a high snow ratio when the snow ratio decreased. Liu et al. [31] raised a hypothesis that a higher snow ratio might lead to annual streamflow increases in the wet year and decreases in the dry year. The snow ratio elasticity coefficient of runoff in the upper reaches of the Kaidu River basin might partially prove this hypothesis. However, the correlation between precipitation and snow ratio was still not fully clear and needed further study.

5. Discussion

5.1. Effects of Different Methods for Distinguishing Precipitation Types on the Results

Han et al. [54] collected the air temperature and precipitation type data from 643 meteorological stations in China from 1961 to 1979 and summarized the critical temperature for distinguishing rainfall and snowfall in the whole country. The result showed that the critical temperature of rainfall and snowfall in northwestern China was about 3.70 °C (empirical temperature). Since we combined the amount of snowfall and sleet, the temperature that distinguished the precipitation types using T_w was equal to the T_max in Equation (1) in this study. Then, we summarized the monthly average T_max. Figure 8a indicates that the T_max was smaller than the empirical temperature in all months.

We used the empirical temperature to distinguish the daily precipitation types and obtained the monthly and annual snowfall and snow ratios calculated by the empirical temperature. As shown in

Table 4. The value of annual precipitation, snowfall, average air temperature, snow ratio, E_P, runoff, and basin parameters during the base period (1960–1995) and change period (1996–2010).

Runoff Period	Precipitation (mm)	Snowfall (mm)	Temperature (°C)	EP (mm)	Snow Ratio	Runoff (mm)	Basin Parameter
1960–1995	339.9	100.0	−2.7	1151.3	0.295	171.4	0.756
1996–2010	401.8	121.5	−1.7	1127.0	0.302	219.0	0.726

Note: The snowfall and snow ratio were calculated by the empirical temperature, and the basin parameters were recalculated after the change of snowfall and snow ratio.

, compared with the result obtained by Equation (1), the snowfall calculated by the empirical temperature in the two periods increased by about 7.00 mm, and the snow ratio increased by about 0.02. The difference was derived from the identification of precipitation type at 0–3.7 °C. This phenomenon usually occurs in the warm season from May to September (Figure 8b). Although the snowfall was different, the changes in snow ratio in the two periods calculated by empirical temperature and Equation (1) were quite close. Moreover, the contribution rates of different factors to runoff change calculated by the two methods had little differences (Figure 9).

Equation (1) considered the altitude, air pressure, and humidity of each grid to determine the precipitation type, which was more comprehensive and had high precision. The empirical temperature was determined by the records from meteorological stations and was more convenient for application. However, the complex topography and the lack of in situ observations might make the records in mountain areas less representative and then affect the results. Therefore, we recommended using Equation (1) to estimate snowfall in northwest China.

5.2. The Physical Response of Runoff to Different Variables

The runoff is greatly influenced by precipitation, and the increase in precipitation in the change period directly enhances the runoff. Padrón et al. [48] found that the parameter n was related to the landscape conditions. The vegetation coverage is the main factor influencing the parameter n [49]. Yang et al. [55] suggested that the increase in the parameter n represented an increase in vegetation coverage, and the increase in vegetation coverage would lead to a decrease in runoff in arid or semi-arid regions. Therefore, the decrease in parameter n in the study basin indicated a reduction in the vegetation coverage from the base period to the change period, which was consistent with the decreasing trend of NDVI from 1998 to 2012 in northwestern China [56]. The degradation of vegetation will result in the reduction of canopy interception and a decrease in transpiration, which is conducive to runoff generation.

The accumulated snowfall in the snow season will extensively melt into liquid water in the melting season. In the early stage of the melting season, the unfrozen ground makes it difficult for the meltwater to infiltrate [57]. Therefore, the increase in air temperature and snowfall during the change period led to the increase in snowmelt, which directly enhanced the surface runoff and total runoff.

5.3. Uncertainties and Limitations

There are still some limitations in this study. The framework used in this study ignored some terms in water balance. Zhu et al. [58] analyzed and compared the two glacier survey data in 1970 and 2009 and found that the numbers of glaciers in the Kaidu River basin decreased by 103, and the reserves decreased by 8.81 km³ in 40 years. Sorg et al. [59] measured ice-snow meltwater in the Tianshan Mountains, and their results showed that the contribution rate of glacier meltwater variation to runoff change was between 3.5% and 67.5%. However, because of lacking long term and high-quality available observed data for glaciers, we could not accurately quantify the process of glacier change and link it to the Budyko framework, which was a limitation in this study and needed future investigation.

Some assumptions adopted by the modified Budyko framework might also cause uncertainties. The sublimation of snow cover in the snow season was ignored in this study. However, the fraction of the sublimation losses in different locations was between 10% and 90% in winter [60]. In a small basin, the snow water could reach confluence quickly with little evaporation loss. In fact, this part of evaporation loss was borne by precipitation, causing the precipitation involved in the evapotranspiration circulation to be lower than the actual value. In addition, we neglected the catchment groundwater storage change in the basin at an interannual scale. These factors might lead to uncertainties in the estimation of available water and available energy.

On the other hand, the analysis framework in this study ignored the relationship between different factors. In addition to the direct effect of different factors on runoff, the coupling and interactions of different related factors also had an indirect effect on runoff [61]. Evidently, this framework could not quantify these indirect effects, which needs further studies to improve the framework.

6. Conclusions

This study analyzed the long-term change in runoff in the upper reaches of the Kaidu River basin from 1960 to 2010 using a modified Budyko framework that considered the snowmelt. Based on the contribution rate and elasticity coefficient of runoff estimated by the modified Budyko framework, the influence of different factors on runoff change was evaluated, and the mechanism between snow ratio change and runoff alteration was discussed. The major findings are as follows:

(1)

From 1960 to 2010, the precipitation, snowfall, air temperature, and runoff exhibited an upward trend in the upper reaches of the Kaidu River basin, while E_P showed a downward trend, but the snowfall ratio did not change significantly. In the two periods divided by 1995, the relative change in runoff was evident.

(2)

From the base period to the change period, the runoff increased by 27.76%. The major factor leading to the increase in runoff was the change in precipitation, whose contribution rate was 81.42%. The contributions of landscape, snow ratio, and E_P change to the runoff alteration accounted for 9.07%, 2.64%, and 1.83%, respectively. The elasticity coefficients of runoff to precipitation, landscape, snow ratio, and E_P in the base period were 1.24, −0.74, 0.22, and −0.24, and that in the change period were 1.23, −0.69, 0.19, and −0.23, respectively.

(3)

Compared with the original Budyko framework, the contribution of precipitation and E_P to the runoff change decreased, and that of landscape increased after considering the snowmelt in the modified Budyko framework. This is because the snowmelt did not participate in the evapotranspiration process, which altered the water-energy condition.

(4)

A larger snow ratio led to a larger runoff elasticity coefficient, suggesting that the variation of runoff increases with the increase in the snow ratio.

The results of this study helped us to better understand the response of runoff to the climate and landscape changes in the upper reaches of the Kaidu River basin. The modified Budyko framework used in this study can be applied to other mountain catchments, and it also needs further improvement to better describe the effect of complex cryosphere processes in cold regions.