Attribution discernment of climate change and human interventions to runoff decline in Huangshui River Basin, China

Water resources are the key factors limiting socio-economic development, ecological protection, and ecological civilization construction in arid and semi-arid areas. High-intensity human water extraction activities, basin landscape pattern changes, and climate change (CC) have seriously disturbed the intensity, pathways, fluxes, and dynamics of the water cycle. This has resulted in increased uncertainty in the variability of water cycle elements and changes in the spatial and temporal patterns and processes of the terrestrial hydrological situation (Huntington 2006). CC and human interventions (HI) are the two major factors affecting the water cycle (Zhai & Tao 2017). Therefore, identifying the mechanisms of runoff effects due to CC and HI is not only the scientific basis for the adaptive management of water resources under environmental changes but also an intrinsic need for eco-hydrological regulation of sustainable development in watershed systems.

Runoff changes and attribution in watersheds under changing environments have become a research hot topic for Chinese and foreign scholars. This is an important topic emphasized by the Panta Rhei of the International Hydrological Decade (2013–2022) (Montanari et al. 2013). Runoff in the northern rivers of China is on a declining trend and suffers from severe droughts and floods, especially in the Yellow River Basin and the Hai River Basin (Feng et al. 2016). Whether the main cause of runoff changes is CC or HI, and how to effectively identify and quantitatively separate the impacts of CC and HI on runoff changes, were the main scientific questions of our study. Attribution analysis of runoff changes methods are as follows: (1) statistical regression analysis (Ma et al. 2008; Xu 2011; Kong et al. 2016); (2) climate elasticity method, total differential method, and complementary method based on the Budyko hypothesis (Xu et al. 2014; Zhou et al. 2016; Wu et al. 2017); (3) hydrological modeling method (Wu et al. 2019; Zhou et al. 2022). Among them, the Budyko framework is extensively adopted because of its low requirement for historical data, its simple calculation, and the availability of basic physical mechanisms, and the method has been proven to be reasonably credible (Wu et al. 2017). The functional analytic formulations of the Budyko hydrothermal coupling framework usually utilize the Choudhury-Yang (C-Y) function (Yang et al. 2008) and the Fu function (Fu 1981). The Budyko hypothesis has confirmed good applicability in the Yellow River Basin, Qinghai-Tibet Plateau region, and the Loess Plateau region in China (Sun et al. 2007; Yang et al. 2015). However, most of the climate elasticity methods can only separate and quantify the contributions of CC and HI to runoff changes at the annual scale, and their application is limited at the smaller spatio-temporal scale. In most studies, the Budyko equation parameters are constant during the baseline period (BP) and perturbation period (PP), and both CC and HI potentially affect equation parameters directly or indirectly (Wang et al. 2018). To address the evolution of watershed landscape parameters over time and their discontinuities, the time-varying Budyko framework has been developed and applied (Jiang et al. 2015; Xu et al. 2021).

The hydrological modeling field has evolved from lumped to distributed models. Lumped models, which treat the watershed as a singular entity, overlook the spatial variability of hydrological phenomena and simplify watershed variables using average values. These models are typically employed for rainfall–runoff simulations. In contrast, distributed models are characterized by robust physical underpinnings. They not only account for geographic heterogeneity but also elucidate the mechanisms of hydrological processes within a watershed. Consequently, distributed models offer superior capabilities over lumped models in quantifying the impacts of CC and HI on runoff variations across diverse spatial and temporal scales. However, their effectiveness is sometimes constrained by challenges related to model spatio-temporal consistency and inherent uncertainties. Soil & Water Assessment Tool (SWAT) is a process-driven and semi-distributed hydrological model with strong physical mechanisms that was developed by the United States Department of Agriculture (USDA) (Arnold et al. 2012). It can better simulate watershed-scale hydrological processes and is widely used for hydrological process simulation and runoff change attribution analysis (Douglas-Mankin et al. 2010; Grusson et al. 2015; Senbeta & Romanowicz 2021). Furthermore, Water and Energy transfer Processes in Large River basin (WEP-L) (Xu et al. 2022), Variable Infiltration Capacity Macro-scale Hydrological mode (VIC) (Yong et al. 2013), and Distributed Time Variant Gain Model (DTVGM) (Ning et al. 2016) have also been applied to analyze the impacts of CC, subsurface change, and socio-economic water withdrawal on runoff changes. Among the aforementioned hydrological models, the SWAT stands out for its versatility in simulating scenarios pertaining to climate changes and anthropogenic impacts. Notably, SWAT remains the most widely utilized model for investigating runoff responses to climatic variations and human activities, as evidenced by numerous studies (Li et al. 2009; Lin et al. 2012; Lin et al. 2015; Wang et al. 2017). Its efficacy in modeling the impact of environmental changes on surface hydrology has been validated across diverse geographical and climatic contexts (Li et al. 2022). Consequently, this study employs SWAT to assess the contributions of CC and HI to the runoff alterations in the Huangshui River Basin (HRB).

