A probabilistic conceptual model to attribute runoff variations to human activity

ABSTRACT

A new probabilistic conceptual model, named the Assessment Model of Human Impact on Runoff Based on Stationarity Hypothesis (AMHIRBSH), was developed to attribute runoff variations to human activities through evaluating the changes in the runoff–runoff relationship (i.e. that between runoff of different watersheds). The AMHIRBSH was then applied to the Baiyangdian drainage basin, North China. After applying the Mann-Kendall test for trend and the Pettitt test for breakpoint, the runoff over the Baiyangdian drainage basin was shown to have had a breakpoint in approximately 1979, which split the time series into a baseline period (pre-1979) and an altered period (post-1979). The runoff was greatly reduced between these two periods. The runoff–runoff relationship of the Wangkuai and Xidayang sub-watersheds was evaluated using a copula function which showed it to have a non-stationarity degree of 35.64%, likely induced by massive afforestation programmes. This method provides an innovative perspective for attributing runoff variations to human activities.

KEYWORDS:

• probabilistic conceptual model
• stationarity hypothesis
• human activity
• runoff variations
• copula function

Editor S. Archfield Associate editor K. Ryberg

1 Introduction

Both climate change and human activities can greatly alter hydrological processes and can subsequently have critical impacts on runoff variations (Tu et al. 2015, Yuan et al. 2016, Zhang et al. 2016, Paul et al. 2018). As human activities intensify, understanding the human contributions to runoff variation is important for developing adaptive water management (Milly et al. 2008, Wu et al. 2017, Liu et al. 2019, Paul et al. 2019). However, this has become a challenge due to the complicated interactions among human activity, climate and hydrology (Wang and Hejazi 2011, Jiang et al. 2015). Human activities can not only alter runoff directly by affecting the hydrological processes but also indirectly by disturbing the climate variables (Wang and Hejazi 2011). In this study, we use the term “human activities” to reflect direct human activities, including water withdrawal, land use and hydrologic constructions.

In recent years, numerous research studies have been carried out to separate the impacts of human activities from those of climate change on runoff variations (Wei and Zhang 2010, Velpuri and Senay 2013, Dey and Mishra 2017, Suttles et al. 2018, Zipper et al. 2018). Three categories of methods have generally been utilized worldwide: empirical statistics, elasticity-based methods and process-based methods. Empirical statistics, such as regression analyses, time-trend analyses and the double-mass curve method, separate human contributions to runoff variations based on the change in statistical relationships between runoff and relevant climatic drivers (e.g. precipitation) at different periods (Wang 2014, Zhao et al. 2014, 2015, Kong et al. 2016, Zhang et al. 2016). Elasticity-based methods, encompassing the nonparametric method and the Budyko framework, assess the runoff variations attributed to climate change and human activities using the elasticities of runoff to changes in precipitation and potential evapotranspiration (Ye et al. 2013, Mwangi et al. 2016, Wu et al. 2017, Liu et al. 2019). Process-based methods use physical hydrological models, such as the simple lumped conceptual daily rainfall–runoff model (SIMHYD) and the Soil & Water Assessment Tool (SWAT), to simulate runoff variations with changed land use or water utilization characteristics (Ahn and Merwade 2014, Xu et al. 2014, Zeng et al. 2015).

There are two main logical ways to estimate human impacts on runoff variations: (1) modelling the change in runoff under varied human activities using physical hydrological models, which is direct but requires time-consuming calibration and validation involving topography, vegetation, climate and hydrology (Wei and Zhang 2010, Bao et al. 2012, Wang et al. 2013, Xu et al. 2014); and (2) estimating human contributions by removing predicted climatic contributions to runoff, which are usually calculated using the climate–runoff relationship (Milly et al. 2008, Zhang et al. 2011). However, the watershed climate–runoff relationship is very complex, as it is non-linear, time-varying and uncertain, especially for extreme climatic variations such as high-intensity storms (Wang et al. 2014), which induces great challenges for estimating climatic and human contributions to runoff variations. Besides runoff data, therefore, a long time series of climatic data such as precipitation and temperature is also mandatory. However, these kinds of climatic data are not always available, especially for mountainous regions without gauging stations.

To overcome the limits discussed above, we put two adjacent sub-watersheds with similar natural backgrounds into one research system, so that the climate factors could be excluded. Then, we develope a new probabilistic conceptual model named the Assessment Model of Human Impact on Runoff Based on Stationarity Hypothesis (AMHIRBSH) to assess the contribution of human activities to runoff variations. This was accomplished by first identifying the baseline period and the altered period using the sequential Mann-Kendall (M-K) test and the Pettitt test, and then measuring human impacts on runoff by analysing the non-stationarity of the runoff–runoff relationship between the two periods, without climate data. This method can provide an innovative perspective on and approach to assessing the contribution of human activities to runoff variations without including the complex climate–runoff relationship. This method is especially useful in cases where climate data are scarce, such as in mountainous regions without gauging stations.

2 Methodology

2.1 Probabilistic conceptual model development

