Impacts of Precipitation Type Variations on Runoff Changes in the Source Regions of the Yangtze and Yellow River Basins in the Past 40 Years

Abstract

Variations of precipitation type can exert substantial impacts on hydrological processes, yet few studies have quantified the impacts of precipitation type variations on runoff changes in high−altitude regions. In this study, we attempted to examine the potential impacts of precipitation type variations induced by the warming climate on the runoff changes of the source regions of the Yangtze River and Yellow River basins from 1979 to 2018, where the mean elevation is over 4000 m. A modified precipitation type identification method using the wet-bulb temperature, and a runoff change attribution method based on a modified Budyko framework has been applied. Results showed that fluctuations of precipitation contributed to the majority of the runoff variations in the source regions of the Yangtze River basin, which accounted for 51.64%. However, the changes of characteristic parameter n, which indicates the impacts of the underlying surface, explained 56.22% of the runoff changes in the source regions of the Yellow River. It was shown that the trend of shifting from snowfall to rainfall due to a warming climate could result in runoff decreasing, which contributed to 24.06% and 11.29% of the runoff changes in the two source regions, comparatively.

1. Introduction

Precipitation (P) has been recognized as one of the most sensitive variables to runoff changes. The amount, type, intensity, and duration of P can vary the hydrological processes to various extents [1,2]. Particularly, the variations of precipitation type (PT) induced by changing climate, which indicate the changes of the proportion of snowfall and rainfall, have been regarded as one of the driving forces of changes in hydrological processes and distribution of water resources [3,4,5,6]. P falling as snow is beneficial to the formation of snow cover and glaciers in cold regions, which alter the seasonal and annual rhythm of river discharge. On the contrary, P falling as rain can accelerate the melting of glaciers and snow cover and generate surface runoff directly. Therefore, it is of great importance to estimate the impacts of PT variations on hydrological processes.

It is estimated that PT changed significantly at the regional scale, which is particularly true in cold regions, e.g., the Qinghai-Tibetan Plateau (QTP). Results from Wang et al. [7] and Bo et al. [8] showed that the snow fraction significantly decreased on the QTP with a mean slope of −1.1% per decade from 1961 to 2016, while strong temporal variation and spatial heterogeneity were presented. The largest decreases of snow fraction occurred in autumn, followed by summer. Significant decreasing trends were found in the middle and east regions of the QTP. Zhu et al. [9], examined the changes of snow/precipitation ratio in seasonal frozen regions and permafrost regions in the QTP, which found that significant trends in these regions existed. Therefore, the temporal and spatial differences in the variations of PT give rise to the necessity of understanding the regional impacts on hydrological changes, particularly runoff changes.

Quantifying the contribution of anthropogenic factors and climate change to runoff changes is the basis for formulating adaptation strategies of water resources utilization against changing climate [10,11]. A number of study have used various statistical methods, hydrological models and paired watershed methods to examine the impacts of climate change and human activities on runoff changes [12,13,14]. However, the limitations of these studies still exist. For example, there always exists a lack of physical basis and adequate observation data to support the implementation of statistical methods [15]. Hydrological models facilitate rigorous physical interpretation and play an important role in exploring the water cycle. However, it has inherent drawbacks in terms of high input data requirements and complex model validation calibration, which increase the uncertainties of simulation results [15,16,17]. Since the paired watershed method is usually suitable for small watersheds, it could hardly be applied in large watersheds [18,19]. Compared with the above methods, the Budyko framework [20] requires only meteorological data as inputs and only one parameter to be calibrated [19]. Meanwhile, the Budyko method takes into account the water and energy constraints in long-term hydrological processes in a concise manner, estimating the contribution of climate variability and underlying surface changes to runoff variability in a more intuitive way, and on a clear physical basis.

Budyko [20] proposed a framework based on the long-term average water balance theory of a watershed that is controlled by the supply and demand of water (P and evaporation) in the atmosphere. Since 1974, a variety of forms of Budyko equations have been proposed based on the Budyko curve [21,22]. Due to the comprehensiveness and effectiveness of Budyko in studying the effects of climate and underlying surface conditions on hydrological processes, the framework has been widely applied in different regions globally, and has proven its excellent accuracy in large basins and at long-term scales [17,23]. Previous studies on the runoff changes attribution using the Budyko framework-based method, particularly in the QTP, mainly focused on separating the influence of climate and human factors on runoff variations [24,25]. However, few of them attempted to quantitatively evaluate the influences of snowfall on runoff variations. Although Zhang et al. [26] have developed a Budyko framework considering snowfall, they did so for the whole of China and only 1% of the alpine source area basins. There are still few studies on the source regions of alpine basins. Particularly, the differences of these effects between different alpine basins need to be further studied.

In this study, a modified PT identification method using the wet-bulb temperature was applied to investigate the spatio-temporal variations of snowfall proportions. Particularly, the improved Budyko framework considering snowfall, and a new division method of complementarity was used to further evaluate the impacts of snowfall proportions on runoff variations. The contribution rates of environmental factors to runoff changes were quantified. Two adjacent alpine river basins, i.e., the source regions of the Yangtze River (SRLR) and the Yellow River (SRYR), were selected as the study areas.

