Analysis of runoff variation and driving mechanism in Huangfuchuan River Basin in the middle reaches of the Yellow River, China

To address the prominent problem of declining runoff in many rivers around the world, studying the law of runoff change and attribution analysis is very important for the planning and management of watershed water resources and has practical significance for solving the imbalance between supply and demand of watershed water resources and maintaining the healthy development of rivers. Three commonly used coupled water-energy balance equations based on Budyko hypothesis are selected to estimate the elasticity coefficient of runoff change to each driving factor, and the contribution rate of different factors to runoff change in the study area is quantified by the total differential method and the complementary method, respectively. The results show that the runoff of Huangfuchuan River basin showed a significant decreasing trend from 1954 to 2015, and the runoff mutation points were 1979 and 1996; in the alteration period I (1979–1996), precipitation was the main factor leading to the runoff reduction in Huangfuchuan River basin, followed by the influence of underlying surface; the contribution rate of underlying surface to runoff alterations ranged from 63.7% to 65.46%; the impact of potential evapotranspiration was slightly smaller. In the alteration period II (1997–2015), the underlying surface played a dominant role in runoff reduction of Huangfuchuan River basin. The contribution rate of the underlying surface to runoff change ranged from 80.21% to 86.34%, followed by precipitation, and the potential evapotranspiration had the least impact. The impact of human activities on the whole watershed increased with the passage of time. The land use change, the overall increase of NDVI (vegetation cover) and the construction of water conservation projects are important reasons for the reduction of runoff in Huangfuchuan River basin.


Introduction
Climate variability and human activities have affected the regional hydrological cycle process to varying degrees (Abbott et al. 2019;Krajewski et al. 2021;Lakshmi et al. 2018;Rodriguez-Iturbe 2000). The runoff of many rivers in the world has decreased significantly under the dual effects of climate change and human activities (Huntington 2006;Silberstein et al. 2012;Zhang et al. 2017). For example, from the late 1970s, the runoff of the Yellow River also showed a significant reduction trend (Gao et al. 2011), resulting in the increasingly tense relationship between supply and demand of water resources in the Yellow River Basin . Studying the law and characteristics of runoff change, and analyzing the attribution and influence mechanism of runoff change can provide basis for the optimal allocation, scientific dispatching and rational utilization of water resources in the basin, which is of great significance to solve the imbalance between supply and demand of water resources in the basin and maintain the healthy development of river (Abera et al. 2019).
In recent years, scholars from various countries have conducted a lot of research on runoff variation and attribution analysis (Guo et al. 2020). Miao et al. (2011) studied the correlation between precipitation and runoff in the base period at seven river gauging stations since 1970 by using the hydrographic method, and clarified the contribution rate of climate change and human activities to the reduction of runoff in the upper, middle and lower reaches of the Yellow River; Pirnia et al. (2019) applied the double cumulative curve method to monitor and analyze the impact of climate variability and human activities changes in stream flow of the Haraz River basin in northern Iran from 1975 to 2010; Patterson et al. (2013) used Budyko curves to ascribe changes in streamflow due to climate and human factors between two time periods in both natural and human-modified basins in the South Atlantic; Hu et al. (2021) assessed the long-term  runoff dynamics of Amu Darya River basin applying the elasticity coefficient method coupled with the water balance and Budyko framework; Santra Mitra et al. (2021) quantified the relative contribution of climate change and land-use change to runoff response in the Kangshabati River basin of West Bengal by using SWAT model; Lv et al. (2022) calculated the effects of climate variation and human activities on streamflow change in a poorly gauged arid inland basin based on SWAT and Budyko-based Mezentsev-Choudhurdy-Yang model driven by multi-source precipitation datasets.
Huangfuchuan River basin is the first-class tributary of the middle reaches of the Yellow River. It is located at the junction of the Loess Plateau and the Ordos Plateau. It belongs to the transition zone of water erosion and wind erosion and is the concentrated source of coarse sediment in the sandy and coarse sediment area of the middle reaches of the Yellow River. The arsenic sandstone in the river basin is exposed in a large area, and the original vegetation is seriously damaged. These characteristics are extremely rare in China and even in the world. In addition, Huangfuchuan River is a typical semi-arid River, and water is the main factor limiting the vegetation restoration and ecological construction of the basin. Huangfuchuan River basin is currently facing the environmental problems of soil erosion and drought and water shortage. Therefore, it is typical and representative to study the impact of climate change and human activities on runoff in this basin, which can also provide data support and theoretical help for local water resources planning and management. It is of great practical significance.
Taking Huangfuchuan River basin in the middle reaches of the Yellow River as the research object, this paper analyzes the runoff evolution from 1954 to 2015 by using linear tendency estimation, Mann-Kendall test, Sen's slope estimation, moving average method, double cumulative curve method and cumulative anomaly method. Three commonly used coupled water-energy balance equations under Budyko hypothesis are selected to analyze the attribution of runoff changes in Huangfuchuan River basin, and the contribution rates of different factors to runoff changes in the study area are quantitatively estimated by the total differential method and the complementary method, respectively, which provides a scientific basis for water resources management and decision-making in the Yellow River basin and has guiding significance for the optimal allocation of water conservation measures in the Loess Plateau.

