Attribution Analysis of Seasonal Runoff in the Source Region of the Yellow River Using Seasonal Budyko Hypothesis

: Previous studies mainly focused on quantifying the contribution rate of different factors on annual runoff variation in the source region of the Yellow River (SRYR), while there are few studies on the seasonal runoff variation. In this study, the monthly water storage and monthly actual evaporation of SRYR were calculated by the monthly ABCD model, and then a seasonal Budyko frame was constructed. Finally, the contribution rate of climatic and anthropic factors on the seasonal runoff variation in Tangnaihai hydrological station were quantitatively calculated. It turned out that: (1) The changing point of runoff data at Tangnaihai hydrological station is 1989. (2) The ABCD monthly hydrological model could well simulate the monthly runoff variation of Tangnaihai hydrological station. (3) Anthropic factors play a major role in runoff change in spring, summer, and winter, while climatic factors play a major role in runoff change in autumn.


Introduction
The source region of the Yellow River (SRYR), with the area of about 121,900 km 2 , accounting for 16.2% of the basin, is the main runoff producing area and water conservation area in the Yellow River Basin. Its annual average runoff is about 248.72 × 10 8 m 3 , accounting for more than 40% of the average annual runoff of the basin [1]. The runoff in SRYR has gradually decreased, which has endangered the ecological security of the basin and harmed its high-quality economic development [1]. In addition, the reduction of water resources also endangered the grain production in the basin.
Climatic and anthropic factors are the two major driving forces of the hydrological cycle change [2][3][4]. Climate change mainly affects runoff through variations in precipitation and temperature [5,6]. The impact of human activities on runoff can be divided into the direct impact of anthropic factors and the indirect impact of underlying surface changes on runoff. The direct impacts include man-made water intake, inter-basin water transfer, and water conservancy project construction [7][8][9][10]. The indirect effects include urbanization and afforestation projects on the underlying surface conditions of the basin [11][12][13][14].
Many researches distinguished the impact of climatic and anthropic factors on annual runoff in SRYR. Jiang  the main factor, with a contribution rate of more than 90% [15]. Wang et al. distinguished the contribution rate of climate change, frozen soil degradation, and land type change on annual runoff change in SRYR from 1965 to 2015. Compared with 1965-1989, the contribution rates of climate change, frozen soil degradation, and land type change on annual runoff change in 1990-2003 were 55.0%, 32.6%, and 12.4%, respectively. Compared with 1990-2003, climate change, frozen soil degradation, and land use change contributed 126.0%, −38.5%, and 12.5% on annual runoff change in 2004-2015 [16]. Wang et al. analyzed the contribution rate of climate and underlying surface changes on annual runoff change.
The results show that an increase of 1% in precipitation will lead to an increase of 1.97% in runoff, and while an increase of 1% in evaporation or underlying surface parameters will would lead to a reduction of 0.97% or 1.18% in runoff. The contribution rate of precipitation change is 24.0%, and that of evaporation is 18.0% [17]. Zhang et al. used the elastic coefficient method to distinguish the impact of climatic and anthropic factors on the annual runoff variation in SRYR from 1965 to 2013. The results showed that the contribution rate of climate change was about 39%-45% and that of anthropic factors was about 55%-61% from 1990 to 1999. Anthropic factors were the leading reason from 2000 to 2013 [18]. Yan et al. (2021) used Budyko-type equations for separating the impacts of climate and vegetation change on annual runoff variation in SRYR from 1961 to 2015 [19]. However, previous studies aimed to quantify the impact of climatic and anthropic factors on annual runoff change, while few attempts have been made to quantitatively assess the impact of climatic and anthropic factors on river runoff changes at seasonal or monthly scales.
The commonly used methods for quantifying the impact of climatic and anthropic factors on runoff change include the Budyko hypothesis and hydrologic model simulation. The Budyko hypothesis takes into account the water and energy balance in the hydrological process and is easy to calculate. The time scale in the Budyko hypothesis is defined as the average value of multi-year meteorological and hydrological data [20]. Therefore, the Budyko hypothesis, based on an assumption that the change of catchment storage water for multi-year water balance is considered as 0, is widely applied to quantitatively calculate the impact of climate change and anthropic factors on long-term annual runoff change, while ignoring their impact on intra-annual runoff distribution [21][22][23][24][25][26]. For multiyear water balance, the change of water storage is usually negligible compared with the average annual precipitation depth. However, monthly or seasonal Budyko frameworks must consider the impact of water storage change on monthly or seasonal runoff change. Therefore, the change of water storage should be regarded as an important part of Budyko's steady-state hypothesis on monthly or seasonal scales. For the hydrologic model simulation method, the impact of a factor (climatic or anthropic) on runoff can be calculated by changing one parameter and fixing other parameters at the same time, and the contribution rate of climatic and anthropic factors to runoff variation can be separated by comparing the differences between observed and simulated runoff during the post-change period and observed runoff in the baseline period [27][28][29]. Common hydrological models include the Soil Water Assessment Tool [30,31], Distributed time-variant gain model [32], Geomorphology-Based Hydrological Model [33], Precipitation Runoff Modeling System [34], ABCD model [35], and Variable Infiltration Capacity model [36] et al. Hydrological model simulation can accurately express the mechanism of runoff evolution and simulate monthly water storage change. However, uncertainties are existed in both model parameters and input data [37][38][39][40][41].
Compared with other hydrological models, the ABCD model has the following advantages: (1) less parameters, high simulation accuracy and easy to optimize; (2) observed soil water content and groundwater data are unnecessary; (3) high adaptability due to unseparated direct runoff and groundwater outflow when calculating runoff; (4) the model can simulate the changing process of soil water content without estimating the parameters of soil aquifer; (5) the ABCD model is calculated by genetic algorithm, which is simple and can be easily implemented in MATLAB [42][43][44][45]. In short, for grasping the intra-annual runoff distribution mechanism of SRYR, this study attempts to quantitatively analyze the impact of climatic and anthropic factors on the seasonal runoff variation in SRYR by following four steps: (1) distinguishing the mutation year of runoff time series data by Mann-Kendall and cumulative anomaly methods; (2) calculating the monthly water storage and monthly actual evaporation of SRYR by the monthly ABCD model; (3) constructing a seasonal Budyko framework of SRYR; and (4) quantifying the contribution rate of climatic and anthropic factors on the runoff variation of Tangnaihai hydrological station in spring, summer, autumn, and winter by the decomposition method.