2. Materials and Methods

2.1. Study Area

The source regions of the Yangtze and Yellow River basins, abbreviated by the SRLR and SRYR, were selected as the study areas, which are located on the QTP with a mean elevation of 4533 m. The SRLR and SRYR, which are usually defined as the regions of the river basin with the outlets at Zhimenda and Tangnaihai hydrological stations, cover an area of 148,000 km² and 122,000 km², respectively (Figure 1). These two regions are regarded as high-altitude river basins, with an elevation ranging from 2569 to 6560 m. However, divergent climate characteristics are found in the two regions. The mean temperature of SRLR and SRYR are 0.16 °C and 0.86 °C, ranging from −1.49 °C in the winter to 2.08 °C in the summer. Nevertheless, the SRLR is relatively drier compared with the SRYR, with a mean annual P of 343.89 mm and 502.72 mm, respectively. It has been recorded that the SRLR and SRYR are getting warmer and wetter (You et al. [26], Zhang et al. [27], Chen et al. [28]), with an increasing rate of the mean annual temperature of 0.53 °C/10a and 0.48 °C/10a, respectively. Meanwhile, the increasing rate of mean annual P of SRLR and SRYR are 23.3 mm/10a and 34.6 mm/10a, respectively (Figure 2).

2.2. Identification of Precipitation Types

In this study, a PT identification method based on wet-bulb temperature (T_w), proposed by Ding et al. [29], was applied to distinguish the proportion of snowfall, sleet, and rainfall of P for the SRLR and SRYR. This method was regarded to have more satisfactory performance than other methods [4,30], e.g., methods based on temperature profile and other atmospheric conditions [31], methods using surface air conditions [32], etc. Particularly, an improvement of 35.42% in the average accuracy of PT identification was found in the QTP compared with T_w [29]. The specific steps of the method are described as follows.

First, the daily T_w was calculated for the date when P occurred:

$${\mathrm{T}}_{\mathrm{w}}={\mathrm{T}}_{\mathrm{a}}-\frac{{\mathrm{e}}_{\mathrm{s}}\left({\mathrm{T}}_{\mathrm{a}}\right)·\left(1-\mathrm{RH}\right)}{0.000643{\mathrm{P}}_{\mathrm{a}}+\Delta }$$

(1)

$${\mathrm{e}}_{\mathrm{s}}\left({\mathrm{T}}_{\mathrm{a}}\right)=6.1078\exp \left(\frac{17.27{\mathrm{T}}_{\mathrm{a}}}{{\mathrm{T}}_{\mathrm{a}}+237.2}\right)$$

(2)

where, ${\mathrm{T}}_{\mathrm{a}}$ is the daily mean temperature, °C, ${\mathrm{e}}_{\mathrm{s}}\left({\mathrm{T}}_{\mathrm{a}}\right)$ is the saturated vapor pressure at ${\mathrm{T}}_{\mathrm{a}}$, kPa, $\mathrm{RH}$ is the relative humidity, the range is 0–1, ${\mathrm{P}}_{\mathrm{a}}$ is the atmospheric pressure, kPa, ∆ is the slope of the saturation water vapor pressure-temperature curve (kPa/°C), $\Delta =4098{\mathrm{e}}_{\mathrm{s}}/{\left({\mathrm{T}}_{\mathrm{a}}+237.2\right)}^2$.

Then, the critical temperatures T_A and T_B are calculated to distinguish between PT: (T_A is Tmin, T_B is Tmax)

