Characteristic Analysis and Uncertainty Assessment of the Joint Distribution of Runoff and Sediment: A Case Study of the Huangfuchuan River Basin, China

Abstract

Exploring the relationship between runoff and sediment elements in a river basin is a prerequisite for realizing the scientific management scheme of runoff and sediment. In this study, six commonly applied probability distributions are utilized to fit the marginal distribution, and three Archimedes copulas are used to fit the joint distribution to build a joint probability distribution model of river runoff and sediment in sandy areas. The synchronous and asynchronous encounter probabilities of runoff and sediment are calculated. The uncertainties of marginal distribution, parameter estimation, and copula function in the process of constructing the joint distribution model framework are analyzed. The results indicate that: (1) The runoff and sediment series from 1954 to 2015 of the Huangfuchuan River basin are divided into three stages by using the cumulative anomaly method and the double mass curve method, and the runoff and sediment in the three stages have strong correlations. In the $T_a$ (1954–1978) and $T_b$ (1979–1996) stages, the optimal joint distribution functions of runoff and sediment are Gumbel, and in the $T_c$ (1997–2015) stage the optimal joint distribution function is Clayton; (2) The synchronous probabilities of runoff and sediment series in the three stages are 69.84%, 84.82%, and 70.72%, respectively, which are much greater than the asynchronous frequencies of abundance and depletion, and this showed that the conditions of runoff and sediment in the river basin are consistent; (3) The joint distribution function is sensitive to the choice of marginal distributions, parameters, and copula functions, and the optimal marginal distribution function, optimal copula function, and the parameters selected by the maximum likelihood estimation method can better fit the runoff-sediment relationship in the river basin and reduce the process uncertainty.

1. Introduction

The interdependence and response processes of runoff and sediment elements in river basins have always been the focus of water science research [1,2,3]. Runoff and sediment are carriers of material transport in river systems and outputs of runoff and sediment production processes in river basins [4,5,6]; their relationship and synergistic evolution play an irreplaceable role in the function of watershed systems [7]. The research of the joint probability of runoff and sediment is a prerequisite for realizing the scientific management scheme of runoff and sediment, which is of great importance in the utilization of water resources, sustainable development of the river basin, flood control, soil and water conservation, and ecological construction [8,9,10,11].

Hydrological events are complex stochastic events coupled with multiple factors, and correlations commonly exist among their internal variables. The use of traditional univariate hydrological frequency analysis is difficult to fully reflect the dependencies among the variables of hydrological events and has limitations [12,13,14]. Traditional multivariate hydrological analysis calculations also mainly tend to adopt linear correlation multivariate distribution models with the same marginal distribution [15,16], limiting the accuracy of their results. Copula function is a breakthrough method for establishing multivariate joint distributions [17], which can connect two or more correlated variables [18] and study the marginal distribution and correlation structure, respectively, without restricting the type of marginal distributions. It has high adaptability and is widely used in hydrological and water resources application fields, such as frequency analysis, risk analysis, stochastic simulation of rainfall, floods, droughts, etc. [19,20,21,22,23] Golian et al. employed the copula method to study the joint probability distribution of rainfall depth and peak flow [24]; Bacchi and Balistrocchi proposed a criterion for deriving flood frequency curves from rainfall and duration distribution, and evaluated the return period of bivariate rainfall events [25]; Zhang et al. utilized copula functions to focus on the relationship between groundwater depth changes and three selected control factors from a probability perspective, and compared the differences between the results of the two-dimensional and three-dimensional copula functions [18]; Requena et al. applied the copula approach to perform multivariate flood frequency analyses, and the uncertainty involved in selecting copula to characterize the dependency structure of short data sequences is analyzed [26]; Chang et al. constructed a two-dimensional copula drought risk model of drought duration and severity, and the drought risk of Weihe River Basin was fully assessed [27]; Li et al. constructed a model of runoff and sediment over two stations in the Anning watershed 2010~2015 using the copula analysis method, analyzing the synchronous and asynchronous probabilities of runoff and sediment from both time and space perspectives [28]; Dodangeh et al. established a joint probabilistic model of the extreme rainfall–runoff, using annual maximum precipitation, corresponding historical and simulated runoff data [29]; Salvadori et al. introduced a framework for multivariate copula for handing multivariate disaster scenarios and evaluated the threat probability of natural disasters [30].

Due to its complex physical mechanisms and numerous influencing factors (including climate, meteorology, topography, geomorphology, underlying surface, etc.,), uncertainties in hydrological processes are prevalent, and the degree of uncertainty is directly related to the complexity of the system [31]. Uncertainties in the framework of multidimensional variable models based on copula theory include the uncertainty in the marginal distribution, the uncertainty in the selection of copula functions, and the uncertainties in sample sampling and parameter estimation [32,33]. The uncertainty in the process of constructing the copula joint distribution model has a significant impact on the accuracy of the conclusion and the correctness of the decision. Comprehensive assessment of uncertainty in joint distribution models has complexity, and the quantitative analysis of uncertainty has rarely been addressed in previous studies, and related research is currently at the initial stage.

The Yellow River has the highest sediment concentration in the world and is the world’s most famous sandy river. In recent years, under the combined influence of climate change and human activities, the runoff–sediment process of the Yellow River has changed significantly [34]. The Huangfuchuan River basin is a typical basin in the concentrated source area of sandy and coarse sediment in the middle reaches of the Yellow River (Heilong interval), which is one of the main production areas of coarse sediment in the Yellow River and the basin with the highest erosion intensity among the tributaries of the Yellow River. Since the 20th century, due to the over-exploitation of coal resources in the Huangfuchuan River basin, the groundwater level has continued to decline, runoff has shown a decreasing trend, and droughts have continued to occur in the arid zone. In addition, arsenic sandstone has been exposed in a large area in the basin, and the original vegetation has been severely damaged. These features are extremely rare in China and even around the world. Under the changing environment, the modeling and analysis of the synergistic evolution of runoff and sediment in the Huangfuchuan River basin can help clarify the runoff-sediment relationship and the change patterns in the Yellow River basin and help guide watershed management planning and soil erosion control. It is the theoretical basis for establishing a runoff-sediment regulation system for the harmonious development of man and nature. It provides the scientific basis and data support for runoff-sediment resource utilization in the region, is theoretical support for the formulation of regional sustainable development countermeasures, and is of practical significance for the ecological environment construction of the Loess Plateau, the healthy management of the Yellow River, and the high-quality development of the Yellow River basin.