Stationarity – the idea that natural systems fluctuate within an unchanging envelope of variability – is a foundational concept that permeates training and practice in water-resource engineering (Milly et al. 2008). It implies that any single variable has a time-invariant marginal probability density function (marginal PDF), and any two related variables have a time-invariant joint probability density function (joint PDF) which depends on both the marginal PDF of the two variables and their specific relationship. Under stationarity, the relationship between two related variables is therefore also inferred to be time invariant, such as the climate-runoff relationship and runoff-sediment relationship for one basin as a system. Then, for two adjacent basins (e.g. basin S₁ and basin S₂), as a system themselves are supposed to have similar hydrological process changes (e.g. the trend, the approximate occurrence times of mutation points) with each other for similar climate changes. Thus, an assumption was made that the runoff–runoff relation of basins S₁ and S₂ is time invariant in natural systems, which we call the “runoff–runoff stationarity” assumption. As it is generally recognized that climate–runoff stationarity could be changed due to anthropogenic effects (Milly et al. 2008), we can separate anthropogenic effects from climatic effect on runoff variations by estimating the change in the climate–runoff relationship. Likewise, the question arises whether or not the runoff–runoff stationarity relationship will also be compromised into non-stationarity by human disturbances. The answer may provide a new perspective for assessing the contribution of human activities on runoff variations.

The previous study by Wang et al. (2015) was utilized to verify whether the runoff–runoff stationarity would be compromised into non-stationarity by human disturbances. Wang et al. (2015) adopted the double mass curve (DMC) method to identify the driving force of streamflow variations for the Chaohe and Baihe river catchments from 1963 to 2011. As the Chaohe and Baihe rivers are twin rivers, they have undergone similar climate conditions and similar streamflow regimes. With the segmented fitting method, the study period was separated into four periods based on different regressive slopes of DMCs over the entire study duration. The periods were defined as follows: Period 1: 1963–1973, the baseline period; Period 2: 1974–1982, the low-impact period fairly close to the baseline; Period 3: 1983–1997, the altered period proven to be dramatically affected by human activities; and Period 4: 1998–2011, the altered period proven to be affected by cumulative human activities.

We evaluated whether the stationarity of runoff–runoff relationships had been changed under human effects among the four periods in the observed records, based on the statistical ellipse model providing a simple and effectivity summary of a bivariate relation (Friendly et al. 2013, Polajzer et al. 2016). As shown in Fig. 1, with 80% confidence ellipses, the coverage area of Period 2 is fairly close to that of Period 1, indicating the presence of stationarity. In contrast, there were only small elliptical areas coincident with Period 1 for Period 3 as well as for Period 4, which indicates that their probability distributions had been changed greatly from the baseline period due to human activities, indicating non-stationarity. Hence, the above analysis clearly showed that stationarity would be compromised by human disturbances within areas such as the Chaobai River basin.

Figure 1. Joint distribution plot of annual runoff in the Baihe and Chaohe rivers in four periods between 1963 and 2011, proving that runoff–runoff stationarity will be compromised into non-stationarity if affected by human disturbances. Period 1: 1963–1973, the baseline period; Period 2: 1974–1982, the low-impact period fairly close to the baseline; Period 3: 1983–1997, the altered period proven to be dramatically affected by human activities; and Period 4: 1998–2011, the altered period proven to be affected by cumulative human activities (Wang et al., 2015)

Following from the above conclusion, a probabilistic conceptual model based on the stationarity hypothesis (i.e. AMHIRBSH) can be established to separate the effects of human activity from those of climate change on runoff variations. As shown in Fig. 2(a), the joint behaviour of R₁ (i.e. the runoff in basin S₁) and R₂ (i.e. runoff in basin S₂) is assumed to be in line with the joint distribution D₁ within the baseline period, and it would probably alter from D₁ to D₂ if affected by human activities in the altered period. While R₁ was given a certain value r₁, the value of R₂ (i.e. r₂) fluctuated within a confidence interval from the minimum value e₂ to maximum value e₁ when subjected to deterministic stationary relationships and stochastic variations. This implies that R₂ has a specific conditional probability density function (conditional PDF) when given R₁, such as in Fig. 2(b). Obviously, the conditional density of R₂ given R₁ (i.e. $f_{R_2|R_1}(r_2|r_1)$) for any sample falling within the confidence interval is greater than that falling outside the confidence interval. In the baseline period, most runoff pairs obeying D₁ distribution would fall within the confidence interval. In the altered period, however, due to anthropogenic effects some runoff pairs would probably fall outside the confidence interval. Obviously, the mean of $f_{R_2|R_1}(r_2|r_1)$ for samples from baseline period should be greater than that from altered period, providing a reasonable way to separate the impacts of human activities from those of climate change on runoff variations by analysing the difference of $f_{R_2|R_1}(r_2|r_1)$ between two periods. The degree of runoff–runoff variations can be expressed by the Human Contribution (HC) index:

(1) $\mathrm{H}\mathrm{C}=1-\frac{\overline{f_{R_2|R_1}\left(r_{2B}|r_{1B}\right)}}{\overline{f_{R_2|R_1}\left(r_{2A}|r_{1A}\right)}}$(1)

Figure 2. The probabilistic conceptual model, Assessment Model of Human Impact on Runoff Based on Stationarity Hypothesis (AMHIRBSH), to separate the impacts of human activities from those of climate change on runoff variations: (a) example of the runoff–runoff relationship (expressed by the joint distribution) having changed from D₁ to D₂, affected by human activities; and (b) example of the conditional probability density function of R₂ given R₁ = r_1.

where $\overline{f_{R_2|R_1}(r_{2B}|r_{1B})}$ is the arithmetic mean of $f_{R_2|R_1}(r_{2B}|r_{1B})$ for samples from the baseline period, while $\overline{f_{R_2|R_1}(r_{2A}|r_{1A})}$ is the arithmetic mean of $f_{R_2|R_1}(r_{2A}|r_{1A})$ for samples from the altered period. Compared to $\overline{f_{R_2|R_1}(r_{2B}|r_{1B})}$, the $\overline{f_{R_2|R_1}(r_{2B}|r_{1B})}$ is smaller and the HC is greater, indicating a stronger human contribution to runoff variations, and vice versa.