$${\mathrm{T}}_{\mathrm{A}}=\{\begin{array}{c}{\mathrm{T}}_0-\Delta \mathrm{S}·\ln \left({\mathrm{e}}^{\frac{\Delta \mathrm{T}}{\Delta \mathrm{S}}}-2·{\mathrm{e}}^{-\frac{\Delta \mathrm{T}}{\Delta \mathrm{S}}}\right),\Delta \mathrm{T}/\Delta \mathrm{S}>\ln 2\\ {\mathrm{T}}_0, \Delta \mathrm{T}/\Delta \mathrm{S}\le \ln 2 \end{array}$$

(3)

$${\mathrm{T}}_{\mathrm{B}}=\{\begin{array}{c}2{\mathrm{T}}_0-{\mathrm{T}}_{\mathrm{A}}, \Delta \mathrm{T}/\Delta \mathrm{S}>\ln 2\\ \hfill {\mathrm{T}}_0, \Delta \mathrm{T}/\Delta \mathrm{S}\le \ln 2\end{array}$$

(4)

where $\Delta \mathrm{T}$ and $\Delta {\mathrm{S}\mathrm{and}\mathrm{T}}_0$ are functions of relative humidity (RH) and elevation (z, km), respectively: $\Delta \mathrm{T}=0.215-0.099·\mathrm{RH}+1.018·{\mathrm{RH}}^2$, $\Delta \mathrm{S}=2.37-1.634·\mathrm{RH}$ and ${\mathrm{T}}_0=-5.87-0.1042·\mathrm{Z}+0.0885·{\mathrm{Z}}^2+16.06·\mathrm{RH}-9.614·{\mathrm{RH}}^2$. $\Delta \mathrm{T}/\Delta \mathrm{S}=\ln 2$, equivalent to RH = 78%. If RH > 78%, the ${\mathrm{T}}_{\mathrm{A}}$ and ${\mathrm{T}}_{\mathrm{B}}$ values are different, otherwise, ${\mathrm{T}}_{\mathrm{A}}$ and ${\mathrm{T}}_{\mathrm{B}}$ are the same. This means that the double threshold method is used when RH > 78% and the single threshold method is used when RH ≤ 78%. The latter means that the probability of sleet occurring at low relative humidity is negligible.

Finally, the magnitude of the wet-bulb temperature is compared with the critical temperature to determine the type of P occurring on that day:

$$\mathrm{PT}=\{\begin{array}{cc}\hfill \mathrm{snow},& {\mathrm{T}}_{\mathrm{w}}\le {\mathrm{T}}_{\mathrm{A}}\hfill \\ \hfill \mathrm{sleet},& {\mathrm{T}}_{\mathrm{A}}<{\mathrm{T}}_{\mathrm{W}}<{\mathrm{T}}_{\mathrm{B}}\hfill \\ \hfill \mathrm{rain},& {\mathrm{T}}_{\mathrm{w}}\ge {\mathrm{T}}_{\mathrm{B}}\hfill \end{array}$$

(5)

2.3. Modified Budyko Framework Considering Snowfall Factor

A modified Buydko framework was proposed by Zhang et al. [26] in this study to quantify the impacts of PT variations on runoff changes in high-altitude regions. A variety of forms of the Budyko framework have been developed in the past decades [22,33]. Particularly, the Budyko equation, developed by Choudhury and Yang [33], Equation (6), has been used to study the relationships between P, evapotranspiration, and runoff at the river basin scale [33].

$$1-\frac{\mathrm{Q}}{\mathrm{P}}={\left[1+{\left(\frac{{\mathrm{ET}}_{\mathrm{p}}}{\mathrm{P}}\right)}^{-\mathrm{n}}\right]}^{-\frac{1}{\mathrm{n}}}$$

(6)

where, P and ET_p represent precipitation, and potential evapotranspiration, respectively, all of which are multiyear averages, and n is a characteristic parameter of the underlying surface.

Regarding the Equation (6), the potential evapotranspiration (ET_p) at each station was estimated using the Penman-Monteith equation recommended by the Food and Agriculture Organization (FAO) [34]. The Ordinary Kriging interpolation method was used to calculate the spatial distribution of precipitation and potential evapotranspiration. Ordinary Kriging interpolation uses the generalized least squares method to unbiased and optimally estimate the above expressions to obtain weight values, the limitation of which is that the interpolation results are affected by spatial stationarity. However, the ordinary Kriging interpolation method can be used for spatial interpolation, and the prediction results and prediction errors can be obtained at the same time, which is convenient for evaluating and analyzing the accuracy of the prediction results, and it has good smoothing effect to reduce the influence of extreme outliers on the overall distribution of the sample.

$${\mathrm{ET}}_{\mathrm{p}}=\frac{0.408\Delta \left({\mathrm{R}}_{\mathrm{n}}-\mathrm{G}\right)+\mathsf{\gamma }\frac{900}{\mathrm{T}+273}{\mathsf{\mu }}_2\left({\mathrm{e}}_{\mathrm{s}}-{\mathrm{e}}_{\mathrm{a}}\right)}{\Delta +\mathsf{\gamma }\left(1+0.34{\mathsf{\mu }}_2\right)}$$

(7)

where ET_p is the potential evapotranspiration (mm day⁻¹), R_n represents the net radiation of the crop surface (MJ m⁻² day⁻¹), G represents the heat flux density of the soil (MJ M⁻² Day⁻¹), T represents the daily mean air temperature (°C), ${\mathsf{\mu }}_2$ is the wind speed at 2 m altitude (m s⁻¹), e_s is the saturated vapor pressure (kPa), e_a is the actual vapor pressure (kPa), e_s − e_a is the insufficient vapor pressure (kPa), Δ is the slope of the vapor pressure (kPa°C⁻¹), and γ is the enthalpy-humidity constant (kPa°C⁻¹).

When the surface temperature rises to a certain level, the snow begins to melt and enters the water cycle process as moisture [35]. Therefore, a variable, snow ratio (SR) was introduced into the Budyko framework to quantify the impacts of PT variations on annual runoff changes. The snow ratio (SR) was calculated by the ratio of the mean annual snowfall to the mean annual P.

$$1-\frac{\mathrm{Q}}{\mathrm{P}}={\left[{\left(1-\mathrm{SR}\right)}^{-{\mathrm{n}}^{\prime }}+{\left(\frac{{\mathrm{ET}}_{\mathrm{p}}}{\mathrm{P}}\right)}^{-{\mathrm{n}}^{\prime }}\right]}^{\frac{1}{-{\mathrm{n}}^{\prime }}}$$

(8)

The sensitivity factors for each of these elements are:

$$\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}=1-\frac{\mathrm{P}-\mathrm{Q}}{\mathrm{P}}\frac{{\left({\mathrm{ET}}_{\mathrm{p}}\right)}^{{\mathrm{n}}^{\prime }}}{{\left[\left(1-\mathrm{SR}\right)\mathrm{P}\right]}^{{\mathrm{n}}^{\prime }}+{\left({\mathrm{ET}}_{\mathrm{p}}\right)}^{{\mathrm{n}}^{\prime }}}$$

(9)

$$\frac{\partial \mathrm{Q}}{\partial {\mathrm{ET}}_{\mathrm{p}}}=-\frac{\mathrm{P}-\mathrm{Q}}{{\mathrm{ET}}_{\mathrm{p}}}\frac{{\left[\left(1-\mathrm{SR}\right)\mathrm{P}\right]}^{{\mathrm{n}}^{\prime }}}{{\left[\left(1-\mathrm{SR}\right)\mathrm{P}\right]}^{{\mathrm{n}}^{\prime }}+{\left({\mathrm{ET}}_{\mathrm{p}}\right)}^{{\mathrm{n}}^{\prime }}}$$

(10)

$$\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}=\frac{\mathrm{P}-\mathrm{Q}}{1-\mathrm{SR}}\frac{{\left({\mathrm{ET}}_{\mathrm{p}}\right)}^{{\mathrm{n}}^{\prime }}}{{\left[\left(1-\mathrm{SR}\right)\mathrm{P}\right]}^{{\mathrm{n}}^{\prime }}+{\left({\mathrm{ET}}_{\mathrm{p}}\right)}^{{\mathrm{n}}^{\prime }}}$$

(11)

$$\frac{\partial \mathrm{Q}}{\partial {\mathrm{n}}^{\prime }}=-\frac{\mathrm{P}-\mathrm{Q}}{{\mathrm{n}}^{\prime }}\left(\frac{\ln \left[{\left(1-\mathrm{SR}\right)}^{{\mathrm{n}}^{\prime }}{\mathrm{P}}^{{\mathrm{n}}^{\prime }}+{\left({\mathrm{ET}}_{\mathrm{p}}\right)}^{{\mathrm{n}}^{\prime }}\right]}{{\mathrm{n}}^{\prime }}-\frac{{\left[\left(1-\mathrm{SR}\right)\mathrm{P}\right]}^{{\mathrm{n}}^{\prime }}\ln \left[\left(1-\mathrm{SR}\right)\mathrm{P}\right]+{\left({\mathrm{ET}}_{\mathrm{p}}\right)}^{{\mathrm{n}}^{\prime }}\ln \left({\mathrm{ET}}_{\mathrm{p}}\right)}{{\left[\left(1-\mathrm{SR}\right)\mathrm{P}\right]}^{{\mathrm{n}}^{\prime }}+{\left({\mathrm{ET}}_{\mathrm{p}}\right)}^{{\mathrm{n}}^{\prime }}}\right)$$

(12)

Equation (9) can be abbreviated as Q = f (P, ET_p, SR, n′), and the runoff elasticity coefficient for each independent variable x can then be calculated as follows:

$${\mathsf{\epsilon }}_{\mathrm{x}}=\frac{{\partial }_{\mathrm{Q}}}{{\partial }_{\mathrm{x}}}\times \frac{\mathrm{x}}{\mathrm{Q}}$$

(13)

where ${\mathsf{\epsilon }}_{\mathrm{x}}$ is the elasticity coefficient, and x represents each influence factor (P, ET_p, SR, n′). A positive elasticity coefficient for a variable indicates that Q increases as the variable increases, while a negative elasticity coefficient indicates that Q decreases as the variable increases.

A complementary method, proposed by Zhou et al. [36], has been applied to examine the relative contributions of environmental factors to runoff changes. Particularly, the impacts of PT variations were quantified by calculating the relative contribution of SR on annual runoff changes, when factors including P, evapotranspiration, and snowfall ratio were considered. The annual runoff changes could be calculated using Equation (14).

$$\begin{array}{cc}\hfill \Delta \mathrm{Q}=\mathsf{\alpha }[{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}\right)}_1\Delta \mathrm{P}& +{\left(\frac{\partial \mathrm{Q}}{\partial {\mathrm{ET}}_{\mathrm{p}}}\right)}_1\Delta {\mathrm{ET}}_{\mathrm{p}}+{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}\right)}_1\Delta \mathrm{SR}+{\mathrm{P}}_2\Delta \left(\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}\right)+{\mathrm{ET}}_{\mathrm{p},2}\Delta \left(\frac{\partial \mathrm{Q}}{\partial {\mathrm{ET}}_{\mathrm{p}}}\right)+{\mathrm{SR}}_2\Delta \left(\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}\right)]\hfill \\ \hfill & +\left(1-\mathsf{\alpha }\right)\left[{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}\right)}_2\Delta \mathrm{P}+{\left(\frac{\partial \mathrm{Q}}{\partial {\mathrm{ET}}_{\mathrm{p}}}\right)}_2\Delta {\mathrm{ET}}_{\mathrm{p}}+{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}\right)}_2\Delta \mathrm{SR}+{\mathrm{P}}_1\Delta \left(\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}\right)+{\mathrm{ET}}_{\mathrm{p},1}\Delta \left(\frac{\partial \mathrm{Q}}{\partial {\mathrm{ET}}_{\mathrm{p}}}\right)+{\mathrm{SR}}_1\Delta \left(\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}\right)\right]\hfill \end{array}$$

