Nonstationary Analysis for Bivariate Distribution of Flood Variables in the Ganjiang River Using Time-Varying Copula

Abstract

Nonstationarity of univariate flood series has been widely studied, while nonstationarity of some multivariate flood series, such as discharge, water stage, and suspended sediment concentrations, has been studied rarely. This paper presents a procedure for using the time-varying copula model to describe the nonstationary dependence structures of two correlated flood variables from the same flood event. In this study, we focus on multivariate flood event consisting of peak discharge (Q), peak water stage (Z) and suspended sediment load (S) during the period of 1964–2013 observed at the Waizhou station in the Ganjiang River, China. The time-varying copula model is employed to analyze bivariate distributions of two flood pairs of (Z-Q) and (Z-S). The main channel elevation (MCE) and the forest coverage rate (FCR) of the basin are introduced as the candidate explanatory variables for modelling the nonstationarities of both marginal distributions and dependence structure of copula. It is found that the marginal distributions for both Z and S are nonstationary, whereas the marginal distribution for Q is stationary. In particular, the mean of Z is related to MCE, and the mean and variance of S are related to FCR. Then, time-varying Frank copula with MCE as the covariate has the best performance in fitting the dependence structures of both Z-Q and Z-S. It is indicated that the dependence relationships are strengthen over time associated with the riverbed down-cutting. Finally, the joint and conditional probabilities of both Z-Q and Z-S obtained from the best fitted bivariate copula indicate that there are obvious nonstationarity of their bivariate distributions. This work is helpful to understand how human activities affect the bivariate flood distribution, and therefore provides supporting information for hydraulic structure designs under the changing environments.

1. Introduction

Flood has been one of most common natural hazards, increasingly posing a significant risk to human life and environment [1]. Flood events can be described in terms of the multivariate characteristic variables of peak discharge, water stage and suspended sediment load, which are all relevant to risk analyses. Univariate frequency analysis for these hydrological variables mentioned above can be found in many literatures [2,3,4,5]. Unlike the common frequency distribution, theoretically derived distributions of flood peak are constructed based on dominant runoff generation mechanisms [6]. Flood frequency curves can also be established by use of a continuous simulation model [7,8]. Yet, the designs of hydraulic structures (e.g., dam spillways, dikes, and river channels), cross drainage structures (e.g., culverts and bridges) and urban drainage systems require not only peak discharge (Q) value but also peak water stage (Z) and suspended sediment load (S). In fact, the multivariate frequency analysis can provide more comprehensive understanding of the flood event than simple univariate analysis [9,10,11,12,13].

Human intervention in river basins (e.g., urbanization), the effect of climatic variability (e.g., Pacific Decadal Oscillation), and increased greenhouse gasses have been suggested to cause the changes in the magnitude and frequency of extreme floods [14]. Therefore, statistical characteristics (e.g., mean and variance) of univariate flood variables as well as their dependence structure could be time-varying under changing environments (climate change and/or human activities) [15]. To better adapt such background, it is essential to identify a reasonable distribution framework to dynamically capture the evolution of statistical characteristics [16]. Under the time-varying moments model, distribution parameters can be expressed as functions of some physical variables including atmospheric and land cover indices [17,18,19]. Thus, it is a feasible way to describe the external effects in the construction of dependence structures. Copulas [20] allow researchers to easily construct the distribution by estimating marginal distributions and their dependence separately [21]. The method could provide a more convenient statistical tool for more flexibility in modeling the marginal distributions and their dependence [22]. In recent years, copula function has become a popular tool for multivariate flood frequency analysis [23,24,25,26,27,28,29].

In this study, we focus on Q, Z and S observed at Waizhou station in the Ganjiang River basin, which has undergone extensive urban and suburban development over the past 50 years. The mean elevation of the main channel dominated by sediment supply deficit or surplus at the bed can exert a direct influence on water stage [30]. The section-cross elevation in the downstream of Ganjiang River displays a change process over the past decades [31]. Forest cover rate is a valuable indicator of forests in a country or region, controlling runoff processes and sediment generation [32]. Forest coverage rate of the Ganjiang River basin has increased from 41% to 69% due to large-scale forest planting of the basin. Given the aforementioned physical correlations, the main channel elevation (MCE) and forest cover rate (FCR) are chosen as explanatory variables for modeling the time variations of Q, Z, and S of the Ganjiang River. Using time-varying copula model with the covariates MCE and FCR, the nonstationary bivariate joint distributions of different flood variables (i.e., Z-Q and Z-S) are constructed to analyze the evolution of the bivariate relationships of both Z-Q and Z-S, including their joint probabilities and conditional probabilities of Q and S given water stage.