2.2 Framework and statistical methods for attributing runoff variations to human activities

Based on the AMHIRBSH, a framework consisting of three steps and corresponding statistical methods is shown in Fig. 3, to attribute runoff variations to human activities.

Figure 3. Framework for attributing human activities to runoff variations by using the probabilistic conceptual model

2.2.1 Identifying the baseline period and altered period

The M-K test (Yue et al. 2002) (see Appendix) was used to investigate the trend (i.e. increase or decrease) and degree of runoff variations. The Pettitt test (Pettitt 1979) (see Appendix) was then used to detect the breakpoint of the time series of runoff variations. Based on the breakpoint detection, the time series was split into two periods: the baseline period (pre-breakpoint) and the altered period (post-breakpoint). The baseline period was assumed to have negligible human activity while the altered period was assumed to have measurable effects from climate change and human activity on runoff variations.

2.2.2 Fitting natural joint distributions of runoff pairs in the baseline period

Copula functions (Nelsen 2006) were used to fit the joint distributions of runoff pairs (i.e. R₁ and R₂) in the baseline period. Copula functions provide a clear view of joint behaviours of multi-variables and allow for intuitive relationship analysis between variables. Over the past several years, copulas have been widely used in hydrology and published in many papers with detailed procedures (Wahl et al. 2015, Wang et al. 2017, 2018).

For two random variables (X and Y) with marginal cumulative distribution functions F(x) and F(y), Sklar’s theorem (Nelsen 2006) states that the joint cumulative distribution function (CDF) can be denoted by a copula function C as

(2) $F_{X,Y}\left(x,y\right)=C\left(F_X\left(x\right),F_Y\left(y\right)\right)$(2)

Thus, the dependence structure between X and Y is well established. Assuming the marginal distributions are continuous with PDFs f(x) and f(y), the bivariate PDF related to the CDF is described as

(3) $f_{X,Y}\left(x,y\right)=f_X\left(x\right)f_Y\left(y\right)c\left(F_X\left(x\right),F_Y\left(y\right)\right)$(3)

where

(4) $c\left(u_1,u_2\right)=\frac{{\mathrm{\partial }}^{2}C\left(\left(u_1,u_2\right)\right)}{\mathrm{\partial }{u}_1\mathrm{\partial }{u}_2}$(4)

The first step in fitting a copula is to construct the marginal distributions (i.e. U and V) by fitting the empirical distributions with the theoretical distributions, such as lognormal, normal, gamma and Weibull distributions. The fitted parameters of marginal distributions were estimated by the maximum likelihood method, and their goodness-of-fit tests were implemented using the Kolmogorov-Smirnov (K-S) test (Marsaglia et al. 2003).

With the identification of marginal probability distributions, their joint distributions can be directly developed by linking U and V through specific copula functions (i.e.$F_{R_1,R_2}\left(r_{1B},r_{2B}\right)=C\left(U,V\right)$). Thus, the runoff–runoff relationships embracing a non-linear relation and a tail dependence between two sub-watersheds in the baseline period could be effectively analysed with the joint distributions.

2.2.3 Attributing runoff variations to human activities by measuring the degree of change of the conditional density from the baseline period to the altered period

Based on the joint distributions developed in Equation (2)(2) $F_{X,Y}\left(x,y\right)=C\left(F_X\left(x\right),F_Y\left(y\right)\right)$(2) , the $f_{R_2|R_1}(r_2|r_1)$ (i.e. the conditional density of R₂ given R₁) of historical samples from the baseline period and the altered period were judged using the conditional PDF as:

(5) $f_{R_2|R_1}\left(r_2|r_1\right)=\frac{f_{R_1,R_2}\left(r_1,r_2\right)}{f_{R_1}\left(r_1\right)}$(5)

where $f_{R_1,R_2}\left(r_1,r_2\right)$ is the joint density of R₁ and R₂, which was estimated from the copula-based joint distribution, and $f_{R_1}\left(r_1\right)$ is the marginal density for R₁, which was estimated from the marginal distribution of runoff in basin S₁.

For two observed runoff sequences $E_B=\left\{\left(r_{11},r_{21}\right),\dots ,\right.$ $\left.\left(r_{1t},r_{2t}\right)\right\}$ in the baseline period and $E_A=\left\{\left(r_{1t+1},r_{2t+1}\right),\dots ,\right.$ $\left.\left(r_{1n},r_{2n}\right)\right\}$ in the altered period, two groups of conditional densities would be obtained as $f_{R_2|R_1}(r_{2B}|r_{1B})$ and $f_{R_2|R_1}(r_{2A}|r_{1A})$. Then the HC index can be calculated following Equation (1)(1) $\mathrm{H}\mathrm{C}=1-\frac{\overline{f_{R_2|R_1}\left(r_{2B}|r_{1B}\right)}}{\overline{f_{R_2|R_1}\left(r_{2A}|r_{1A}\right)}}$(1) to analyse the human impacts on runoff variations.

3 Study area and data acquisition

3.1 Study area