(14)

where, α is a weighting factor taking values in the range [0, 1] and following the recommendation of [36], using α = 0.5 is quite close to the best estimate of the magnitude of the climate effect and basin effect, with subscripts 1 and 2 denoting two different stages considering the changes of climatic factors, respectively. Hence, the relative contributions of P, ETp, SR, and n’ to the changes in Q are:

$$\mathrm{C}\_\left(\mathrm{P}\right)=\mathsf{\alpha }\left[{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}\right)}_1\Delta \mathrm{P}\right]+\left(1-\mathsf{\alpha }\right)\left[{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}\right)}_2\Delta \mathrm{P}\right]$$

(15)

$$\mathrm{C}\_\left({\mathrm{ET}}_{\mathrm{p}}\right)=\mathsf{\alpha }\left[{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}\right)}_1\Delta {\mathrm{ET}}_{\mathrm{p}}\right]+\left(1-\mathsf{\alpha }\right)\left[{\left(\frac{\partial \mathrm{Q}}{\partial {\mathrm{ET}}_{\mathrm{p}}}\right)}_2\Delta {\mathrm{ET}}_{\mathrm{p}}\right]$$

(16)

$$\mathrm{C}\_\left(\mathrm{SR}\right)=\mathsf{\alpha }\left[{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}\right)}_1\Delta \mathrm{SR}\right]+\left(1-\mathsf{\alpha }\right)\left[{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}\right)}_2\Delta \mathrm{SR}\right]$$

