Evaluation of Return Period and Risk in Bivariate Non-Stationary Flood Frequency Analysis

: The concept of a traditional return period has long been questioned in non-stationary studies, and the risk of failure was recommended to evaluate the design events in ﬂood modeling. However, few studies have been done in terms of multivariate cases. To investigate the impact of non-stationarity on the streamﬂow series, the Yichang station in the Yangtze River was taken as a case study. A time varying copula model was constructed for bivariate modeling of ﬂood peak and 7-day ﬂood volume, and the non-stationary return period and risk of failure were applied to compare the results between stationary and non-stationary models. The results demonstrated that the streamﬂow series at the Yichang station showed signiﬁcant non-stationary properties. The ﬂood peak and volume series presented decreasing trends in their location parameters and the dependence structure between them also weakened over time. The conclusions of the bivariate non-stationary return period and risk of failure were different depending on the design ﬂood event. In the event that both ﬂood peak and volume are exceeding, the ﬂood risk is smaller with the non-stationary model, which is a joint effect of the time varying marginal distribution and copula function. While in the event that either ﬂood peak or volume exceed, the effect of non-stationary properties is almost negligible. As for the design values, the non-stationary model is characterized by a higher ﬂood peak and lower ﬂood volume. These conclusions may be helpful in long-term decision making in the Yangtze River basin under non-stationary conditions.


Introduction
A flood is a natural hazard that may cause great damages to human civilization, and it is studied worldwide using hydrologic as well as hydraulic models [1,2].Flood frequency analysis, a commonly used method in hydrological modeling, has been applied in the calculation of various design values, such as dams, bridges, culverts and many other hydraulic constructions [3,4].However, the flood frequency analysis is based on some assumptions in terms of statistics, one of which is the well-known stationary assumption.As climate change and human activities intensify, the stationary assumption has faced more and more challenges [5].Much research has been done to study the non-stationarity in a hydrological series.For example, Zhang et al. [6] used three methods to investigate the effect of dams and lakes on the streamflow series in the Yangtze River basin.Ishak et al. [7] used a regional Mann-Kendall test to explore the non-stationary properties in annual maximum flood series for 491 catchments in Australia.O'Brien and Burn [8] used a non-stationary index-flood technique to estimate the extreme quantiles of annual maximum streamflow.Most of the studies supported the proposition by Milly et al. [5] and concluded that a non-stationary modeling framework may provide valuable information for policy implementation and hazard prevention.
The return period is an important indicator of risk for a design flood event; however, in the context of non-stationary, it may conceal the actual meaning of underlying mechanism [9] and even lead to misconceptions.To tackle such problems, its concept has been reformulated and extended to accommodate the non-stationary conditions.Salas and Obeysekera [10] presented a review of how to extend the concept of the return period into a non-stationary framework.Furthermore, Serinaldi et al. [11] suggested that the exceeding probability and risk of failure can provide more coherent tools for risk assessment.The routinely used methods were revived in terms of non-stationary conditions and it has been widely accepted in the latest literature.For example, Du et al. [12] used the non-stationary return period by employing meteorological covariates to investigate the low flow in the Wei River.Yan et al. [13] compared four methods of non-stationary hydrologic design strategies and highlighted the significance of considering a project's design life period.
Until now, most of the non-stationary studies were focused on univariate cases; however, the dependence structure between different hydrological variables may also have changed under the changing environment.Jiang et al. [14] proposed a time varying copula model to explore how the dependence structure between river flows changed with the construction of upstream reservoirs on the Hanjiang River.Furthermore, Ahn and Palmer [15] used the non-stationary copula to evaluate the non-stationarity in low flow characteristics and to forecast bivariate frequency in the future.Sarhadi et al. [16] developed a dynamic Bayesian copula to model the time varying dependence structure between drought characteristics under climate change in Tehran, Iran.Yet, few studies have been done on the evaluation of the return period and risk under non-stationary conditions in multivariate cases.Hence, there are two objectives for this article-the first is to evaluate the effect of the non-stationarity on marginal distributions of flood characteristics and the dependence structure between them in the Yangtze River, and the second is to investigate how these effects change the design values of flood events when evaluated by bivariate non-stationary return periods and risk of failure.