Baiyangdian drainage basin is located at 113°39′–116°12′E and 38°3′–40°4′N (Fig. 4). The basin encompasses an area of 31 200 km², of which 53% occupies the western mountains and hills, with elevation of over 1000 m a.s.l., and 47% occupies the eastern alluvial plains, with elevation ranging from 10 to 100 m a.s.l. Within a semi-arid temperate continental monsoon climate, the basin receives an average annual precipitation of 548 mm, mostly in the summer (70–80%). There are nine tributaries into the Baiyangdian Lake, which is known as the largest wetland in North China, playing an important role in maintaining regional ecosystem functions such as retaining water and sustaining the biodiversity of the area.

Figure 4. Location of the Baiyangdian drainage basin, North China

In recent decades, however, the Baiyangdian drainage basin has undergone a substantial decrease in runoff. Since the 1950s, the annual average temperature in the Baiyangdian drainage basin has increased by approximately 0.4–1.2°C/10 years, while the annual average precipitation has declined by 22.91 mm/10 years (Gao et al. 2009). More than 150 reservoirs were built along the nine tributaries after 1958. Affected by such climate change and human activities, the upper streamflow and downstream water inflow into the Baiyangdian Lake have both dramatically decreased, resulting in a decreased water level and increased dried-up frequency of Baiyangdian Lake. Baiyangdian Lake dried up only once during the period 1919–1965, while it dried up frequently during 1965–2002 (Hu et al. 2012). The drought of Baiyangdian Lake was generally recognized as a cause of the continuous decline in ecological functions followed by various eco-environmental problems (Hu et al. 2012, Li et al. 2014). Therefore, identifying the role of human activities in runoff variations in Baiyangdian drainage basin is very important for water resource management and ecological restoration.

The rivers Tang and Zhulong are two of the main tributaries into Baiyangdian Lake. Xidayang sub-watershed (4420 km²) and Wangkuai sub-watershed (3770 km²), respectively located in the upper Tang River catchment and Zhulong River catchment, are critical areas for making this attribution. The two sub-watersheds are close to each other and thus have similar climate conditions. Xidayang sub-watershed has an annual average temperature of 12.1°C and an annual average precipitation of 530–660 mm/year, while Wangkuai sub-watershed has an annual average temperature of 12.5°C and an annual average precipitation of 536–745 mm/year.

3.2 Data acquisition

The total annual runoff of three hydrologic gauging stations located in the Upper Sha River, the Tang River and the Juma River, i.e. Fuping, Daomaguan and Zijingguan stations, respectively (shown in Fig. 4), and the annual inflow into Baiyangdian Lake were used to analyse the trends and breakpoints of runoff over the Baiyangdian drainage basin. Then the annual runoff in the Xidayang and Wangkuai sub-watersheds was used to attribute runoff variations to human activities in the Baiyangdian drainage basin. These data were collected or calculated from the Hydrological Yearbook of the Haihe River basin (Ministry of Water Resources of the P.R.China 2012).

The data processes of breakpoint detection, copula function establishment and conditional density calculation were carried out using R software, and the corresponding packages are described in the Appendix.

4 Results

4.1 Hydrological trends and breakpoints

In previous studies, the breakpoint of runoff in Baiyangdian drainage basin was found to occur in approximately 1980 (Liu et al. 2011, Hu et al. 2012). Here, the Pettitt test was utilized to identify breakpoints of runoff for Baiyangdian drainage basin and three typical stations (Fuping, Daomaguan and Zijingguan). The results showed that the four long-term runoff series all had breakpoints at a 5% significance level (see Table 1), and three of them occurred in approximately 1979, when abrupt change was most likely to occur, except for Daomaguan Station where it occurred in approximately 1982.

Table 1. Pettitt test results and mean of runoff series in the baseline and altered periods

Hydrologic station	Change point	p value	Long-term mean runoff (10^8 m^3)
			Baseline period	Altered period
Baiyangdian Lake	1979	0.0031**	15.46	2.73
Fuping Station	1979	0.0001**	3.83	1.88
Daomaguan Station	1982	0.0007**	3.07	1.71
Zijingguan Station	1979	0.0138*	2.94	1.76

* p < 0.05; ** p < 0.01.

As the Pettitt test has been found to be sensitive to monotonic trends and might determine a breakpoint in what is really a series with an up (or down) gradual trend (Busuioc and von Storch 1996, Ryberg et al. 2020), monotonic trends across the Baiyangdian watershed were also analysed with the M-K test. As shown in Fig. 5, the observed annual runoff levels of the three typical hydrologic gauging stations during 1960 to 2002 all display a significant (5% significance level) decreasing trend. In contrast, when the M-K test was used for two periods (i.e. 1960–1979 and 1980–2002), the annual runoff series all showed insignificant decreasing trends, in Baiyangdian drainage basin and the three stations. This indicates that the Pettitt test was not influenced by monotonic trends.

Figure 5. Trends (1960–2002) of annual runoff in typical hydrologic gauging stations in the Baiyangdian drainage basin. The grey band shows the 95% confidence interval of the LOESS smoothing (in R software) for runoff

Consequently, in this study, the runoff series over the basin was recognized as experiencing an abrupt change in 1979. Thus, the runoff series was divided into two periods, the baseline period (i.e. pre-1979) and the altered period (i.e. post-1979). Compared to the baseline period, the average annual runoff in the altered period decreased by 81.4%, 44.5%, 38.7% and 34.3% for Baiyangdian drainage basin, Fuping Station, Daomaguan Station and Zijingguan Station, respectively. Further analysis is necessary to determine the contribution of human activities to this reduction in runoff.