Study Area
SRYR is located between 32 • N~36 • N and 95 • E~104 • E ( Figure 1). The topography of the basin is uneven. The altitude ranges from 2600 to 6200 m, with an area of about 121,900 km 2 , accounting for 16.2% of the basin. It is an important runoff producing area in the Yellow River Basin. Its annual average runoff is about 248.72 × 10 8 m 3 , accounting for more than 40% of the average annual runoff of the basin.
Land 2021, 10, x FOR PEER REVIEW 3 of 5 model can simulate the changing process of soil water content without estimating the parameters of soil aquifer; (5) the ABCD model is calculated by genetic algorithm, which is simple and can be easily implemented in MATLAB [42][43][44][45]. In short, for grasping the intra-annual runoff distribution mechanism of SRYR, this study attempts to quantitatively analyze the impact of climatic and anthropic factors on the seasonal runoff variation in SRYR by following four steps: (1) distinguishing the mutation year of runoff time series data by Mann-Kendall and cumulative anomaly methods; (2) calculating the monthly water storage and monthly actual evaporation of SRYR by the monthly ABCD model; (3) constructing a seasonal Budyko framework of SRYR; and (4) quantifying the contribution rate of climatic and anthropic factors on the runoff variation of Tangnaihai hydrological station in spring, summer, autumn, and winter by the decomposition method.

Study Area
SRYR is located between 32° N~36° N and 95° E~104° E ( Figure 1). The topography of the basin is uneven. The altitude ranges from 2600 to 6200 m, with an area of about 121,900 km 2 , accounting for 16.2% of the basin. It is an important runoff producing area in the Yellow River Basin. Its annual average runoff is about 248.72 × 10 8 m 3 , accounting for more than 40% of the average annual runoff of the basin.