Study Area and Data
The Yangtze River, 6300 km long, originates from the Tibetan Plateau and ends in the Chongming Island in Shanghai (shown in Figure 1).The minimum annual rainfall of this basin is about 270 mm in the west and the maximum is about 1900 mm in the Southeast [17].Because of the monsoon climate, more than 70% of its annual precipitation occurs in summer.The Yangtze River basin is prone to huge flood events and it has caused great damages to social prosperity and loss of life, and the flood in 1998 led to an economic loss of more than $20 billion [18].To this end, many measures have been taken in recent years.The Three Gorges Project was built in the Yangtze River and it has been completed in 2003, the main function of which is flood control.Besides, the streamflow is also sensitive to climate changes [19], which makes it more necessary and important to reassess the flood conditions and forecast future risks in this region.
According to its geographical and climatic characteristics, the Yangtze River is subdivided into three sections, where Yichang is the dividing point between the upper and middle reaches.The Yichang station, located 44 km downstream of the Three Gorges Reservoir, is an important monitoring station for hydrological regimes in the upper Yangtze River.In this study, daily averaged streamflow records at the Yichang station from 1949 to 2012 were collected for research, and the flood peak and 7-day flood volume series were extracted from the data using the block maximum approach.The flood peak (P), is defined as the discharge of the annual maximum flood; the 7-day flood volume (V), is defined as the integral of the annual maximum 7-day flood and time, it does not necessarily include the day of the annual maximum flood.There are no missing values in the records for this period.According to its geographical and climatic characteristics, the Yangtze River is subdivided into three sections, where Yichang is the dividing point between the upper and middle reaches.The Yichang station, located 44 km downstream of the Three Gorges Reservoir, is an important monitoring station for hydrological regimes in the upper Yangtze River.In this study, daily averaged streamflow records at the Yichang station from 1949 to 2012 were collected for research, and the flood peak and 7-day flood volume series were extracted from the data using the block maximum approach.The flood peak (P), is defined as the discharge of the annual maximum flood; the 7-day flood volume (V), is defined as the integral of the annual maximum 7-day flood and time, it does not necessarily include the day of the annual maximum flood.There are no missing values in the records for this period.

Univariate Trend Test
The Mann-Kendall (M-K) trend test [20,21] is a rank-based method which is widely adopted in hydro-meteorological studies for detecting the trend in a time series.For a hydrological variable x with n years of records, the test is based on the statistic value S defined by the following formula:

Univariate Trend Test
The Mann-Kendall (M-K) trend test [20,21] is a rank-based method which is widely adopted in hydro-meteorological studies for detecting the trend in a time series.For a hydrological variable x with n years of records, the test is based on the statistic value S defined by the following formula: Under a null hypothesis of no trend, the statistical S is approximately normally distributed with zero mean and its variance: The flood peak and volume series were examined with a two-sided test given a 0.05 significance level in this study.

Multivariate Trend Test
The Covariance Sum Test (CST), introduced by Hirsch and Slack [22], was adopted as an extension of the univariate M-K test for detecting trends in multivariate correlated data in this study.In multivariate cases, the univariate statistic S (u) (u = 1, 2, . . ., d) for each component of the variables is calculated first.S = (S (1) , S (2) , . . ., S (d) ) is a d-dimensional vector with a covariance matrix C M = (c u,v ) u,v=1, . . . . . .d , in which c u,v = cov(M (u) ,M (v) ).The covariate terms are calculated with the following equations [23]: Water 2019, 11, 79 4 of 14 where The multivariate M-K test is based on the statistic S m , given as: Under a null hypothesis, the statistic S m is asymptotically normal with zero mean and its variance is given as: For a given significance level, the application of the CST is similar to that of the univariate M-K test.