4.2 Joint distribution of runoff pairs in Xidayang and Wangkuai sub-watersheds

4.2.1 Marginal probability distribution of runoff pairs

The historical annual runoff datasets for Wangkuai and Xidayang sub-watersheds were denoted WR and XR, respectively, for convenience. The marginal distributions of WR and XR were selected from a pool of four distributions commonly used in hydrology: lognormal, normal, gamma and Weibull. The parameters estimated using the maximum-likelihood method and the corresponding goodness-of-fit test results are presented in Table 2.

Table 2. Goodness-of-fit tests for the marginal distributions of runoff pairs. The underlined values are the best distributions of runoff in the specified site. AIC: Akaike information criterion; D: Kolmogorov’s goodness-of-fit statistic

Sub-watershed	Distribution	Parameter		AIC	D*	p value
Wangkuai	Lognorm	1.53	0.93	142.55	0.10	0.94
	Norm	6.58	4.94	148.78	0.16	0.56
	Gamma	1.58	0.24	139.78	0.10	0.97
	Weibull	1.33	7.17	139.67	0.10	0.94
Xidayang	Lognorm	1.48	0.64	122.07	0.076	0.99
	Norm	5.36	3.42	131.12	0.157	0.54
	Gamma	2.71	0.50	122.34	0.08	0.99
	Weibull	1.68	6.04	123.71	0.10	0.95

*H₀: The true distribution function is equal to the hypothesized distribution function.

The Akaike information criterion (AIC; Honjo et al. 1994) was used to assess the goodness-of-fit of candidate distributions. Lower AIC values indicate better fit for a given set of data; given this, the best marginal distribution for WR is the Weibull distribution. The Kolmogorov-Smirnov (K-S) test was used for the goodness-of-fit test. The K-S statistic is 0.1, with a p value of 0.94, which indicates insufficient evidence to reject the null hypothesis at α = 0.1. Therefore, the marginal distribution of WR can be assumed to follow the Weibull distribution. Likewise, the marginal distribution of XR can be assumed to follow the lognormal distribution, due to the minimum AIC value and K-S test of not rejecting the null hypothesis with a p value of 0.99.

4.2.2 Joint probability distribution of runoff pairs

A strong dependence between WR and XR was first verified (with a Pearson correlation coefficient of 0.94, a Kendall correlation coefficient of 0.72 and a Spearman correlation coefficient of 0.89), which demonstrated the potential of using copula functions. Five well-known and commonly used copula functions in hydrology – two elliptical copula families (the Gaussian copula and Student’s copula) and three Archimedean copula families (the Clayton copula, the Frank copula and the Gumbel copula) (Wang et al. 2017, Jhong et al. 2018, Pandey et al. 2018) – were selected as candidates to fit the joint probability distributions of runoff pairs in the baseline period. The results of parameters fitted using the maximum-likelihood method and goodness-of-fit test following Kojadinovic et al. (2011) are shown in Table 3. Among the five copulas, the Gumbel copula (specializing in upper tail dependence) was recognized, due to its small AIC value, as the best model to effectively present the dependence structure of runoff pairs. The joint probability distribution is given as:

(6) $C\left(u,v\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left\{-{\left[{\left(-lnu\right)}^{4.25}+{\left(-lnv\right)}^{4.25}\right]}^{\frac{1}{4.25}}\right\}$(6)

Table 3. Goodness-of-fit tests for joint probability distributions of the runoff pairs. AIC: Akaike information criterion; Sn: von Cramér’s goodness-of-fit statistic

Copula	Parameter(s)		AIC	Sn*	p value
Gaussian	0.93		−40.42	0.02	0.78
Student’s t	0.90	1.11	−40.73	0.03	0.89
Clayton	5.07		−39.51	0.03	0.70
Gumbel	4.25		−42.55	0.02	0.97
Frank	12.55		−34.55	0.03	0.37

*H₀: The dependence structure of a multivariate distribution is well represented by the specific parametric family of copulas.

where u represents the marginal CDFs of the Weibull distribution with a shape parameter of 1.33 and a scale parameter of 7.17, while ν represents the marginal CDFs of the lognormal distribution with a mean of 1.48 and a standard deviation of 0.64.

The cumulative probability distributions are presented in Fig. 6(a) and Fig. 6(b). It can be seen that there was a probability of approximately 90% of runoff occurring in a range from 1 to 18 × 10⁸ m³ for Wangkuai sub-watershed and from 2 to 12 × 10⁸ m³ for Xidayang sub-watershed. The cumulative probability values of runoff pairs are dependent on one another. For example, the maximum value of the cumulative probability was only approximately 40% when the runoff in Wangkuai sub-watershed was 4.5 × 10⁸ m³ or the runoff in Xidayang sub-watershed was 3.8 × 10⁸ m³. According to Fig. 6(c) and Fig. 6(d), the probability density at co-low or co-high values of runoff pairs was relatively higher than that at low-high values, indicating a positive correlation of runoff in Wangkuai and Xidayang sub-watersheds.

Figure 6. Joint probability distributions for annual runoff in the Wangkuai sub-watershed and the Xidayang sub-watershed. (a) and (b) are the cumulative probability distributions; (c) and (d) are the probability densities of distributions

4.3 Attributing runoff variations to human activities

With the joint distribution function developed in Section 4.2, the conditional densities of XR given WR (i.e. $f_{R_2|R_1}(\mathrm{X}\mathrm{R}|\mathrm{W}\mathrm{R})$) in the baseline period and altered period were calculated. Figure 7 shows that the maximum and minimum values of $f_{R_2|R_1}(\mathrm{X}\mathrm{R}|\mathrm{W}\mathrm{R})$ in the altered period (i.e. $f_{R_2|R_1}(\mathrm{X}{\mathrm{R}}_{\mathrm{A}}|\mathrm{W}{\mathrm{R}}_{\mathrm{A}})$) were both less than those in the baseline period (i.e.$f_{R_2|R_1}(\mathrm{X}{\mathrm{R}}_{\mathrm{B}}|\mathrm{W}{\mathrm{R}}_{\mathrm{B}})$), and that $f_{R_2|R_1}(\mathrm{X}{\mathrm{R}}_{\mathrm{A}}|\mathrm{W}{\mathrm{R}}_{\mathrm{A}})$ occurring at a maximum frequency was close to 0.05, while $f_{R_2|R_1}(\mathrm{X}{\mathrm{R}}_{\mathrm{B}}|\mathrm{W}{\mathrm{R}}_{\mathrm{B}})$was close to 0.35. The one-tailed Student’s test (Ruxton and Neuhäuser 2010) was applied to test the null hypothesis that the $f_{R_2|R_1}(\mathrm{X}{\mathrm{R}}_{\mathrm{B}}|\mathrm{W}{\mathrm{R}}_{\mathrm{B}})$ had a mean less than $f_{R_2|R_1}(\mathrm{X}{\mathrm{R}}_{\mathrm{A}}|\mathrm{W}{\mathrm{R}}_{\mathrm{A}})$. The results show values of t = 1.7929 and p = 0.04015, which is strong evidence against the null hypothesis. In other words, with 95% confidence, the runoff variations were significantly affected by human activities.

Figure 7. Scatter plots and violin plots of the conditional densities of XR (annual runoff in Xidayang sub-watershed) given WR (annual runoff in Wangkuai sub-watershed) in the baseline and altered periods. The two layers are: (i) mapped points (scatter plot) of the conditional densities of XR given WR (i.e. $f_{R_2|R_1}(\mathrm{X}\mathrm{R}|\mathrm{W}\mathrm{R})$), and (ii) the compact display (violin plot) of the entire distribution of $f_{R_2|R_1}\left(\mathrm{X}\mathrm{R}|\mathrm{W}\mathrm{R}\right);$the median, interquartile range and lower/upper adjacent values are also given

To further measure the degree of human contribution to runoff variations, the HC index was calculated following Equation (1)(1) $\mathrm{H}\mathrm{C}=1-\frac{\overline{f_{R_2|R_1}\left(r_{2B}|r_{1B}\right)}}{\overline{f_{R_2|R_1}\left(r_{2A}|r_{1A}\right)}}$(1) . It equalled 0.3564 when the mean of $f_{R_2|R_1}(\mathrm{X}\mathrm{R}|\mathrm{W}\mathrm{R})$ changed from 0.2823 in the baseline period to 0.1817 in the altered period. This indicates a non-stationarity degree of 35.64% in the runoff–runoff relationship before and after 1979, suggesting that human activities contributed to runoff variations.

In addition, according to the probability theory (Fitelson and Hájek 2017), generally $f_{A|B}(a|b)\ne f_{B|A}(b|a)$, indicating that different results may be introduced with varied antecedents. Therefore, we additionally verified the case of $f_{R_2|R_1}(\mathrm{W}\mathrm{R}|\mathrm{X}\mathrm{R})$ to make the conclusion more robust. The results led to the same conclusion of human activities significantly contributing to runoff variations, with t = 1.987 and p = 0.02744. The non-stationarity degree of the runoff–runoff relationship was 45.85%.

5 Discussion

5.1 Model validation

A commonly used method, the double mass curve (DMC) method, was used to verify whether runoff variations in Baiyangdian watershed were significantly affected by human activities. The results of segmented linear fitting of cumulative precipitation and runoff (Fig. 8) show that the ratio of runoff to precipitation followed a notable downward trend from the baseline period to the altered period in both Wangkuai and Xidayang sub-watersheds. This indicates that the runoff was influenced by human activities, which is consistent with the conclusion obtained by the AMHIRBSH method.

Figure 8. Double mass curve (DMC) of annual precipitation and runoff for the period 1960–2002 in (a) Wangkuai sub-watershed and (b) Xidayang sub-watershed, showing the linear regression of cumulative precipitation and runoff in the baseline period (blue) and altered period (orange)

5.2 Land-use change and its impact on runoff decline

Identifying the impacts of specific human activity on runoff decline is important for water resource management. According to previous studies (Wang et al. 2009, Ma et al. 2010, Xu et al. 2014), land-use change is probably the main human activity that can be attributed to runoff decline in North China, where massive afforestation programmes were launched in the beginning of the 1980s. Fu (2010) identified the land-use changes in five sub-basins of Baiyangdian drainage basin, and their results showed a decrease in farmland and an increase in forest and grassland since 1980. Using Landsat Multispectral Scanner (MSS) data and Landsat Thematic Mapper (TM) data, we identified the temporal variations of the normalized difference vegetation index (NDVI) obtained every May within Xidayang and Wangkuai sub-watersheds, respectively (see Fig. 9). Since 1981, the NDVI in Xidayang sub-watershed showed a positive trend of 0.0057/year, and that in Wangkuai sub-watershed also showed a positive trend, of 0.0080/year, indicating the forest biomass had increased.

Figure 9. Variations of the normalized difference vegetation index (NDVI) in the Xidayang and Wangkuai sub-watersheds

To verify whether such afforestation had increased evapotranspiration and thus resulted in decreased runoff in the upper reservoir catchment, we applied the Budyko framework (Budyko 1974) to estimate the effect of increasing forest biomass on evapotranspiration, following Jaramillo et al. (2018). To eliminate infrastructure influence on runoff, hydrometeorological data from Daomaguan Station, located in the Upper Xidayang Reservoir catchment, and Fuping Station, located in the Upper Wangkuai Reservoir catchment, being far from irrigation areas and uncontrolled by any great reservoir, were used to identify possible forest-related effects on evapotranspiration. The meteorological data for precipitation and temperature were downloaded from the Climatic Research Unit Database (CRU TS 2.1) of the Consortium of Spatial Information of the Consultative Group on International Agricultural Research (CGIAR-CSI; Mitchell and Jones 2005). Based on previous studies (Budyko 1974, Zhang et al. 2001, Jaramillo et al. 2013, 2018), the actual evapotranspiration (AET) is estimated using a water balance method, as shown in Equation (7)(7) $\mathrm{A}\mathrm{E}\mathrm{T}=P-Q$(7) , assuming there is no annual water storage change, and the climate-driven AET (i.e. $\mathrm{A}\mathrm{E}{\mathrm{T}}_{\mathrm{c}}$) is estimated using Equation (8)(8) $\mathrm{A}\mathrm{E}{\mathrm{T}}_{\mathrm{c}}=P\sqrt{\left(1-\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{P}\mathrm{E}\mathrm{T}/P\right)\right)\times \frac{\mathrm{P}\mathrm{E}\mathrm{T}}{P}\times \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\frac{P}{\mathrm{P}\mathrm{E}\mathrm{T}}\right)}$(8) . As such, the residual effect on evaporative ratio between the baseline and altered periods, $\mathrm{\Delta }{\left(\mathrm{A}\mathrm{E}\mathrm{T}/P\right)}_{\mathrm{r}}$, is calculated using Equation (9)(9) $\mathrm{\Delta }{\left(\mathrm{A}\mathrm{E}\mathrm{T}/P\right)}_{\mathrm{r}}=\mathrm{\Delta }\left(\mathrm{A}\mathrm{E}\mathrm{T}/P\right)-\mathrm{\Delta }\left(\frac{\mathrm{A}\mathrm{E}{\mathrm{T}}_{\mathrm{c}}}{P}\right)$(9) .