In this paper, the Huangfuchuan River basin is taken as the research object, the joint probability distribution model of runoff and sediment in a sandy area is constructed, and the uncertainties in the process of constructing the copula joint distribution function are analyzed. Firstly, the cumulative anomaly method, the double mass curve (DMC) method, and the Mann-Kendall test are used to analyze the runoff and sediment discharge series from 1954 to 2015, determine the mutation points, and divide them into three stages; the copula function is introduced to construct the joint distribution model of runoff and sediment in each stage, calculate the synchronous and asynchronous encounter probabilities of runoff and sediment, and analyze the joint return period and the co-occurrence return period. Finally, the uncertainties of the marginal distribution, parameters estimation, and the copula function are analyzed.

2. Materials and Methods

2.1. Study Area and Data Sources

Huangfuchuan is located on the upper right bank of the middle reaches of the Yellow River between Hekou town and the Longmen area. It is a first-class tributary of the Yellow River, which originates in the eastern part of the Ordos Plateau, the transition zone between the Loess Plateau and the desert grassland, flowing through the Jungar Banner of Inner Mongolia Autonomous Region, and flows into the Yellow River in Fugu County, Shanxi Province. As shown in Figure 1, the mainstream has a length of 137 km and a basin area of 3246 km². The basin is located inland, influenced by the southeast monsoon, and has a continental monsoon climate with an average precipitation of 350–450 mm per year, and the distribution is extremely uneven during the year, with more than 80% of annual precipitation concentrated in June to September. Huangfuchuan belongs to a sandy and coarse sediment river, with an average annual sediment discharge of 61.2 million tons and an average annual sediment modulus of 19.1 million tons/(km²·a), which is one of the most serious zones of soil erosion in the Loess Plateau area and transports a large amount of coarse sediment to the Yellow River. The Huangfuchuan hydrologic station is close to the basin outlet, controlling a basin area of 3199 km². The main tasks of the station include precipitation, water level, flow, sediment testing, etc. It is a national basic hydrological station.

The hydrological data used in this study were from the Yellow River Conservancy Commission of the Ministry of Water Resources and the Yellow River Basin Hydrological Yearbook, and the measured runoff and sediment data from 1954 to 2015 were selected from Huangfuchuan Hydrological Station, the outlet station of Huangfuchuan River basin. We collected data from 12 rainfall stations in the Huangfuchuan River basin and used the interpolation method of selecting neighboring stations to interpolate the missing data. The Thiessen polygons module of Arcgis 10.3 software was used to obtain basin surface precipitation data.

2.2. Mutation Analysis

2.2.1. The Cumulative Anomaly Method

The cumulative anomaly method is basic and widely used in mutation detection [35,36]. This method cumulates the departure value between each datum and the average of the time series, and according to the change of cumulative value slope and extreme points, the mutation point is detected to judge the time of time series mutation. For the time series $X_i\left(X_1,X_2,\cdots ,X_n\right)$, cumulative anomalies $C_t$ can be expressed as:

$$C_t={\displaystyle \sum _{i=1}^t(X_i-\overline{X})}\hspace{1em}(t=1,2,3,\cdots ,n)$$

(1)

where $X_i$ is the discrete values of the data series, and $\overline{X}$ is the average of the data series.

2.2.2. The Double Mass Curve Method

The double mass curve (DMC) method is a method based on linear regression assessment to investigate long-term trends and detect change points of hydro-meteorological time series [37]. The basic principle of this method is to accumulate the two variables under investigation step by step for the same length of time, take the two accumulated values as horizontal and vertical coordinates, respectively, and observe the change in slope. If the ratio between the two variables remains constant, the relationship is linear over the same time. Conversely, it indicates that the original relationship between the two variables has changed. The time corresponding to the point where the slope of the straight line deviates significantly is the year when the variable begins to undergo significant changes. Namely:

$${\displaystyle \sum S=a{\displaystyle \sum R+b}}$$

(2)

where S and R denote sediment discharge and runoff, respectively, and a and b are coefficients.

2.3. Trend Analysis

2.3.1. The M-K Trend Test

The Mann-Kendall trend test [38,39] is a non-parametric test method that does not require samples to follow a certain distribution. It is often used for analyzing the trend of time series changes in hydrological and meteorological elements. In the M-K test, the calculation formula for the test statistic S is:

$$S={\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=i+1}^n\mathrm{sgn}\left(X_j-X_i\right)}}$$

(3)

where n is the length of the time series, and $X_i$ and $X_j$ represent the values at the $i_{\mathrm{th}}$ and $j_{\mathrm{th}}$ time points, respectively $\left(j>i\right)$. Among them:

$$\mathrm{sgn}\left(X_j-X_i\right)=\left\{\begin{array}{cc}+1\hfill & if \left(X_j-X_i\right)>0\hfill \\ 0\hfill & if \left(X_j-X_i\right)=0\hfill \\ -1\hfill & if \left(X_j-X_i\right)<0\hfill \end{array}\right.$$

(4)

$$Var\left(S\right)=\frac{1}{18}\left[n\left(n-1\right)\left(2n+5\right)-{\displaystyle \sum _{p=1}^q{t}_p\left(t_p-1\right)\left(2t_p+5\right)}\right]$$

(5)

where q is the number of the tied values, and $t_p$ is the number of ties for the $p_{\mathrm{th}}$ value.

The calculation formula for Mann-Kendall trend test statistics is:

$$Z_{S\_MK}=\left\{\begin{array}{cc}\frac{S-1}{\sqrt{Var\left(S\right)}}\hfill & if S>0\hfill \\ 0\hfill & if S=0\hfill \\ \frac{S+1}{\sqrt{Var\left(S\right)}}\hfill & if S=0\hfill \end{array}\right.$$

(6)

In a two-tailed test, if $\left|Z_{s\_MK}\right|\ge Z_{1-\alpha /2}$, this indicates a significant trend in the sequence at significance level $\alpha$.

2.3.2. The TFPW—MK Test

In trend analysis, in order to eliminate serial autocorrelation, it is usually necessary to perform pre-whitening processing. The Mann-Kendall test with trend-free pre-whitening (TFPW-MK) [40] is a pre-whitening treatment based on M-K, which effectively reduces the influence of series autocorrelation and ensures the accuracy of the test. The operation steps are as follows:

(1)

Construct the trend-free series Y_t:

$$Y_t=X_t-\beta t$$

(7)

where $X_t$ is the series value at the year t, and $\beta$ is the trend slope.

(2)

Calculate the first-order autocorrelation coefficient of the series Y_t, and perform a significance test. When $r_1$ is not significant, a detrended series is used for trend detection. On the contrary, the detrended series is pre-whitened, and the residual series is recorded as $Z_t$.

$$r_1=\frac{\frac{1}{n-1}{\displaystyle \sum _{i=1}^{n-1}\left(x_i-E\left(x_i\right)\right)\left(x_{i+1}-E\left(x_i\right)\right)}}{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(x_i-E\left(x_i\right)\right)}^2}}$$