Data Sources
The observed monthly runoff data in the period 1967-2016 in Tangnaihai station were obtained from the Yellow River Water Conservancy Commission ( Table 1).
The 25 weather stations' data from 1967 to 2016 in and around SRYR were obtained from China Meteorological Administration (Table 1). First, the daily precipitation and potential evaporation of 25 meteorological stations were obtained by distribution calculation, and then the monthly precipitation and potential evaporation of 25 meteorological stations were calculated. The monthly precipitation and potential evaporation with a resolution of 1 × 1 km were interpolated from 25 meteorological stations' data in and around the SRYR.

Data Sources
The observed monthly runoff data in the period 1967-2016 in Tangnaihai station were obtained from the Yellow River Water Conservancy Commission (Table 1). The 25 weather stations' data from 1967 to 2016 in and around SRYR were obtained from China Meteorological Administration (Table 1). First, the daily precipitation and potential evaporation of 25 meteorological stations were obtained by distribution calculation, and then the monthly precipitation and potential evaporation of 25 meteorological stations were calculated. The monthly precipitation and potential evaporation with a resolution of 1 × 1 km were interpolated from 25 meteorological stations' data in and around the SRYR.

Mutation Analysis Methods
Mann-Kendall mutation analysis has been widely used to test abrupt point of hydrometeorological elements due to its advantages of no interference from a few outliers and simple calculation [46][47][48][49]. The cumulative anomaly method is a nonlinear statistical method which can directly judge the data change trend from the curve. Its calculation process is simple, and it is also widely used to test the abrupt change points of hydrometeorological elements [50,51].
For distinguishing the mutation year of runoff time series data as accurately as possible, this paper uses the Mann-Kendall and cumulative anomaly methods to analyze abrupt points of the runoff time series data in Tangnaihai hydrological station, and then by combining the results of the two methods, the runoff time series data of Tangnaihai hydrological station were divided into several periods.

Monthly ABCD Hydrological Model
The monthly ABCD model proposed by Thomas in 1981 is a simple and efficient four parameter conceptual hydrological model. The main input data of the model are monthly rainfall data, monthly potential evapotranspiration data, and monthly runoff depth observation data, and the output data are monthly runoff depth simulated data, monthly soil water content data, monthly groundwater content data, and monthly actual evaporation data. The basic principle of the monthly ABCD model is the water balance principle.
The water balance equation in the soil aquifer can be expressed as: where P t is the monthly rainfall (mm), ET t is the actual monthly evapotranspiration (mm), DR t is the surface direct runoff (mm), GR t represents groundwater recharge (mm), and S t and S t1 represent the soil moisture content (mm) of the current month and the last month. The probable evapotranspiration (Y t ) is the maximum water amount leaving the basin in the form of evapotranspiration. The effective water (W t ) is the sum of Y t and soil aquifer discharge (including surface runoff and groundwater recharge).
The nonlinear functional relationship between Y t and W t is shown in Equation (4): where parameter a is the probability of runoff formation before the soil water content is fully saturated and parameter b is the upper limit of unsaturated aquifer water storage capacity. The range of parameter a is (0, 1) and the range of parameter b is (0, 1000). The curve of possible evapotranspiration and effective water quantity in the "ABCD" model is shown in Figure 2 With the increase of W t , Y t will be more and more close to b. Y t will remain unchanged when it reaches b, so b can be considered as the upper limit of unsaturated aquifer water storage. When Y t reaches b, the excess rainfall will be converted into surface runoff and groundwater recharge. The curve of possible evapotranspiration and effective water quantity in the "ABCD" model is shown in Figure 2. Yt ≤ Wt, Y't (0) = 1, Y't (∞) = 0. With the increase of Wt, Yt will be more and more close to b. Yt will remain unchanged when it reaches b, so b can be considered as the upper limit of unsaturated aquifer water storage. When Yt reaches b, the excess rainfall will be converted into surface runoff and groundwater recharge. The ratio between the decrease rate of soil water content (S) caused by evapotranspiration and the potential evapotranspiration is assumed to be St/b in the "ABCD" model.
where PETt is the potential evapotranspiration (mm). PETt can be calculated using the Penman-Monteith formula.
The water balance equation in the groundwater aquifer can be expressed as: where GDt is underground runoff, GRt is groundwater recharge, and Gt and Gt1 are the groundwater storage in the current month and the last month, respectively. GDt and GRt can be expressed as follows: =d t t GD G (10) where parameter c is the proportion of groundwater recharge from the soil aquifer and parameter d is the rate of groundwater outflow. The range of parameter c is (0, 1) and the range of parameter d is (0, 1). The Nash efficiency coefficient (NS) and root mean square error (RMSE) and relative error (RE) are selected as the evaluation indexes in the calibration and validation of hydrological model parameters. The ratio between the decrease rate of soil water content (S) caused by evapotranspiration and the potential evapotranspiration is assumed to be S t /b in the "ABCD" model.
where PET t is the potential evapotranspiration (mm). PET t can be calculated using the Penman-Monteith formula.
The water balance equation in the groundwater aquifer can be expressed as: where GD t is underground runoff, GR t is groundwater recharge, and G t and G t1 are the groundwater storage in the current month and the last month, respectively. GD t and GR t can be expressed as follows: where parameter c is the proportion of groundwater recharge from the soil aquifer and parameter d is the rate of groundwater outflow. The range of parameter c is (0, 1) and the range of parameter d is (0, 1). The Nash efficiency coefficient (NS) and root mean square error (RMSE) and relative error (RE) are selected as the evaluation indexes in the calibration and validation of hydrological model parameters.
Land 2021, 10, 542 6 of 14 In the formulas, Q i represents the runoff simulation value at time i, P i represents the runoff observation value at time i, and P represents the average value of runoff observation values.