Time Varying Copula
Under a changing environment, both the marginal distribution of hydrological variables and the dependence structure between them might be changing and the non-stationary copula is a perfect solution to this problem.In bivariate cases, based on the conventional copula theory, the time varying copula is defined as follows [14]: where F 1 ( ) and F 2 ( ) are marginal distributions for the variables Y 1 and Y 2 , C( ) is the copula function; θ t 1 and θ t 2 are parameters for the marginal distributions, and θ t c is the parameter for the copula function, and the superscript t means the parameters are time varying.The procedure for how to construct a time varying copula model is presented in Figure 2.

Marginal Distribution
To investigate the flood characteristics, various probability distributions have been recommended and studied in flood modeling researches [24].Six widely used distributions, including Gamma, Weibull, Gumbel, Logistic, Lognormal and Pearson type III distributions, were selected as candidates in this study for flood peak and 7-day flood volume series.Parameters of the distributions were estimated by the Generalized Additive Models in Location, Scale and Shape (GAMLSS) [25], which has been widely used in non-stationary modeling recently [26,27].In a GAMLSS framework, the observation of a variable y is assumed to follow a distribution with a probability density function (PDF) f (y|θ), and the parameter vector θ T is described as a function of explanatory variables and random effects.In case there are no additive terms, the distribution parameters can be denoted by a monotonic link function g k ( ) as: where X k denote the explanatory variables, β k are polynomial coefficients, p is the number of parameters, q is the degree of the polynomial, and q = 1 in this study to avoid over-parameterization.The location and scale parameters were related to the time covariate, including either one of them is time varying or both of them as time varying, and the shape parameter was treated as a constant.The stationary model was obtained by assuming that the parameters are independent of the explanatory variables.The optimal distribution was determined by the Akaike Information Criterion (AIC) [28] and the goodness-of-fit test was performed by the worm plot [29].PDFs of the marginal distributions and their link functions are listed in Table 1.

Copula Function
Among all copula families, the Archimedean copula is the most commonly used in hydrological modeling for its simple structure and good applicability [30].In this study, four Archimedean copulas, namely the Clayton, Gumbel-Hougaard (GH), Frank, and Joe copulas, were adopted for the bivariate modeling of flood peak and volume series.
Similar to the marginal distributions, the parameter θ of the copula is linked to the time covariate as follows: where β 0 and β 1 are polynomial coefficients and t is the time.A brief description of these copulas and their link functions are presented in Table 2.
Water 2019, 11, x FOR PEER REVIEW 5 of 15

Marginal Distribution
To investigate the flood characteristics, various probability distributions have been recommended and studied in flood modeling researches [24].Six widely used distributions, including Gamma, Weibull, Gumbel, Logistic, Lognormal and Pearson type III distributions, were Note: µ is the location parameter, σ is the scale parameter and ε is the shape parameter.
Table 2.The description of the Archimedean copulas and their link functions.
Generally, there are three combinations of marginal distributions and copula functions for a time varying copula model, and the final model is determined according to the trend test results.The temporal variations among different copula models are captured by comparing the AIC values and the stationary model is taken as a benchmark.The Inference Function for Margins (IFM) method [31] was applied for the estimation of parameters.