Given that both the Budyko framework and hydrological models such as SWAT have their respective strengths and limitations, there are uncertainties in the attribution results of runoff changes in the same study area obtained by different scholars (Miao et al. 2021). Thus, cross-validation of multiple methods can help us to improve the reliability of hydrological process simulation and runoff change attribution. A single method such as statistical analysis or climate elasticity model has been employed in previous investigations by Chinese researchers (Zhang et al. 2014; Du et al. 2022) to evaluate the contributions of CC and HI to runoff changes in the HRB, but the comparison with the results of other methods is still lacking and the description of hydro-physical processes is weak. Moreover, the diagnosis of runoff breakpoints is inconsistent and the evaluation results of runoff change attribution are different due to noise interference of runoff time series and the influence of the study time scale. To enhance the credibility of the results, it is necessary to validate the runoff breakpoints using multiple methods and quantify runoff change attribution coupling Budyko framework and SWAT hydrological model.

The HRB, located in eastern Qinghai Province, China, serves as a crucial ecological barrier and is the political, economic, and cultural hub of the province. Confronted with the dual pressures of CC and HI, the basin's rivers predominantly exhibit a declining trend in annual runoff. This trend exacerbates the already critical contradiction between water supply and demand in the HRB, positioning water as the primary limiting factor for sustainable development. Moreover, the implementation of the Lanzhou-Xining Urban Agglomeration Development Planning and Ecological Protection and High-quality Development in the Yellow River Basin amplifies the necessity for integrated water resources management (IWRM). Accurately elucidating the runoff response and its driving forces in the changing environmental context, and discerning the attribution of runoff changes, hold substantial scientific importance. These insights are crucial for adaptive water resources management and the strategic planning of water network projects in the HRB. Therefore, the objectives of our study were as follows: (1) examining runoff change trend and diagnosing runoff breakpoints with multiple methods in HRB; (2) quantifying climate elasticity and human activity elasticity of runoff change in HRB using the Budyko framework; (3) differentiating the impacts of CC and HI on runoff changes using the climate elasticity method, complementary method, and SWAT hydrological model; (4) identifying the major factors and drivers of anthropogenic disturbances affecting runoff changes.

The study area is the HRB (100°42′–103°01′E, 36°02′–37°28′N), a major tributary of the upper reaches of the Yellow River, covering a total area of about 16,383 km². The HRB originates in Haiyan County of Qinghai Province and is located in the transition zone between the Tibetan and Loess Plateaus, an area of arid and semi-arid continental climate. The average annual precipitation is 300–500 mm, and the average annual temperature is 0.6–7.9 °C. The HRB has a population of about 3.39 million (2018) and constitutes about 60% of the total economy of the province and about 56% of the arable land of the province. HRB's average annual runoff is 2.16 billion m³. Per capita water resources are about 1/3 of that of the national water resources, and per capita water consumption is 2/3 of that of the national average water resources. The utilization rate is close to 40%, and part of the utilization rate of branch ditch surface water is greater than 50%. The resource and engineering water shortage problems are prominent in HRB. Minhe hydrological station (MHS) is the control station for the HRB.