(17)

$$\mathrm{C}\_\left({\mathrm{n}}^{\prime }\right)=\mathsf{\alpha }\left[{\mathrm{P}}_2\Delta \left(\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}\right)+{\mathrm{ET}}_{\mathrm{p},2}\Delta \left(\frac{\partial \mathrm{Q}}{\partial {\mathrm{ET}}_{\mathrm{p}}}\right)+{\mathrm{SR}}_2\Delta \left(\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}\right)\right]+\left(1-\mathsf{\alpha }\right)\left[{\left(\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}\right)}_2\Delta \mathrm{SR}+{\mathrm{P}}_1\Delta \left(\frac{\partial \mathrm{Q}}{\partial \mathrm{P}}\right)+{\mathrm{ET}}_{\mathrm{p},1}\Delta \left(\frac{\partial \mathrm{Q}}{\partial {\mathrm{ET}}_{\mathrm{p}}}\right)+{\mathrm{SR}}_1\Delta \left(\frac{\partial \mathrm{Q}}{\partial \mathrm{SR}}\right)\right]$$

(18)

2.4. Statistical Analysis

The Pettitt nonparametric test, which has been widely used in hydrologic and climatic fields [37,38,39], is applied to identify points of abrupt change point of the annual runoff series [40]. The state variables of the Pettitt nonparametric test were calculated using the following equations.

$${\mathrm{U}}_{\mathrm{k},\mathrm{N}}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{k}}\mathrm{sgn}\left({\mathrm{X}}_{\mathrm{i}}-{\mathrm{X}}_{\mathrm{j}}\right)$$

(19)

$$\mathrm{K}=\underset{1\le \mathrm{k}\le \mathrm{N}}{\max }\left|{\mathrm{U}}_{\mathrm{k},\mathrm{N}}\right|$$

(20)

$$\mathrm{P}\cong 2\exp [-6{\mathrm{K}}^2/({\mathrm{N}}^3+{\mathrm{N}}^2\left)\right]$$

(21)

where, X_i and X_j are the sequence values of time i and j, respectively, N is the total number of elements in the series, and the breakpoint K is determined by the maximum value of U_K,N, and the significance level of this study was set at 0.05.

2.5. Data