Non-Stationary Return Period
The return period is essential in flood frequency analysis and usually used in risk quantification and evaluation [32].In general, there are two interpretations of the return period, one of which is the expected waiting time (EWT) until the first event occurs.For a known design flood event (x = 0), X is defined as the waiting time until the event reoccurs for the first time, then the probability for X = x can be expressed using the following equation [10]: where P i is the exceeding probability of a design event for each year.The EWT return period T of the design event can thus be expressed as: Water 2019, 11, 79 7 of 14 Another interpretation of the return period is that the expected number of events (ENE) in T years is one [33].For a period of T years, N is the number of the design event, the expected value of which is calculated by: where I is an indicator variable.The ENE return period T is obtained when the solution to Equation ( 13) is equal to 1, as follows: In terms of bivariate frequency analysis, the design event is always divided into two cases, namely case AND (Y 1 ≥ y 1 and Y 2 ≥ y 2 ) and case OR (Y 1 ≥ y 1 or Y 2 ≥ y 2 ), which can be defined by two exceeding probabilities using copula theory as follows: Under the stationary assumption, the exceeding probabilities P i for each year are exactly the same so that both the two interpretations lead to the same result, T = 1/P i , which can be divided into two cases, specifically, using the following equations [34]: It is obvious that the stationary return period defined in Equations ( 17) and ( 18) do not provide more information than only the exceeding probabilities.Under non-stationary assumptions, since the exceeding probabilities for each year are different, Equations ( 17) and (18) are no longer valid, and the EWT and ENE definitions can lead to different return periods.However, the definition of EWT includes calculation of the exceeding probabilities for an infinite period, which is hard to predict and has a risk of not converging.Besides, the rate of convergence of EWT may be influenced by the selected distribution to some extent [13].Thus, the EWT definition is not taken for comparison in this study.

Risk of Failure and Reliability
The risk of failure, defined as the probability that a critical event occurs at least once in its design life, has become increasingly popular in non-stationary risk analysis.It provides a different perspective for evaluating flood risks from the return period.For n years of design life, the probability R that at least a failure was observed with exceeding probabilities P i in the ith year is given as: where P i can be either P AND or P OR given above in bivariate cases.Correspondingly, the reliability, r, is defined as the probability that no design flood event will happen in the next n years, given as: Water 2019, 11, 79 8 of 14

The Trend Test
The flood peak P and volume V series were extracted from the data and examined first.Given α = 0.05, the threshold values for significant trends are ±1.96.According to the univariate M-K test, the statistic result for the flood peak is −2.06 and it is −2.05 for the flood volume, indicating a significant decreasing trend for both of the two series.Moreover, according to the multivariate M-K test, the statistic result is −2.40, indicating a significant bivariate decreasing trend for the two variables.The results demonstrated that the stationary assumption may no longer be valid for the streamflow series at the Yichang station-not only the time varying marginal distributions, but also a time varying copula is needed in the flood frequency analysis in this region.Thus, we have two models for research in the subsequent sections, which are: Model 1, the stationary model with constant parameters; and model 2, the non-stationary model with both the marginal distribution and copula function parameters varying with time.

Marginal Distribution
The six distributions were tested to fit the marginal distributions of the flood peak and volume series, and the best-fitted one was selected according to minimum AIC values.For the flood peak series, the Weibull was the best-fitted distribution under the stationary assumption and the Gamma was the best-fitted distribution under the non-stationary assumption.A decreasing trend was detected in the location parameter while the scale parameter showed an increasing trend.On the other hand, the Logistic distribution was the most appropriate distribution for the flood volume series under both stationary and non-stationary assumptions, and a decreasing trend was detected in the location parameter while the scale parameter is constant.The optimal distributions and their parameters are presented in Table 3.The time varying distributions perform better than the stationary ones according to AIC values.The worm plots for time varying marginal distributions were presented in Figure 3.All points are within the 95% confidence interval and the selected distributions have enough fitting quality.

Copula Modeling
The dependence structure between flood peak and volume series was explored first in the multivariate frequency analysis.Using the three widely used relationship coefficients, the Pearson's ρ is 0.90, Kendall's tau is 0.75, and the Spearman's rho is 0.91, indicating that the two variables are

Copula Modeling
The dependence structure between flood peak and volume series was explored first in the multivariate frequency analysis.Using the three widely used relationship coefficients, the Pearson's ρ is 0.90, Kendall's tau is 0.75, and the Spearman's rho is 0.91, indicating that the two variables are strongly correlated.Table 4 summarizes the parameters and goodness-of-fit results for the two models, and it can be found that all the copula functions of model 2 detected a decreasing trend of θ, which means the correlation relationship between the two variables weakened over time.Almost all the AIC values of non-stationary models are smaller than the stationary model when using the same copula function, except for the Clayton copula.The Frank copula of model 2 is the best fit for bivariate modeling of flood peak and volume among all scenarios.The relationship between empirical and theoretical probabilities for the Frank copula under stationary and non-stationary conditions was plotted in Figure 4.All points are located near the 1:1 line, indicating satisfactory fitting performance.