Decomposition Method Based on Seasonal Budyko Framework
Chen et al. put forward the definition of effective precipitation, and then the applicability of the Budyko framework was extended from the annual scale to the seasonal scale [52]. The effective precipitation is equal to the precipitation minus the change value of water storage. The seasonal Budyko framework must consider the change of water storage, but there is no measured value between the change of water storage and the actual evapotranspiration. In this paper, the monthly ABCD model was applied for simulating the monthly water storage change and actual evapotranspiration, which were superimposed on the season to construct the seasonal Budyko framework.
The formula of the seasonal Budyko framework is shown as the following: P represents precipitation, E represents actual evaporation, E p represents potential evaporation, ∆S represents the change of storage water, ϕ represents the corresponding lower bounds of aridity indices, and ω represents the underlying surface parameter.
Wang et al. put forward a decomposition method based on the Budyko hypothesis [53] and this method is simple in structure, easy to use, and does not need a lot of hydrometeorological data. Sun (18) where ΔR represents the difference between the observed runoff in the pre-change period and the observed runoff in the post-change period, and R c  and R h  represent the contribution rate of climatic and anthropic factors on runoff variation, respectively.

Mutation Analysis
In this paper, the Mann-Kendall and cumulative anomaly method were applied for distinguishing the mutation year of runoff time series data in Tangnaihai (Figures 4 and 5).  The runoff change caused by anthropic factors (∆R h ) and caused by climate change (∆R c ) can be calculated by: where ∆R represents the difference between the observed runoff in the pre-change period and the observed runoff in the post-change period, and η R c and η R h represent the contribution rate of climatic and anthropic factors on runoff variation, respectively.

Mutation Analysis
In this paper, the Mann-Kendall and cumulative anomaly method were applied for distinguishing the mutation year of runoff time series data in Tangnaihai (Figures 4 and 5).
tribution rate of climatic and anthropic factors on runoff variation, respectively.