(7) $\mathrm{A}\mathrm{E}\mathrm{T}=P-Q$(7)

(8) $\mathrm{A}\mathrm{E}{\mathrm{T}}_{\mathrm{c}}=P\sqrt{\left(1-\mathrm{e}\mathrm{x}\mathrm{p}\left(-\mathrm{P}\mathrm{E}\mathrm{T}/P\right)\right)\times \frac{\mathrm{P}\mathrm{E}\mathrm{T}}{P}\times \mathrm{t}\mathrm{a}\mathrm{n}\mathrm{h}\left(\frac{P}{\mathrm{P}\mathrm{E}\mathrm{T}}\right)}$(8)

(9) $\mathrm{\Delta }{\left(\mathrm{A}\mathrm{E}\mathrm{T}/P\right)}_{\mathrm{r}}=\mathrm{\Delta }\left(\mathrm{A}\mathrm{E}\mathrm{T}/P\right)-\mathrm{\Delta }\left(\frac{\mathrm{A}\mathrm{E}{\mathrm{T}}_{\mathrm{c}}}{P}\right)$(9)

where P is precipitation, Q is runoff, PET is potential evapotranspiration. A positive $\mathrm{\Delta }{\left(\mathrm{A}\mathrm{E}\mathrm{T}/P\right)}_{\mathrm{r}}$ can be an indicator of the effect of afforestation on the AET (Jaramillo et al. 2018). The results showed that the residual effect on the evaporative ratio (Equation 9(9) $\mathrm{\Delta }{\left(\mathrm{A}\mathrm{E}\mathrm{T}/P\right)}_{\mathrm{r}}=\mathrm{\Delta }\left(\mathrm{A}\mathrm{E}\mathrm{T}/P\right)-\mathrm{\Delta }\left(\frac{\mathrm{A}\mathrm{E}{\mathrm{T}}_{\mathrm{c}}}{P}\right)$(9) ) was 0.07664 and 0.02864 in Daomaguan and Fuping, respectively, indicating the increased vegetation as a result of afforestation programmes contributed to increased evapotranspiration and thus to decreased runoff.

5.3 Applications and implications of AMHIRBSH

The AMHIRBSH was used to correlate human activities with runoff variations through testing the non-stationarity of the runoff–runoff relationship. Although it cannot directly quantify the specific contribution of human activities, such as land use and hydrologic constructions, on runoff variations, without further analysis as described in Section 5.2, it can be used as an alternative tool to detect the degree of influence of human activities on runoff without considering climate change. This benefits those areas where climate data are difficult to obtain. In addition, the current climate data in some areas are mainly obtained by interpolation algorithms based on neighbouring meteorological data, which is possibly inconsistent with actual changes. In such a case, AMHIRBSH can be used as a supplement to other methods to ensure more accurate results.

More importantly, AMHIRBSH provides us with a new perspective for identifying the impact of human activities on runoff or other factors with respect to water resources or ecology, by analysing changes in the relationship among natural internal elements within the surface system. This is conducive to revealing the human-induced changes within the entire system (including variables and their relationships), rather than just a single variable, which may be more instructive for sustainable water resources management.