The Non-Stationary Return Period
Based on the selected copula and estimated parameters for each model, the bivariate nonstationary return periods can be estimated by substituting Equations ( 15) and ( 16) to Equation ( 14).It should be noted that the non-stationary return period was defined for each year rather than the whole series.Following the Chinese flood guidelines [35], the unique combination of design flood peak and volume is determined by setting them to be of the same univariate frequency, which was a uniform procedure for similar calculations of design values in past studies [36].Specifically, for the year 2012, the univariate exceeding probability against the bivariate non-stationary return period was presented in Figure 5.In case AND from Figure 5a, for a given return period, the univariate exceeding probability of model 2 is bigger than that of model 1, implying that the non-stationary model is safer than the stationary model under the same design standard.On the other hand, the non-stationarities  Based on the selected copula and estimated parameters for each model, the bivariate non-stationary return periods can be estimated by substituting Equations ( 15) and ( 16) to Equation (14).It should be noted that the non-stationary return period was defined for each year rather than the whole series.Following the Chinese flood guidelines [35], the unique combination of design flood peak and volume is determined by setting them to be of the same univariate frequency, which was a uniform procedure for similar calculations of design values in past studies [36].Specifically, for the Water 2019, 11, 79 10 of 14 year 2012, the univariate exceeding probability against the bivariate non-stationary return period was presented in Figure 5.In case AND from Figure 5a, for a given return period, the univariate exceeding probability of model 2 is bigger than that of model 1, implying that the non-stationary model is safer than the stationary model under the same design standard.On the other hand, the non-stationarities have a very small influence on the return period of case OR, as shown in Figure 5b.The result was further evaluated by the design values.For example, the univariate exceeding probabilities under 10 years of the bivariate return period in case OR are both 0.065 for the two models, while the designed flood peaks (m 3 /s) are 62,202 and 63,127, respectively, and the designed flood volumes (10 6 m 3 /s) are 33,054 and 31,165, respectively.Model 2 is characterized by a higher flood peak and lower flood volume compared with model 1, which may be attributed to the time varying marginal distributions.

Risk of Failure and Reliability
Given the initial flood peak (63,059 m 3 /s) and volume (33,738 × 10 6 m 3 /s) with 20 years of univariate return period under the stationary assumption, the risk of failure of the design flood events for model 1 and model 2 in the next 50 years are calculated and plotted in Figure 6. Figure 6a shows the risk of failure for case AND.Under the stationary assumption, the risk of failure R increased rapidly with the design life n, as shown by the red line.However, due to non-stationarities, the flood risk is smaller than the stationary model, and the rate of growth decreased as presented by the blue line.On the other hand, Figure 6b shows the risk of failure for case OR, the conclusions of which are completely different from that of case AND.In total, the R values of case OR are nearly the same between model 1 and model 2, indicating that the non-stationary properties have a very small influence on the risk of failure for the flood event case OR. Figure 7 presents the results of reliability for the two models, which are similar to those of Figure 6.The results suggested that the hydrologic facilities controlled by design flood event case OR are more stable under the changing environment in the study region, but the risks are much bigger than that of case AND.