Mutation Analysis
In this paper, the Mann-Kendall and cumulative anomaly method were applied for distinguishing the mutation year of runoff time series data in Tangnaihai (Figures 4 and 5).  As can be seen from Figure 4, the UF and UB curve intersected in 1969, 1989, 2009, and 2015, and all of them are within the range of the 0.05 significance level. As can be seen from Figure 5, the change trend of the cumulative anomaly value of runoff in the Tangnaihai station changed in 1989. Combining the results of two mutation analysis methods, As can be seen from Figure 4, the UF and UB curve intersected in 1969, 1989, 2009, and 2015, and all of them are within the range of the 0.05 significance level. As can be seen from Figure 5, the change trend of the cumulative anomaly value of runoff in the Tangnaihai station changed in 1989. Combining the results of two mutation analysis methods, the change point of annual runoff series data in Tangnaihai hydrological station from 1967 to 2016 is 1989.

Parameter Calibration and Verification of Monthly ABCD Model
According to the analysis results in Section 4.1, the period 1967-2016 was segmented into two periods: the pre-change period (T1: 1967-1989) and the post-change period (T2: 1990-2016). The T1 period and the T2 period are both divided into the calibration period and the validation period, respectively, and the specific time division is shown in Tables 2 and 3.  Table 3. Monthly runoff simulation results in the post-change period using the monthly ABCD model. The ABCD model was calculated by genetic algorithm in MATLAB; the value of parameter a, b, c, and d in the calibration period (1967. 1-1978.6) in the T1 period were calculated ( Table 2). The NS of the calibration period (1967. 1-1978.6) in the T1 period is 0.80, the RMSE is 5.11, and the RE is −3.62% (Table 2). Then, the value of parameter a, b, c, and d in the calibration period (1967. 1-1978.6) in the T1 period were used to simulate the monthly simulated runoff depth data in the validation period (1978. 7-1989.12). The NS of the validation period (1978. 7-1989.12) in the T1 period is 0.82, the RMSE is 5.84, and the RE is 3.37% ( Table 2).

Hydrological
The ABCD model was calculated by genetic algorithm in MATLAB; the value of parameter a, b, c, and d in the calibration period (1990. 1-2003.6) in the T2 period were calculated ( Table 3). The NS of the calibration period (1990. 1-2003.6) in the T2 period is 0.79, the RMSE is 3.76, and the RE is 0.31% (Table 3). Then, the values of parameter a, b, c, and d in the calibration period in the T2 period were used to simulate the monthly simulated runoff depth data in the validation period (2003.7-2016.12). The NS of the validation period (2003. 7-2016.12) in the T2 period is 0.80, the RMSE is 4.59, and the RE is −1.65% (Table 3).
Therefore, the monthly ABCD model constructed in this paper can well simulate the monthly runoff changes of Tangnaihai hydrological stations in T1  and T2 , and was applied for simulating monthly water storage change and actual evaporation in SRYR. The monthly runoff simulation results of Tangnaihai hydrological station in T1  and T2 (1990-2016) are displayed in Figures 6 and 7. Table 3. Monthly runoff simulation results in the post-change period using the monthly ABCD model.

Parameter Fitting of Budyko Curve
According to the distribution of seasonal precipitation in SRYR, the annual runoff data are segmented into four seasons: spring, summer, autumn, and winter. Before applying the vertical decomposition method to quantitatively calculate the contribution rate of climatic and anthropic factors on runoff variations, it was necessary to fit the parameter of Budyko curves in the four seasons in the pre-change period .
The parameter fitting results of the Budyko curve in the pre-change period of four seasons are displayed in Figure 8.

Parameter Fitting of Budyko Curve
According to the distribution of seasonal precipitation in SRYR, the annual runoff data are segmented into four seasons: spring, summer, autumn, and winter. Before applying the vertical decomposition method to quantitatively calculate the contribution rate of climatic and anthropic factors on runoff variations, it was necessary to fit the parameter of Budyko curves in the four seasons in the pre-change period .
The parameter fitting results of the Budyko curve in the pre-change period of four seasons are displayed in Figure 8.