Meteorological data, including the T, P, relative humidity, sunshine hours, wind speed, and air pressure, were obtained from The China Meteorological Data Service Center (http://www.cma.gov.cn/, accessed on 18 March 2021). Particularly, the China Surface Climate Data Daily Value Dataset (V3.0) was selected due to its high quality and adequate amount of data in this study, in which meteorological data from 1979 to 2018 were extracted. In this study, 29 stations that fall in or surrounding the two sources regions were used, where the basic information of the stations was listed in Supplementary Material (Table S1). Monthly discharge data of Zhimenda and Tangnaihai stations from 1979 to 2018 were derived from the hydrologic data yearbook, which represents the outlets of the source regions of the Yangtze and Yellow river basin, respectively. Basic information about the hydrological stations were demonstrated in the Supplementary Materials. Additionally, the digital elevation model (DEM) data with a spatial resolution of 90 m was obtained from the SRTM (Shuttle Radar Topography Mission), which was derived from the Geospatial Data Cloud (http://www.gscloud.cn/, accessed on 5 February 2022).

3. Results

3.1. Variations of the Precipitation Type of the SRLR and SRYR

The PT identification results showed that the average annual snowfall amount of the SRLR and SRYR increased at a rate of 5.62 mm/10a and 2.43 mm/10a, respectively, with a mean value of 113.07 mm and 113.83 mm from 1979 to 2018 (Figure 3). However, the proportions of the PTs in the two source regions have experienced significant changes in the past four decades. The proportions of snowfall, indicated by SR, showed a decreasing trend with a rate of 2.50%/10a and 0.24%/10a for the SRLR and SRYR, respectively (Figure 4). The proportions of snowfall in P of the SRLR and SRYR in the 1980s were 42.75% and 27.82% respectively. However, these proportions decreased to 28.85% and 16.66% in the 2010s.

Strong spatial heterogeneity was presented in terms of SR and SR trends in the SRLR and SRYR (Figure 5). Results showed that SR in SRLR is higher in the northern and eastern regions and lower in the southern regions, while SR in SRYR is higher in the western regions and lower in the northern, eastern, and southern regions. There was a significant decreasing trend in SR in the northern and central regions of SRLR and a significant decreasing trend in SR in the western and northern regions of SRYR. In addition, we divided the watersheds at 500 m intervals with reference to the studies of Liu et al. [41] and Guo et al. [42] in the Tibetan Plateau region and found that the variations of SR were found to be elevation-dependent (Figure 6). Results showed that the SR trend of SRLR was stable below 4500 m and fluctuated greatly above 4500 m, while the SR trend of SRYR was more stable below 3500 m and fluctuated in the range above 3500 m. The mean values of SR trend also showed a weak upward trend with the increasing of altitude, with an SR trend of −0.83 for SRLR below 4500 m and an SR trend of −0.20 for regions above 4500 m, and an SR trend of −0.15 for SRYR below 3500 and an SR trend of −0.14 for regions above 3500 m. This trend may be related to the altitude dependence of warming on the Tibetan Plateau, that is, warming is significantly higher at higher altitudes than at lower altitudes, and the temperature rise at high altitudes tends to increase with increasing altitude.

3.2. Variations and Attributions of Runoff Changes of the SRLR and SRYR

This study analyzed the evolutionary characteristics of runoff based on measured annual runoff data from the study basin. The results showed that the long-term average runoff depths for SRLR and SRYR are 99.95 mm and 163.40 mm, respectively. According to the 5-year moving average curve, the SRLR runoff depth increased significantly at a rate of 7.41 mm/10a, while the SRYR runoff decreased at a rate of 7.64 mm/year (Figure 7). The Pettitt method was used to analyze the runoff time series data of the Zhimenda and Tangnaihai hydrological station from 1979 to 2018. Figure 8 displays the trend and catastrophe analysis results of runoff in both basins, with SRLR runoff depth mutating in 2004 and SRYR runoff depth mutating in 1990. Therefore, the runoff series can be divided into two periods, the baseline period and the change period, according to the year of the abrupt change. The average runoff depths for the SRLR baseline period (1979–2004) and the change period (2005–2018) were 90.58 mm and 117.35 mm. A comparison of the two periods indicated a more significant increase in SRLR runoff. Similarly, the average runoff depths for the SRYR baseline period (1979–1990) and the change period (1991–2018) were 190.31 mm and 151.87 mm, indicating a decrease of 38.44 mm.

It was shown that P was the dominant factor affecting SRLR runoff (Figure 9), with a contribution of 51.64%, and n′ has the largest contribution to the reduction of runoff in the SRYR, accounting for 56.22%. ET_p has the least influence on the runoff in the two study basins, with contributions of 2.61% and 1.26%. The annual decrease in SR contributed −24.06% to SRLR runoff and −11.29% to SRYR runoff. The different roles played by climatic factors in the attribution of runoff changes are consistent with the results of their trend analysis.

The elasticity coefficients of runoff to P, SR, ET_p, and n′ were calculated by the Budyko formula. According to

Table 1. Characteristics of hydrometeorological number elements and elasticity coefficients for SRLR and SRYR.

							Elasticity Coefficient
Watershed	Period	ETp/mm	Q/mm	P/mm	SR	n′	βET	βp	βSR	βn′
SRLR	1979–2018	701.5	99.95	343.89	0.36	1.32	−2.01	1.43	1.15	−0.86
	1979–2004	703.13	90.58	323.20	0.38	1.25	−2.12	1.44	1.29	−0.95
	2005–2018	698.48	117.35	383.32	0.34	1.15	−2.23	1.53	1.14	−1.07
SRYR	1979–2018	678.54	163.40	558.27	0.24	1.00	−1.49	1.93	0.47	−1.61
	1979–1990	681.30	190.31	517.31	0.27	0.99	−2.31	1.62	0.85	−1.13
	1991–2018	677.36	151.87	575.82	0.23	1.04	−2.84	2.10	0.84	−1.81

, the differences in hydrological variables and the estimated watershed parameter n′ between the two stages were compared, and the absolute value of the elasticity coefficient of ET_p to runoff was the largest. For SRLR, a 1% increase in P would result in a 1.43–1.53% increase in runoff, and for SRYR, a 1% increase in P would result in a 1.62–2.10% increase in runoff. The elasticity coefficients for runoff to ET_p and n′ were both negative, with a 1% increase in ET_p leading to a 2.01–2.23% decrease in runoff and a 1% increase in n′ leading to a 0.86–1.07% decrease in the runoff for SRLR. For SRYR, a 1% increase in ET_p would result in a 1.49–2.84% decrease in the runoff, and a 1% increase in n’ would result in a 1.13–1.81% decrease in the runoff. There was a positive correlation between SR and runoff in both catchments, with a 1% increase in SR in SRLR leading to a 1.14–1.29% increase in runoff and a 1% increase in SR in SRYR leading to a 0.47–0.85% increase in runoff. Comparing the elastic coefficients of all the factors in the base period and the change period, the absolute values of all the factors except SR increased, indicating that the sensitivity of these factors increased. The elasticity coefficient of SR decreased in both basins, indicating that the sensitivity of runoff to SR decreased.

Based on the calculation of the elasticity coefficients, the complementary method was used to distinguish the influence degree of each factor on the runoff changes.

Table 2. Changes in runoff due to changes in factors estimated by Budyko (mm).

	ETp	P	SR	n′	$\Delta \mathbf{Q}\mathbf{s}\mathbf{i}\mathbf{m}$	$\Delta \mathbf{Q}\mathbf{o}\mathbf{b}\mathbf{s}$	Error
SRLR	1.31	25.88	−12.06	10.87	26.00	26.77	2.88%
SRYR	1.35	33.57	−12.14	−60.44	−37.66	−38.45	2.05%

summarized the changes in runoff caused by each factor and it can be seen that the simulated runoff variation $\Delta \mathrm{Qsim}$ has an error of less than 5% compared to the measured runoff $\Delta \mathrm{Qobs}$, indicating that the Budyko equation considering SR performed well in the attribution analysis of runoff. As shown in Figure 9, P is the dominant factor affecting SRLR runoff.

It should be noted that the effect of snowfall on runoff is important, and the contribution rate of SR to runoff were negative in both basins. The study found that the phase shift from snow to rain significantly reduced the average flow. In addition to the effects of P and ET_P, the changes in runoff were also largely attributed to the changes in n′. Overall, from 1979 to 2018, underlying surface changes contributed positively to both SRLR and SRYR runoff.

4. Discussion

4.1. Precipitation Type Identification in High-Altitude Regions

P in the SRLR and SRYR has increased significantly over the last four decades, with a maximum rate of increase of 34.6 mm/10a. Previous studies have shown that annual P has increased in most regions of the Tibetan Plateau since the 1960s [43]. Yan et al. [44] found that the annual mean temperature increase at over 2000 m on the Tibetan Plateau was greater than the global rate of increase, indicating that the surface temperatures rose faster at higher altitudes. The significant rising of temperature could play a key role in affecting the snowfall and its proportions in precipitation [45]. Deng et al. [46] showed that 59.1% of stations on the Tibetan Plateau showed an increasing trend in snowfall, which was consistent with the results found in this study. However, the impacts of warming on snowfall may be reflected in the proportions to P, i.e., the SR, as can be seen in Figure 4. The results showed that the SR of SRLR and SRYR decreased at 2.50%/10a and 0.24%/10a, respectively. It has been found that more proportions of precipitation fall as rain than that as snow in cold regions, e.g., the QTP [7], the western United States [47] and New England [48], making SR decline a global problem. The results were similar to those of Wang et al. [49] who found that SR declined at a rate of −0.164/10a throughout the year except for November on the QTP. Spatially, the SR declined more significantly in the southwestern and central regions of the QTP. It can been seen from the results that the rate of SR decline in the SRLR region in this study was ten times greater than that of the SRYR, which is located in the central part of the Tibetan Plateau (Figure 1), and thus its SR decline rate was faster than that of the SRYR.

Results also showed that the changes of SR were closely related to altitude, with SRLR and SRYR taking 4500 m and 3500 m as the critical thresholds, respectively, below which SR trends are more stable, and SR trends increase with altitude. Deng et al. [46] showed that the effect of altitude on snowfall trends is also significant. Li et al. [50] found that P occurs as snowfall in areas above 3500 m in Xinjiang. Guo et al. [51] in the Tianshan region found that SR decreased significantly below 2500 m in altitude between 1961 and 2010, and at high altitudes above 3500 m The decrease in SR at high altitudes above 3500 m was small.

4.2. Impacts of Precipitation Type Variations on Runoff Changes at River Basin Scale

Reliable attribution of runoff changes is fundamental to our understanding of the hydrological cycle, and essential for decision makers to manage water resources sustainably [52]. Many scholars have analyzed the contribution of SRLR and SRYR to the long-term variation in runoff using a variety of methods. For instance, Jia et al. [53] studied the runoff characteristics of SRLR from 1964 to 2018 and found that the contribution of P to the changes of runoff was 67.5%. Zheng et al. [54] found that land use and land cover changes dominated the effect of SRYR from 1990–2000, contributing about 70% of the runoff changes. Snowfall is an important source of irrigation and drinking water in alpine regions. Due to the different trends in the phase transformation of P, the sources of runoff recharge change significantly, including the types of possible recharge sources and their respective proportions. The effects of rainfall or snowfall on runoff could not therefore be considered simply. The results of this study, which took snowfall into account, show (Figure 9) that SR contributes 24.06% and 11.29% to SRLR and SRYR runoff variations, respectively, in the context of the Budyko hypothesis, indicating that SR has a significant influence on runoff, which has often been overlooked in previous studies. Changes in SR may lead to future regional water shortages [55], and studies in the Alpine basin have shown that declines in SR can lead to summer water shortages or flood disasters [56]. The decline in SR will lead to a delay in the start of the snow season and an earlier end of the season, and will also affect the accumulation and melting processes of glaciers. Thus influenced by changes in PTs, more uncertainty is expected in runoff predictions for snowmelt and P recharge rivers, and changes in water flow and water resources will be more complex. Therefore, clarifying the mechanisms influencing SR changes and the sensitivity of runoff to SR changes will be important for future water resources management policies.

4.3. Impacts of Underlying Surface Conditions on Runoff Changes in High-Altitude Regions

Many studies have shown that parameter n’ is a function of human activities, soil properties, topography, and land-cover vegetation [57,58,59,60]. According to Zhang et al. [61], the SRYR grassland area has decreased by approximately 10% since 1990, and the sand area has increased by approximately 4%. Liu et al. [62] found an increase in grassland and wetland area in SRLR, and a significant decrease in grassland cover in SRYR over the nearly 30-year period from the late 1970s to 2012, with a more dramatic overall change in SRYR. The increase in vegetation area may intercept more P and will increase actual evapotranspiration, negatively impacting runoff [63]. The QTP has a wide distribution of cryosphere, complex topography, and ecosystems, and climate warming has caused permafrost degradation and glacier retreat in the study area [64,65]. Jin et al. [66] found that from 1980 to 2016, the multi-year permafrost area in the SRYR shrank by about 2000 km² and the active layer thickness expanded at a rate of 35 cm/10a. Shi et al. [67] showed that the active layer thickness (ALT) in the SRLR increased at a rate of 4.2 cm/10a from 1982 to 2015. These processes can have a considerable impacts on hydrological processes. The results of this paper show that the main factor of SRYR runoff reduction is parameter n′, which may be due to the drastic changes of the underlying surface and the obvious degradation of frozen soil. The comparison showed that the rate of permafrost degradation is much greater in the SRYR than in the SRLR, and that the overall magnitude of land use changes is more dramatic, which may lead to the parameter n’ being the dominant factor in runoff changes in the SRYR rather than climate.

5. Conclusions

In this study, a wet-bulb temperature-based PT identification method was presented to detect the changes of the proportions of snowfall and rainfall from 1979 to 2018. Afterward, a modified Budyko framework considering the snowfall ratio was applied to examine the impacts of environmental factors on annual runoff changes. The main conclusions are presented as follows.

(1)

The proportions of snowfall decreased at 2.50%/10a and 0.24%/10a for the source regions of the Yangtze and Yellow River basins, respectively, while the snowfall amount increased due to the increase of P amount.

(2)

Strong spatial heterogeneity of the SR variations has been found in the two regions. Particularly, SR trends for SRLR and SRYR are more stable below the 4500 m and 3500 m thresholds, respectively, and more volatile above the thresholds.

(3)

The impacts of PT variations on runoff changes in the SRLR were relatively larger than that in the SRYR, with a relative contribution of −24.06% and −11.29%, respectively. Particularly, the impacts of underlying surfaces on annual runoff changes were not negligible in the SRYR, since a significant increasing trend was found in the vegetation recovery in this region. Apart from that, the variations of PT and cryospheric conditions will exert substantial impacts on the inter-annual and intra-annual runoff changes.