(8)

$$Z_t=Y_t-r_1{Y}_{t-1}$$

(9)

(3)

The residual sequence is reorganized with the trend term to form a trend-free pre-whitening series $X_t^{*}$:

$$X_t^{\ast }=Z_t-\beta t.$$

(10)

(4)

Perform an M-K trend test on series $X_t^{*}$.

2.4. Marginal Distribution Model

Normal (NOR), gamma (GAM), exponential (EXP), Weibull (WBL), generalized extreme value (GEV), and log-normal (LOGN) are the commonly used types of probability distributions. The above six univariate distribution functions are selected to fit the runoff and sediment discharge of Huangfuchuan hydrologic station as their marginal distribution functions. The maximum likelihood method is used to estimate the parameters, and the Kolmogorov-Smirnov (K-S) test is used to test whether the samples obeyed a certain marginal function, and the principle of the minimum deviation square (OLS) and Akaike information criterion (AIC) [41] are used to evaluate the goodness-of-fit of the marginal function.

(1) Kolmogorov-Smirnov (K-S) Test: The K-S test is a test based on the cumulative distribution function to test whether two empirical distributions or an empirical distribution and a theoretical distribution are different. The hypothesis test problem is that the overall sample obeys a particular distribution (H₀) and the overall sample distribution does not obey a particular distribution (H₁).

(2) Ordinary least square (OLS): The OLS is the minimum deviation square criterion. The OLS is used to measure the deviation between the theoretical and empirical values, and the smaller the value the better the curve fit, which is calculated as follows:

$$OLS=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=1}^n{\left(P_{ei}-P_i\right)}^2}}$$

(11)

where $p_i$ and $p_{ei}$ are theoretical and empirical frequencies, respectively.

(3) Akaike information criterion (AIC): AIC is one of the criteria used to evaluate the merit of the model. The smaller the AIC value, the better the curve fit, and its formula is as follows:

$$AIC=n\ln (MSE)+2m$$

(12)

$$MSE=\frac{1}{n-m}{\displaystyle \sum _{i=1}^n{\left(p_{ei}-p_i\right)}^2}$$

(13)

where MSE is the root-mean-square deviation, n is the sample size, and m is the number of model parameters.

2.5. Copula Joint Distribution Model

The copula function was proposed by Sklar in 1959 [42]. It is a multi-dimensional function defined on the interval [0, 1], connecting the joint distribution and the respective marginal distribution.

Sklar’s theorem is the theoretical basis of the copula function. For two random variables X and Y, there are two-dimensional joint distribution functions $F_{X,Y}\left(x,y\right)$, and the marginal distribution functions of the variables X and Y are $R_X\left(x\right)$ and $S_Y\left(y\right)$, respectively. Then there exists a copula function C, such that:

$$F_{X,Y}\left(x,y\right)=C\left(R_X\left(x\right),S_Y\left(y\right)\right)\triangleq C(u,v)$$

(14)

If the marginal distribution functions $R_X\left(x\right)$ and $S_Y\left(y\right)$ are continuous functions, the copula function C is unique.

There are many classifications of copula functions, among which Archimedean copulas are widely used in the field of hydrology because of their simple form and strong adaptability. In this paper, the three most commonly used two-dimensional Archimedean copula functions are selected to construct the two-dimensional joint distribution, and their cumulative distribution functions are shown in

Table 1. Common Archimedean copulas.

Copulas	$\mathit{C}\left(\mathit{u},\mathit{v}\right)$	Parameters
Gumbel	$\exp \left[-{\left({\left(-\ln u\right)}^{\theta }+{\left(-\ln v\right)}^{\theta }\right)}^{1/\theta }\right]$	θ ≥ 0
Clayton	${\left(u^{-\theta }+v^{-\theta }-1\right)}^{-1/\theta }$	θ ≠ 0
Frank	$-\frac{1}{\theta }\ln \left[1+\left(\frac{\left(e^{-\theta u}-1\right)\left(e^{-\theta v}-1\right)}{e^{-\theta }-1}\right)\right]$	θ ≥ 1

. The parameters θ of Copula functions are estimated by the maximum likelihood estimation method, and the OLS criterion, AIC information criterion, and Bayesian information criterion are used for testing the goodness of fit to determine the optimal copula functions.

Bayesian information criterion (BIC):

$$BIC=n\ln \left(MSE\right)+m\ln \left(n\right)$$

(15)

where n is the sample size, and m is the number of model parameters. The smaller the BIC value, the better the fitting effect.

2.6. Runoff-Sediment Encounter Probability Analysis

Using pf = 37.5% and pk = 62.5% as the frequency criteria for the division of runoff and sediment rich or poor status, respectively, the hydrological probability method is used to divide the rich-poor of runoff and sediment, and the encounter combination between the two can be divided into the following nine scenarios shown in

Table 2. Encounter combination of rich-poor runoff and sediment.

		Runoff
		Rich	Normal	Poor
Sediment	Rich	$P_1=p\left(R\ge r_{pf},S\ge s_{pf}\right)$	$P_2=p\left(r_{pk}<R<r_{pf},S\ge s_{pf}\right)$	$P_3=p\left(R\le r_{pk},S\ge s_{pf}\right)$
	Normal	$P_4=p\left(R\ge r_{pf},s_{pk}<S<s_{pf}\right)$	$P_5=p\left(r_{pk}<R<r_{pf},s_{pk}<S<s_{pf}\right)$	$P_6=p\left(R\le r_{pk},s_{pk}<S<s_{pf}\right)$
	Poor	$P_7=p\left(R\ge r_{pf},S\le s_{pk}\right)$	$P_8=p\left(r_{pk}<R<r_{pf},S\le s_{pk}\right)$	$P_9=p\left(R\le r_{pk},S\le s_{pk}\right)$

Note: R represents runoff; S represents sediment discharge; $r_{pf}$ ($s_{pf}$) represents the rich runoff (sediment) determined by the frequency of 37.5%; $r_{pk}$ ($s_{pk}$) represents the poor runoff (sediment) determined by the frequency of 62.5%.

