A New Assessment of Hydrological Change in the Source Region of the Yellow River

Hydrological responses to climate change are a widely concerning question, particularly for the source region of the Yellow River (SRYR), which is sensitive to climate change and is widely underlain by frozen ground. In considering climate change impacts on catchment properties, the traditional separation approach based on the Budyko framework was modified to identify and quantify the climatic causes of discharge changes. On the basis of the decomposition method, the traditional separation method and the modified separation method were used to analyse the discharge change in the SRYR. Using the observed annual maximum frozen depth (MFD) to indicate the frozen ground level, the impacts of frozen-ground degradation on the discharge change were further considered using the modified separation method. Our results show that the traditional separation approach underestimated climate-induced discharge change; over the past half-century, the discharge change in the SRYR has been primarily controlled by climate change. Increasing air temperature is generally a negative force on discharge generation; however, it also causes frozen ground to degrade—a positive factor for discharge generation. Such conflicting effects enhance the uncertainty in assessments of hydrological responses to climate change in the sub-basins with widely distributed permafrost.


Introduction
Over the past five decades, the climate has changed dramatically at the global scale [1], resulting in long-term impacts on global rivers [2][3][4][5].The response of hydrological basins in the Tibetan Plateau as a result of climate change is a major concern for the management of water resources, particularly in China.The water balance of headwater basins in the Tibetan Plateau is influenced by climate change at different levels [6][7][8].Originating from the Tibetan Plateau, the Yellow River is the second longest river in China.Over the past five decades, the temperature has increased rapidly in the region, accompanied by significant changes in the river's streamflow that were primarily exhibited as a sharp decrease for many sub-basins during the 1990s [9,10].
The source region of the Yellow River (SRYR) is situated in the eastern Tibetan plateau in a large basin above the Tangnaihai hydrological station (Figure 1), which contributes about 35% of the total annual runoff but has an area that is only 16% (121,972 km 2 ) of the total Yellow River Basin [11].In the 1990s, streamflow from the SRYR significantly decreased with increasing cases of flow cutoff [10,12,13].The causes of the decreased flow within the SRYR and other sub-basins of the Yellow River have been extensively investigated in the literature.However, there are disagreements as to the causes Water 2018, 10, 877 2 of 18 of the observed decreases in streamflow.Much research based on data analysis and process-based hydrological modelling has suggested that climate change was the primary cause of the decrease in streamflow [10,12,[14][15][16][17][18][19].For instance, using the non-parametric Mann-Kendall test and moving t-test, Zheng [12] analysed the streamflow change in the SRYR and concluded that the decrease in streamflow was mainly caused by wet-season rainfall decrease.By analysing discharge data and meteorological data, Hu [16] found that a decrease in precipitation and an increase in temperature resulted in the streamflow decrease in the SRYR.Using a process-based land-surface hydrological model (VIC), to simulate the streamflow change in the SRYR, Cuo [14] and Meng [19] indicated that climate change (changes in precipitation and evapotranspiration) was the primary cause for the streamflow change.Other research based on the Budyko framework proposed that the decrease in discharge was mainly caused by changes in catchment properties [20,21].For example, Zheng [20] found that the decrease in streamflow in the 1990s was a result of climate change and land-use changes, where the land-use changes were more important and contributed about 70% of the total change in streamflow.Zhao [21] suggested that 65% of the streamflow change in the upper catchment of the Jimai hydrological station could be attributed to the effects of human activities.This conclusion is potentially at odds with the very low population density in the area above the Tangnaihai hydrological station (about 6 people/km 2 ; 2003 census data) and in the area above the Huangheyan station (0.34 people/km 2 ) [15].From 1990 to 2000, the change in land use in the SRYR was generally less than 5%, although a few sites exhibited 5-15% change [22].Therefore, it is not easy to attribute the negligible impact of human activities and land-use changes to more than 50% of the change in the discharge of the SRYR.
A poorly understood factor of the SRYR is the hydrological impacts of changes in the frozen ground (permafrost and seasonally frozen ground).Permafrost is defined as a subsurface region in which the temperature is below 0 • C for two or more consecutive years [23].Frozen ground in the SRYR has undergone significant degradation since the 1980s [24], which has been manifested as a shrinking permafrost area, a decreasing maximum frozen depth (MFD), and a reducing freezing duration of the seasonally frozen ground [25,26].
In considering frozen-ground degradation, recently Wang [27] used a climate elasticity method combined with a sequential average method to reassess the discharge change in the SRYR and concluded that the 32.6% discharge decrease in the 1990s was caused by frozen-ground degradation while the 12.4% discharge decrease was caused by land-use changes.These results seem inconsistent with existing studies indicating that permafrost degradation could enlarge baseflow [28][29][30][31][32][33] and increase the runoff discharge [34,35].Nevertheless, Wang [27] showed the potential of considering frozen-ground degradation in the Budyko framework to remedy the overestimation of land-use changes in a previous study on the hydrological changes in the SRYR.Meanwhile, the catchment-specific parameter in the Budyko framework could also be dependent on the climate condition [36][37][38][39]; however, this effect was not accounted for in most of the investigations on the SRYR [20,21,27].
As described above, previous assessments of hydrological response in the SRYR on the basis of the Budyko framework were inconsistent with conclusions obtained by other methods, and climate change impacts on the catchment-specific parameter were sparsely considered in previous studies.Frozen ground is a climate-sensitive landscape that is widely distributed in the SRYR, but its hydrological impacts on the SRYR are still disputed.Accordingly, this study aims to (1) modify the traditional separation method in considering both climate change and frozen-ground degradation impacts on the catchment-specific parameter, (2) reassess the causes of long-term discharge changes in the SRYR, and (3) provide a better understanding of hydrological responses to climate change and frozen-ground degradation in the SRYR.A catchment-specific parameter, on the basis of the Budyko framework, was applied to represent the properties of the drainage basins.Correlation analysis was performed to reveal the relationships between the time-varying catchment parameter and potential drivers of hydrologic change, such as climatic forcing and frozen-ground degradation.

Catchments and Sub-Basins
The Yellow River originates from the Tibetan plateau, flows across the north of China, and eventually empties into the Baohai Sea.The SRYR is situated in the northeastern Tibetan plateau, the region between 95 • 30 E-103 • 30 E and 32 • N-36 • 40 N, which is a transition zone of seasonally frozen ground and permafrost (Figure 1).Three hydrological stations, Jimai, Maqu, and Tangnaihai, monitor the mainstream of the SRYR from upstream to downstream.The SRYR was divided into three sub-basins in this study according to the three hydrological stations: the Jimai (JM) Basin, the Jimai-Maqu (MQ) Basin, and the Maqu-Tangnaihai (TNH) Basin, from upstream to downstream, respectively.Two major lakes, Gyaring and Ngoring, are located within the JM sub-basin.The JM and TNH sub-basins are relatively dry, with a high elevation and cold climate, whereas the MQ sub-basin, the major runoff generation area of the SRYR, has a relatively humid and warm climate and contributes 51% of the runoff to the SRYR [19].
The distribution of frozen ground in the 1980s was obtained from the Map of Snow, Ice, and Frozen Ground in China (1988), published by the Environmental and Ecological Science Data Center for West China, National Natural Science Foundation of China [40].The classifications of the frozen ground were further simplified into plateau permafrost, alpine permafrost, and seasonally frozen ground.As shown in Figure 1, permafrost occupies 80% and 46% of the basin areas of JM and TNH, respectively, whereas in the MQ, only 23% of the basin area is occupied by permafrost.

Data Processing
Daily meteorological data from 1961 to 2013 at 18 meteorological stations inside and close to the SRYR were collected from the China Meteorological Administration (CMA), including precipitation, air temperature, wind speed at 10 m height, relative humidity, and sunshine duration.The FAO Penman-Monteith equation [41] was employed to calculate potential evapotranspiration.The annual Water 2018, 10, 877 4 of 18 precipitation intensity (I) was defined as the ratio of annual precipitation to the annual number of raining days.The annual precipitation (P) (mm), annual precipitation intensity (I) (mm/d), annual mean temperature (T) ( • C), and annual potential evapotranspiration (E 0 ) (mm) from 1961 to 2013 were derived for the meteorological stations.Then, the inverse-distance weighted technique (IDW) was applied to gridded data of P, E 0 , T, and I at a resolution of 0.01 Monthly streamflow data from 1961 to 2013 at the Jimai, Maqu and Tangnaihai hydrological stations were obtained from the Yellow River Conservancy Commission (YYCC).Firstly, these monthly streamflow data were translated into annual runoff volume data for each hydrological station.Annual runoff volume data for the JM and TNH stations were directly used to represent annual runoff volume generated in JM and the SRYR, respectively.The MQ-JM difference and TNH-MQ difference in annual runoff volume were used to represent the individual runoffs in MQ and TNH, respectively.Finally, the obtained runoff volume data for each basin were divided by their respective basin area to obtain the annual runoff depth, which was defined as discharge in following analyses.
Observed daily frozen depth data for 10 meteorological stations (Figure 1) were collected from the CMA.Firstly, the daily frozen depth was used to calculate the monthly mean frozen depth for each year; then the maximum monthly mean frozen depth value for each year was defined as the annual MFD.Daily frozen depth data were available for 1961 to 2004, but some missing data existed, and the length of data varied for different stations.The variation patterns of the annual MFD and annual mean air temperature (MAT) (T) for the three sub-basins are shown in Figure 2.
Water 2018, 10, x FOR PEER REVIEW 4 of 18 mean temperature (T) (°C), and annual potential evapotranspiration (E0) (mm) from 1961 to 2013 were derived for the meteorological stations.Then, the inverse-distance weighted technique (IDW) was applied to gridded data of P, E0, T, and I at a resolution of 0.01° × 0.01°.Monthly streamflow data from 1961 to 2013 at the Jimai, Maqu and Tangnaihai hydrological stations were obtained from the Yellow River Conservancy Commission (YYCC).Firstly, these monthly streamflow data were translated into annual runoff volume data for each hydrological station.Annual runoff volume data for the JM and TNH stations were directly used to represent annual runoff volume generated in JM and the SRYR, respectively.The MQ-JM difference and TNH-MQ difference in annual runoff volume were used to represent the individual runoffs in MQ and TNH, respectively.Finally, the obtained runoff volume data for each basin were divided by their respective basin area to obtain the annual runoff depth, which was defined as discharge in following analyses.
Observed daily frozen depth data for 10 meteorological stations (Figure 1) were collected from the CMA.Firstly, the daily frozen depth was used to calculate the monthly mean frozen depth for each year; then the maximum monthly mean frozen depth value for each year was defined as the annual MFD.Daily frozen depth data were available for 1961 to 2004, but some missing data existed, and the length of data varied for different stations.The variation patterns of the annual MFD and annual mean air temperature (MAT) (T) for the three sub-basins are shown in Figure 2. As shown in Figure 2, the MFD values for JM and TNH varied from 120 to 300 cm; however, MFD values for MQ were generally lower than 120 cm.As air temperature increased, the MFD exhibited a decreasing trend at most of the observation sites, particularly in the period after 1980.The As shown in Figure 2, the MFD values for JM and TNH varied from 120 to 300 cm; however, MFD values for MQ were generally lower than 120 cm.As air temperature increased, the MFD exhibited a decreasing trend at most of the observation sites, particularly in the period after 1980.The MFDs at each station in and around the same sub-basin showed a similar trend; thus the station with the longest record was selected to account for the representative MFD data in each sub-basin.Therefore, M5, M10, and M14 were selected for JM, MQ, and TNH, respectively.In the MQ sub-basin, M17 also had the longest record but showed little variation compared to the other sites, because it was located at alpine permafrost.Therefore, M17 was not selected to represent the MFD of MQ.The locations of M5, M10, and M14 are shown in Figure 1.The station M14 was not located in the MQ sub-basin, but it exhibited the same variation trend in MFD as that observed at stations in the MQ basin (M12, M13, and M16) with shorter records.

Methods
To assess the hydrological response of the SRYR and its three sub-basins to climate change and frozen-ground degradation, a new analysis method was used in this study.First, the Budyko framework was applied to investigate the shift in mean annual water balance in different periods.Second, using a decomposition method, the discharge change was preliminarily partitioned into two components on the basis of the traditional separation approach (Figure 3a): the change along a Budyko curve (∆Q 0 ) and the change induced by catchment-specific parameter change (∆Q w ).Finally, the impacts of climate change and the degrading frozen ground on ∆Q w were further identified with the modified separation approach.Accordingly, the total climatic contribution to the discharge change (∆Q c ) was identified (Figure 3b).MFDs at each station in and around the same sub-basin showed a similar trend; thus the station with the longest record was selected to account for the representative MFD data in each sub-basin.Therefore, M5, M10, and M14 were selected for JM, MQ, and TNH, respectively.In the MQ sub-basin, M17 also had the longest record but showed little variation compared to the other sites, because it was located at alpine permafrost.Therefore, M17 was not selected to represent the MFD of MQ.The locations of M5, M10, and M14 are shown in Figure 1.The station M14 was not located in the MQ sub-basin, but it exhibited the same variation trend in MFD as that observed at stations in the MQ basin (M12, M13, and M16) with shorter records.

Methods
To assess the hydrological response of the SRYR and its three sub-basins to climate change and frozen-ground degradation, a new analysis method was used in this study.First, the Budyko framework was applied to investigate the shift in mean annual water balance in different periods.Second, using a decomposition method, the discharge change was preliminarily partitioned into two components on the basis of the traditional separation approach (Figure 3a): the change along a Budyko curve (ΔQ0) and the change induced by catchment-specific parameter change (ΔQw).Finally, the impacts of climate change and the degrading frozen ground on ΔQw were further identified with the modified separation approach.Accordingly, the total climatic contribution to the discharge change (ΔQc) was identified (Figure 3b).

Budyko Framework
Evapotranspiration in a catchment can be calculated from the water balance equation as follows: where the annual evapotranspiration, E; the annual precipitation, P; and the annual discharge, Q, are calculated as the average flux on the basin, and ΔS is the increment in storage during an accounting period.Over a long period (usually more than 10 years), the change in storage (ΔS) is negligible for the mean annual water balance, such that E = P − Q can be applied.Budyko [42] demonstrated that the mean annual evapotranspiration ratio (E/P) is mainly controlled by the mean aridity index (E0/P, where E0 is the annual potential evapotranspiration).This hypothesis is the underlying principle of the Budyko framework when building a straightforward estimation formula for E/P versus E0/P.Non-parametric and one-parameter formulas have been proposed in the literature; however, the one-parameter formulas are more practical in comparison to the non-parametric formulas, as the catchment-specific parameter can represent the variable

Budyko Framework
Evapotranspiration in a catchment can be calculated from the water balance equation as follows: where the annual evapotranspiration, E; the annual precipitation, P; and the annual discharge, Q, are calculated as the average flux on the basin, and ∆S is the increment in storage during an accounting period.Over a long period (usually more than 10 years), the change in storage (∆S) is negligible for the mean annual water balance, such that E = P − Q can be applied.Budyko [42] demonstrated that the mean annual evapotranspiration ratio (E/P) is mainly controlled by the mean aridity index (E 0 /P, where E 0 is the annual potential evapotranspiration).This hypothesis is the underlying principle of the Budyko framework when building a straightforward estimation formula for E/P versus E 0 /P.Non-parametric and one-parameter formulas have been proposed in the literature; however, the one-parameter formulas are more practical in comparison to the non-parametric formulas, as the catchment-specific parameter can represent the variable catchment properties [43].Accordingly, the long-term mean evapotranspiration ratio can be calculated as follows: where w is the catchment-specific parameter and F() denotes a one-parameter formula.The long-term mean discharge can then be calculated as where f () is the function derived from F().The parameter w can be inversely determined with observation data of Q, E 0 , and P for a long period using the mean annual water balance or for year-by-year patterns using a moving average method [39].

The Traditional Separation Approach
In the traditional separation approach, the discharge change between two periods, indicated as 3a, is primarily separated into two parts on the basis of the Budyko framework: As shown in Figure 3a, Q 0 and Q 1 are basin-scale discharges with respect to different periods that are represented by points located on different Budyko curves.∆Q 0 is the discharge change from Q 0 to Q a caused by the change in E 0 /P along the same Budyko curve.∆Q w is the discharge change from Q a to Q 1 caused by the change in specific parameter (∆w) when the E 0 /P value remains the same.
There are three common used methods to calculate ∆Q 0 and ∆Q w : the climate elasticity method [44][45][46], the sensitivity method [47], and the decomposition method [48].As compared by Wang and Hejazi [48], these separation methods obtained similar results when they were applied for the same catchments.However, the decomposition method eliminates the discretization errors caused by the calculation of partial derivatives in the other two methods [39,49].By the decomposition method [48], the different parts of the change in discharge can be estimated as In previous studies, using the traditional separation method, the discharge change caused by the specific parameter change (∆Q w ) was directly attributed to land-use changes or human activity [20,21].However, for explicitly explaining the discharge change in the SRYR, further attribution analysis was imperative.

Further Attribution Analysis Based on the Modified Separation Approach
As indicated in Figure 3a, the discharge change caused by ∆w is defined as ∆Q w .The further attribution analysis of ∆Q w or ∆w in previous studies generally contained two steps: firstly, obtaining time-varying simples (∆Q w or w) by the moving average method [39] or the sequential average method [27]; secondly, determining potential explanatory variables from stepwise regression analysis [39] or a curve-fitting method [27].In this study, the 11 year moving average method was Water 2018, 10, 877 7 of 18 used to obtain the time-varying specific parameter w, and the stepwise regression method was used to determine the explanatory variables.
Previous studies have reported that the catchment-specific parameter depends on several catchment properties (the relative infiltration capacity, relative soil water storage, etc.) that can be influenced by climate change, such as by changes in the precipitation intensity (I), evapotranspiration (E 0 ), and temperature (T) [36,37,39].To identify the climatic impacts on the catchment properties, we applied an indirect method that was similar to the approach used by Jiang [39].The w value was assumed to be linearly correlated with potential explanatory variables by the following: where β 0 , β 1 , β 2 , and β 3 are coefficients; T is the annual MAT; I represents the characteristics of precipitation intensity; and E 0 is the annual potential evapotranspiration.Each of the variables was normalized by the absolute mean annual value, indicated by subscript m (I m , E m , and T m ).
The β i (i = 0, 1, 2, . . . ) values were determined through stepwise regression analysis.In each step of the stepwise regression analysis, a variable was considered for addition to, or subtraction from, the set of explanatory variables on the basis of a sequence of F-tests, and the variance inflation factor (VIF) was used to diagnose the collinearity.The selected model had the highest R 2 value and contained variables with a VIF of <10.A high VIF value indicates the presence of multicollinearity; a VIF value of 10 is often adopted to diagnose multicollinearity.In practice, a VIF of <10 indicates ignorable multicollinearity among such variables [50][51][52][53].
According to Equation ( 7), the climatic change in the catchment-specific parameter, ∆w c , can be estimated as where to Equation (7), a modified formula is proposed in this study to check the impact of change in the frozen ground: where D and β 4 are the normalized MFDs used to indicate the status of the frozen ground and the relevant coefficient, respectively.In most of the existing studies based on decomposition methods, ∆Q 0 was considered as the impact of climate change on discharge, whereas ∆Q w was attributed to the change in land use or human activities.However, in this study, we propose an additional climatic discharge change (∆Q wc ) that links with ∆w c ; it is also induced by variations in climate conditions and relevant properties (MFD).This additional climatic discharge change is estimated in a similar way to Equation (6) (Figure 3b): Accordingly, we define the total climatic change in discharge as the following summation: Because of the very low population density (2003 census data [15]) and a slight land-use change [22] in the SRYR, human activities were not directly accounted for in this study, but the impact of land-use change on discharge could then be approximately evaluated by ∆Q l = ∆Q − ∆Q c (Figure 3b).
In this study, the normal climatic impacts on the discharge change in the SRYR and its three sub-basins were assessed using Equation (8), and the impact of frozen-ground degradation on the three sub-basins with different frozen-ground distributions were particularly analysed using Equation (9).The comparison between the results of Equations ( 8) and ( 9) indicated the impacts of the frozen-ground degradation, as presented in Sections 4.3 and 4.4.

Selection of the Budyko Formula
Fu's equation [54,55] was selected in this study.Substituting Fu's equation into Equation (3), we obtain To inversely obtain the w value from the observation data, we simply used E = P -Q but included an uncertainty of the impact of the ignored change in water storage.To minimise the effect of the change in storage, the 11 year moving average method was applied [39].

Climate and Discharge Changes
The long-term discharge changes of the SRYR and its three sub-basins are shown in Figure 4; the variation patterns generally can be divided into three periods: the pre-change period from 1961 to 1990 and two post-change periods-a low-flow period from 1991 to 2002 and a recent period from 2003 to 2013.In this section, the climate and discharge changes of the 1991-2002 and 2003-2013 periods are investigated in comparison with the 1961-1990 period, which was the same as the climate pre-change period defined by IPCC [1].Climate change showed temporal and spatial differences in the SRYR (Table 1).For the period 1991-2002, P decreased (3.0-6.6%) and E 0 slightly increased (1.24-1.50%) in the whole study area.The largest change in P (−6.56%) and E 0 (1.54%) occurred in MQ.In the 21st century (2003-2013), P increased (5.52-11.05%)across the whole study area, with a sharp increase in the JM (11.05%) and TNH (6.25%) sub-basins.Potential evapotranspiration clearly increased in the 21st century, and ∆E 0 varied from 2.89% to 4.41%, with the largest increases in E 0 (4.40%) and P (11.05%) occurring in MQ and JM, respectively.
Water 2018, 10, x FOR PEER REVIEW 8 of 18 To inversely obtain the w value from the observation data, we simply used E = P -Q but included an uncertainty of the impact of the ignored change in water storage.To minimise the effect of the change in storage, the 11 year moving average method was applied [39].

Climate and Discharge Changes
The long-term discharge changes of the SRYR and its three sub-basins are shown in Figure 4; the variation patterns generally can be divided into three periods: the pre-change period from 1961 to 1990 and two post-change periods-a low-flow period from 1991 to 2002 and a recent period from 2003 to 2013.In this section, the climate and discharge changes of the 1991-2002 and 2003-2013 periods are investigated in comparison with the 1961-1990 period, which was the same as the climate pre-change period defined by IPCC [1].Climate change showed temporal and spatial differences in the SRYR (Table 1).For the period 1991-2002, P decreased (3.0-6.6%) and E0 slightly increased (1.24-1.50%) in the whole study area.The largest change in P (−6.56%) and E0 (1.54%) occurred in MQ.In the 21st century (2003-2013), P increased (5.52-11.05%)across the whole study area, with a sharp increase in the JM (11.05%) and TNH (6.25%) sub-basins.Potential evapotranspiration clearly increased in the 21st century, and ΔE0 varied from 2.89% to 4.41%, with the largest increases in E0 (4.40%) and P (11.05%) occurring in MQ and JM, respectively.As shown in Figure 4 and Table 1, discharge in the SRYR and its three sub-basins decreased by more than 20% in 1991-2002 but showed a recovery trend in 2003-2013, particularly in JM.The discharge changes in the two periods 1991-2002 and 2003-2013 were further partitioned using the decomposition method, and the results are listed in Table 1.It is indicated that both ΔQ0 and ΔQw were significant.In general, ΔQw contributed more than half of the total absolute change, indicating that the change in the catchment parameter plays a more important role in the hydrological behaviour of the SRYR, which is similar to the findings of previous research for the other catchments [20,21].A relatively different result was found for the 21st century for the JM basin, where ΔQ was positive but  As shown in Figure 4 and Table 1, discharge in the SRYR and its three sub-basins decreased by more than 20% in 1991-2002 but showed a recovery trend in 2003-2013, particularly in JM.The discharge changes in the two periods 1991-2002 and 2003-2013 were further partitioned using the decomposition method, and the results are listed in Table 1.It is indicated that both ∆Q 0 and ∆Q w were significant.In general, ∆Q w contributed more than half of the total absolute change, indicating that the change in the catchment parameter plays a more important role in the hydrological behaviour of the SRYR, which is similar to the findings of previous research for the other catchments [20,21].A relatively different result was found for the 21st century for the JM basin, where ∆Q was positive but ∆Q w was negative, and the normal direct climatic change, ∆Q 0 , became a dominant component in ∆Q.
If ∆Q w is completely explained as the human-induced change in the discharge, we would obtain an unreasonable conclusion on the basis of the results in Table 1: human activities were able to cause significant negative impact on the hydrological processes in areas with a very low population density.This indicates that something is missing from the analysis; ∆w should not be completely referred to as the change in land use.Therefore, a further analysis is necessary to extract the impact of the climate change on ∆Q w .

Time-Varying w Value
The shift paths of the basin-scale water balance were plotted in the Budyko space (Figure 5a).Each point in the Budyko space could be located to a given Budyko curve by changing the w value for Fu's equation.Typical Budyko curves for w = 2.4, 2.0, and 1.7 are also shown in Figure 5a for comparison.The gray area in Figure 5a indicates a zone confined by the Budyko curves with the maximum and minimum w values.Using the 11 year moving average values of P, E 0 , and Q, the time-varying w value was obtained using Fu's equation.As shown in Figure 5b, the time-varying values of w in the three sub-basins and the SRYR significantly increased during 1980-2000.In the 21st century, this increasing trend continued for MQ but stopped and transitioned to a declining trend for JM, TNH, and the SRYR.MQ had the lowest dryness index (E 0 /P), but its w value ranged between the w values of JM and TNH.The data points generally showed a positive relationship between the evapotranspiration ratio (E/P) and E 0 /P, with slopes that were steeper than the slope of traditional Budyko curves with a constant parameter, indicating that ∆Q w plays a more important role than ∆Q 0 in the discharge change (∆Q) of the SRYR and the three sub-basins (∆Q w and ∆Q 0 refer to Figure 3a).century, this increasing trend continued for MQ but stopped and transitioned to a declining trend for JM, TNH, and the SRYR.MQ had the lowest dryness index (E0/P), but its w value ranged between the w values of JM and TNH.The data points generally showed a positive relationship between the evapotranspiration ratio (E/P) and E0/P, with slopes that were steeper than the slope of traditional Budyko curves with a constant parameter, indicating that ΔQw plays a more important role than ΔQ0 in the discharge change (ΔQ) of the SRYR and the three sub-basins (ΔQw and ΔQ0 refer to Figure 3a).

Results of Using Equation (8)
In this study, we used the new method (Section 3.3) to further investigate ∆Q w with the regression formula (Equation ( 8)).The relationship between the moving average w values and the climatic forces was analysed with Equation (7).Four optimal regression models (with a significance level of 0.001) were chosen through covariate analysis (Table 2).For the sub-basin JM, the air temperature (T) had a more essential role than the potential evapotranspiration (E 0 ) in controlling the w value, which contrasts with MQ.In TNH and the SRYR, both T and E 0 were significant factors that correlated to the w value.The precipitation intensity (I) was a significant factor both in the SRYR and the three sub-basins.The selected variables in the respective model were not the same, but the explanatory variables contained in each model could be generally divided into two groups: water-supply-and energy-supply-related.Water-supply-related explanatory variables included I, and the energy-supply-related explanatory variables included T or E 0 .According to the Budyko hypothesis, the mean annual water balance in a basin is a synthesis of water supply from precipitation and energy-supply-controlled loss of water by evapotranspiration.Both water-supply-and energy-supply-related variables were included in each formula (Table 2), indicating the reasonability of the obtained results.To assess the impact of the climate change on discharge at the basin scale, ∆Q wc and ∆Q c were estimated with Equations ( 10) and (11), respectively.The change was related to the moving average data of the first window of 1961-1971.Results are shown in Figure 6, where ∆Q, ∆Q c , ∆Q 0 , and ∆Q wc are presented as the percentages of the mean annual discharge during 1961-1971 (substituted by 1966 in Figure 6).As indicated, ∆Q c was close to ∆Q for the three sub-basins and the SRYR, demonstrating that climate change is the primary cause of the discharge change in the SRYR.The absolute values of ∆Q wc were generally higher than those of ∆Q 0 , indicating that the change in the catchment-specific parameter plays a more important role than the direct impacts of the changes in the precipitation and potential evapotranspiration.In JM, TNH, and the SRYR, as shown in Figure 6a,c, the absolute value of ∆Q wc was more than twice the absolute value of ∆Q 0 between 1988 and 2006.The contributions of the catchment-specific parameter changes obtained by the decomposition methods(∆Q w ) and shown in Table 1 are similar to those of the modified method(∆Q−∆Q 0 ≈ ∆Q wc ), which indicates that the modified method can obtain time-varying results without expense of accuracy.

Results of Using Equation (9)
To assess different frozen-ground distribution impacts in the three sub-basins, normalised MFD change series ( ) of the three sub-basins were added as candidate explanatory variables in Equation (7), and further stepwise regression analyses were carried out for JM, MQ, and TNH.The results showed that improved the regression model of the two sub-basins JM and TNH, where permafrost is widely distributed, as follows: These regression formulas could increase the correlation coefficient to R 2 = 0.92 and R 2 = 0.93 for JM and TNH, respectively.However, was not selected by the stepwise regression analysis in MQ, which indicates that the impact of frozen-ground degradation on MQ was negligible.For the SRYR, w was more significantly correlated with the MFD of TNH rather than with the mean MFD of the whole area; this may be linked with the fact that TNH is close to the outlet of the whole area.By considering the MFD of TNH in the SRYR stepwise regression analysis, the regression coefficient slightly increased from R 2 = 0.94 to R 2 = 0.95, indicating that the impact of frozen-ground degradation also existed for the discharge change in the whole SRYR.
The coefficients of were positive and close to 1.0, indicating that the MFD was a positive factor for the catchment parameter.Equations ( 13) and ( 14) also showed that the w value was positively correlated with the air temperature, T, and negatively correlated with the precipitation intensity, I, in both sub-basins and was positively correlated with E0 in TNH.These results for the relationship between the catchment parameter and the normal climatic forces were consistent with those shown in Table 2, but the values of the coefficients were different.
The change in discharge derived from land-use changes was then estimated by ΔQl = ΔQ − ΔQc from the new results of ΔQc.As shown in Figure 7, the ΔQl values were generally less than 5% of the total change and appeared to show a random oscillation, which indicates that the land-use-change

Results of Using Equation
To assess different frozen-ground distribution impacts in the three sub-basins, normalised MFD change series (D) of the three sub-basins were added as candidate explanatory variables in Equation ( 7), and further stepwise regression analyses were carried out for JM, MQ, and TNH.The results showed that D improved the regression model of the two sub-basins JM and TNH, where permafrost is widely distributed, as follows: These regression formulas could increase the correlation coefficient to R 2 = 0.92 and R 2 = 0.93 for JM and TNH, respectively.However, D was not selected by the stepwise regression analysis in MQ, which indicates that the impact of frozen-ground degradation on MQ was negligible.For the SRYR, w was more significantly correlated with the MFD of TNH rather than with the mean MFD of the whole area; this may be linked with the fact that TNH is close to the outlet of the whole area.By considering the MFD of TNH in the SRYR stepwise regression analysis, the regression coefficient slightly increased from R 2 = 0.94 to R 2 = 0.95, indicating that the impact of frozen-ground degradation also existed for the discharge change in the whole SRYR.
The coefficients of D were positive and close to 1.0, indicating that the MFD was a positive factor for the catchment parameter.Equations ( 13) and ( 14) also showed that the w value was positively correlated with the air temperature, T, and negatively correlated with the precipitation intensity, I, in both sub-basins and was positively correlated with E 0 in TNH.These results for the relationship Water 2018, 10, 877 12 of 18 between the catchment parameter and the normal climatic forces were consistent with those shown in Table 2, but the values of the coefficients were different.
The change in discharge derived from land-use changes was then estimated by ∆Q l = ∆Q − ∆Q c from the new results of ∆Q c .As shown in Figure 7, the ∆Q l values were generally less than 5% of the total change and appeared to show a random oscillation, which indicates that the land-use-change impacts on discharge were negligible.

The Role of Temperature in Discharge Change
In this study, a modified method based on the Budyko framework was used to investigate hydrological responses to climate change in the SRYR, by Equations ( 9) and (11).It was found that temperature (T), precipitation intensity (I), and evapotranspiration (E0) were controls on the discharge change in the three sub-basins by altering the catchment-specific parameter.In the past five decades, the temperature has increased in the Yellow River Basin [56][57][58][59]; the effect of temperature changes on discharge is also highlighted in the largest tributary of the Yellow River: the Weihe River Basin [39].This indicates that air temperature is an important factor in altering catchment discharge by disturbing the true evaporation [39,56].
Wang [27] used the freezing index (DDF) as an indicator of the frozen-ground change to analyse the discharge change in the SRYR and concluded that a 32.6% discharge decrease in the 1990s was caused by frozen-ground degradation.The DDF is defined as the absolute value of accumulated negative air temperature in a hydrological year (from September 1 to August 31 in the SRYR).In that study, the DDF was considered as the only explanatory variable for change in ΔQw.However, one should be aware of the fact that the DDF is significantly correlated with annual MAT; attributing DDF impacts on the discharge change to frozen-ground degradation may cause a misleading assessment of the pure hydrological impact of frozen-ground degradation in the SRYR.
The temperature effect on the discharge change should be carefully considered when the freezing index (DDF) is used in analysis.Accordingly, we further analysed the relationship between the air temperature and ΔQw by nonlinear regression as shown in Figure 8.In the three sub-basins and the SRYR, the air temperature was significantly correlated with ΔQw, which was similar to the relationship between the DDF and ΔQw obtained by Wang [27].However, ΔQw was not monotonically decreasing with increasing temperature, particularly in JM, TNH, and the SRYR.This indicates that using only temperature change could not sufficiently explain the change in ΔQw.In this study, we selected three climate factors (T, E0, and I) as candidate explanatory variables in the stepwise regression analysis and obtained four optimal models (Table 2).These optimal models explained the ΔQw change in the SRYR and the three sub-basins well (Figure 6).Because T and E0 were positive controls for the w value, they acted as negative controls for ΔQw.However, precipitation intensity (I) was a negative factor for w change and thus a positive factor for ΔQw change; hence the nonmonotonic relationship between ΔQw and T was caused by increased precipitation.

The Role of Temperature in Discharge Change
In this study, a modified method based on the Budyko framework was used to investigate hydrological responses to climate change in the SRYR, by Equations ( 9) and (11).It was found that temperature (T), precipitation intensity (I), and evapotranspiration (E 0 ) were controls on the discharge change in the three sub-basins by altering the catchment-specific parameter.In the past five decades, the temperature has significantly increased in the Yellow River Basin [56][57][58][59]; the effect of temperature changes on discharge is also highlighted in the largest tributary of the Yellow River: the Weihe River Basin [39].This indicates that air temperature is an important factor in altering catchment discharge by disturbing the true evaporation [39,56].
Wang [27] used the freezing index (DDF) as an indicator of the frozen-ground change to analyse the discharge change in the SRYR and concluded that a 32.6% discharge decrease in the 1990s was caused by frozen-ground degradation.The DDF is defined as the absolute value of accumulated negative air temperature in a hydrological year (from September 1 to August 31 in the SRYR).In that study, the DDF was considered as the only explanatory variable for change in ∆Q w .However, one should be aware of the fact that the DDF is significantly correlated with annual MAT; attributing DDF impacts on the discharge change to frozen-ground degradation may cause a misleading assessment of the pure hydrological impact of frozen-ground degradation in the SRYR.
The temperature effect on the discharge change should be carefully considered when the freezing index (DDF) is used in analysis.Accordingly, we further analysed the relationship between the air temperature and ∆Q w by nonlinear regression as shown in Figure 8.In the three sub-basins and the SRYR, the air temperature was significantly correlated with ∆Q w , which was similar to the relationship between the DDF and ∆Q w obtained by Wang [27].However, ∆Q w was not monotonically decreasing with increasing temperature, particularly in JM, TNH, and the SRYR.This indicates that using only temperature change could not sufficiently explain the change in ∆Q w .In this study, we selected three climate factors (T, E 0 , and I) as candidate explanatory variables in the stepwise regression analysis and obtained four optimal models (Table 2).These optimal models explained the ∆Q w change in the SRYR and the three sub-basins well (Figure 6).Because T and E 0 were positive controls for the w value, they acted as negative controls for ∆Q w .However, precipitation intensity (I) was a negative factor for w change and thus a positive factor for ∆Q w change; hence the non-monotonic relationship between ∆Q w and T was caused by increased precipitation.
Water 2018, 10, x; doi: FOR PEER REVIEW www.mdpi.com/journal/waterwas significantly correlated with ΔQw, which was similar to the relationship between the DDF and ΔQw obtained by Wang [27].However, ΔQw was not monotonically decreasing with increasing temperature, particularly in JM, TNH, and the SRYR.This indicates that using only temperature change could not sufficiently explain the change in ΔQw.In this study, we selected three climate factors (T, E0, and I) as candidate explanatory variables in the stepwise regression analysis and obtained four optimal models (Table 2).These optimal models explained the ΔQw change in the SRYR and the three sub-basins well (Figure 6).Because T and E0 were positive controls for the w value, they acted as negative controls for ΔQw.However, precipitation intensity (I) was a negative factor for w

The Role of Frozen-Ground Degradation in Discharge Change
To analyse the impact of frozen-ground degradation on the discharge change, the observed MFD was used to indicate the frozen-ground change in this study.As shown in Figure 2, the MFD has generally decreased since the 1980s, whereas the air temperature has increased.To analyse this correlation, we compared the MFD data for 10 stations with the annual MAT and freezing index (DDF).The locations of these stations

The Role of Frozen-Ground Degradation in Discharge Change
To analyse the impact of frozen-ground degradation on the discharge change, the observed MFD was used to indicate the frozen-ground change in this study.As shown in Figure 2, the MFD has generally decreased since the 1980s, whereas the air temperature has increased.To analyse this correlation, we compared the MFD data for 10 stations with the annual MAT and freezing index (DDF).The locations of these stations are shown in Figure 1, and the correlation analysis results are shown in Table 3.The MFD was positively correlated with the DDF and negatively correlated with the MAT.In general, the significance of the correlation between the MFD and DDF was higher than the correlation between the MFD and MAT, except for the M8 station.This indicates that the decrease in MFD was mainly controlled by the increase in the air temperature for the cold seasons.However, the R 2 correlation coefficients were not high, indicating that frozen ground in the region may be influenced by other factors, such as the soil temperature and moisture, snowpack, and canopy [60][61][62][63].Long-term observed MFD data series may be more accurate in indicating frozen-ground degradation than the freezing index.If the negative correlation between the MFD and MAT is significant, it would influence the meaning of the regression formulas for the catchment parameter, w.According to Equations ( 13) and ( 15), w can be determined by where w f are the components related to the other factors and β 3 and β 4 are the coefficients for the MAT and MFD (both are positive), respectively.Further, we can describe the MFD-MAT correlation by where ζ is the absolute slope of the correlation line and η is the residual resulting from other controls.Substituting Equation ( 16) into Equation ( 15) yields Because β 3 , β 4 , and ζ are positive coefficients, then β 3 − β 4 ζ must be less than β 3 .Accordingly, the MFD-MAT correlation may reduce the positive sensitivity (when ζ is less than β 3 /β 4 ) or cause negative sensitivity (when ζ is greater than β 3 /β 4 ) of the catchment parameter to the change in air temperature.This may explain why the coefficient of T in Table 2 was less than that of Equations ( 13) and (14).If β 3 ≈ β 4 ζ, then the role of the MAT cannot be explicitly included in Equation ( 15).This may be the reason why T is not included in the formulas of the sub-basin MQ in Table 2, although the effects of air temperature should also exist in the basin.It is speculated that in MQ, the w value is also positively correlated with the MFD but is offset by the effect of air temperature.
As indicated in Equations ( 13) and (14), this positive relationship between w and MFD implies that the degradation of frozen ground (decreasing MFD) will cause negative changes in the w and E/P values and thus increase the discharge in JM and TNH, which are widely underlain by permafrost.In general, increasing air temperature increases evapotranspiration and decreases the discharge in catchments [39,56].However, the degradation of permafrost due to increasing air temperature would increase the discharge, which partially counteracts the negative impact of the increasing air temperature.This effect leads to a more complex response of streamflow in the cold basins with permafrost than that in warm basins without permafrost.
Why does the decreasing MFD act as a positive factor for the discharge of basins in the SRYR?This is still an unsolved problem that was not resolved using only the Budyko framework.However, these results are consistent with the degradation of permafrost due to increasing air temperature being able to increase the discharge of groundwater [28][29][30][31][32][33]; the decrease in the MFD is a major factor for the increase in baseflow in the Qilian Mountain [64], and melt ice within permafrost and increasing hydrological connectivity following permafrost thaw will increase runoff discharge [34,35].
The changes in the catchment-specific parameters (w) in the SRYR and the three sub-basins were mainly caused by changes in the water-supply-related variable, I, and energy-supply-related variables, T and E 0 .In addition, the catchment-specific parameters of JM and TNH were correlated with the MFD of frozen ground, which indicates that the impacts of the frozen-ground degradation on hydrological responses in the permafrost dominated the sub-basins.(4) Increasing air temperature is generally a negative force for discharge.However, in the SRYR, it also causes permafrost degradation that can act as a positive factor for discharge.Such conflicting effects enhance the uncertainty in assessments of the hydrological response to climate change in the sub-basins with widely underlain permafrost.

Figure 1 .
Figure 1.Relative location of the three sub-basins (Jimai (JM), Jimai-Maqu (MQ), and Maqu-Tangnaihai (TNH)), distribution of frozen ground, and meteorological and hydrological stations (the red circles indicate stations that have the longest frozen-ground record in their respective sub-basin).

Water 2018 ,
10, x FOR PEER REVIEW 5 of 18

Figure 3 .
Figure 3. Partition of discharge change in the Budyko framework.(a) The traditional separation approach: discharge change (ΔQ) is primarily separated into two components, ΔQ0 and ΔQw.(b) The modified separation approach: ΔQw is further separated into ΔQwc and ΔQl in considering the impacts of climate change and frozen-ground degradation on the catchment-specific parameter.

Figure 3 .
Figure 3. Partition of discharge change in the Budyko framework.(a) The traditional separation approach: discharge change (∆Q) is primarily separated into two components, ∆Q 0 and ∆Q w .(b) The modified separation approach: ∆Q w is further separated into ∆Q wc and ∆Q l in considering the impacts of climate change and frozen-ground degradation on the catchment-specific parameter.
(= I/I m ), E 0 (= E 0 /E m ), and T(= T/T m ) are the normalized values of I, E 0 , and T, respectively.Using the MFD to indicate frozen ground change and adding normalized MFD change series (D)

Figure 4 .
Figure 4. Time series of annual discharge, Q, in the three sub-basins and the source region of the Yellow River (SRYR).The dashed, red, and blue lines indicate the mean values in different periods, 11 year moving average value, and annual mean value of discharge, respectively.

Figure 4 .
Figure 4. Time series of annual discharge, Q, in the three sub-basins and the source region of the Yellow River (SRYR).The dashed, red, and blue lines indicate the mean values in different periods, 11 year moving average value, and annual mean value of discharge, respectively.

Figure 5 .
Figure 5. Plots of the moving average water balance in the Budyko space: (a) time-varying values of the catchment parameter (w) in Fu's formula, and (b) for the three sub-basins.In (a), 2.4, 2.0, and 1.7 refer to the w value for the three types of curves.

Figure 5 .
Figure 5. Plots of the moving average water balance in the Budyko space: (a) time-varying values of the catchment parameter (w) in Fu's formula, and (b) for the three sub-basins.In (a), 2.4, 2.0, and 1.7 refer to the w value for the three types of curves.

Water 2018 , 18 Figure 6 .
Figure 6.Decomposition results of the climatic change in discharge (ΔQ0, ΔQwc, and ΔQc) in comparison with the total change in discharge (ΔQ) for the three sub-basins and the source region of the Yellow River (SRYR).(a-d) The results of using Equation (8) for Jimai (JM), Maqu (MQ), Tangnaihai (TNH), and the SRYR.

Figure 6 .
Figure 6.Decomposition results of the climatic change discharge (∆Q 0 , ∆Q wc , and ∆Q c ) in comparison with the total change in discharge (∆Q) for the three sub-basins and the source region of the Yellow River (SRYR).(a-d) The results of using Equation (8) for Jimai (JM), Maqu (MQ), Tangnaihai (TNH), and the SRYR.

Water 2018 , 18 Figure 7 .
Figure 7. Time series of ΔQl in the three sub-basins and the source region of the Yellow River (SRYR).

Figure 7 .
Figure 7. Time series of ∆Q l in the three sub-basins and the source region of the Yellow River (SRYR).

Figure 8 .
Figure 8. Relationship between temperature and ΔQw in the three sub-basins and the source region of the Yellow River (SRYR).

Figure 8 .
Figure 8. Relationship between temperature and ∆Q w in the three sub-basins and the source region of the Yellow River (SRYR).

Table 1 .
Results of decomposition method.

Table 2 .
Results of w covariate analysis for the period between 1961 and 2013.

Table 3 .
Correlations between maximum frozen depth (MFD) and DDF at weather stations, compared to correlations between MFD and mean air temperature (MAT).