This paper is organized as follows. The study region and the available data used in this study are described in the next section. Then, the methodologies including marginal distribution and copula model with time-varying parameters are presented. Based on deriving the best fitted marginal distributions and copulas, the results about the time-variation of joint and conditional probability for both Z-Q and Z-S are calculated. Finally, the main conclusions of the case study are summarized.

2. Study Region and Data

The Ganjiang River, which is the largest tributary of the Poyang Lake basin, is situated in Jiangxi Province (in Southeastern China) with a drainage area of 83,500 km² (Figure 1). This study area belongs to hilly region with the main land use type of woodland [33]. The Ganjiang River basin is located in the subtropical humid monsoon climate zone with distinct seasonal variations, where the annual mean precipitation is about 1680 mm. The heaviest rainfall occurs in the main flood season from April and June, and often lasts 15 to 20 days due to monsoon and typhoon rainstorms [34].

The catchment area upstream of the Waizhou station (115°50′ E, 28°38′ N) accounts for about 97% of the total basin area of the Ganjiang River, and covers an area of 80,948 km² (Figure 1). The hydrological data used in this study include daily river discharge, daily water stage, daily suspended sediment rate and yearly section-cross elevation at the hydrological station from 1964 to 2013. These data are provided by the Hydrological Bureau of Jiangxi province (http://www.jxsw.cn/). To reflect the capacity of suspended sediment transport in a flood event, suspended sediment load that measures the absolution amount of sediment appears to be more reasonable than suspended sediment rate, since the latter is a relative value for sediment. The daily suspended sediment load (S) is calculated by multiplication between the daily discharge and the corresponding suspended sediment rate. Then, the hydrological series of Q, Z and S are extracted from the same flood event in each year with the maximum discharge value criterion. Besides, the annual forest cover rates of the region are obtained from the book of China Compendium of Statistics 1949–2008 [35] and Jiangxi provincial statistical yearbook from 1983 to 2014 (http://www.jxstj.gov.cn/).

3. Methodology

The time-varying copula model, in which both marginal distributions parameters and copula parameter are expressed as functions of explanatory variables, is constructed to describe the time variations of both marginal distribution and dependence structure of bivariate flood variables. Then, joint probability and conditional probability that are derived from time-varying bivariate copula are illustrated to present nonstationarity of bivariate flood variables over time.

3.1. Calculation of Explanatory Variables

Change in the main channel elevation of alluvial rivers is a normal result about the hydraulic adjustments to adapt the variations of discharge and sediment load. In this study, the observed cross-section elevation data for each year are tabulated and graphed as the profiles of the river cross-section. Then, the main channel elevations of 20 points are extracted from the cross-section digitized graph from left to right bank [36]. To weaken the influence of some distortion/variability points, the point elevations are divided into four quartiles and the 2nd quartile value is considered as the main channel mean elevation (MCE) as follow:

$$MC{E}^t=Me(H^t)$$

(1)

where $Me(·)$ refers to the median of vector value $H^t=(h_1^t,{\mathrm{h}}_2^t,\dots ,{\mathrm{h}}_{20}^t)$, which represents the bed elevation value of the selected point.

The forest landscape structure usually contains of the six forest cover classes (i.e., coniferous forest, broadleaf forest, bamboo forest, mixed forest, economic forest, and shrubbery land). Each forest cover with crown cover percent (cp) bigger than 20% has good effects on soil and water conservation [37]. Here, the criterion with this threshold value used in national forest definition is employed to calculate forest area. Thus, the forest cover rate (FCR) is the percentage of six forest cover area to total land area in a region and can be expressed as follow:

$$FC{R}^t={\displaystyle \sum _{i=1}^6\frac{A_i^t(cp>0.2)}{A_L}}\times 100\%$$

(2)

where $A_i^t(cp>0.2)$ (i = 1,2, …, 6) is the six forest cover classes area with crown cover percent bigger than 20%. A_L stands for the total land area of the study region.

3.2. Marginal Distribution with Time-Varying Parameters

To construct the dependence structure of bivariate hydrological variables by copulas, marginal distribution of each variable should be determined firstly. In this study, five probability distributions, including four two-parameter distributions (i.e., Lognormal, Weibull, Logistic, and Gamma) and one three-parameter distributions (Pearson type III distribution) are selected as candidate marginal distributions for Q, Z and S. These distributions have been widely applied in flood frequency analysis [38,39]. The marginal distribution of a flood variable denoted by $Y$ can be specified through a parametric cumulative distribution function (CDF) $F_Y\left(y|\mu ,\sigma ,\nu \right)$, where $\mu$, $\sigma$ and $\nu$ represent location, scale and shape parameters, respectively, and are denoted by the vector $\mathsf{\theta }=(\mu ,\sigma ,\nu )$ [40].

The Generalized Additive Models for Location Scale and Shape (GAMLSS) introduced by Rigby and Stasinopoulos [40] has been widely employed in nonstationary hydrological frequency analysis for its flexibility [15,18,19,41]. If candidate marginal distribution function $F_Y\left(y^t|{\mathsf{\theta }}^t\right)$ is chosen to fit the distribution of the variable $y^t$ at any time $t$, parameters of marginal distribution can be expressed as a function of explanatory variables as follow:

$$g({\theta }^t)={\alpha }_0+{\alpha }_{1}MC{E}^t+{\alpha }_{2}FC{A}^t$$

(3)

where $g(·)$ represents the monotonic link function, which depends on the domain of statistical parameter, i.e., if the domain of the distribution parameter ${\theta }^t$ is ${\theta }^t\in ℝ$, the link function is $g\left({\theta }^t\right)={\theta }^t$, or if ${\theta }^t>0$, $g\left({\theta }^t\right)=\ln \left({\theta }^t\right)$. ${\theta }^t$ represents one of marginal distribution parameters (${\mu }_{}^t,{\sigma }_{}^t,{\nu }_{}^t$). ${\alpha }_0$, ${\alpha }_1$, and ${\alpha }_2$ are the GAMLSS parameters. $MC{E}^t$ and $FC{A}^t$ represent two candidate explanatory variables. In practice, the shape parameter ${\nu }^t$ (if it exists) is often treated as constant since it is quite unstable and difficult to estimate [42], whereas time variations are considered only in location parameter ${\mu }^t$ and scale parameter ${\sigma }^t$. What’s more, if the parameters are independent from the explanatory variables, it would be a stationary distribution with constant parameters.

The best fitted marginal distribution is determined by the Corrected Akaike Information Criterion (AICc) [43], which is derived from the likelihood function with a penalty determined by the number of model parameters. In addition, because of the potential drawbacks in the quality of the fitting of AICc [44], the goodness-of-fit (GOF), which describes how well the selected distribution fits a set of observations, for candidate distributions is assessed by the Kolmogorov–Smirnov (KS) test and Root Mean Square Error (RMSE). Besides, visual assessment of the residual plot (QQ-plot) [45] is used to examine the best fitted marginal distribution.

3.3. Time-Varying Bivariate Copula Model

In practice, the implementation of the time-varying copula model could be divided into two steps: fitting the time-varying marginal distribution of each variable firstly, and then estimating the time-varying dependence structure of the copula. In other words, both the parameters of marginal distribution and copula dependence parameter could be treated as time variation in building the time-varying copula function. According to the definition of the copula [20], time-varying bivariate copula function $H(·)$ for the hydrological variable pairs of $(Z_{}^t,Q_{}^t)$ and $(Z_{}^t,S_{}^t)$ at any time $t$ can be expressed as follows:

$$\begin{gathered} H_{Z,Q}(z_{}^t,q_{}^t)=C(F_Z(z_{}^t|{\theta }_z^t),F_Q(q_{}^t|{\theta }_q^t)|{\theta }_{zq}^t) \\ {H}_{Z,S}(z_{}^t,s_{}^t)=C(F_Z(z_{}^t|{\theta }_z^t),F_S(s_{}^t|{\theta }_s^t)|{\theta }_{zs}^t) \end{gathered}$$

(4)

where $C(·)$ represents bivariate copula function with time-varying dependence parameter ${\theta }_{zq}^t$ or ${\theta }_{zs}^t$. $F_Z(·)$, $F_Q(·)$ and $F_S(·)$ represent marginal cumulative distribution functions with corresponding time-varying parameter vectors ${\theta }_z^t=({\mu }_z^t,{\sigma }_z^t,{\nu }_z)$, ${\theta }_q^t=({\mu }_q^t,{\sigma }_q^t,{\nu }_q)$, and ${\theta }_s^t=({\mu }_s^t,{\sigma }_s^t,{\nu }_s)$, respectively. The joint distribution can be constructed by three Archimedean copula functions (i.e., Clayton, Gumbel–Hougaard, and Frank, as shown in

Table 1. The applied time-varying bivariate copulas in this study.

Copula	Cumulative Distribution Function with Time-Varying Parameters	Parameters
Clayton	$C(u,\nu |{\theta }_{}^t)={\left({(u)}_{}^{-{\theta }_{}^t}+{(\nu )}_{}^{-{\theta }_{}^t}-1\right)}^{-1/{\theta }_{}^t}$	${\theta }_{}^t>0$
Gumbel–Hougaard	$C(u,\nu |{\theta }_{}^t)=\exp \left(-{\left({\left(-\ln u\right)}^{{\theta }_{}^t}+{\left(-\ln \nu \right)}^{{\theta }_{}^t}\right)}^{1/{\theta }_{}^t}\right)$	${\theta }_{}^t>1$
Frank	$C(u,\nu |{\theta }_{}^t)=-\frac{1}{{\theta }_{}^t}\ln \left(1+\left(\exp \left(-u{\theta }_{}^t\right)-1\right)\times \left(\exp \left(-\nu {\theta }_{}^t\right)-1\right)/\left(\exp \left(-{\theta }_{}^t\right)-1\right)\right)$	${\theta }_{}^t\ne 0$

).

Because of the impacts of external forces on the Ganjiang River basin, the dependence structure of both Z-Q and Z-S could be nonstationary. Similar to the formula expression in Equation (3), the copula dependence parameter could be expressed as a function of the two explanatory variables (FCR and MCE) to reflect the nonstationarity of dependence structure. The totally four scenarios for copula dependence parameter in this paper are listed as follows:

$$g_c\left({\theta }_c^t\right)={\beta }_0$$

(5)

$$g_c\left({\theta }_c^t\right)={\beta }_0+{\beta }_{1}FC{R}^t$$

(6)

$$g_c\left({\theta }_c^t\right)={\beta }_0+{\beta }_{1}MC{E}^t$$

(7)

$$g_c\left({\theta }_c^t\right)={\beta }_0+{\beta }_{1}FC{R}^t+{\beta }_{2}MC{E}^t$$

(8)

where $g_c(·)$ depends on the domain of copula dependence parameter ${\theta }_c^t$, i.e., if ${\theta }_c^t\in ℝ$, the link function is $g_c\left({\theta }_c^t\right)={\theta }_c^t$ or if ${\theta }_c^t>0$, $g_c\left({\theta }_c^t\right)=\ln \left({\theta }_c^t\right)$. ${\beta }_0$, ${\beta }_1$ and ${\beta }_2$ are the model parameters.

The dependence parameter of Archimedean copula family can be estimated by Inference Function for Margins method (IFM) [46]. The IFM method estimates the parameters through maximization of the log-likelihood function of a copula. The most fitted copula function with time-varying dependence parameter is determined by AICc criterion. The Cramér-von Mises test (CM) [47] and RMSE are used to test the GOF of the copula functions. In addition, similar with the test of the univariate distribution, QQ plot is also used to examine the best fitted copula based on the joint probability derived by the selected copula.

3.4. Joint and Conditional Probability under Nonstationary Framework

Joint probability gives the probability that each component falls in any particular range or discrete set of values specified for these variables. Referring to the multivariate frequency analysis, the joint probability of two variables usually includes three cases in general (i.e., KEN, OR, AND). Take AND for the further analysis in this study, AND ($\wedge$) means the space in which all the variables exceed corresponding values simultaneously [48]. Under the time-varying copula framework, the joint probability $P_{ZQ}^{\wedge }$ and $P_{ZS}^{\wedge }$ can be calculated based on the best fitted copula function as follows:

$$\begin{array}{l}{P}_{ZQ}^{\wedge }=P\left(Z\ge z_{}^{\ast }\wedge Q\ge q_{}^{\ast }\right)\\ \hspace{1em}\hspace{1em}=1-F_Z(z_{}^{\ast }|{\mathsf{\theta }}_z^t)-F_Q(q_{}^{\ast }|{\mathsf{\theta }}_q^t)+C\left(F_Z(z_{}^{\ast }|{\theta }_z^t),F_Q(q_{}^{\ast }|{\theta }_q^t)|{\theta }_{zq}^t\right)\\ P_{ZS}^{\wedge }=P\left(Z\ge z_{}^{\ast }\wedge S\ge s_{}^{\ast }\right)\\ \hspace{1em}\hspace{1em}=1-F_Z(z_{}^{\ast }|{\mathsf{\theta }}_z^t)-F_S(s_{}^{\ast }|{\mathsf{\theta }}_s^t)+C\left(F_Z(z_{}^{\ast }|{\theta }_z^t),F_S(s_{}^{\ast }|{\theta }_s^t)|{\theta }_{zs}^t\right)\end{array}$$

(9)

where $z^{\ast }$, $q^{\ast }$, and $s^{\ast }$ are the threshold values of Z, Q and S, respectively. In fact, a given joint probability can correspond infinite data couples, among which there exist the most likely combination with the largest probability density [49]. The most likely combinations conditioned on the given joint probability k can be expressed as follows:

$$\begin{gathered} \left(z,q\right)=\underset{H(z,q)=k}{\arg \max }c\left(F_Z(z|{\mathsf{\theta }}_z^t),F_Q(q|{\mathsf{\theta }}_q^t)|{\theta }_{ZQ}^t\right)·f_Z\left(z|{\mathsf{\theta }}_z^t\right)·f_Q\left(q|{\mathsf{\theta }}_q^t\right) \\ \left(z,s\right)=\underset{H(z,s)=k}{\arg \max }c\left(F_Z(z|{\mathsf{\theta }}_z^t),F_S(s|{\mathsf{\theta }}_s^t)|{\theta }_{ZS}^t\right)·f_Z\left(z|{\mathsf{\theta }}_z^t\right)·f_S\left(s|{\mathsf{\theta }}_s^t\right) \end{gathered}$$

(10)

where $c(·)$ represents copula density function. $f_Z(·)$, $f_Q(·)$ and $f_S(·)$ represent the marginal distribution density functions of Z, Q and S, respectively. $\underset{H(·)=k}{\arg \max }$ stands for argument of the maxima, which is the set of points, $\left(z,q\right)$ or $\left(z,s\right)$, for which the function $c(·)$ attains the function’s largest value conditioned on $H(·)=k$.

To display the time variation of flood variables with given water stage, the conditional probability of peak discharge or suspended sediment load could be derived from the joint probability as Equation (9). The conditional probability based on nonstationary bivariate copula can be expressed as follows:

$$\begin{array}{l}{P}_{}^{Q|Z}=P\left(Q\ge q_{}^{\ast }|Z\ge z_{}^{\ast }\right)=1-\frac{P\left(Z\ge z_{}^{\ast },Q\ge q_{}^{\ast }\right)}{P\left(Z\ge z_{}^{\ast }\right)}\\ \hspace{1em}\hspace{1em}=\frac{1-F_Z\left(z_{}^{\ast }|{\theta }_z^t\right)-F_Q\left(q_{}^{\ast }|{\theta }_q^t\right)-C\left(F_Z\left(z_{}^{\ast }|{\theta }_z^t\right),F_Q\left(q_{}^{\ast }|{\theta }_q^t\right)|{\theta }_{ZQ}^t\right)}{1-F_Z\left(z_{}^{\ast }|{\theta }_z^t\right)}\\ P_{}^{S|Z}=P\left(S\ge s_{}^{\ast }|Z\ge z_{}^{\ast }\right)=1-\frac{P\left(Z\ge z_{}^{\ast },S\ge s_{}^{\ast }\right)}{P\left(Z\ge z_{}^{\ast }\right)}\\ \hspace{1em}\hspace{1em}=\frac{1-F_Z\left(z_{}^{\ast }|{\theta }_z^t\right)-F_S\left(s_{}^{\ast }|{\theta }_s^t\right)-C\left(F_Z\left(z_{}^{\ast }|{\theta }_z^t\right),F_S\left(s_{}^{\ast }|{\theta }_s^t\right)|{\theta }_{ZS}^t\right)}{1-F_Z\left(z_{}^{\ast }|{\theta }_z^t\right)}\end{array}$$

(11)

4. Results

The temporal trends of peak discharge, peak water stage and suspended sediment load are investigated for selection of explanatory variables. Then, the marginal distributions of three flood variables are described by time-varying distributions, and time-varying copulas are applied to construct joint distributions of both Z-Q and Z-S. Finally, joint probability and conditional probability of bivariate flood variables are estimated to display time variation of their joint distributions over time. The details of the results are provided in the following sub-sections.

4.1. Temporal Trend Analysis

The nonstationarities of the flood series Q, Z and S at Waizhou station are examined by three trend analysis methods, including the MK test, Spearman test and Kendall test. Results of the trend analysis for Q, Z, and S are shown in

Table 2. Trend analysis of the three flood peak series at the Waizhou station during 1964–2013.

Series	Annual Mean		MK	Spearman	Kendall
	1964–2003	2004–2013
$Q$	12.01	10.47	−1.33	−0.16	−0.12
$Z$	23.25	20.96	−2.92 **	−0.41 **	−0.29 **
$S$	418.23	123.90	−5.33 **	−0.70 **	−0.52 **

Note: Q: peak discharge (10³ m³/s); Z: peak water stage (m); S: suspended sediment load (10³ kg/d). –: delineates negative trends; * and ** delineate significant trend at 0.05, 0.01 significance level, respectively.

. It is seen that the series of Z and S have undergone significant decreasing trends at the 1% significance level during 1964–2013. Peak discharge Q has presented some degree of negative trends, but cannot pass the trend tests at the 5% significance level. The change of forest cover does not have strong effects on flood peak discharge during large storm events, especially for large basins [50].

The suspended sediment load has undergone significant decreasing trend, particularly after 2003 because of projects implementation about afforestation and conservation measures [51]. The river elevation has undergone significant decrease due to river sand mining [31], which has intensified obviously since 2003. In order to display the sharp changes of Q, Z, and S in the period of last decades, the year of 2003 is chosen as a separated point from the inherent physical cause-effect connection. The annual mean values of Z and S are about 23.25 m and 418.22 × 10³ kg/day from 1964 to 2003, and then decrease to 20.96 m and 123.90 × 10³ kg/day for the last decade. Through the aforementioned analysis, it is reasonable to conclude that both Z and S display significant nonstationarity during the period from 1964 to 2013, whereas Q is stationary. The trend identification of Q, Z and S are consistent with some previous research conclusions [52,53].

4.2. Nonstationary Marginal Distributions

According to the cause-effect analysis, MCE and FCR could be two potential physical driving forces for the nonstationarity of Z and S. MCE series (Figure 2a) at Waizhou station have presented an obvious decreasing trend from 1964 to 2013, especially during last ten years. FCR series (Figure 2b) of Ganjiang River basin displays a slight decreasing trend before 1980s, but it expands rapidly from 1980s to 2010s due to artificial afforestation and forest conservation in Jiangxi province. Figure 2c,d illustrate significant positive correlation between Z and MCE, and significant negative correlation between S and FCR, respectively. Considering their inherent physical connection as well as statistical correlation, MCE and FCR are separately used as explanatory variables of the nonstationarity of Z and S.

Since annual peak discharge displays stationarity during 1964–2013, Q is fitted by five candidate distributions with all parameters treated as constants. In fitting Z and S series by nonstationary models, there are three main possible situations for each candidate distribution: (1) only the location parameter is time-varying, (2) only the scale parameter is time-varying and (3) both the location and scale parameters are time-varying. The final distribution model for each candidate is determined from the three models above by the selection criterion of AICc (i.e., model with the minimum AICc is the best). For the three flood variables Q, Z and S, the estimated parameters and the results of the GOF test of all candidate marginal distributions are summarized in

Table 3. Parameters and goodness-of-fit of the candidate marginal distributions fitted to Q, Z, and S at the Waizhou station during 1964–2013, respectively.

Variable	Distribution	Estimated Parameters	AICc	KS-Test
				Statistic	p-Value
Q	LNO	m = 9.309, σ = 0.346	970.57	0.083	0.881
	WEI	μ = 13090, σ = 3.154	973.81	0.101	0.683
	LOG	μ = 11451, σ = 2305	976.57	0.108	0.602
	GAM	μ = 11703, σ = 0.339	970.48	0.097	0.737
	PⅢ	μ = 11694, σ = 0.350, γ = 0.474	971.89	0.080	0.910
Z	LNO	μ = exp (1.081 + 0.004MCEt)	177.31	0.084	0.843
		σ = 0.059
	WEI	μ = exp (3.005 + 0.011MCEt)	181.46	0.101	0.646
		σ = exp (2.123 + 0.061MCEt)
	LOG	μ = 18.461 + 0.326MCEt	179.09	0.081	0.873
		σ = 0.782
	GAM	μ = exp (2.946 + 0.014MCEt)	177.24	0.089	0.793
		σ = 0.059
	PⅢ	μ = exp(2.936+0.014MCEt)	178.67	0.093	0.745
		Σ = 0.060, γ = 0.251
S	LNO	μ = exp (2.150 − 0.886FCRt)	646.40	0.097	0.701
		σ = exp (-1.332 + 1.378FCRt)
	WEI	μ = exp (8.072 − 4.544FCRt)	648.35	0.136	0.285
		σ = exp (1.476 − 1.448FCRt)
	LOG	μ=950.510 − 12.693FCRt	657.23	0.131	0.332
		σ = exp (6.137 − 3.412FCRt)
	GAM	μ = exp (7.943 − 4.530FCRt)	645.95	0.113	0.515
		σ = exp (−1.336 + 1.291FCRt)
	PⅢ	μ = exp (8.098 − 4.854FCRt)	646.05	0.084	0.845
		σ = 0.537, γ = 0.702

Note: LNO, WEI, LOG, GAM and PⅢ are the abbreviations of Lognormal, Weibull, Logistic, Gamma and Pearson type III distribution, respectively. MCE and FCR stand for the main channel mean elevation (m) and forest coverage rate (%), respectively. $\mu$, $\sigma$ and $\nu$ represent location, scale and shape parameters of marginal distribution, respectively. The best appropriate distribution marked with bold fonts.

. The P-value of the KS test was simulated using the Monte Carlo method. All the distributions pass the KS test at the 0.01 significance level. According to AICc, the most appropriate marginal distributions of Q, Z and S are Gamma distribution. The location parameter $\mu$ of the Gamma distribution for describing Z is positively related to MCE, whereas the scale parameter $\sigma$ is constant. Meanwhile, both location and scale parameters of the Gamma distribution for describing S are positively related to FCR. The QQ plot (Figure 3a) of the best fitted distribution indicates that these selected distributions have a quite good fitting quality.

Due to the balance of riverbed erosion and deposition during 1964–1994, the mean value of Z at Waizhou station is about 23.26 m. In the period of 1995–2013, the value of MCE has underwent a sharp decrease from 12.74 to 5.90 m, thus the variable value for Z has dropped obviously as well. Similarly, the mean and coefficient of variation (Cv) of the suspended sediment load at Waizhou station are 469.81 × 10³ kg/day and 0.47 during 1964–1994. After 1995, the mean and Cv of S are getting smaller significantly, because massive afforestation activities have been implemented in the study basin with growth rate of forest cover at 1.17% year by year.

4.3. Nonstationary Dependence of Bivariate Flood Variables

Bivariate copulas under stationary and nonstationary conditions are constructed based on the estimated marginal distributions. In modelling time-varying copula, selection of the explanatory variables (i.e., MCE and FCR) is determined by the two corresponding marginal variables. In detail, about the link function of copula parameter, MCE is selected as covariate expressed in Equation (6) for Z-Q, while MCE and/or FCR are chosen as covariates described in Equations (7) and (8) for Z-S.

The results of the estimated parameters and GOF are summarized in

Table 4. Parameters and goodness-of-fit of the candidate bivariate copulas fitted to Z-Q and Z-S at the Waizhou station during 1964–2013, respectively.

	Copula	Parameter (θ)	AICc RMSE		CM-Test
					Statistic	p-Value
Z-Q	Clayton	5.702	−103.91	0.031	0.053	0.457
		exp (1.811 − 0.005MCEt)	−101.93	0.031	0.053	0.458
	GH	4.008	−94.76	0.039	0.056	0.426
		exp (1.743 − 0.026MCEt)	−93.33	0.039	0.055	0.433
	Frank	17.683	−102.79	0.035	0.055	0.434
		41.713 − 1.747MCEt	−104.25	0.035	0.054	0.448
Z-S	Clayton	1.250	−22.70	0.040	0.048	0.538
		exp (−1.760 + 0.063FCRt)	−25.11	0.039	0.038	0.716
		exp (1.739 − 0.114MCEt)	−24.01	0.040	0.044	0.596
		exp (−5.025 + 0.118FCRt − 0.113MCEt)	−23.48	0.040	0.052	0.479
	GH	1.796	−28.04	0.038	0.050	0.506
		exp (0.350 + 0.007FCRt)	−26.25	0.039	0.047	0.541
		exp (1.195 − 0.045MCEt)	−27.40	0.038	0.047	0.556
		exp (3.563 − 0.039FCRt − 0.128MCEt)	−26.76	0.039	0.057	0.410
	Frank	5.319	−28.77	0.032	0.035	0.759
		−1.528 + 0.226FCRt	−28.79	0.032	0.029	0.855
		16.169 − 0.782MCEt	−30.29	0.031	0.028	0.870
		24.186 − 0.132FCRt − 1.072MCEt	−28.45	0.032	0.030	0.840

. The P-value of the CM test is simulated using the Monte Carlo method. All the applied copula functions pass the CM test at the significance level of 0.01. Then, RMSE and AICc, which claim that the model with smaller values is the better, are used as selected criteria [54]. Performances of the candidate copulas are not different obviously from RMSE because of a little difference between them. But the optimum time-varying copulas perform better than stationary ones for Z-Q and Z-S in terms of AICc. Frank is found to be the more appropriate bivariate copula for both Z-Q and Z-S, and parameters ${\theta }_{ZQ}^t,{\theta }_{ZS}^t$ expressed in Equations (12) and (13), respectively, as below:

$${\theta }_{ZQ}^t=41.713-1.747MC{E}^t$$

(12)

$${\theta }_{ZS}^t=16.169-0.782MC{E}^t$$

(13)

Figure 3b shows the QQ plot of the two bivariate copulas above. It displays a good agreement between empirical distribution and theoretical distribution. In addition, graphical GOF of the selected Frank copulas for both Z-Q and Z-S are shown in Figure 4. The simulation series of Q, Z, and S, which are 30 scatters per year from 1964 to 2013, are generated by Monte Carlo method. Data are transformed to the real space by use of the corresponding marginal distributions. This graphical GOF also demonstrates that the selected time-varying copula functions have satisfactory fitting performances.

Therefore, evolutions of the correlation parameters (i.e., ${\theta }_{ZQ}^t$, ${\theta }_{ZS}^t$) present the overall upward trend significantly from 1964 to 2013, which can demonstrate that the dependence structures of Z-Q and Z-S are nonstationary. These results demonstrate that riverbed down-cutting instead of forest coverage is the main external effect for both Z-Q and Z-S at Waizhou station. It is possible to analyze the multivariate flood frequency in the future through the prediction of explanatory variables at Waizhou station based on their relationships [55,56]. Furthermore, it can give supporting information for the flood control design of hydraulic structures under the changing environments [14,39].

4.4. Temporal Variation in Joint and Conditional Probabilities

Contours for various joint probabilities denoted in Figure 5 are calculated by using the selected copulas with parameters in the years of 1970, 1990 and 2010. With the same probability, Z maintains generally higher and almost equivalent value before 1990, whereas moves downward greatly from 1990 to 2010 due to the decrease in mean of the peak water stage. Similarly, the value of S falls rapidly in last ten years due to the significant decrease in mean and Cv of S after 1995. For example, the contours for 0.2 in 2010 for both Z and S are lower than lines for 0.7 in 1970 and 1990, which indicates that the nonstationarity play a great influence on joint probability in the last decade. Upon closer inspection, due to the strengthening in dependence for both Z-Q and Z-S, the corner of contour with the same joint probability becomes angular, especially for the lower probability.

As shown in Figure 6, the temporal variations of joint probability as Equation (9) are assessed for each time step from 1964 to 2013. Meanwhile, the magenta dots on the probability contours represent the corresponding the most likely combination derived from Equation (10). The probability contours cover a broad range, with marginal values ranging from 22.12 m to 25.06 m for Z and 721.78 × 10³ kg/day to 151.58 × 10³ kg/day for S. It can be seen that the suspended sediment load is more sensitive to flood event than the peak water stage from comparison of their value ranges. For the most likely combination, the peak discharge at Waizhou station in the Ganjiang River displays a slight increasing trend from 15.97 × 10³ m³/s to 16.41 × 10³ m³/s over time, while the values of peak water stage and suspended sediment load have obviously decreased, especially in the last decade. As shown in Figure 6b,d, it is indicated that Z and S in the most likely combination are obviously affected by riverbed down-cutting and the change of forest cover, respectively. Moreover, since the dependence structure of Z-Q is nonstationary as Equation (12), Q in the most likely combination is impacted by riverbed down-cutting as well.

The conditional probability as Equation (11) is calculated by the best fitted copula function for both Z-Q and Z-S. In practice, flood control planners and managers usually focus on the water stage rather than discharge and suspended sediment load in severe flood flow. Given the warning water stage (23.5 m) at Waizhou station, the time variations of conditional probability ($P^{Q|Z},P^{S|Z}$) from 1964 to 2013 are presented in Figure 7. The result can be of great important and noteworthy for spillway design and flood control. It can be seen that $P^{Q|Z}$ for the fixed discharge has increased over time, while $P^{S|Z}$ has decreased during last fifty years. That’s to say, even if Z maintains the same value at Waizhou station, the fact should attract more attention that probabilities of inundation and riverbed erosion are getting higher from 1964 to 2013.

What is more, the proposed method also allows us to obtain more information concerning the conditional probabilities under various given water stages. Figure 8 has shown the results, where MCE and FCR are assumed to be 6 m and 70%, respectively. The conditional probabilities become bigger for both Q and S with the higher peak water stage. When values of Q and S become bigger, the differences of conditional probabilities are getting smaller under various given water stages, especially for $P^{S|Z}$ (Figure 8b). In addition, the conditional probability $P^{Q|Z}$ is more sensitive than $P^{S|Z}$ on condition of various water stages at Waizhou station. Those results can give us quantitative information about temporal variation of flood variables under the changing environments in the Ganjiang River basin.

5. Conclusions

This study employs time-varying copula model to investigate the evolution of the relationships of Z-Q and Z-S depending on main channel elevation (MCE) and forest coverage ratio (FCR) in the Ganjiang River basin during 1964–2013. The main conclusions are presented as follows:

• It is obvious that both the mean and variance of S have significantly decreased, while only the mean has reduced for Z, particularly in the recent decades. Furthermore, Gamma distribution with location parameter expressed as a function of MCE is best fitted distribution for Z, and Gamma with parameters of location and scale expressed as functions of FCA is for S, while the best fitted distribution of Q is the Gamma with constant parameters.
• It is found that the most fitted bivariate copulas for both Z-Q and Z-S are Frank copula, the parameters of which are expressed as the function of MCE. Therefore, riverbed down-cutting at Waizhou station plays the dominant role in strengthening dependences of both Z-Q and Z-S from 1964 to 2013.
• The results of joint probability and conditional probability show that the corner of contour lines enhanced more greatly due to the strengthening dependences over time, especially for the lower probability. In addition, it can be seen that values of Z and S fall rapidly in the last ten years due to the decreasing mean of these two variables.