Risk of Failure and Reliability
Given the initial flood peak (63,059 m 3 /s) and volume (33,738 × 10 6 m 3 /s) with 20 years of univariate return period under the stationary assumption, the risk of failure of the design flood events for model 1 and model 2 in the next 50 years are calculated and plotted in Figure 6. Figure 6a shows the risk of failure for case AND.Under the stationary assumption, the risk of failure R increased rapidly with the design life n, as shown by the red line.However, due to non-stationarities, the flood risk is smaller than the stationary model, and the rate of growth decreased as presented by the blue line.On the other hand, Figure 6b shows the risk of failure for case OR, the conclusions of which are completely different from that of case AND.In total, the R values of case OR are nearly the same between model 1 and model 2, indicating that the non-stationary properties have a very small influence on the risk of failure for the flood event case OR. Figure 7 presents the results of reliability for the two models, which are similar to those of Figure 6.The results suggested that the hydrologic facilities controlled by design flood event case OR are more stable under the changing environment in the study region, but the risks are much bigger than that of case AND.
completely different from that of case AND.In total, the R values of case OR are nearly the same between model 1 and model 2, indicating that the non-stationary properties have a very small influence on the risk of failure for the flood event case OR. Figure 7 presents the results of reliability for the two models, which are similar to those of Figure 6.The results suggested that the hydrologic facilities controlled by design flood event case OR are more stable under the changing environment in the study region, but the risks are much bigger than that of case AND.

Discussion
The multivariate analysis is more complicated than the univariate analysis in terms of different combinations of hydrological variables and various design events.In past literature on time varying copulas, the assessment of dynamic risks is a difficult task.Jiang et al. [14] used the joint return period for bivariate non-stationary modeling of low flows in the Hanjiang River, but the definition is similar to that under the stationary assumption.Zhang et al. [37] abandoned the return period and used the exceeding probabilities instead for non-stationary flood frequency analysis in the Wangkuai Reservoir catchment.The bivariate non-stationary return period and risk of failure used in this study provided an explicit measure of risk with a robust theoretical basis, thus filling up the current knowledge gap.The performances of case AND and case OR can be completely different in this study, and more attention can be paid to different kinds of return periods for multivariate flood events, such as Kendall's return period, the structure-based return period, and so on.However, we should also keep in mind that T AND and T OR and any other kind of return period are only different descriptions of the underlying mechanisms of design structures and must be used accordingly.It is important to figure out what type of event is used in the practical application and how to evaluate the effect of non-stationarity on flood risk under this specific condition.

Conclusions
In summary, this paper built a non-stationary copula model for bivariate flood frequency analysis in the Yangtze River basin, where the streamflow records from 1949 to 2012 at the Yichang station were used as a practical application.Flood peak is defined as the annual maximum flooding,

Discussion
The multivariate analysis is more complicated than the univariate analysis in terms of different combinations of hydrological variables and various design events.In past literature on time varying copulas, the assessment of dynamic risks is a difficult task.Jiang et al. [14] used the joint return period for bivariate non-stationary modeling of low flows in the Hanjiang River, but the definition is similar to that under the stationary assumption.Zhang et al. [37] abandoned the return period and used the exceeding probabilities instead for non-stationary flood frequency analysis in the Wangkuai Reservoir catchment.The bivariate non-stationary return period and risk of failure used in this study provided an explicit measure of risk with a robust theoretical basis, thus filling up the current knowledge gap.The performances of case AND and case OR can be completely different in this study, and more attention can be paid to different kinds of return periods for multivariate flood events, such as Kendall's return period, the structure-based return period, and so on.However, we should also keep in mind that T AND and T OR and any other kind of return period are only different descriptions of the underlying mechanisms of design structures and must be used accordingly.It is important to figure out what type of event is used in the practical application and how to evaluate the effect of non-stationarity on flood risk under this specific condition.