Study area
Huangfuchuan River basin ( Fig. 1) (110°18'-112°12' E, 39°12'-39°54'N) is located in the upper part of the middle reaches of the Yellow River. It is located in the transition zone between Ordos Plateau and Loess Plateau. It flows through Jungar banner of Inner Mongolia Autonomous Region. The main stream is 137 km long and covers an area of 3246 km 2 . Huangfuchuan belongs to the first-class tributary of the Yellow River, which is composed of Shilichangchuan, Ganchabangou, Hushigou, Geqiugou and other main tributaries. The two main tributaries of Nalinchuan and Shilichangchuan are evenly distributed on both banks of the basin and flow into Huangfuchuan in an approximate parallel, which is a parallel stream. Huangfuchuan River basin is located in the inland arid and semi-arid area. Affected by the southeast monsoon, it belongs to continental monsoon climate, with an annual average temperature of 9. 1 ℃; precipitation varies greatly between and within years. The average annual precipitation is 350~450 mm, which is concentrated from June to September, accounting for more than 80% of the annual precipitation. The temperature and precipitation decrease from southeast to northwest in spatial distribution. Huangfuchuan is a river with many and coarse sands, with an average annual sediment discharge of 61.2 million tons. The Huangfuchuan hydrologic station is adjacent to the outlet and controls approximately 98.5% of the river basin (3199 km 2 ). Since 1954, daily scale hydrological observation data can be obtained in this river basin (Shi et al. 2015).