.

2.7. Return Period Analysis

The return period refers to the average recurrence interval between the occurrences of an observation greater than or equal to 1 time during the statistical period of the data. The return period is an important indicator in water engineering. There are two cases of binary return period: one is called the joint return period, which means the return period when the design value of either of the two variables is exceeded ($X\ge x$ or $Y\ge y$), and here it is expressed by $T^{or}$; the second is called the co-occurrence return period, which means the return period when the design values of two variables are exceeded at the same time ($X\ge x$ and $Y\ge y$), and it is expressed by $T^{and}$.

$$\text{The joint return period}: T^{or}=\frac{1}{1-C\left(u,v\right)}$$

(16)

$$\text{The co-occurrence return period}: T^{and}=\frac{1}{1-u-v+C\left(u,v\right)}$$

(17)

3. Results

3.1. Runoff and Sediment Characteristics

Runoff and sediment discharge are important hydrological information to ensure the ecological security of rivers. To understand the changes between runoff and sediment in the Huangfuchuan watershed over 62 years, the interannual variation of annual precipitation, corresponding annual runoff, and annual sediment in the Huangfuchuan watershed from 1954 to 2015 is plotted, as shown in Figure 2. It can be seen from the figure that the interannual trends of runoff and sediment are basically consistent.

The cumulative anomaly method is used to test the mutation points of the runoff series and sediment series of the Huangfuchuan River basin from 1954 to 2015, and the cumulative anomaly curves of annual runoff and annual sediment discharge in the Huangfuchuan River basin are shown in Figure 3. Observing the trend of cumulative anomaly curves, it can be found that from 1954 to 1978, the cumulative anomaly curves of annual runoff and annual sediment discharge show both fluctuating upward trends, indicating that the runoff and sediment in the Huangfuchuan River basin were in a relatively abundant state during this period. From 1979 to 1996, the cumulative anomaly curves of annual runoff and annual sediment discharge fluctuated but are relatively flat overall, indicating that the Huangfuchuan River basin was in a relatively medium state during this period. From 1997 to 2015, the cumulative anomaly curves of annual runoff and annual sediment discharge show a significant decreasing trend, indicating that the runoff and sediment in the Huangfuchuan River basin were in a relatively low period.

To further determine the year of the mutation, the double mass curve method is adopted to conduct the test of abrupt changes on both the annual runoff series and annual sediment series of the Huangfuchuan River basin, as shown in Figure 4. It can be found that the runoff-sediment accumulation curve of the Huangfuchuan River basin significantly shifted between 1979 and 1996. The period before the shift is generally regarded as the period of no or less disturbance by human activities, i.e., the base period. After the base period, the slope of the accumulation curve changes again, and the slope of the curve decreases, indicating that there is a downward trend of sediment discharge relative to the runoff, so the periods of 1979–1996 and 1997–2015 can be classified as human activity periods. The determination coefficient R² of the fitted relational equation for each period is above 0.97, which is a good fit.

According to the above two mutation tests, the results of variability analysis can be obtained. The years 1979 and 1996 are the mutation points of the runoff and sediment series in the Huangfuchuan River basin, and the runoff and sediment series in the Huangfuchuan River basin are divided into three stages: 1954–1978 ($T_a$), 1979–1996 ($T_b$), and 1997–2015 ($T_c$).

Table 3. Statistical characteristics of runoff and sediment at each stage.

Statistic	${\mathit{T}}_{\mathit{a}}$		${\mathit{T}}_{\mathit{b}}$		${\mathit{T}}_{\mathit{c}}$		1954~2015
	Runoff	Sediment	Runoff	Sediment	Runoff	Sediment	Runoff	Sediment
$Z_{s\_MK}$	−0.9576	−0.4437	−1.1363	−0.8333	−2.0991 **	−2.7289 ***	−5.7825 ***	−5.3209 ***
$Z_{S\_TFPW}$	−0.3721	−0.0758	−0.5355	−0.1236	−1.8939 *	−2.6514 ***	−5.4948 ***	−5.1339 ***

Note: *, **, and *** represent significance at 90%, 95%, and 99% confidence levels corresponding to critical z values of 1.65, 1.96, and 2.58, respectively.

shows the results of trend tests for runoff and sediment in the Huangfuchuan River basin at these three stages. It can be seen from the table that there are differences in the trends of annual runoff and sediment discharge at the three stages of division. The $Z_{s\_MK}$ value of the M-K trend test for runoff at the stage $T_c$ is −2.0991, $\left|Z_{s\_MK}\right|$ > 1.96, has passed the significance test with a confidence level of 95%. The $Z_{s\_MK}$ value of the M-K test for sediment at the stage $T_c$ is −2.7289, $\left|Z_{s\_MK}\right|$ > 2.58, has passed the significance test of 99% confidence level. The statistics of the TFPW_MK test are slightly smaller than those of the M-K test, but the significance and trend are basically consistent. This indicates that the runoff and sediment in the Huangfuchuan River basin show a significant decreasing trend in the $T_c$ stage (1997–2015). In other periods, the runoff and sediment show a decreasing trend, but the trends are not significant.

The main statistical characteristic values of runoff and sediment discharge in these three stages are given in

Table 4. Statistical characteristics of runoff and sediment in each stage.

Stage	Characteristic Variables	Mean	Maximum	Minimum	Standard Deviation	${\mathit{C}}_{\mathit{v}}$	$\mathit{S}\mathit{K}$
$T_a$	Runoff ($MCM$)	185.437	508	41.080	115.247	0.621	1.534
	Sediment ($MT$)	58.006	171.071	5.224	41.639	0.718	1.153
$T_b$	Runoff ($MCM$)	134.775	435.940	25.110	96.965	0.719	1.699
	Sediment ($MT$)	43.726	147.507	5.155	39.608	0.906	1.337
$T_c$	Runoff ($MCM$)	37.676	103.449	0	34.475	0.915	0.886
	Sediment ($MT$)	8.667	29.076	0	9.702	1.119	1.058

. As can be seen from Table 4, the average runoffs are 185.437, 134.775, and 37.676 million cubic meters, respectively, and the average sediments are 58.006, 43.726, and 8.667 million tons, respectively, at the three stages. The runoff and sediment in each stage show a decreasing trend as a whole with the continuation of time, which is consistent with the conclusion obtained by the cumulative anomaly method (Figure 3). The coefficient of variation $C_v$ represents the degree of dispersion of the data. From the overall perspective, the dispersion degree of sediment discharge is greater than that of runoff, and the dispersion degree of both runoff and sediment discharge tends to increase with the continuation of time at each stage. The skewness coefficient SK describes the degree of distribution deviation from symmetry, and the skewness of both runoff and sediment shows a trend of first increasing and then decreasing in all three stages.