Conclusions
In summary, this paper built a non-stationary copula model for bivariate flood frequency analysis in the Yangtze River basin, where the streamflow records from 1949 to 2012 at the Yichang station were used as a practical application.Flood peak is defined as the annual maximum flooding, while 7-day flood volume is the integral of annual maximum 7-day flood and time.The non-stationary return period and risk of failure were evaluated for flood risk assessment.The main conclusions are listed as follows.
(1) Both univariate and multivariate M-K trend tests were used to examine the temporal variations between flood characteristics, and the results revealed severe non-stationary properties, indicating that the stationary assumption is invalid, and a non-stationary copula model is needed for hydrological modeling in this region.(2) The flood peak and volume series were fitted by 6 marginal distributions under stationary and non-stationary situations, and the appropriate distributions vary correspondingly.The Weibull and Logistic distributions were selected for flood peak and volume under the stationary assumption, while the Gamma and Logistic distributions were selected under the non-stationary assumption.The flood peak and volume series presented decreasing trends in their location parameters.Four copula functions were applied to investigate the dependence structure between flood variables, and the Frank copula was selected as the most appropriate.The non-stationary model performs better than the copula with constant parameters according to AIC values.The copula parameters of the non-stationary model detected a decreasing trend, which means the dependence structure between flood variables also weakened over time.(3) The joint non-stationary return period was calculated for comparison, and the results varied with different design flood events.Compared with the stationary model, in case AND, the flood risk of non-stationary models decreased as a result of the decreased parameters of both marginal distributions and copula functions.While in case OR, the effect of the non-stationary properties is almost negligible.As for the design values, the non-stationary model is characterized by a higher flood peak and lower flood volume.The bivariate risk of failure of 50 years of design life was also estimated, and the conclusions were similar to those obtained by the return period.The effect of non-stationarity is more inclined to case AND rather than case OR, and the non-stationary model is safer than the stationary model in case AND.
The streamflow series at the Yichang station showed significant non-stationary properties according to our research, and the non-stationary copula model can better describe the temporal variations.The bivariate non-stationary return period and risk of failure may throw light on the extension of conventional concepts and serve as valuable tools for flood risk assessment in a changing environment.Time was used as a covariate in this study for both the marginal distribution and copula function.However, it is not recommended in non-stationary modeling because of its lack of physical meanings.Moreover, the time varying trend obtained through past data may not be consistent after 50 or 100 years, and the reliability of predicted results needs to be carefully considered.Some physical-based variables, such as the climate index, can be taken into account in future studies.

Figure 1 .
Figure 1.A map of the Yangtze River basin.

Figure 1 .
Figure 1.A map of the Yangtze River basin.

Figure 2 .
Figure 2. The procedure for constructing the time varying copula model.

Figure 2 .
Figure 2. The procedure for constructing the time varying copula model.

Figure 3 .
Figure 3. Worm plots for the time varying marginal distribution of (a) the flood peak series and (b) the flood volume series.

Figure 3 .
Figure 3. Worm plots for the time varying marginal distribution of (a) the flood peak series and (b) the flood volume series.

Figure 4 .
Figure 4.A comparison of empirical and theoretical probabilities for the Frank copula (a) under the stationary condition and (b) under non-stationary condition.

Figure 4 .
Figure 4.A comparison of empirical and theoretical probabilities for the Frank copula (a) under the stationary condition and (b) under non-stationary condition.

Water 2019 , 15 Figure 5 .
Figure 5.The univariate exceeding probability against the non-stationary return period for (a) case AND and (b) case OR.

Figure 5 .
Figure 5.The univariate exceeding probability against the non-stationary return period for (a) case AND and (b) case OR.

Figure 6 .
Figure 6.Bivariate risk of failure as a function of design life for (a) case AND and (b) case OR.Figure 6. Bivariate risk of failure as a function of design life for (a) case AND and (b) case OR.

Figure 6 . 15 Figure 7 .
Figure 6.Bivariate risk of failure as a function of design life for (a) case AND and (b) case OR.Figure 6. Bivariate risk of failure as a function of design life for (a) case AND and (b) case OR.Water 2019, 11, x FOR PEER REVIEW 12 of 15

Figure 7 .
Figure 7. Reliability as a function of design life for (a) case AND, and (b) case OR.

Table 1 .
Probability density functions of the marginal distributions and their link functions.

Table 3 .
The optimal marginal distributions and their parameters.
Water 2019, 11, x FOR PEER REVIEW 9 of 15

Table 3 .
The optimal marginal distributions and their parameters.

Table 4 .
The parameters and goodness-of-fit tests for candidate copulas.