6 Conclusions

This study proposed a stationarity hypothesis in which the runoff–runoff relationship of two adjacent sub-watersheds was time invariant within a natural system. Statistical ellipse analysis of 49 years of runoff records observed in the Miyun Reservoir basin showed that a stationary system could be changed to non-stationary if affected by human activities. Based on this, a probabilistic conceptual model (AMHIRBSH) was established to correlate human activities with runoff variations. This was accomplished by analysing the non-stationarity degree (quantified by the HC index) of the runoff–runoff relationship between a baseline period and an altered period, without climate data. This method provides a new perspective for correlating human activities with runoff variations, especially relevant to mountainous areas that lack climate data.

Our application of AMHIRBSH in the Baiyangdian drainage basin indicated a significant reduction in runoff beginning in the 1980s. The degree of non-stationarity of the runoff–runoff relationship was 35.64% between the baseline and altered periods, suggesting that there was a human contribution to the decrease in runoff. With further analysis of changes in evapotranspiration, the specific human activity was identified as being probably related to the massive afforestation programmes conducted in North China since the 1980s, which resulted in increased evapotranspiration and decreased runoff.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Key Research and Development Program of China [2017YFC0404505]; Major Science and Technology Program for Water Pollution Control and Treatment of China [2018ZX0711001]; National Natural Science Foundation of China [51679008,52070024]; Innovative Research Group of the National Natural Science Foundation of China [51721093]; Svenska Forskningsrådet Formas [942-2015-740]; Swedish National Space Agency [180/18].

Appendix

Mann-Kendall test

The trend assessment of hydrometeorological variables (i.e. precipitation, temperature and runoff) was carried out using the popular nonparametric Mann-Kendall (M-K) test (Yue et al. 2002). For time series samples $X=x_1,x_2,\dots ,x_n$, the null hypothesis is that the samples are independent and identically distributed, and the alternative hypothesis of a two-sided test is that the samples are not distributed identically. The M-K test statistic S obeying normal distribution is given by:

(A1) $S=\sum _{k=1}^{n-1}\sum _{j=k+1}^n\mathrm{s}\mathrm{g}\mathrm{n}\left(x_j-x_k\right)$(A1)

The mean and variance of S are given by:

(A2) $E\left[S\right]=0$(A2)

(A3) $\mathrm{v}\mathrm{a}\mathrm{r}\left[S\right]=\frac{n\left(n-1\right)\left(2n+5\right)-{\sum }_{m=1}^n{t}_{m}m\left(m-1\right)\left(2m+5\right)}{18}$(A3)

where t_m is the number of ties of extant m. Next, the standardized test statistic Z is computed by:

(A4) $Z=\left\{\begin{array}{cc}\frac{S-1}{\sqrt{\mathrm{v}\mathrm{a}\mathrm{r}\left(S\right)}},& S>0\\ 0,& S=0\\ \frac{S+1}{\sqrt{\mathrm{v}\mathrm{a}\mathrm{r}\left(S\right)}},& S<0\end{array}\right.$(A4)

The null hypothesis will be rejected when |Z| ≥ Z_1-α/2 is satisfied with a significance level of α, which indicates a significant increasing or decreasing trend for the time series. A positive Z value indicates an increasing trend, while a negative Z value shows a negative trend.

Pettitt test

The breakpoint in the annual runoff time series between the baseline period and the altered period was defined using the nonparametric Pettitt test (Pettitt 1979). For a sequence of random variables divided into two groups represented by x₁, …, x_t and x_t+1, …, x_n, the null hypothesis is that they are from the same population with no changing point, whereas the alternative hypothesis is that each group has its corresponding common distribution function, and the breakpoint is defined at t. The test statistic is given by:

(A5) $U_{t,n}=\sum _{i=1}^t\sum _{j=t+1}^n\mathrm{s}\mathrm{g}\mathrm{n}\left(x_i-x_j\right)\mathrm{f}\mathrm{o}\mathrm{r}t=1,\dots ,n$(A5)

The test statistic counts the number of times the value of a member of the first group exceeds the value of a member of the second group. The most probable breakpoint is found by the statistic k(t):

(A6) $k\left(t\right)=\underset{1<t<n}{max}\left|U_{t,n}\right|$(A6)

and the associated probabilities used in significance testing are given as

(A7) $p\cong 2exp\left\{-\frac{6{\left(k_n\right)}^2}{\left(n^3+n^2\right)}\right\}$(A7)

The null hypothesis will be rejected if the value of p is smaller than that of the specific significance level.

Data process in R software

The M-K test was implemented using “tfpwmk” in the R package “modifiedmk,” and the breakpoint detection was implemented using “pettitt.test” in the R package “trend.” Two R packages were used to build a copula function. The first package, “fitdistrplus,” helps to fit a parametric univariate distribution to a given dataset, and in this way the marginal distributions of runoff pairs were obtained. In the second package, “copula,” the “mvdc” function helps to construct bivariate distributions given their margins. Next, the probability density functions and cumulative distribution functions of the bivariate distributions can be easily attained using the functions “dmvdc” and “pmvdc,” respectively. The function “gofCopula” helps to carry out the goodness-of-fit test. In addition, the function “selectedCopula” in package “VineCopula” can also be used for a fast selection of the best copula. These R packages are freely available at the Comprehensive R Archive Network (cran.r-project.org).