The intra-annual distributions of the mean values of runoff and sediment in three stages of the Huangfuchuan River basin are shown in Figure 5. The runoff and sediment are very unevenly distributed within the year, mainly concentrated in the flood season (June to September), and the sediment discharge in each stage of the flood season accounts for 98.11%, 99.19%, and 99.80% of the annual sediment discharge, respectively, and the runoff in each stage of the flood season accounts for 78.18%, 85.47%, and 95.20% of the annual runoff, respectively. The unevenness of the sediment discharge is greater than that of the runoff.

Copula is a joint probability distribution function constructed by connecting the marginal distribution function of related variables. Therefore, the correlation measurement between variables is an important prerequisite. An important step is to draw a correlation analysis chart between runoff and sediment, as shown in Figure 6. The linear fitting correlation coefficient R² = 0.8449 between runoff and sediment. The vast majority of the correlation analysis point distances between runoff and sediment transport in the Huangfuchuan River Basin are densely distributed near the linear fitting line. The 95% prediction band (estimated to be the area where 95% of the data points belong) relatively evenly includes the linear fitting line, and the point distances of each era are distributed on both sides of the linear fitting line. This indicates that there is no significant systematic deviation in the relationship between water and sediment, and there is a strong correlation between runoff and sediment transport at the Huangfuchuan River basin from 1954 to 2015.

Further, the Pearson linear correlation coefficient γ, the Kendall rank correlation coefficient τ, and the Spearman rank correlation coefficient ρ are adopted to analyze the correlation between the two indicators of runoff and sediment discharge, and the closer the absolute value of the correlation coefficient is to 1, the stronger the correlation between the two variables. According to the calculation results in

Table 5. Correlation coefficients between the runoff and sediment.

Correlation Coefficient	${\mathit{T}}_{\mathit{a}}$	${\mathit{T}}_{\mathit{b}}$	${\mathit{T}}_{\mathit{c}}$
γ	0.844	0.970	0.924
τ	0.673	0.895	0.801
ρ	0.850	0.973	0.919

Note: γ denotes the Pearson linear correlation coefficient, τ represents the Kendall rank correlation coefficient, and ρ indicates the Spearman rank correlation coefficient.

, there is a very strong correlation between the runoff and sediment discharge in the three stages. Among them, the Pearson linear correlation coefficients exceed 0.84, the Kendall rank correlation coefficients are over 0.67, and the Spearman rank correlation coefficients surpass 0.85 in the three stages. The basic requirements of the copula function for variables are met.

3.2. Marginal Distributions and Copula Joint Distributions

3.2.1. Selection of Marginal Distributions

Six commonly used theoretical distribution functions, namely, the normal (NOR) distribution function, the gamma (GAM) distribution function, the exponential (EXP) distribution function, the Weibull (WBL) distribution function, the generalized extreme value (GEV) distribution function, the and lognormal (LOGN) distribution function, are used to fit the annual runoff and annual sediment data of the Huangfuchuan River basin, as shown in

Table 6. Parameters estimation and goodness of fit of the marginal distribution.

Margin	Parameters and Tests	Stage
		${\mathit{T}}_{\mathit{a}}$		${\mathit{T}}_{\mathit{b}}$		${\mathit{T}}_{\mathit{c}}$
		Runoff	Sediment	Runoff	Sediment	Runoff	Sediment
Normal	μ	1.854	5800.600	1.348	4372.600	0.377	866.706
	σ	1.176	4249.700	0.998	4075.600	0.354	996.833
	K-S	0	0	0	0	0	0
	OLS	0.123	0.077	0.089	0.097	0.096	0.121
	AIC	−100.778	−123.977	−82.909	−79.885	−84.932	−76.369
Gamma	${\mathsf{\alpha }}_1$	3.194	1.800	2.151	1.200	0.538	0.300
	${\mathsf{\alpha }}_2$	0.581	3257.600	0.627	3609.500	0.700	2867.300
	K-S	0	0	0	0	0	0
	OLS	0.082	0.041	0.051	0.055	0.0600	0.087
	AIC	−120.950	−155.961	−103.278	−100.280	−103.138	−88.927
Exponential	μ	1.8544	5800.600	1.348	4372.600	0.377	866.706
	K-S	1	0	0	0	0	0
	OLS	0.136	0.071	0.083	0.039	0.047	0.142
	AIC	−95.876	−128.155	−85.526	−112.704	−112.415	−70.276
Weibull	${\mathsf{\alpha }}_1$	2.099	6386.900	1.501	4549.100	0.326	514.919
	${\mathsf{\alpha }}_2$	1.753	1.400	1.494	1.100	0.703	0.426
	K-S	0	0	0	0	0	0
	OLS	0.093	0.043	0.054	0.053	0.051	0.085
	AIC	−114.895	−153.718	−101.288	−101.493	−108.963	−89.825
Generalized Extreme Value	$\kappa$	0.216	0.100	0.213	0.800	0.488	1.256
	α	0.666	2850.600	0.591	1712.900	0.187	265.187
	ξ	1.294	3724.400	0.865	1787.800	0.165	167.711
	K-S	0	0	0	0	0	0
	OLS	0.054	0.043	0.048	0.062	0.043	0.067
	AIC	−142.254	−152.881	−105.071	−95.941	−115.704	−98.983
Lognormal	μ	0.453	8.359	0.048	7.917	−2.143	4.486
	σ	0.580	0.886	0.757	1.067	3.123	5.306
	K-S	0	0	0	0	0	1
	OLS	0.062	0.052	0.049	0.053	0.139	0.161
	AIC	−135.023	−143.847	−104.694	−101.602	−71.133	−65.324

Note: Bold indicates the selected Marginal distribution function.

.