The daily precipitation, temperature, relative humidity, wind speed, and sunshine hours were derived from the China National Meteorological Information Center (https://data.cma.cn). Potential evaporation was calculated using the Penman–Monteith equation (Allen et al. 1998). Missing meteorological data were acquired by linear interpolation using data from neighboring stations. Batch extraction and processing of daily meteorological data were implemented in Python.

Observed runoff data of the hydrological stations were obtained from the Qinghai Hydrology and Water Resources Forecasting Center. The multi-year average rainfall at representative stations in HRB was derived from the Qinghai Provincial Hydrology Handbook. The average precipitation and potential evapotranspiration of the HRB were calculated by the Tyson polygon method (Bayraktar et al. 2005; Wang et al. 2008). The annual NDVI dataset was acquired from the Data Center for Resources and Environmental Sciences of the Chinese Academy of Sciences (https://www.resdc.cn). Vector data such as river system, government residence, and administrative region boundaries were derived from the National Geomatics Center of China (https://www.ngcc.cn/ngcc). Socio-economic data were extracted from Qinghai Statistical Yearbook and Xining National Economic and Social Development Statistical Bulletin. Water resources development and utilization data were obtained from the Qinghai Water Resources Bulletin and Water Network Planning of the HRB.

For the analysis of long-term trends in the annual runoff time series spanning from 1959 to 2014, several statistical methods were utilized. These included linear regression, the non-parametric Mann–Kendall test (Woo & Shin 2022), the Spearman test, and the sliding average method (Lou et al. 2019). To identify breakpoints in the runoff data, we employed the non-parametric Pettitt mutation test, the sliding t-test, and the cumulative anomaly curve method. The entire runoff sequence was segmented into two distinct periods: the BP and the PP, demarcated by the identified abrupt change point.

SWAT is a multi-process model that integrates hydrology, ecology, agriculture, and water quality. We utilized the SWAT model to discretize the HRB into 85 sub-basins (SUBs) and 1,480 hydrological response units (HRUs). A preheat period of 3 years was set to simulate the monthly runoff of the study area from 1962 to 2014. The potential evapotranspiration (expressed with ET₀) was calculated through the FAO-56 Penman–Monteith equation in the SWAT model. Surface runoff was computed using Soil Conservation Service Curve Number (SCS-CN) and the intra-channel confluence process was calculated using the Muskingum principle. Twenty-eight parameters were selected to calibrate and validate the SWAT model using the Sequential Uncertainty Fitting Ver-2 (SUFI-2) algorithm, which was built and introduced into the SWAT-CUP software. The model was set up for an initial run of 500 times. The larger the t-Stat determined by the SUFI-2 algorithm, the more sensitive the parameter is. Besides, the smaller the p-value, the more important the parameter is. We used the parameters with p < 0.1 as key sensitive parameters and tuned the parameters using a combination of automatic and manual calibration until the Nash–Sutcliffe efficiency (NSE) no longer could be improved (Abbaspour et al. 2015). We assessed the SWAT model performance by using four indicators (Abou Rafee et al. 2019), including determination coefficient (R²), Nash–Sutcliffe efficiency (NSE), percentage bias (PBIAS), and root mean square error (RMSE) to observations of standard deviation ratio (RSR) (Ren et al. 2021).

The results of the analysis of runoff elasticity changes of MHS in HRB are presented in Table 2. The potential evapotranspiration (ET₀) and aridity index (φ) increased by 1.93 and 4.24%, respectively, and the precipitation (P) and runoff (R) decreased by 2.27 and 16.67%, respectively, during the PP in comparison to the BP. The watershed landscape parameters (n) calculated by the Fu and C-Y functions all showed an increase in PP compared with BP, indicating that the influence of human activities tended to be enhanced. Based on the multi-year average hydro-meteorological data from 1959 to 2014, the elasticity coefficients calculated using the C-Y equation were 2.19 for P, −1.19 for ET₀, and −1.52 for n. The elasticity coefficients demonstrated that a 1% increase in P, ET₀, and n gave rise to a 2.19% increase, a 1.19% decrease, and a 1.52% decrease in the runoff, respectively. Obviously, the sensitivity of runoff to P was higher than ET₀ and basin landscape changes in HRB.

Table 2

Analysis of hydro-meteorological factors and runoff elasticity changes of MHS in HRB

Parameters .		Baseline period .	Perturbation period .	Simulation period .	Change amount .
		1959–1990 .	1991–2014 .	1959–2014 .
Hydrothermal parameters	/mm	511.9	500.3	504.5	−11.6
	/mm	846.1	862.4	854.2	16.2
	/mm	110.4	92.0	102.5	−18.4
		1.65	1.72	1.69	0.07
Fu function	ω	2.27	2.38	2.30	0.12
		2.10	2.23	2.14	0.12
		−1.10	−1.23	−1.14	−0.12
		−2.12	−2.26	−2.19	−0.15
C-Y function	n	1.55	1.67	1.59	0.12
		2.14	2.28	2.19	0.13
		−1.14	−1.28	−1.19	−0.13
		−1.46	−1.60	−1.52	−0.14

Parameters .		Baseline period .	Perturbation period .	Simulation period .	Change amount .
		1959–1990 .	1991–2014 .	1959–2014 .
Hydrothermal parameters	/mm	511.9	500.3	504.5	−11.6
	/mm	846.1	862.4	854.2	16.2
	/mm	110.4	92.0	102.5	−18.4
		1.65	1.72	1.69	0.07
Fu function	ω	2.27	2.38	2.30	0.12
		2.10	2.23	2.14	0.12
		−1.10	−1.23	−1.14	−0.12
		−2.12	−2.26	−2.19	−0.15
C-Y function	n	1.55	1.67	1.59	0.12
		2.14	2.28	2.19	0.13
		−1.14	−1.28	−1.19	−0.13
		−1.46	−1.60	−1.52	−0.14

View Large

Table 3 presents the results of the attribution analysis of runoff changes under the Budyko framework. The calculations using the Fu equation suggested that the contribution of CC to runoff reduction was 39.85% (P contributed 27.51% and ET₀ contributed 12.34%) and the contribution of HI to runoff attenuation was 60.15%. Using the C-Y equation, it was found that the contribution of CC to the decrease of runoff was 38.96% (P contributed 27.04% and ET₀ contributed 11.92%) and the contribution of HI to the reduction of runoff was 61.04%. The complementary method usually used α = 0.5 as the best estimate for the weighting factor, which calculated the contribution of CC to runoff decrease to be 39.92% (P contributed 27.47% and ET₀ contributed 12.45%) and the contribution of HI to runoff reduction was 60.08%.

Table 3

Results of attribution analysis of runoff changes using the Budyko framework

Budyko function and method .		Runoff changes (mm) .				Contribution rate of runoff change (%) .
		.	.	.	.	.	.	.	.
Fu function		−5.15	−2.31	−11.27	−18.73	27.51	12.34	39.85	60.15
C-Y function		−5.05	−2.23	−11.40	−18.67	27.04	11.92	38.96	61.04
Complementary method		−5.10	−2.31	−11.16	−18.58	27.47	12.45	39.92	60.08
		−4.85	−2.21	−11.52	−18.58	26.10	11.89	37.98	62.02
		−5.36	−2.42	−10.80	−18.58	28.84	13.02	41.86	58.14

Budyko function and method .		Runoff changes (mm) .				Contribution rate of runoff change (%) .
		.	.	.	.	.	.	.	.
Fu function		−5.15	−2.31	−11.27	−18.73	27.51	12.34	39.85	60.15
C-Y function		−5.05	−2.23	−11.40	−18.67	27.04	11.92	38.96	61.04
Complementary method		−5.10	−2.31	−11.16	−18.58	27.47	12.45	39.92	60.08
		−4.85	−2.21	−11.52	−18.58	26.10	11.89	37.98	62.02
		−5.36	−2.42	−10.80	−18.58	28.84	13.02	41.86	58.14

Note: The watershed parameter (n) was estimated using the C-Y equation in the complementary method.

View Large

The estimates of runoff changes attributed to the above three methods were similar, and HI played a primary role in runoff changes. The complementary method can determine the upper and lower limits of CC and HI impacts on runoff using weighting factors α = 0 and α = 1. The thresholds for CC contribution to runoff reduction were 37.98–41.86% and for HI contribution to runoff decrease were 58.14–62.02%. The attribution results of runoff changes identified by the Fu equation and C-Y equation were within the thresholds estimated by the complementary method. The runoff reductions calculated by the Fu equation, C-Y equation, and complementary method were 18.73, 18.67, and 18.58 mm, respectively, and the observed runoff reduction was 18.4 mm. The relative errors between the estimated and measured values were 1.79, 1.45, and 0.98%, respectively, reflecting that the calculation accuracy of the Budyko hypothesis was high. Hence, the results of the attribution of runoff changes based on the Budyko framework were reliable.

The attribution results of runoff changes calculated by the SWAT model in HRB are presented in Table 5. The contribution of CC to runoff changes was 38.72% and that of HI was 61.28%. Hence, the impact of HI on runoff changes was more obvious than that of CC.

Table 5

Runoff changes and attribution discernment based on SWAT simulation in HRB

Average flow (m3/s) .				Flow changes (m3/s) .			Contribution rate (%) .
.	.	.	.	.	.	.	.	.
53.48	51.64	44.13	48.02	−9.35	−3.62	−5.73	38.72	61.28

Average flow (m3/s) .				Flow changes (m3/s) .			Contribution rate (%) .
.	.	.	.	.	.	.	.	.
53.48	51.64	44.13	48.02	−9.35	−3.62	−5.73	38.72	61.28

View Large

We determined the impact of CC and HI on runoff using the runoff elasticity method based on the Budyko hypothesis, complementary method, and SWAT hydrological model, respectively. The Budyko hypothesis determined the contribution of CC to runoff changes to be 37.98–41.86% and that of HI to be 58.14–62.02%. And SWAT hydrological model estimated the contribution of CC and HI to runoff changes to be 38.72 and 61.28%, respectively. Although the Budyko framework and SWAT model were substantially different in terms of computational principles and time scales, attribution results of runoff changes were close. Zhang et al. (2014) utilized the change rate of the accumulation slope method to determine the contribution of CC to runoff reduction in the 1966–2010 periods and found that it was 35.46% in HRB, and that of HI was 64.54%. Du et al. (2022) investigated the attribution of runoff changes on different time scales based on the Fu function and concluded that landscape had the most significant impact on runoff, followed by precipitation, and potential evapotranspiration had the least impact. The results of existing studies were generally consistent with our conclusion that human activities were the primary factor that contributed to the decrease in the runoff. It was worth noting that the watershed landscape parameter (n) in the Budyko framework was influenced by both CC and HI, and that CC could indirectly affect runoff changes by influencing the watershed landscape, which may result in a small contribution of CC to runoff changes under the Budyko framework (Miao et al. 2021).

The effects of CC and HI on runoff are different at different spatial and temporal scales, the interactions between the two are complex, and it is difficult to separate them completely. Uncertainty in the results of the attribution of runoff changes in this study arose principally from the following aspects. (1) The division of the BP and PP was determined by the detection accuracy of runoff time series breakpoint, and the spatio-temporal harmony and spatio-temporal scale effects of the hydrological model had significant effects on the results. (2) The physical meaning and mathematical expression of the parameter n(ω) in the Budyko equation is still unclear, and the parameters affected by the water conservancy project and watershed landscape pattern dynamic time-varying are not sufficiently represented. (3) The role of CC and HI is interactive, and the two factors are not completely independent. (4) Geographic divergence and temporal variability of hydro-meteorological parameters and watershed landscape patterns in the SWAT hydrological model give rise to differences in eco-hydrological processes at different spatial and temporal scales. The primary sources of uncertainty in the SWAT model are associated with input data, model parameters, and the structure of calibration tools. These factors collectively contribute to unavoidable deviations in model simulations. More importantly, the uncertainty in contribution to runoff changes due to the ‘equifinality for different parameters’ of the hydrological model cannot be ignored.

We used the Budyko hypothesis and SWAT hydrological model to quantitatively differentiate the attribution of CC and HI to runoff changes, which is important for adaptive water resources management in HRB. The main findings of the study are as follows:

• (1)

  The runoff in HRB on the time scale of 1959–2014 showed a decreasing trend, and the average decrease rate of runoff was 0.212 m³/s·a⁻¹ (0.44 mm/a). The significant runoff breakpoint was 1990, and the study period was divided into BP (1959–1990) and PP (1991–2014). The average runoff in PP had decreased by 18.4 mm compared with the BP.

• (2)

  The Budyko framework determined the contribution of CC to runoff changes in HRB to be 37.98–41.86% and that of HI to be 58.14–62.02%. The SWAT hydrological model estimated the contribution of CC and HI to runoff decline to be 38.72 and 61.28%, respectively. The study concluded that HI was the primary factor in runoff decline, and direct human activities such as economic and social water withdrawals and water conservancy project construction had a substantial impact on runoff changes.

• (3)

  The main reasons for runoff decline driven by HI were: (a) the increase in water consumption for domestic, industrial, and agricultural use to support economic and social development; (b) increased watershed storage and disruption of natural runoff processes owing to reservoir construction; (c) increase in green water flux and green water storage because of reforestation. However, how to effectively and quantitatively distinguish the effects of different anthropogenic disturbances (e.g., Land use/cover change (LUCC), dam operation, and water extraction) on hydrological processes requires further investigation and research.

The authors would like to thank the reviewers and editors for their insightful comments on improving the paper. Besides, the authors would like to express their gratitude to EditSprings (https://www.editsprings.cn) for the expert linguistic services provided. This work was supported by the Qinghai Minzu University Planning Project [grant number: 23GH14] and the Applied Basic Research Project of Qinghai Province in China [grant number: 2024-ZJ-773].

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.