Parameter Fitting of Budyko Curve
According to the distribution of seasonal precipitation in SRYR, the annual runoff data are segmented into four seasons: spring, summer, autumn, and winter. Before applying the vertical decomposition method to quantitatively calculate the contribution rate of climatic and anthropic factors on runoff variations, it was necessary to fit the parameter of Budyko curves in the four seasons in the pre-change period .
The parameter fitting results of the Budyko curve in the pre-change period of four seasons are displayed in Figure 8. The values of ϕ in the pre-change period are 0.46 in spring, 0.27 in summer, 0.23 in autumn, and 0.52 in winter, respectively. The values of ω in the pre-change period are 1.08 in spring, 1.40 in summer, 1.14 in autumn, and 1.04 in winter, respectively. The values of ϕ and ω reflect the effects of vegetation, soil, and topography on the evaporation and water storage in different seasons in SRYR [52]. In addition, most of the purple dots are above the blue dots, indicating that anthropic factors have a reducing impact on runoff of the four seasons in the post-change period (1990-2015).

Quantitative Analysis of Seasonal Runoff by Climatic and Anthropic Factors
The eigenvalues of climatic and hydrological variables in the pre-change period  and the post-change period (1990-2016) are displayed in Table 4. Compared with the pre-change period , precipitation in spring, summer, autumn, and winter of the post-change period (1990-2016) changed by 1.96, −4.86, −3.91, and 1.93 mm, respectively; the actual evaporation increased by 8.69 mm in spring, 15.51 mm in summer, 8.01 mm in autumn, and 3.34 mm in winter; the runoff depth decreased by 4.12 mm in spring, 12.18 mm in summer, 18.9 mm in autumn, and 1.05 mm in winter, respectively. Water storage changed by −2.61 mm in spring, −8.19 mm in summer, 6.98 mm in autumn, and −0.36 mm in winter, respectively.  Table 5 displays the contribution rate of climatic and anthropic factors on seasonal runoff variation in the post-change period. Compared with the pre-change period , the runoff in the post-change period (1990-2016) decreased by 1.22 mm in spring, 6.06 mm in summer, 15.01 mm in autumn, and 0.43 mm in winter due to climatic factors, accounting for 29.72%, 49.75%, 79.45%, and 40.89%, respectively. The runoff variation caused by anthropic factors in the T2 period is −2.90 mm in spring, −6.12 mm in summer, −3.88 mm in autumn, and −0.62 mm in winter, accounting for 70.28%, 50.25%, 20.55%, and 59.11%, respectively. In summary, anthropic factors played a major role in runoff change in spring, summer, and winter, while climatic factors played a major role in runoff change in autumn.

Conclusions and Discussion
In this paper, based on the monthly ABCD hydrological model and the seasonal Budyko framework, the contribution rate of climatic and anthropic factors on the seasonal runoff variation in SRYR was quantitatively calculated. It turned out that the ABCD monthly hydrological model could well simulate the monthly runoff variation of Tangnaihai hydrological station in T1  and T2 (1990-2016) periods. The contribution rate of climatic and anthropic factors on the runoff variation in spring, summer, autumn, and winter are 29.72% and 70.28%, 49.75% and 50.25%, 79.45% and 20.55%, and 40.89% and 59.11%, respectively.
The main reason for the higher contribution rate of anthropic factors in spring is that the spring crops (wheat and rapeseed) require a large amount of water, which extracted from the river basin for irrigation. In summer, the runoff is relatively large, and Huangheyuan hydropower station has a great influence on the regulation of river runoff, especially when the hydropower station is in full load operation for power generation and preventing flood, resulting in the decrease of the runoff at Tangnaihai in summer. In winter, the runoff is relatively small; reservoirs and hydropower stations often regulate and control water levels for power generation, resulting in the decrease of the runoff at Tangnaihai. In autumn, the water level can meet the basic requirements of reservoirs and hydropower stations; its regulation mode is "release as much as possible", which indicates that the impact of anthropic factors on the runoff change at Tangnaihai in autumn is relatively small.
Although the data used in this study were carefully controlled, there were still some uncertainties. Precipitation and potential evapotranspiration data were interpolated from low-density station data from national meteorological stations in and around the SRYR. Although the theoretical basis of the monthly ABCD hydrological model is rigorous, there are still some uncertainties when using mathematical formulas to describe the process of runoff generation. Additionally, the monthly ABCD hydrological model is a lumped hydrologic model, which cannot show the influence of the uneven distribution of the underlying surface conditions and precipitation on the runoff formation process.