The parameter estimates corresponding to these six marginal distribution functions in the three stages are listed in Table 6, and the K-S test is used to fit the marginal functions. A K-S test result of 0 indicates that the distribution is satisfied by the K-S test. Except for the exponential distribution of runoff in stage $T_a$ and the lognormal distribution of sediment in stage $T_c$, the six theoretical distributions can be used to describe the annual runoff and annual sediment discharge in all other stages. The results of OLS and AIC criteria tests are also listed in the table. The test criteria are: The smaller the OLS and AIC values are, the better the fit is. It can be seen from the table that the runoff in stage $T_a$ obeys the generalized extreme value distribution and the sediment obeys the gamma distribution; the runoff in stage $T_b$ obeys the generalized extreme value distribution and the sediment obeys the exponential distribution; the runoff in stage $T_c$ obeys the generalized extreme value distribution and the sediment obeys the generalized extreme value distribution. The comparison of the cumulative probability distribution of the six theoretical distribution functions and the cumulative probability distribution of the empirical distribution functions for runoff and sediment in the three stages are shown in Figure 7. The closer the six theoretical distributions are to the empirical distribution, the better the fitting effect is.

3.2.2. Choice of Copula Joint Distribution

The three copula functions are used to construct the joint distributions of runoff and sediment, using the principles of the OLS, AIC, and BIC, that is, the smaller the value of OLS, AIC, and BIC, the better the fit is. It can be seen from

Table 7. Parameters estimation and fitting evaluation of the copula function.

Copula	Parameters and Tests	Stage
		${\mathit{T}}_{\mathit{a}}$	${\mathit{T}}_{\mathit{b}}$	${\mathit{T}}_{\mathit{c}}$
Gumbel	${\mathsf{\theta }}_1$	2.824	6.062	3.674
	OLS	0.026	0.039	0.061
	AIC	−178.881	−112.607	−102.371
	BIC	−176.444	−110.826	−100.482
Clayton	${\mathsf{\theta }}_2$	1.984	3.575	3.638
	OLS	0.044	0.058	0.057
	AIC	−151.894	−98.698	−105.086
	BIC	−149.457	−96.917	−103.197
Frank	${\mathsf{\theta }}_3$	10.937	23.927	11.485
	OLS	0.027	0.040	0.061
	AIC	−175.828	−112.037	−102.050
	BIC	−173.390	−110.256	−100.161

Note: Bold indicates the selected copula function.

that the optimal joint distribution copula functions of runoff and sediment in the $T_a$ and $T_b$ stages are both Gumbel, and the optimal joint distribution copula function of runoff and sediment in the stage $T_c$ is Clayton.

The following figure shows that the empirical cumulative frequencies of three Archimedean copula functions in the three stages are compared with the theoretical cumulative frequency consistency. By calculating R² and according to the principle of maximum R² the copula function can be determined as the optimal function. It can be seen from the Figure 8 that all the three copula functions of the three phases fall near the 45° line, and the correlation coefficients are relatively high (all are greater than 0.9593). The Gumbel copula’s fit is the highest for the stage $T_a$ and the stage $T_b$. The Clayton copula’s fit is the highest for the stage $T_c$. These are consistent with the conclusions of OLS, AIC, and BIC tests. That is, Gumbel copula is selected as the optimal copula function to connect the optimal marginal distribution functions of runoff and sediment in the Huangfuchuan River basin in the stages $T_a$ and $T_b$, and the Clayton copula is selected as the optimal copula function to connect the optimal marginal distribution functions of runoff and sediment in the Huangfuchuan River basin in the stage $T_c$.

3.3. Runoff-Sediment Encounter Probability Analysis

The copula joint distribution models of runoff and sediment determined above are applied to study the rich-poor encounter combinations of runoff and sediment at each stage in the Huangfuchuan River basin. The frequency method with $pf$ = 37.5% and $pk$ = 62.5% as the division frequency is used to classify the runoff and sediment at each stage into three states of rich, normal, and poor, and nine different encounter combinations, respectively, to analyze the changing characteristics of runoff–sediment combination and the encounter probabilities.

It can be seen from

Table 8. The encounter probabilities of runoff and sediment.

Stage	Synchronous Probability (%)				Asynchronous Probability (%)
	Rich–Rich	Normal–Normal	Poor–Poor	Total	Rich–Normal	Rich–Poor	Normal–Rich	Normal–Poor	Poor–Rich	Poor–Normal	Total
$T_a$	29.84	11.46	28.54	69.84	6.12	1.54	6.12	7.42	1.54	7.42	30.16
$T_b$	34.04	17.48	33.30	84.82	3.39	0.07	3.39	4.13	0.07	4.13	15.18
$T_c$	28.02	11.57	31.13	70.72	8.25	1.22	8.25	5.17	1.22	5.17	29.28

that: ① The synchronous probabilities of runoff and sediment in different periods are greater than the asynchronous probabilities, which indicates that the runoff and sediment in each stage before and after the mutation points are closely related, with a strong correlation. Among the synchronous probabilities in the three stages, the encounter probabilities of annual rich runoff and rich sediment are the largest, the probabilities of poor runoff and poor sediment are slightly smaller than the probabilities of same rich, and the probabilities of annual normal runoff and normal sediment are the least, which indicates that the synchronization with the same rich is the strongest in the runoff-sediment relationship. ② Due to the symmetry of the copula function, in the asynchronous probabilities, the encounter probabilities of rich runoff and normal sediment are equal to those of normal runoff and rich sediment, rich runoff and poor sediment are equal to those of poor runoff and rich sediment, normal runoff and poor sediment are equal to those of poor runoff and normal sediment. Among them, the encounter probabilities of the combinations of runoff–rich and sediment–poor, as well as runoff–poor and sediment–rich, are the smallest among the six asynchronous probabilities, which indicates that the correlation between the runoff and sediment series is strong, and the probability of occurrence of the two diametrically opposite situations is the smallest. ③ The asynchronous probabilities in the two human activity periods ($T_b$ and $T_c$) are 84.82% and 70.72%, respectively, which are increased compared with the base period ($T_a$). This is because the marginal distribution of runoff and sediment variables, the impact mechanism of runoff and sediment, and the relationship between runoff and sediment may have changed to a certain extent under the combined influence of climate change and human activities during the human activity periods. The calculation of the runoff–sediment encounter probability in the Huangfuchuan River basin can further quantify the frequency of runoff and sediment disasters in the Huangfuchuan River basin, which in turn can provide decision-making and guidance for the scientific management of runoff-sediment projects in the Huangfuchuan River basin.

3.4. Return Period Analysis