Data sources
The measured runoff data from 1954 to 2015 of Huangfu hydrological station, the outlet station of Huangfuchuan River basin, are selected as the hydrological data, and the data are from the hydrological yearbook of the Yellow River Basin; the daily meteorological data of Hequ station in the same period are selected as the meteorological data, and the data are from China Meteorological data network (http:// data. cma. cn); there are 13 rainfall stations in the basin, the interpolation method of selecting adjacent stations is used to interpolate the missing measurement data, and the Thiessen polygons module of the Arc-gis10.3 software is used to finally obtain the basin surface Page 3 of 16 234 precipitation data; the Penman-Monteith formula modified by FAO is used to calculate the potential evapotranspiration. The calculation formula is as follows: where ET 0 is the potential evaporation (mm/d); ∆ is the slope of the saturation water pressure curve (kPa/℃); R n is the net surface radiation (MJm -2 d -1 ); G is the soil heat flux (MJm -2 d -1 ); γ is the dry and wettable constant (kPa/℃); U 2 is the wind speed at 2m (m/s); T is the average daily air temperature (℃); e s is the saturation water pressure (kPa); e a is the actual water pressure (kPa); (e s -e a ) is the vapor pressure deficit (kPa). (1)

Trend and mutation analysis
Linear tendency estimation method, Mann-Kendall trend test method, Sen's slope estimation method and moving average method were adopted to analyze the trend of runoff in Huangfuchuan River basin from 1954 to 2015. The linear tendency estimation method reflects the trend of hydrological factors over time by establishing a linear regression equation with unary. The Mann-Kendall nonparametric statistical method recommended and widely used by the World Meteorological Organization (WMO) can effectively distinguish whether a natural process is in natural fluctuation or has a definite change trend and determine the significance of the change trend of hydrological time series (Das and Banerjee 2021;Tunde et al. 2022;Zhang et al. 2019). Sen's slope estimation method, also known as Theil-Sen estimator, is a nonparametric method for estimating the trend slope in data samples. It has the advantages of high calculation efficiency and insensitive to outliers. It is often used in the trend analysis of long-time series data. The moving average method is equivalent to low-pass filtering, filtering out the highfrequency components of the sequence and highlighting the long-term variation trend (Assal et al. 2021;Caloiero and Aristodemo 2021).
Mann-Kendall test, double-mass curve (DMC) method and cumulative anomaly method were used to synthetically diagnose the change points of runoff in Huangfuchuan River basin from 1954 to 2015 to reduce omissions and errors from only using a single method. Mann-Kendall mutation test is a nonparametric statistical test method (Xing et al. 2018). Its advantage is that it does not need the sample to follow any distribution and is not disturbed by a few outliers. It is more suitable for type variables and sequential variables. It can clarify the starting time of mutation and point out the mutation area. Therefore, it has become a widely used mutation test method. The double-mass curve is a method used to investigate the consistency, long-term trend and detect mutation points of hydrometeorological time series. Its basic principle is that two variables are gradually accumulated according to the same time length, one variable is the abscissa and the other variable is the ordinate, so as to describe the trend changes of the two (Rozbeh et al. 2019). If the proportion between the two variables remains unchanged, it will show a linear relationship within the same time. If the slope changes, it indicates that the original relationship between the two variables has changed. The point where the linear slope deviates significantly correspond to the year when the y-axis variable begins to change significantly. The cumulative anomaly method compares the discrete data of the data series with the average value of the data series. When the discrete data are larger than the average value of the series, the cumulative anomaly increases and the cumulative anomaly curve shows an upward trend. On the contrary, the curve shows a downward trend. According to the slope change, obvious mutation points can be detected and the time of mutation of the data series can be judged (Liu et al. 2020).

Budyko hypothesis
Budyko theory (Du et al. 2016;Greve et al. 2016;Lhomme and Moussa 2016;Malagò et al. 2018;Moussa and Lhomme 2016) is based on the water balance equation, and the expression is: where R is the annual average runoff, mm; P is the annual average precipitation, mm; ET is the annual average actual evapotranspiration, mm; and △S is the variation of water storage in the basin, mm.
Budyko believes that the multi-year evapotranspiration rate of a basin in a long-time scale is determined by the potential evapotranspiration and precipitation. Therefore, Budyko hypothesis can be simply expressed as: where ET 0 is the multi-year average potential evapotranspiration of the basin, mm.

Coupled water-energy balance equations based on budyko hypothesis
Under certain climate and underlying surface conditions, the long-term water-energy balance relationship of the basin can be expressed by the equation describing the relationship between precipitation, actual evapotranspiration and potential evapotranspiration, that is, the coupled water-energy balance equations. Based on Budyko hypothesis, many scholars have proposed different forms of modified Budyko models. Among them, Fu (1981); Choudhury-Yang (Choudhury 1999;Yang et al. 2008); and Zhang et al. (2001), the three single parameter Budyko equations are widely used. Their descriptions are shown in Table 1: Table 1 Coupled water-energy balance equations n represents the underlying surface parameters of the basin, which is mainly related to the underlying surface characteristics such as terrain, soil and vegetation characteristics; ω indicates the plant-available water efficiency, which is mainly related to the vegetation characteristics of the underlying surface Serial No.
Proposer Equations

Contribution rate of climate and underlying surface to runoff change
(1)Total differential method According to the mutation point, the study period is divided into the base period and the alteration period. Assuming that the annual average runoff depth in the base period and alteration period is R 1 and R 2 , respectively, the change of annual runoff from the base period to the alteration period can be expressed by the difference dR of annual average runoff depth in the two periods, that is: The runoff change caused by factors can be estimated by the product of factor change and its partial derivative. The calculated runoff change can be obtained by calculating the contribution of each factor to runoff change, respectively: where dR * is the calculated runoff variation, mm; The expression of the above formula in the form of total differential is the total differential method. Ignoring the higherorder term of Taylor expansion here, we can get: where ΔR P 、ΔR ET 0 、ΔR n∕ are estimated precipitation, potential evapotranspiration and runoff variations caused by underlying surface; ΔP、ΔET 0 and Δn∕ are the changes of annual average precipitation, potential evapotranspiration and underling surface parameters from the base period to the alteration period.
The contribution rate of each impact factor to runoff change is: where x represents P、ET 0 、 or n∕ , C x represents the contribution rate of each impact factor to runoff change.
(2)Complementary method Zhou et al. (2016) assumed that precipitation (P) and potential evapotranspiration (ET 0 ) are relatively independent and proposed that the elastic coefficient of runoff to precipitation and to potential evapotranspiration have a complementary relationship, that is: Due to the uncertainty of climate and underlying surface conditions, and considering the complementary relationship between precipitation and potential evapotranspiration elastic coefficients, the complementary method of introducing weight factor α to divide runoff variation is proposed. The expression is as follows: where ΔR is the variation of runoff, mm; α is the weight factor, 0 ≤ α ≤ 1, subscripts 1 and 2 represent the base period and the alteration period before and after the mutation point.
Based on the above formula, the runoff change can be divided as follows: where ΔR P is the contribution of precipitation change to runoff change; ΔR ET 0 is the contribution of potential evapotranspiration change to runoff change; ΔR n is the contribution of underlying surface change to runoff change.
The contribution rate of each impact factor to runoff change is: where x represents P、ET 0 、or n∕ , C x represents the contribution rate of each impact factor to runoff change.

Trend test
In order to study the variation and trend of runoff in Huangfuchuan area in 62 years, the annual runoff series was analyzed. Figure 2 shows the annual precipitation, corresponding annual runoff depth, linear trend and 5a moving average curve of Huangfuchuan River basin from 1954 to 2015. From 1954 to 2015, the annual average runoff depth of Huangfuchuan River basin was 38.65mm. Through linear trend analysis, it was found that the runoff of Huangfuchuan River basin showed a downward trend, and the reduction rate of runoff depth was −1.1541mm/a. The Z value of M-K trend test is −5.78 ( |Z| > 2.576) , which has passed the significance test of 99% confidence, indicating that the runoff of Huangfuchuan River basin has shown a significant decreasing trend in recent 60 years. Sen's slope estimation β=− 0.0292, indicating that the runoff of Huangfuchuan River basin shows a downward trend. The runoff of Huangfuchuan River basin is taken as 5a for moving average, which weakens the impact of sequence high-frequency oscillation, that is, the influence of runoff in special years on the variation trend of runoff is weakened. To sum up, the runoff of Huangfuchuan River basin shows an obvious downward trend. The conclusion of trend analysis is shown in Table 2.

Mutation test
M-K mutation test is used to test the mutation of Huangfuchuan runoff series from 1954 to 2015, as shown in Fig. 3. The unique intersection of the runoff statistics of the two curves of UF and UB is at the 0.01 significance level, which indicates that the runoff of Huangfuchuan River basin may have changed abruptly in 1996. In order to further determine the mutation year, the double-mass curve method is used to carry out the mutation test on the annual runoff series of Huangfuchuan, as shown in Fig. 4. From the figure, it can be found that the precipitation-runoff cumulative curve of Huangfuchuan River basin showed abrupt alterations in 1979 and 1996. The stage before the deviation of precipitation-runoff double-mass curve is generally regarded as the period when the runoff is not disturbed or less disturbed by human activities, that is,  the base period. The slope of the cumulative curve after the base period has changed to a certain extent and the slope of the curve has become smaller, indicating that the runoff has a decreasing trend relative to precipitation. The years 1979-1996 and 1997-2015 are divided into alteration periods. The decision coefficient R 2 of the fitting relationship in each stage is more than 0.96, and the fitting degree is good. Finally, the mutation point of Huangfuchuan annual runoff series is verified by the cumulative anomaly method, as shown in Fig. 5. From the change trend of runoff cumulative anomaly curve, it can be seen that the runoff change of Huangfuchuan can be divided into three stages: from 1954 to 1978, the cumulative anomaly curve showed a fluctuating upward trend, indicating that Huangfuchuan River basin was in a relatively wet period at this stage; the cumulative anomaly curve from 1979 to 1996 fluctuated slightly, but tended to be flat as a whole, indicating that Huangfuchuan River basin was in the medium period at this stage; the cumulative anomaly curve from 1997 to 2015 showed a significant downward trend, indicating that Huangfuchuan River basin is in a dry period at this stage.   19961979, 19961979, 1996Conclusion 1979, 1996 234 Page 8 of 16 The conclusion of mutation analysis is shown in Table 3.

Elastic analysis of runoff on climatic elements and underlying surface
According to the results of mutation analysis, the study period sequence is divided into base period , which are not affected by human activities or less affected by human activities, alteration period I (1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996) and alteration period II (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015). According to the collected and calculated annual average precipitation, annual average runoff depth and annual average potential evapotranspiration in the base period and alteration period, calculate the corresponding underlying surface parameters in each period, and calculate the elasticity coefficient corresponding to the base period and alteration period. The results are shown in Table 4. As can be seen from Table 4, among the hydro-meteorological factors affecting the runoff change in Huangfuchuan River basin, the precipitation, runoff, evaporation and runoff coefficient (R/P) in the two alteration periods (1979-1996 and 1997-2015) decreased significantly compared with the base period , and the largest decrease was the annual runoff, with the decline of the two alteration periods relative to the base period reaching 27.32% and 79.68%, respectively; the second is the runoff coefficient, with a decrease of 17.06% and 76.26%, respectively; the precipitation decreased by 12.37% and 14.41%, respectively; The potential evapotranspiration decreased by 13.58% and 9.77%, respectively. Except potential evaporation, the variation amplitude of other factors in alteration period II was greater than that in alteration period I. The basin drought index (ET 0 /P) decreased slightly (1.38%) in the alteration period I and rebounded again (5.42%) in the alteration period II. For the underlying surface parameter n/ω characterizing the impact of human activities, the alteration periods of each equation are significantly increased compared with the base period, and the increase in alteration period II is huge. To sum up, the change trend of relevant factors in each equation is basically the same. Both climate factors and human activity factors tend to indicate that there is a trend of aridity in the basin.
The precipitation elasticity coefficient calculated by each equation is positive, while the potential evapotranspiration elasticity coefficient and underlying surface elasticity coefficient are negative. The results show that the runoff depth of Huangfuchuan is positively correlated with precipitation and negatively correlated with underlying surface changes and potential evapotranspiration, which is consistent with the actual hydrological law. The absolute value of elasticity coefficient reflects the sensitivity of runoff to various influencing factors. Take the elasticity coefficient obtained by Fu equation in the base period as an example, P 、 ET 0 and n are 2.11, − 1.11 and − 3.13, respectively, indicating that when P, ET 0 and n increase by 1%, the change of runoff depth are increased by 2.11%, decreased by 1.11% and decreased by 3.13%, respectively. Therefore, Huangfuchuan runoff is the most sensitive to the change of underlying surface and the least sensitive to potential evapotranspiration.

Attribution analysis of runoff change
The total differential method and complementary method are used to analyze the attribution of runoff change in Huangfuchuan River basin. The two methods are compared. Table 5 shows the quantitative impact of various factors on runoff change when the attribution analysis of total differential method is adopted. The impact of various coupled water-energy balance equations on runoff change is basically the same.

Total differential method
In the alteration period I, the range of runoff change caused by precipitation reduction calculated by the three equations is − 13.91~ − 15.51mm, and the contribution rate is 91.21%~97.22%; the potential evapotranspiration leads to the increase range of Huangfuchuan runoff is 8.37~10mm, and the contribution rate is -54.91%~ − 62.68%; the range of runoff reduction caused by the underlying surface is − 9.71~ − 10.45mm, and the contribution rate is 63.70%~65.46%. Therefore, in the alteration period I, precipitation is the main factor leading to the reduction of runoff in Huangfuchuan, followed by the influence of  underlying surface, and the impact of potential evapotranspiration is slightly smaller.
In the alteration period II, the range of runoff variation caused by precipitation reduction calculated by the three equations is − 10.61~ − 12.44mm, and the contribution rate is 22.20%~32.73%; the increase range of runoff in Huangfuchuan River basin caused by potential evapotranspiration is 4.08~4.92mm, and the contribution rate is − 8.54%~ − 12.94%; the range of runoff reduction of Huangfuchuan River basin caused by the underlying surface is − 30.49~ − 41.27mm, and the contribution rate is 80.21% ~ 86.34%. Therefore, in the alteration period II, the underlying surface is the leading factor conducing to the reduction of runoff in Huangfuchuan, followed by the influence of precipitation, and the potential evapotranspiration is the least.
The total differential method is one of the most widely used methods in the attribution analysis of hydrological cycle changes in the basin. However, due to the neglect of the high-order term of Taylor expansions in the calculation, there are certain errors. In this calculation, in the alteration period I, the relative errors between the runoff changes calculated by the three equations and the actual observed values are: − 2.30%, − 2.28% and 2.25%, respectively; the relative errors between the runoff variation calculated by the three the equations and actual observed values in the alteration period II are 5.02%, 4.51% and − 16.49%, respectively. The relative errors of alteration period I are less than those of alteration period II. Table 6 shows the quantitative impact of various factors on runoff variation when the attribution analysis of complementary method is adopted. The impact of various coupled water-energy balance equations on runoff change is basically the same.

Complementary method
In the alteration period I, the range of runoff change caused by precipitation reduction calculated by the three equations is − 14.08~ − 15.52mm, and the contribution rate is 90.08%~99.49%; the potential evapotranspiration leads to the increase range of Huangfuchuan runoff is 8.42~10mm, and the contribution rate is − 53.89%~ − 64.13%; the Table 5 Contribution rate of runoff variation in Huangfuchuan River basin calculated by total differential method ΔR P 、ΔR ET 0 、ΔR n∕ in the table represent the runoff change caused by precipitation, potential evapotranspiration and human activities, respectively; η P 、η ET 0 、η n∕ represent the contribution rate of precipitation, potential evapotranspiration and underlying surface to runoff change, respectively

Budyko
Base period Alteration period Runoff changes caused by various factors/(mm) Relative error (%) Contribution rate/(%) Fu 1954Fu -1978Fu 1979Fu -1996 range of runoff reduction caused by the underlying surface is − 9.97~ − 10.08mm, and the contribution rate is 63.81%~64.64%. Therefore, in the alteration period I, precipitation is the main factor leading to the reduction of runoff in Huangfuchuan, followed by the influence of underlying surface, and the impact of potential evapotranspiration is slightly smaller.
In the alteration period II, the range of runoff variation caused by precipitation reduction calculated by the three equations is − 11.75~ − 12.44mm, and the contribution rate is 25.85%~27.31%; the increase range of runoff in Huangfuchuan River basin caused by potential evapotranspiration is 4.44~4.92mm, and the contribution rate is − 9.76%~ − 10.80%; the range of runoff reduction of Huangfuchuan River basin caused by the underlying surface is − 38.04~ − 38.12mm, and the contribution rate is 83.49%~83.91%. Therefore, in the alteration period II, the underlying surface is the leading factor conducing to the reduction of runoff in Huangfuchuan, followed by the influence of precipitation, and the potential evapotranspiration is the least.
The complementary method can accurately divide the runoff change into precipitation influence amount, potential evapotranspiration influence amount and underlying surface change influence amount. There is basically no error, and the α = 0.5 is selected in this study for calculation. Compared with the total differentiation method, this method has a wide range of applicability and more accurate results. Therefore, the complementary method is a better choice in calculating the attribution analysis of changes in the hydrological cycle of the watershed.
In the complementary method, the runoff changes estimated by the three Budyko equations are very close to the actual runoff variation. In this calculation, the relative errors between the runoff changes calculated by the Budyko equations in the alteration period I and the actual observed values are 0.14%, 0.11% and 0.08%, respectively; the relative errors between the runoff variation calculated by the Budyko equations and the actual observed values in the alteration period II are 0.2%, 0.03% and 0.09%, respectively. It can be seen that the calculation accuracy of the three selected Budyko equations is very high. In addition, combining the results of total differential method and complementary method, it can be seen that in the alteration period I, precipitation is the main factor affecting the reduction of runoff, and the contribution rate of underlying surface to runoff change ranges from 63.70% to 65 46%. In the alteration period II, the underlying surface is the leading factor conducing to the reduction of runoff. The range of the contribution rate of the underlying surface to the runoff change is 80.21%~86.34%, that is, the contribution rate calculated by the three Budyko equations has little difference, which can also confirm each other that the three Budyko equations selected are relatively reliable in the attribution analysis of runoff change.
Comparing the attribution analysis results of runoff changes in the two stages of the alteration period, it can be seen that the impact of underlying surface on runoff change has increased significantly in the alteration period II (1997II ( -2015 compared with the alteration period I (1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996), indicating that the drastic underlying surface changes in recent 20 years have had a greater impact on the water cycle and available water of the river basin.
Some studies on runoff change have been carried out in Huangfuchuan River basin. For example, Zhou, et al. (2014) applied principal components regression method to quantify the relation of runoff to climate change (aridity index) and anthropogenic activities. It was concluded that during 1979-1998, the contribution rate of human activities to runoff reduction was 60%, of which soil conservation measures accounted for 28% and irrigation accounted for 31%. Wang, et al. (2012) proposed a method of slope change ratio of accumulative quantity (SCRAQ) to calculate the contributions of the precipitation and human activities. According to this method, the relative contribution rates of precipitation and human activities to the runoff change in the Huangfuchuan River basin were 36. 43% and 63.57% during 1980 to 1997, and 16.81% and 83.19% from 1998 to 2008, respectively. The above related studies are consistent with the conclusions of this research.

Verification of contribution rate results by double-mass curve method
The double-mass curve method can not only be used to test the mutation points of hydro-meteorological time series, but also be used to quantitatively calculate the contribution rate of climate change and human activities to runoff change. If the slope of the double-mass curve changes, it indicates that the underlying surface of the basin has changed, and human activities have an impact on runoff. The double-mass curve of runoff and precipitation in Huangfuchuan River basin is shown in Fig. 4. There are two inflection points in 1979 and 1996. The period before the inflection points is set as the base period that is not affected by human activities or less affected by human activities , and the periods after the inflection points are set as the alteration periods. The fitting slope of the curve in the alteration periods is less than the base period, indicating that human activities reduce the runoff of the basin.
The results of runoff change contribution rate in Huangfuchuan River basin by double-mass curve method are shown in Table 7. In the alteration period I, the contribution rate of climate change to the reduction of Huangfuchuan runoff is 29.08%, and the contribution rate of human activities to the reduction of Huangfuchuan runoff is 70.92%. In the alteration period II, the contribution rate of climate change to the reduction of Huangfuchuan runoff is 12.99%, and the contribution rate of human activities to the reduction of Huangfuchuan runoff is 87.01%. Compared with the coupled water-energy balance method based on Budyko hypothesis, the contribution rates of human activities are slightly larger, mainly because these two models are quite different in model assumptions and hydrological processes representation. However, the conclusions from both models conclude that runoff decline in Huangfuchuan River basin was dominated by human activities. Moreover, it is observed that the impact of human activities on runoff change is greater in alteration period II than in alteration period I; that is, the impacts of human activities on the whole basin are enhanced with time.

Impact mechanism of human activities on runoff reduction in Huangfuchuan River basin
The causes of runoff change in Huangfuchuan River basin are complex and multifactorial. Climate change and underlying surface change are the two most important aspects. According to the attribution analysis of runoff change in Huangfuchuan River basin, the change of underlying surface is the leading factor of runoff reduction. The change of underlying surface conditions is mainly reflected in land use change, vegetation cover change and the construction of water conservation projects (Ahn and Merwade 2017; Kazemi et al. 2019;Yonaba et al. 2021). In order to clarify the specific driving effect of underlying surface on runoff reduction, these aspects need to be studied.

Land use change analysis
This study analyzes the land use change in Huangfuchuan River basin in recent decades, including the change of land use quantity and land use type. Different land use types have different runoff yield and runoff concentration mechanisms, and their contribution to surface runoff concentration is different. According to the classification standard of land use status in China (GB/T21010-2007), land types can be divided into cropland, forestland, grassland, water bodies, urban & construction land and unused land. The land use types of Huangfuchuan River basin are shown in Fig. 6.

Fig. 6
Land use types of Huangfuchuan River basin in typical years the cropland and unused land showed a decreasing trend, and the forest land and grassland remained basically unchanged. In the alteration period II (1995~2015), cropland and grassland showed a decreasing trend, while forest land, water bodies, urban & construction land and unused land showed an increasing trend. Figure 7 is a histogram showing the percent stacked column of land use in Huangfuchuan River basin in different periods. It can be seen from the chart that the main land use type in Huangfuchuan basin is grassland. Generally speaking, the land use type in Huangfuchuan River basin is relatively stable and has little fluctuation. Table 9 shows the land use type transition matrix of Huangfuchuan River basin from 1978 to 2015. It can be seen from the figure that the main land use changes occurred as the transformation from cropland to forest and grassland. This is directly related to China's attention to ecological restoration in the West and the gradual deepening of water conservation projects. Especially in 1999, Huangfuchuan River basin began a largescale project of returning cropland to forest and grassland, which played a role in soil and water conservation in a certain sense. However, the growth trend of urban & construction land, the development of urbanization and the increase of urban construction land area have changed the vegetation types and damaged the ecological balance, which is also the reason for the reduction of runoff in Huangfuchuan River basin.

NDVI change analysis
Normalized difference vegetation index (NDVI) is an index that can effectively reflect the characteristics of surface vegetation cover. This index is widely used in the analysis and research of spatio-temporal changes of vegetation cover. The spatial distribution data set of China's annual vegetation   (Fig. 9). It can be clearly seen from the figure that since 2000, the NDVI value of the basin has shown an increasing trend, the ecological construction of the basin has gradually spread from the side near the Yellow River to the whole basin, and the NDVI has shown a growth trend of the whole basin. Generally speaking, the greater the vegetation cover, the greater the rainfall infiltration and the less the runoff. Therefore, the increase of vegetation cover is one of the main reasons for the decrease of runoff in Huangfuchuan River basin.

Water conservancy engineering measures analysis
The check dam construction will have a profound effect on the movement of surface channel water flow. Compared with the situation without dam, the peak flow of the basin with dam will be reduced, and the basin with dam has an obvious regulating-storing function on the runoff process. As shown in A large number of studies show that the causes of runoff change are usually attributed to climate change, land use change, vegetation cover change and water conservancy engineering measures. Although the main influencing factors are these aspects, combined with the diagnosis results of the mutation points in the early stage, the causes of runoff  variation in different stages are discussed, respectively, that is, the changes of these influencing factors are analyzed in each stage, and specific targeted conclusions are drawn.

Uncertainty analysis
There is great uncertainty in the attribution analysis of runoff change based on Budyko hypothesis, and its uncertain factors are mainly reflected in the following aspects: (1) the definition of underlying surface parameters in Budyko equations. The underlying surface parameters should be defined as the land cover that emphasizes the natural attributes of the land and the results of human activities, which is the comprehensive embodiment of the underlying surface characteristics of a region; that is, when attributing human activities, the underlying surface parameters represent the comprehensive effects other than precipitation and evaporation factors. The complexity of hydrological process and the compound effect of various soil and water conservation measures on the Loess Plateau make it difficult to accurately separate the attribution analysis of human activities to runoff change. Vegetation restoration, construction of terraces, check dams and other water conservation engineering measures will cause changes in the underlying surface characteristic parameters. At present, the impact of various water conservancy measures on the underlying surface parameters cannot be clearly quantitatively analyzed.
(2) The mathematical assumption that each influence factor is independent in the elasticity coefficient method. The premise of the elasticity coefficient method based on Budyko hypothesis is that each variable is independent of each other. For example, when calculating the partial derivative of the multi-year average potential evapotranspiration, it is required that the precipitation and underlying surface parameter changes are independent of each other. However, in fact, the underlying surface and the climate system are an organic whole, and various factors always interact with each other, there is no complete independence. Therefore, the model will produce some errors. (3) Budyko theory ignores the variation of water storage in the basin. In the water balance equation based on Budyko hypothesis, the change of water storage △S in the basin is 0 on the multi-year average scale, which obviously ignores the retention of runoff in the basin by water conservation engineering measures such as check dams and various water conservancy projects such as reservoirs. Therefore, the quantitative study on the contribution rate of various water conservation measures to runoff change is still a problem that needs to be further solved in the attribution analysis of water and sediment change in the future, which is of great significance to the control of soil and water loss in the Loess Plateau.
(2) The elasticity coefficient of runoff on each influencing factor is calculated according to coupled water-energy balance equations based on Budyko hypothesis. The results show that Huangfuchuan runoff is the most sensitive to the change of underlying surface. The results of runoff change attribution analysis show that in the alteration period I (1979)(1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996), precipitation was the main factor leading to the reduction of runoff in Huangfuchuan River basin, followed by the influence of underlying surface, and the contribution rate of underlying surface to runoff change ranges from 63.7% to 65.46%; the impact of potential evapotranspiration was slightly smaller. In the alteration period II (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015), the underlying surface was the main factor leading to the reduction of Huangfuchuan runoff. The contribution rate of the underlying surface to runoff change ranged from 80.21% to 86.34%, followed by precipitation, and the potential evapotranspiration had the least impact. Compared with alteration period I, the impact of underlying surface on runoff change increased significantly, that is, the impact of human activities on the whole basin increased with the passage of time.
(3) In the attribution analysis of runoff variation, the complementary method has the characteristics of wide applicability and more accuracy than the total differential method. Therefore, the complementary method is a better choice in calculating the attribution analysis of watershed hydrological cycle change.
(4) Human activities have changed the land use and vegetation cover of river basin, changed the types of land use and destroyed the ecological balance, thus affecting the hydrological process and altering canopy interception, infiltration and evapotranspiration, so as to change the runoff. In addition, the construction of water and soil conservation projects, especially the implementation of the policy of returning cropland to forest and grassland after 1999, has achieved remarkable ecological benefits, which is also an important reason for the reduction of runoff.