By constructing the two-dimensional copula function of runoff and sediment, the joint return period and the co-occurrence return period of the two-dimensional joint distribution of runoff and sediment are calculated using Equations (8) and (9), and the corresponding return period contours are drawn in Figure 9. It can be seen from the figure that the constructed model can obtain the return periods of different magnitudes of annual runoff and annual sediment discharge encounters combinations in the Huangfuchuan River basin. On this basis, the typical years of runoff and sediment can be selected, and the same frequency amplification method can be used to obtain the runoff and sediment ratio system with different designed runoff hydrographs and designed sediment discharge hydrographs with the same return period, which can provide references for different runoff–sediment conditions selection and generalization in the numerical simulation of the river, hydraulic model testing, and planning and dispatching of downstream cascade hubs.

4. Uncertainty Analysis

In multivariate hydrological analysis, the uncertainty of the results is mainly due to the following reasons: (1) The representativeness of the sample. The difficulty of obtaining the measured sample data and the insufficiency of the sample data length make the sample prone to some deviation in representing the population. (2) The selection of the model function. The differences in the selection of the model function and the determination of the function parameters can cause great uncertainty in the results.

4.1. Uncertainty Analysis of Margin Distributions

The multivariate joint distribution is usually constructed by specifying the univariate marginal distribution and copula function. Determination of the discrete or continuous marginal distribution is the first step of the standard two-step procedure, after which the copula parameters and functions are determined based on the marginal distribution determination [43]. Therefore, it is of great significance to determine a reasonable marginal distribution for the determination of the final joint distribution. The runoff and sediment discharge at each stage of the Huangfuchuan River basin are fitted with six marginal distribution functions. When the return period of runoff and sediment discharge is 0~100 years, the inverse functions of these six marginal distribution functions can be obtained, and the eigenvalues of runoff and sediment discharge corresponding to various return periods can be obtained and compared with those of the empirical distribution (Figure 9). It can be seen from Figure 10 that there are significant differences between different marginal distribution functions. In stage $T_a$, the return periods of runoff obtained from the generalized extreme value distribution and the lognormal distribution are the closest to those obtained from the empirical distribution, and both pass the K-S test, with the smallest values of the OLS and AIC criteria for the generalized extreme value distribution; the return periods of sediment discharge obtained from the generalized extreme value distribution, the Weibull distribution, and the gamma distribution are the closest to those obtained from the empirical distribution, and they all pass the K-S test, with the smallest values of the OLS and AIC criteria for the gamma distribution. In stage $T_b$, the return period of runoff obtained from the generalized extreme value distribution is the closest to that obtained from the empirical distribution, and the return period of sediment discharge obtained from the exponential distribution is the closest to that obtained from the empirical distribution, and both distributions pass the K-S test; their OLS and AIC criterion values are also the smallest. In stage $T_c$, the return periods obtained from the generalized extreme value distribution for both runoff and sediment are the closest to those obtained from the empirical distribution, and both distributions pass the K-S test; their OLS and AIC criterion values are also the smallest. These validation results indicate that runoff and sediment discharge fitting different marginal distribution functions differ significantly and that the univariate characteristic return periods are sensitive to the selection of marginal distribution functions.

4.2. Uncertainty Analysis of Parameters

In multivariate hydrological analysis, the selection of function parameters determines the correctness of the function and the accuracy of the results. The Monte Carlo method is a commonly used method for analyzing uncertainty problems based on probability statistics. To address the uncertainty of the selected model parameters, the Monte Carlo method [44,45] is chosen in this paper to resample the samples several times to assess the impact of parameter uncertainty on the model. The specific method is based on the joint distribution model of the copula function generated at each stage, selecting binary samples with the same length as the observation values, and applying the Monte Carlo method to simulate, with the simulation times set to 1000. The maximum likelihood estimation method is used to estimate the parameters of the marginal distribution and the runoff–sediment joint distribution model based on the simulated samples. The box plot of the simulated parameter values is shown in Figure 11. In the stage $T_a$, the runoff marginal function is the generalized extreme value (GEV) function, whose parameters are shape parameter κ, scale parameter α, and location parameter ξ; the sediment marginal function is the gamma (GAM) function, whose parameters are shape parameter ${\alpha }_1$, scale parameter ${\alpha }_2$; and the copula function is the Gumbel function with parameter θ. In the stage $T_b$, the runoff marginal function is the generalized extreme value (GEV) function, whose parameters are shape parameter κ, scale parameter α, and location parameter ξ; the sediment marginal function is the exponential (EXP) function, with parameter μ, and the copula function is the Gumbel function with parameter θ. In the stage $T_c$, the marginal functions of both runoff and sediment are the generalized extreme value (GEV) functions with parameters κ, α, and ξ, and the copula function is a Clayton function with parameter θ. It can be seen from Figure 10 that the parameters of the marginal distributions and joint distributions based on the observation values locate in the range of the simulation parameters, which are close to the mean position. Therefore, it can be concluded that the parameters selected in determining the marginal distributions and copula functions are reasonable.

4.3. Uncertainty Analysis of Copula Functions

Every copula has its own specific dependence structure [46]. The difference in the selection of the copula function has a great impact on the calculation results. Selecting a reasonable copula function is the key to constructing the joint distribution function [47]. The aforementioned copula functions are used to construct the joint distribution functions of runoff and sediment, and the differences between the three copula functions are analyzed. Since both runoff and sediment obey six marginal distributions, each year corresponds to 108 joint distribution functions, and the co-occurrence return period and joint return period of runoff and sediment corresponding to the 108 joint distributions can be calculated. For example, the year 1961 is selected for analysis, which is in the stage $T_a$. The annual runoff is 345.160 million m³ and the annual sediment discharge is 91.037 million tons in that year. Three copula functions are used to construct the joint distribution function of runoff and sediment, and the runoff and sediment discharge are also subject to the six marginal distributions mentioned above. When the Gumbel copula function is obeyed, the range of the co-occurrence return period is 5.06~13.64 years; when the Clayton copula function is obeyed, the range of the co-occurrence return period is 6.38–30.88 years; when the Frank copula function is obeyed, the range of the co-occurrence return period is 5.66–19.14 years. From the variation range of the co-occurrence return period, it can be seen that the corresponding results of the three Copula functions are quite different. The optimal marginal distribution functions of runoff and sediment are generalized extreme value (GEV) distribution and gamma (GAM) distribution, and the optimal copula function is the Gumbel copula in the stage $T_a$. The co-occurrence return period of runoff and sediment corresponding to the optimal marginal distribution and the optimal copula function is 10.89 years, while the empirical co-occurrence period of runoff and sediment in 1961 is 10.16 years, which are close to each other, indicating that the selection of the copula function is correct. Therefore, it is significant to select the appropriate copula function for constructing the optimal joint distribution function.

5. Discussions

Many studies have shown that climate change and human activities are important factors that cause changes in runoff and sediment, as well as changes in the runoff–sediment relationship. Therefore, it is necessary to study the response of runoff and sediment to climate and human activities. Many meteorological factors are closely related to changes in runoff and sediment. Precipitation and temperature are selected as representatives of meteorological factors, which act on hydrological cycles and directly or indirectly affect the process of runoff and sediment. It is important to calculate the Pearson correlation coefficient between runoff (sediment) and selected meteorological factors and analyze the impact of climate change (precipitation and temperature) on runoff and sediment changes and runoff-sediment relationships in the Huangfuchuan River basin. The calculation results are shown in

Table 9. Correlation coefficient between meteorological factors and runoff (sediment).

Item	${\mathit{T}}_{\mathit{a}}$		${\mathit{T}}_{\mathit{b}}$		${\mathit{T}}_{\mathit{c}}$
	Precipitation	Temperature	Precipitation	Temperature	Precipitation	Temperature
Runoff	0.859 **	−0.137	0.773 **	−0.008	0.619 **	−0.173
Sediment	0.773 **	−0.235	0.745 **	−0.051	0.498 *	−0.065

Note: * indicates passing the 0.05 significance level test, and ** indicates passing the 0.01 significance level test.

.

As shown in Table 9, there are significant differences in the impact factors on runoff and sediment at each stage, but there are still certain regularities. Among them, precipitation in each stage has a strong correlation with runoff and sediment, but the correlation gradually weakens over time, and the correlation between runoff (sediment) and temperature in each stage is relatively small. In other words, the dynamic change of runoff and sediment is mainly driven by precipitation, and temperature is also one of the indispensable reasons for the change of runoff and sediment. In recent decades, the climate in the Huangfuchuan River basin has shown a relatively obvious trend of warming and drying, which is one of the important reasons for the reduction of runoff and sediment in the Huangfuchuan River basin.

Human activities usually have an impact on the underlying surface of a watershed through methods such as soil and water conservation, urban construction, and domestic production water use. The changes in underlying surface conditions are mainly reflected in land use methods, vegetation cover changes, and water conservancy engineering measures. The construction of check dams is one of the important measures for soil erosion control in the Huangfuchuan River basin. The check dam project can effectively regulate runoff and play a role in reducing flood and sediment. Since the 1950s, soil and water conservation management has been carried out in the Huangfuchuan River basin, with the construction of check dams. However, due to the use of traditional farming methods and the small proportion of remaining storage capacity of check dams before 2000, most of them have already been silted and become ineffective. Therefore, the control speed is relatively slow, and the effect is not significant. Compared with the base period, the synchronous frequency of runoff and sediment increased, while the asynchronous frequency decreased in the stage $T_b$. The frequency of rich runoff and rich sediment, normal runoff and normal sediment, poor runoff and poor sediment all increased, indicating that although soil and water conservation measures at this stage have played a role in reducing runoff and sediment to a certain extent (Table 4), they have not made a contribution to regulating the coupling relationship between water and sediment.

Vegetation can intercept precipitation, increase infiltration, and play a role in reducing flood peaks and erosion. In 1999, the Huangfuchuan River basin started a large-scale project of returning cropland to forests and grassland. The normalized vegetation index (NDVI) showed an upward trend, with an increase rate of 0.095/10a. After 2000, the construction of check dams was further strengthened. By 2015, the number of check dams in the Huangfuchuan River basin had increased to 886, with a total storage capacity of 493 million m³, and the proportion of dam-controlled area in the basin reached 70%. Therefore, compared with the stage $T_b$, the synchronous frequency of runoff-sediment (especially the frequency of the same rich and same normal) in stage $T_c$ significantly decreased, while the asynchronous frequency of runoff–sediment increased during the same period, indicating that soil and water conservation measures in this stage not only reduced runoff and sediment, but also effectively reduced the sediment content of rich runoff, reduced co-occurrence probability, and adjusted the probability of extreme runoff and sediment events.

6. Conclusions

In this paper, the Huangfuchuan River basin in the middle reaches of the Yellow River is taken as the research object, and a joint distribution model of runoff and sediment is constructed by using the copula theory. The characteristics of rich–poor encounter combinations of runoff and sediment are analyzed, the encounter probabilities under each scenario are identified, and the uncertainty factors in the process of constructing the Copula joint distribution function are discussed. The main conclusions are given as follows:

(1)

Based on the cumulative anomaly method and the double mass curve method, 1979 and 1996 are identified as mutation points of runoff and sediment series in the Huangfuchuan River basin, and the runoff and sediment series in the Huangfuchuan River basin are divided into three stages. The runoff and sediment in all three stages show a strong correlation.

(2)

Based on the results of K-S, OLS, and AIC criteria tests, the runoff in the Huangfuchuan River basin obeys the GEV distribution and the sediment obeys the GAM distribution in the stage $T_a$; the runoff obeys the GEV distribution and the sediment obeys the EXP distribution in the stage $T_b$; the runoff and sediment obey the GEV distribution in the stage $T_c$. Based on the results of OLS, AIC, and BIC criteria tests, the Gumbel copula is selected as the optimal joint distribution function of runoff and sediment in the stages $T_a$ and $T_b$, and the Clayton copula is selected as the optimal joint distribution function of runoff and sediment in the stage $T_c$.

(3)

The synchronous probabilities in the three stages are 69.84%, 84.82%, and 70.72%, respectively, which are much greater than the asynchronous probabilities, indicating a strong consistency of runoff and sediment in the River basin. The synchronous probabilities in the stages $T_b$ and $T_c$ have increased compared with those in stage $T_a$, indicating that under the combined influence of climate change and human activities, the impact mechanism of runoff–sediment has become more complex and the relationship between runoff and sediment may have changed to a certain extent.

(4)

The contour maps of the joint return period and the co-occurrence return period of the two-dimensional joint distribution of runoff and sediment discharge are drawn. On this basis, the typical years of runoff and sediment can be selected, and the same frequency amplification method can be used to obtain the runoff and sediment ratio system with different designed runoff hydrographs and designed sediment discharge hydrographs with the same return period.

(5)

The uncertainties and sensitivities of the joint distribution function construction process are analyzed, showing that the optimal marginal distribution functions, the optimal copula functions, and the parameters selected by the great likelihood estimation method can better fit the runoff–sediment relationship in the river basin and reduce the process uncertainties.