Research article 22 Mar 2019
Research article  22 Mar 2019
Multivariate hydrologic design methods under nonstationary conditions and application to engineering practice
 ^{1}School of Environmental Studies, China University of Geosciences, Wuhan 430074, China
 ^{2}State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China
 ^{3}College of Water Conservancy and Hydropower, Hebei University of Engineering, Handan 056002, China
 ^{4}Guangxi Water Resources Management Center, Nanning 530023, China
 ^{5}Department of Geosciences, University of Oslo, P.O. Box 1047 Blindern, 0316 Oslo, Norway
 ^{1}School of Environmental Studies, China University of Geosciences, Wuhan 430074, China
 ^{2}State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China
 ^{3}College of Water Conservancy and Hydropower, Hebei University of Engineering, Handan 056002, China
 ^{4}Guangxi Water Resources Management Center, Nanning 530023, China
 ^{5}Department of Geosciences, University of Oslo, P.O. Box 1047 Blindern, 0316 Oslo, Norway
Correspondence: Cong Jiang (jiangcong@cug.edu.cn)
Hide author detailsCorrespondence: Cong Jiang (jiangcong@cug.edu.cn)
Multivariate hydrologic design under stationary conditions is traditionally performed through the use of the design criterion of the return period, which is theoretically equal to the average interarrival time of flood events divided by the exceedance probability of the design flood event. Under nonstationary conditions, the exceedance probability of a given multivariate flood event varies over time. This suggests that the traditional returnperiod concept cannot apply to engineering practice under nonstationary conditions, since by such a definition, a given multivariate flood event would correspond to a timevarying return period. In this paper, average annual reliability (AAR) was employed as the criterion for multivariate design rather than the return period to ensure that a given multivariate flood event corresponded to a unique design level under nonstationary conditions. The multivariate hydrologic design conditioned on the given AAR was estimated from the nonstationary multivariate flood distribution constructed by a dynamic Cvine copula, allowing for timevarying marginal distributions and a timevarying dependence structure. Both the mostlikely design event and confidence interval for the multivariate hydrologic design conditioned on the given AAR were identified to provide supporting information for designers. The multivariate flood series from the Xijiang River, China, were chosen as a case study. The results indicated that both the marginal distributions and dependence structure of the multivariate flood series were nonstationary due to the driving forces of urbanization and reservoir regulation. The nonstationarities of both the marginal distributions and dependence structure were found to affect the outcome of the multivariate hydrologic design.
A complete flood event or a flood hydrograph contains multiple feature variables, such as flood peak and flood volume, which can be associated with the safety of hydraulic structures (Salvadori et al., 2004, 2007, 2011; Xiao et al., 2009; Xiong et al., 2015; Loveridge et al., 2017; Shafaei et al., 2017). For example, the water level of a reservoir is controlled not only by flood peak flow but also by flood volume (Salvadori et al., 2011). Therefore, multivariate hydrologic design, which takes into account multiple flood characteristics as well as their dependence, provides a more rational design strategy for hydraulic structures compared to univariate hydrologic design (Zheng et al., 2013, 2014; Balistrocchi and Bacchi, 2017).
Multivariate hydrologic design under stationary conditions has been widely investigated, and the design criterion is usually quantified by the return period, similar to univariate hydrologic design. Under the definition of the average recurrence interval between flood events equaling or exceeding a given threshold (Chow, 1964), the return period of a given flood event under stationary conditions theoretically equals the average interarrival time between flood events divided by the exceedance probability (Salvadori et al., 2011). On the other hand, the exceedance probability of a univariate flood event is usually uniquely defined without ambiguity, whereas the exceedance probability of a multivariate flood event could have multiple definitions (Salvadori and De Michele, 2004; Salvadori et al., 2011; Vandenberghe et al., 2011). To date, at least five kinds of different exceedance probabilities for a multivariate flood event have been defined: (1) the OR case in which at least one of the flood features exceeds the prescribed threshold, (2) the AND case in which all flood features exceed the prescribed thresholds, (3) the Kendall case in which the univariate representation transformed from Kendall's distribution function exceeds the prescribed threshold, (4) the survival Kendall case in which the univariate representation transformed from survival Kendall's distribution function exceeds the prescribed threshold, and (5) the structural case in which the univariate representation transformed from a structure function exceeds the prescribed threshold (Favre et al., 2004; Salvadori and De Michele, 2004, 2010; Salvadori et al., 2007, 2013, 2015, 2016; Vandenberghe et al., 2011; Requena et al., 2013; Zheng et al., 2014).
Due to climate change as well as certain anthropogenic driving forces (Milly et al., 2008), the nonstationarities of both univariate and multivariate flood series have been widely reported (Xiong and Guo, 2004; Villarini et al., 2009; Vogel et al., 2011; López and Francés, 2013; Bender et al., 2014; Xiong et al., 2015; Blöschl et al., 2017; Kundzewicz et al., 2018). The multivariate flood distribution exhibits more complex nonstationarity behaviours than the univariate distribution, including nonstationarities of individual margins and the dependence structure between the margins (Quessy et al., 2013; Bender et al., 2014; Xiong et al., 2015; Kwon et al., 2016; Sarhadi et al., 2016; Qi and Liu, 2017; Vezzoli et al., 2017; Bracken et al., 2018; Salvadori et al., 2018). Both nonstationarities of the margins and dependence structure could impact the multivariate hydrologic design. Under nonstationary conditions, the exceedance probability p of a given flood event varies from year to year; thus, the return period, calculated as the average interarrival time between two successive flood events divided by p, is no longer a constant (Salas and Obeysekera, 2014; Jiang et al., 2015a; Kwon et al., 2016; Sarhadi et al., 2016; Yan et al., 2017). As a result, a given flood event would correspond to a timevarying and nonunique return period. Consequently, the traditional returnperiodbased method for estimating hydrologic design may no longer be applicable to engineering practice under nonstationary conditions (Salas and Obeysekera, 2014).
Although increasing attention has been focused on the hydrologic designs under nonstationary conditions in recent years, the focus has mainly been on univariate designs (Obeysekera and Salas, 2014, 2016; Read and Vogel, 2016). To overcome the limitation of the traditional return period under nonstationary conditions, the concept of the return period has been revisited. Salas and Obeysekera (2014) extended two concepts of the return period into a nonstationary framework, defined as the expected waiting time (EWT) for an exceedance to occur (Olsen et al., 1998), and the time period that results in the expected number of exceedances (ENE) equal to 1 over this period (Parey et al., 2010).
Risk and reliability are both important measurements for assessing hydrologic designs (Vogel, 1987; Read and Vogel, 2015). Besides redefinitions of the return period, some riskbased or reliabilitybased metrics have been introduced as the hydrologic design criteria under nonstationary conditions (Rosner et al., 2014). Rootzén and Katz (2013) proposed the concept of the design life level (DLL) to quantify hydrologic risk in a nonstationary climate during the entire design life period of hydraulic structures. Read and Vogel (2015) introduced the concept of average annual reliability (AAR) to estimate the hydrologic designs under nonstationary conditions. Liang et al. (2016) defined the equivalent reliability (ER) to estimate the design flood under nonstationary conditions by linking the DLL to the return period. Salvadori et al. (2018) associated hydrologic designs with both given life times and failure probabilities to calculate bivariate design values under nonstationarity. These design criteria assess the risk or reliability of hydraulic structures associated with the flood distribution during the entire design life period, rather than for a single year. For a given design life period, these criteria can always yield a unique risk or reliability; therefore, they are applicable to the hydrologic designs under both stationary and nonstationary conditions.
Under the multivariate framework, a given design level would correspond to an infinite number of possible hydrologic design events (Hawkes, 2008; Kew et al., 2013; Zheng et al., 2015, 2017); however, these design events are generally not equivalent because their joint probability density values (i.e. likelihood) usually differ (Salvadori et al., 2011; Volpi and Fiori, 2012; Li et al., 2017; Yin et al., 2017). In engineering practice, it should be necessary to determine a typical design event as representative of a specific design level. For example, in Chinese engineering practice, a unique design flood hydrograph corresponding to a given design level is usually required to determine the scale of hydraulic structures (Yin et al., 2017). The mostlikely design event, which theoretically has the largest joint probability density (likelihood) among all possible design events (Salvadori et al., 2011), appears to be the best representative candidate. Besides the mostlikely design event, it is also necessary to identify the confidence interval for an infinite possible number of design events to provide a finite design range for designers (Volpi and Fiori, 2012; Yin et al., 2017). The mostlikely design event and confidence interval for the bivariate hydrologic design under stationary conditions have been identified (Salvadori et al., 2011; Volpi and Fiori, 2012; Li et al., 2017; Yin et al., 2017; Salvadori et al., 2018); however, very few studies have focused on hydrologic designs with higher dimensions under nonstationary conditions.
Therefore, the objective of the present study was to address the issue of multivariate hydrologic design applying to engineering practice under nonstationary conditions, which is achieved through the following steps. First, the nonstationary multivariate flood distribution was constructed using a dynamic canonical vine (Cvine) copula (Aas et al., 2009), which was able to capture the nonstationarities of both marginal distributions and the dependence structure. The design criterion for the multivariate flood event was then quantified according to AAR rather than the traditional return period, since a given multivariate flood event would correspond to a unique AAR under both stationary and nonstationary conditions (Read and Vogel, 2015; Yan et al., 2017). The multivariate hydrologic design for any given AAR was estimated from the nonstationary multivariate flood distribution.
The aforementioned methods for the multivariate hydrologic design under nonstationary conditions were applied to the Xijiang River, China. The fourdimensional (4D) multivariate flood series, including the annual maximum daily discharge, annual maximum 3day flood volume, annual maximum 7day flood volume and annual maximum 15day flood volume of the Xijiang River were chosen as the case study data because they constitute the variables used for deriving the design flood hydrograph for hydraulic structures. It has been found that the natural flood processes of this river have been significantly altered by urbanization and reservoir regulation (Xu et al., 2014), but these two factors have not yet been taken into account in multivariate hydrologic design.
The next section of the present paper describes the study area and data. Section 3 presents the methods developed in this paper. The results of the case study are provided in Sect. 4. Finally, the conclusion and remarks are provided in Sect. 5.
The multivariate flood series of the Xijiang River, South China (see Fig. 1), were selected as a case study to illustrate the multivariate hydrologic design methods under nonstationary conditions. The drainage area of the Xijiang River basin (XRB) is 353 120 km^{2}, with a river length of 2214 km. The basin falls within a humid subtropical monsoon climate region, with the flood season extending from May to October; therefore, floods have always been a serious natural hazard within the basin.
The calculation of design floods in China involving the derivation of flood hydrographs for hydraulic structures requires not only flood peak but also flood volumes with different durations, such as 3, 7, 15 and 30 days (Ministry of Water Resources of People's Republic of China, 1996; Xiao et al., 2009; Xiong et al., 2015; Li et al., 2017). For a large catchment such as the XRB, the duration of a flood process is usually longer than 10 days. Therefore, the annual maximum daily discharge (Q_{1}), annual maximum 3day flood volume (V_{3}), annual maximum 7day flood volume (V_{7}) and annual maximum 15day flood volume (V_{15}) of the Xijiang River were defined as the multivariate flood series $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ for deriving design flood hydrographs. The flood data were from 1951 to 2012 and observed at the Dahuangjiangkou gauge located at the main stream of the Xijiang River, draining a total catchment area of 294 669 km^{2}, approximately 83 % of the total area of the XRB.
Rapid urbanization over recent decades has resulted in increasing river regulation projects built in the XRB, such as artificial levees for protecting urban areas from river flooding. As a result, flood flow has become increasingly constrained to the channel rather than overflowing to the floodplain, resulting in an increase in the observed river flood flow (Xu et al., 2014). For the purpose of flood control and hydropower generation, it is hard to find a river which is not impacted by reservoirs, particularly in rapidly developing China. Reservoir regulation has become an increasingly significant factor affecting flood processes of the XRB and should be seriously considered within downstream flood risk analysis and hydrologic design, particularly after 2006, when two reservoirs with considerable flood control capacities were put into operation. These are the Longtan and Baise reservoirs, with flood control capacities of 5 × 10^{9} m^{3} and 1.64×10^{9} m^{3} and catchment areas of 98 500 and 9600 km^{2}, respectively. Climate change will likely result in flood nonstationarity by altering climatic conditions of the basin. Climatic conditions dominating flood processes in the XRB, such as extreme precipitation, appear to have been stationary over the past decades (Yang et al., 2010). Therefore, the current study introduced only urbanization and reservoir regulation as the potential driving forces of flood nonstationarity and ignored the effect of climate change.
The effect of urbanization on flood processes was quantified using the urban population (Pop). Given the unavailability of urban population data at the basin scale and the fact that the vast majority of cities in the XRB are distributed in Guangxi province, we used urban population data for Guangxi province to represent those of the XRB. The annual urban population data during the observation period were obtained from the China Compendium of Statistics 1949–2008 (Department of Comprehensive Statistics of National Bureau of Statistics, 2010) and the website of the National Bureau of Statistics of China (http://www.stats.gov.cn/tjsj/ndsj/, last access: 20 March 2019). The present study assumed the design life period for hydraulic structures to be from 2013 to 2100. The urban population over the design life period was estimated based on the predicted growth rate of China's urban population reported by He (2014). The reservoir index (RI), which depends on the catchment area and flood controlling capacities of reservoirs, was used to quantify the effects of reservoir regulation on flood processes (López and Francés, 2013). As shown in Table 1, two reservoirs with flood control functions have been completed during the observation period from 1951 to 2012, and a further two are planned for operation during the design life period. Figure 2 illustrates the evolution of the urban population and reservoir index during both the observation and design life periods.
The present study included the following methods: (1) the nonstationary multivariate flood distribution based on a dynamic Cvine copula, allowing for both timevarying marginal distributions and a timevarying dependence structure, and (2) estimation of the multivariate hydrologic design associated with AAR under nonstationary conditions. To correspond to the case study in this paper, the multivariate flood series consisting of the annual maximum daily discharge (Q_{1}), annual maximum 3day flood volume (V_{3}), annual maximum 7day flood volume (V_{7}) and annual maximum 15day flood volume (V_{15}) were chosen to illustrate the multivariate design methods under nonstationary conditions. It is worth noting that the proposed methods can be extended to other multivariate flood series, such as those consisting of flood peak, flood volume and flood duration.
3.1 Probability distribution of the nonstationary multivariate flood series
According to Sklar's theorem (Sklar, 1959), the probability distribution of the 4D flood series $\left({Q}_{\mathrm{1}}\phantom{\rule{0.125em}{0ex}},{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ at time t measured by years ($t=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{2},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}n$, and n is the length of the flood series) can be formulated through a copula C(⋅) as follows:
where ${F}_{\mathrm{1}}\left({q}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}\right)$, ${F}_{\mathrm{3}}\left({v}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t}\right)$, ${F}_{\mathrm{7}}\left({v}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t}\right)$ and ${F}_{\mathrm{15}}\left({v}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}\right)$ denote the marginal distributions for Q_{1}, V_{3}, V_{7} and V_{15}, respectively; u_{1, t}, u_{3, t}, u_{7, t} and u_{15, t} are the marginal probabilities of Q_{1}, V_{3}, V_{7} and V_{15}, respectively; θ_{1, t}, θ_{3, t}, θ_{7, t} and θ_{15, t} are the corresponding distribution parameters; and θ_{c, t} stands for the copula parameter vector, which describes the strength of the dependence structure. ${\mathit{\theta}}_{t}=\left({\mathit{\theta}}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t},{\mathit{\theta}}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t},{\mathit{\theta}}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t},{\mathit{\theta}}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t},{\mathit{\theta}}_{c,\phantom{\rule{0.125em}{0ex}}t}\right)$ is the parameter vector of the entire multivariate distribution, including the marginal distribution parameters as well as the copula parameters.
According to the multivariate distribution of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ defined by Eq. (1), the corresponding density function can be written as:
where ${f}_{\mathrm{1}}\left({q}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}\right)$, ${f}_{\mathrm{3}}\left({v}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t}\right)$, ${f}_{\mathrm{7}}\left({v}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t}\right)$ and ${f}_{\mathrm{15}}\left({v}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}\right)$ are the density functions of the marginal distributions for Q_{1}, V_{3}, V_{7} and V_{15}, respectively, and c(⋅) denotes the density function of copula C(⋅). As shown by Eq. (2), the multivariate distribution of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ can be separated into two modules, including the marginal distributions, i.e. ${f}_{\mathrm{1}}\left({q}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}\right)$, ${f}_{\mathrm{3}}\left({v}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t}\right)$, ${f}_{\mathrm{7}}\left({v}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t}\right)$ and ${f}_{\mathrm{15}}\left({v}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}\right)$, as well as the dependence structure expressed by the copula density function $c\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{c,\phantom{\rule{0.125em}{0ex}}t}\right)$. Under nonstationary conditions, both the margins and dependence structure of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ can vary over time t.
3.1.1 Nonstationary marginal distributions based on the timevarying moment model
The timevarying moment model that expresses the distribution parameters or moments as functions of time or some other explanatory variable or variables have been widely employed to capture the nonstationarities of univariate flood series (Strupczewski et al., 2001; Villarini et al., 2009). In this study, the nonstationary marginal distributions of the multivariate flood series $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ were constructed by the timevarying moment model.
Based on cause–effect analysis, the flood processes of the XRB were found to mainly be impacted by urbanization and reservoir operation. The reservoir index RI and urban population Pop were therefore used as potential covariates for marginal distribution parameters, including the location parameter μ, scale parameter σ and shape parameter ν (if any). In this study, both linear and exponential functions were considered to build the relationships between distribution parameters and covariates (Strupczewski et al., 2001; Vogel et al., 2011; Salas and Obeysekera, 2014; Jiang et al., 2015; Sarhadi et al., 2016; Read and Vogel, 2016; Yan et al., 2017). Taking the location parameter for illustration, the candidate functions of μ were generally formulated as follows.
Here α_{0}, α_{1} and α_{2} are model parameters estimated using the maximum likelihood estimate (MLE) method (Strupczewski et al., 2001). As above, the linear expression in Eq. (3) gives an additive model which suggests that the effects of the covariates RI and Pop on μ are independent, while the exponential expression defines a multiplicative model which is able to take into account the possible interaction between the covariates RI and Pop. It is important to note that Eq. (3) defines four specific nonstationary models: the first one is the most complex nonstationary model where it is assumed that both RI and Pop are the driving factors of marginal distributions, the second and third models illustrate that the marginal nonstationarity is linked only to RI and Pop, respectively, and the final one represents the simplest and stationary model, which does not contain any covariates.
Five probability distributions widely used in flood frequency analysis, namely Pearson type III (PIII), generalized extreme value (GEV), gamma, Weibull and lognormal distributions, were employed as the candidate distributions for margins (Villarini et al., 2009; Yan et al., 2017). The goodness of fit (GoF) of the probability distributions was examined by the Kolmogorov–Smirnov (KS) test with a significance level of 0.05 (Frank and Massey, 1951). The p value of the KS test was simulated by the Monte Carlo method. The relative fitting qualities of the timevarying moment models were assessed by the corrected Akaike information criterion (AICc; Hurvich and Tsai, 1989), which is stricter than the Akaike information criterion (AIC; Akaike, 1974). The best model featured with the smallest AICc value was chosen to describe the marginal distributions from the nonstationary models as expressed by Eq. (3).
3.1.2 Nonstationary dependence structure based on the dynamic Cvine copula
After estimating the marginal distributions, the nonstationary dependence structure of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ as formulated by the copula density function $c\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{c,\phantom{\rule{0.125em}{0ex}}t}\right)$ was constructed. Given that most applied copula functions are for bivariate random variables, $c\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{c,\phantom{\rule{0.125em}{0ex}}t}\right)$ cannot be directly expressed as a specific copula function. The pair copula method has been proven to be powerful for the construction of the distribution of multivariate random variables through the decomposition of the multivariate probability density into a series of bivariate copulas (Aas et al., 2009; Xiong et al., 2015; Shafaei et al., 2017). Therefore this study constructed the dependence structure of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ using the pair copula method.
Numerous pair copula decomposition forms for a multivariate distribution are available, among which two kinds of decompositions with regular vine structures prevail in practice, namely the canonical vine (C vine) and the drawable vine (D vine; Aas et al., 2009). It is known that flood peak (e.g. Q_{1}) is the dominant feature quantifying a flood event as well as being the key factor in hydrologic design (Ministry of Water Resources of People's Republic of China, 1996). The C vine is more suitable when there is a key variable governing multivariate dependence (Aas et al., 2009). In this case, the C vine was employed to construct the joint distribution of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$, with Q_{1} elected as the key variable. Thus, the density function $c\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}t}{\mathit{\theta}}_{c,\phantom{\rule{0.125em}{0ex}}t}\right)$ can be decomposed into six bivariate pair copulas as follows:
where ${\mathit{\theta}}_{c,\phantom{\rule{0.125em}{0ex}}t}=\left({\mathit{\theta}}_{\mathrm{13},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{\mathit{\theta}}_{\mathrm{17},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{\mathit{\theta}}_{\mathrm{115},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{\mathit{\theta}}_{\mathrm{37}\mathrm{1},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{\mathit{\theta}}_{\mathrm{315}\mathrm{1},\phantom{\rule{0.125em}{0ex}}t},\phantom{\rule{0.125em}{0ex}}{\mathit{\theta}}_{\mathrm{715}\mathrm{13},\phantom{\rule{0.125em}{0ex}}t}\right)$ is the parameter vector in the Cvine copula, and
Figure 3 shows the schematic decomposition of the 4D Cvine copula as expressed by Eq. (4). It is evident that the hierarchical structure of the 4D Cvine copula contains three trees and six edges. The first tree (T1) includes three bivariate pair copulas, i.e. ${c}_{\mathrm{13}}\left(\cdot \left{\mathit{\theta}}_{\mathrm{13},\phantom{\rule{0.125em}{0ex}}t}\right.\right)$, ${c}_{\mathrm{17}}\left(\cdot \left{\mathit{\theta}}_{\mathrm{17},\phantom{\rule{0.125em}{0ex}}t}\right.\right)$ and ${c}_{\mathrm{115}}\left(\cdot \left{\mathit{\theta}}_{\mathrm{115},\phantom{\rule{0.125em}{0ex}}t}\right.\right)$, which directly act on the marginal probabilities and describe the bivariate dependencies between the key variable Q_{1} and the other three variables, i.e. V_{3}, V_{7} and V_{15}. The second tree (T2) includes two bivariate pair copulas ${c}_{\mathrm{37}\mathrm{1}}\left(\cdot \left{\mathit{\theta}}_{\mathrm{37}\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}\right.\right)$ and ${c}_{\mathrm{315}\mathrm{1}}\left(\cdot \left{\mathit{\theta}}_{\mathrm{315}\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}\right.\right)$, which act on the conditional distribution functions with u_{1, t} as the conditioning variable. Finally, the third tree (T3) includes only one bivariate pair copula ${c}_{\mathrm{715}\mathrm{13}}\left(\cdot \left{\mathit{\theta}}_{\mathrm{715}\mathrm{13},\phantom{\rule{0.125em}{0ex}}t}\right.\right)$ acting on conditional distribution functions with both u_{1, t} and u_{3, t} as the conditioning variables.
In flood frequency analysis, the upper tail of the flood distribution deserves more attention because it allows the quantification of risks of the more serious flood events. The Gumbel–Hougaard copula, an extremevalue copula widely used in hydrology, accounts for the upper tail dependence and is wellsuited to the dependence structure of a multivariate flood distribution (Salvadori et al., 2007; Zhang and Singh, 2007; Xiong et al., 2015). Consequently, the present study employed the bivariate Gumbel–Hougaard copula to construct the dynamic Cvine copula formulated by Eq. (4). The bivariate Gumbel–Hougaard copula is expressed as follows:
where u and v are the bivariate marginal probabilities and θ_{c} is the single parameter measuring the dependence strength.
Similar to the nonstationary marginal distributions, the nonstationarity of the dependence structure of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ was characterized by the time variations of the copula parameters in T1, T2 and T3. Both linear and exponential functions were considered to characterize the timevarying copula parameters and generally formulated as follows.
Here β_{0}, β_{1} and β_{2} are model parameters estimated using the MLE method (Aas et al., 2009). Here, the exponential expression in Eq. (7) was written as the sum of 1 and an exponential function of the covariates so that the domain range of the copula parameter θ_{c} can be satisfied under any condition. To make it easy in parameter estimation, the model parameters for each pair copula were separately estimated. The model parameters for θ_{13, t}, θ_{17, t} and θ_{115, t} in T1 were first estimated, and those for the remaining copula parameters ${\mathit{\theta}}_{\mathrm{37}\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}$, ${\mathit{\theta}}_{\mathrm{315}\mathrm{1},\phantom{\rule{0.125em}{0ex}}t}$ and ${\mathit{\theta}}_{\mathrm{715}\mathrm{13},\phantom{\rule{0.125em}{0ex}}t}$ in T2 and T3 were then estimated in sequence. It is worth noting that these parameters can be also simultaneously estimated. These two methods could result in possible difference in parameter estimation.
The available GoF tests for vine copulas are very limited, with the probability integral transform (PIT) test appearing to be reliable (Aas et al., 2009). Under a null hypothesis of the multivariate flood variables $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ following a given Cvine copula, the PIT converts the dependent flood variables into a new set of variables that are independent and uniformly distributed on [0, 1]^{4}. The GoF of vine copulas can be obtained through determining whether the resulting variables are independent and uniform in [0, 1]. For more details of the PIT test, readers are referred to Aas et al. (2009). The best nonstationary model for each bivariate pair copula in Eq. (4) was chosen from the nonstationary models generally expressed by Eq. (7) in terms of the AICc value.
3.2 Multivariate hydrologic design under nonstationary conditions
3.2.1 Average annual reliability for multivariate flood events
The AAR introduced by Read and Vogel (2015) was calculated using the arithmeticaverage method, thereby taking into account the reliability of each year with the same weighting factor. A safer design strategy should pay more attention to worse (i.e. lower) annual reliability; however, the arithmeticaverage AAR is not capable of this function. The present study employed the geometricaverage method to calculate AAR, which is dominated more by the minimum than arithmetic average and is theoretically able to yield safer design values. The geometricaverage AAR is also equivalent to the metrics of the DLL (Rootzén and Katz, 2013) and ER (Liang et al., 2016; Yan et al., 2017).
Denoting $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$ as a given multivariate flood event, its exceedance probability p_{t}, which is the occurrence probability of a more dangerous multivariate event than $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$ in a specific hazard scenario, would vary from year to year under nonstationary conditions. The AAR for $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$ was calculated by the geometricaverage method as follows:
where T_{1} and T_{2} stand for the beginning year and ending year of the operation of an assumed hydraulic structure, respectively, ${T}_{\mathrm{2}}{T}_{\mathrm{1}}+\mathrm{1}$ is the length of the design life period of the hydraulic structure, and 1−p_{t} measures the annual reliability of the given multivariate flood event $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$ at time t.
3.2.2 Exceedance probabilities of multivariate flood events
The present study characterized AAR by considering three widely used definitions of the exceedance probabilities of the multivariate flood event $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$, i.e. the OR, AND and Kendall cases (Salvadori and De Michele, 2004, 2010; Favre et al., 2004; Salvadori et al., 2007, 2016; Vandenberghe et al., 2011). The OR case for $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$ defines the case under which at least one of the flood features exceeds the prescribed threshold. The exceedance probability in the OR case at time t was denoted as ${p}_{t}^{\mathrm{or}}$ and was calculated by the following:
where “∨” stands for the OR operator and F(⋅θ_{t}) is defined in Eq. (1).
The AND case for $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$ defines the case under which all of the flood features exceed the prescribed thresholds, and the corresponding exceedance probability ${p}_{t}^{\mathrm{and}}$ at time t was expressed as follows:
where “∧” is the AND operator and f(⋅θ_{t}) is defined in Eq. (2).
Under the Kendall case, the multivariate flood event $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$ was first transformed into a univariate representation via Kendall's distribution function K_{C}(⋅) as follows:
where ${\mathit{\rho}}_{t}=F\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\left{\mathit{\theta}}_{t}\right.\right)$ is the probability level corresponding to the given flood event $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$. The exceedance probability ${p}_{t}^{\mathrm{ken}}$ in the Kendall case at time t was expressed as follows:
For general multivariate cases, the exceedance probabilities ${p}_{t}^{\mathrm{or}}$, ${p}_{t}^{\mathrm{and}}$ and ${p}_{t}^{\mathrm{ken}}$ could have no analytical solutions but can be numerically estimated through the Monte Carlo method (Niederreiter, 1978; Salvadori et al., 2011, 2013). The AAR in the OR, AND and Kendall cases can be calculated by replacing the exceedance probability p_{t} in Eq. (8) by ${p}_{t}^{\mathrm{or}}$, ${p}_{t}^{\mathrm{and}}$ and ${p}_{t}^{\mathrm{ken}}$, respectively.
3.2.3 Mostlikely design event and confidence interval for multivariate hydrologic design
The methods identifying both the mostlikely design event, denoted by $\left({z}_{{Q}_{\mathrm{1}}}^{*},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}}^{*},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}}^{*},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}^{*}\right)$, and the confidence interval for the multivariate hydrologic design $\left({z}_{{Q}_{\mathrm{1}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}\right)$ given AAR=η are introduced below. The average annual probability density, denoted by g(⋅), of $\left({z}_{{Q}_{\mathrm{1}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}\right)$ over the entire design life period from T_{1} to T_{2}, was expressed as follows:
The probability distribution function for AAR≤η can be written as
By denoting the density function of Φ(η) as ϕ(η), the probability density of $\left({z}_{{Q}_{\mathrm{1}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}\right)$ conditioned on AAR=η can be expressed as follows:
The mostlikely design event conditioned on AAR=η was theoretically written as
Unfortunately, the analytical solutions of both the mostlikely design event $\left({z}_{{Q}_{\mathrm{1}}}^{*},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}}^{*},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}}^{*},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}^{*}\right)$ and confidence interval are unavailable but can be approximately estimated through the Monte Carlo simulation method. First, the design events with sample size N conditioned on AAR=η were generated. These design events were then sorted in descending order of their multivariate probability densities, denoted by
where $Nc=N\cdot {p}_{\mathrm{c}}$, and p_{c} is the critical probability level for the confidence interval. Thus, the approximate solution for $\left({z}_{{Q}_{\mathrm{1}}}^{*},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}}^{*},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}}^{*},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}^{*}\right)$ is $\left({z}_{{Q}_{\mathrm{1}}}^{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}}^{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}}^{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}^{\mathrm{1}}\right)$. The lower boundary for the confidence interval was expressed as follows:
The upper boundary for the confidence interval was estimated by the following:
3.2.4 Derivation of design flood hydrographs
In China, the design flood hydrographs for hydraulic structures are usually derived from the design flood events set against a benchmark flood hydrograph, which is chosen from the observed flood processes (Ministry of Water Resources of People's Republic of China, 1996; Xiao et al., 2009, Yin et al., 2017). For example, suppose that a flood hydrograph consists of the features of annual maximum daily discharge, 3day flood volume, 7day flood volume and 15day flood volume. The four features of the benchmark flood hydrograph are denoted by ${Q}_{\mathrm{1}}^{B}$, ${V}_{\mathrm{3}}^{B}$, ${V}_{\mathrm{7}}^{B}$ and ${V}_{\mathrm{15}}^{B}$, respectively. The design flood hydrograph corresponding to the multivariate hydrologic design realization $\left({z}_{{Q}_{\mathrm{1}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}\right)$ can be derived by multiplying the benchmark flood hydrograph by different amplifiers, given as described below.
The amplifier K_{1} for the annual maximum daily discharge was calculated by
The amplifier K_{3−1} for the 3day flood volume except for the annual maximum daily discharge was calculated by
where V(⋅) is the operator transforming daily discharge into flood volume. The amplifier K_{7−3} for the 7day flood volume except for the 3day flood volume was calculated by
Finally, the amplifier K_{15−7} for the 15day flood volume except for the 7day flood volume was calculated by
4.1 Nonstationary analysis for marginal distributions
The timevarying moment model was employed to perform nonstationary analysis for each marginal distribution of the multivariate flood series $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ of the Xijiang River. In general, the candidate distributions for all margins passed the GoF test at the 0.05 significance level. The chosen models featured with the smallest AICc values were shown in Table 2. The results indicated that the GEV distribution provided the best fit for the annual maximum daily discharge series Q_{1}, whereas the Gamma distribution was chosen as the theoretical distribution for the flood volume series V_{3}, V_{7} and V_{15}. All estimated model parameters were found to be statistically significant at the 0.05 level. The 95 % uncertainty intervals for the estimated parameters were calculated by the parametric bootstrap method (Kyselý, 2009). In accordance with the modelling results, it can be seen that the location parameters μ for all flood series were nonstationary, while the scale and shape parameters were stationary. Through an exponential function, the location parameters μ referring to the means of the flood series were generally positively related to the urban population Pop, whereas they were negatively related to the reservoir index RI. This finding revealed the opposite roles played by urbanization and reservoir regulation on the flood processes of the XRB. In particular, more artificial levees were required to protect urban areas from flooding by constraining the flood flow to river channels, which resulted in increasing the river channel flood flow. The reservoirs played an active role in flood control by reducing the flood discharge downstream.
The relationships between μ and covariates were built by the exponential function in Eq. (3). α_{1} and α_{2} are the parameters related to urban population (Pop) and reservoir index (RI), respectively. The symbols ^{a}, ^{b} and ^{c} denote that the estimated model parameters are significant at the levels of 0.1, 0.05 and 0.1, respectively. The numbers in brackets are the 95 % uncertainty interval. p_KS stands for the p value of the KS test for marginal distributions.
More specific to each margin of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$, the location parameters μ of the three shortduration flood series, i.e. Q_{1}, V_{3} and V_{7}, were positively linked to Pop, whereas the RI was the driving factor reducing μ for all flood series, including Q_{1}, V_{3}, V_{7} and V_{15}. Owing to the difference in covariate selections, the shortduration flood series, including Q_{1}, V_{3} and V_{7}, displayed asynchronous nonstationary behaviours with the longduration flood series V_{15} occurring in the observation period of 1951–2012. As shown in Fig. 4, Q_{1}, V_{3} and V_{7} presented significantly increasing trends during 1951–2005, particularly since the 1980s, marking the beginning of a period of rapid urbanization in China. V_{15} tended to follow a stationary process during 1951–2005. After the two flood control reservoirs were put into operation in 2006, all flood series, including Q_{1}, V_{3}, V_{7} and V_{15}, exhibited a sharp decline.
The relationships between copula parameters and covariates were built by the exponential function in Eq. (7). β_{1} and β_{2} are the parameters related to urban population (Pop) and reservoir index (RI), respectively. The symbols ^{a}, ^{b} and ^{c} denote that the estimated model parameters are significant at the levels of 0.1, 0.05 and 0.01, respectively. The numbers in brackets are the 95 % uncertainty interval.
The predicted marginal distributions for $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ during the design life period from 2013 to 2100 were estimated using the timevarying moment model by replacing the observed covariates for μ with those predicted. Figure 4 also shows that the mean values of Q_{1}, V_{3} and V_{7} during the design life period increased with the growth of the urban population, following which they decreased sharply in 2023 after a larger reservoir named Datengxia is expected to be put into operation. After 2023, with no more reservoirs planned, the predicted mean values of Q_{1}, V_{3} and V_{7} would be expected to reach their peaks in the mid21th century followed by a slight declining trend because of a shrinking urban population. Since V_{15} was only related to RI, V_{15} would show an abrupt decline in 2023 due to the regulation of the Datengxia reservoir. In general, the predicted nonstationary marginal distributions for Q_{1} and V_{3} during 2013–2100 were roughly approximate to the marginal distributions under the assumption of stationarity, whereas the predicted nonstationary marginal distributions for V_{7} and V_{15} exhibited smaller mean values than those of the stationary distributions.
4.2 Nonstationary dependence structure for $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$
After estimating the nonstationary marginal distributions for $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$, the multivariate dependence structure was constructed by the dynamic Cvine copula with Q_{1} elected as the key variable. Figure 5 illustrates significant correlations between the flood peak Q_{1} and the flood volumes (i.e. V_{3}, V_{7} and V_{15}). Table 3 shows the estimation results of the dynamic Cvine copula. The PIT test for the nonstationary dependence structure of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ suggested a satisfactory fitting effect, and most estimated parameters were statistically significant at the 0.05 level. The results indicated that the copula parameter θ_{115} for pair (Q_{1}, V_{15}) was found to be nonstationary and expressed as an exponential function of both the urban population Pop and reservoir index RI, whereas other copula parameters indicated stationary dependences. It was seen that the margin of Q_{1} displayed asynchronous nonstationarity behaviours with V_{15} (see Table 2 and Fig. 4). Therefore, the dependence nonstationarity of the pair (Q_{1}, V_{15}) could possibly be attributed to the asynchronous marginal nonstationarities.
According to the regression function, θ_{115} was negatively related to Pop, whereas it was positively related to RI. In other words, growing urbanization weakened the multivariate flood dependence, whereas reservoir regulation played an opposite role, enhancing the dependence. This finding indicated that human activities, including urbanization and reservoir regulation, not only changed the statistical characteristics of the marginal distributions of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ but also affected the dependence of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$. Figure 6 shows the time variations of θ_{115} during the observation period of 1951–2012 as well as during the design life period of 2013–2100. Due to reservoir regulation, θ_{115} presented two obvious upward change points in both 2006 and 2023. Besides this, θ_{115} also exhibited an obvious decreasing trend with urban population growth from 1951 to the mid21th century, followed by a slight increasing trend due to a shrinking urban population. During the design life period, the predicted nonstationary θ_{115} suggests a weaker dependence for $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ than the dependence under the stationary assumption, since it is usually smaller than the stationary estimation.
In addition, the changepoint detection method based on the Cramér–von Mises statistic (Bücher et al., 2014) was employed to detect possible nonstationarities in both the marginal distributions and dependence of the multivariate flood series $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$. Readers are referred to Bücher et al. (2014) and Kojadinovic (2017) for specific steps to implement the changepoint detection. The results indicated that neither the marginal distributions nor dependence displayed change points at the 0.05 significance level (see Table 4), whereas the previous analysis suggested nonstationary margins and dependence due to the joint effects of urbanization and reservoir regulation. These aforementioned inconsistencies could be attributed to the opposite roles of urbanization and reservoir regulation on shifting of the multivariate flood distribution, with urbanization generally enlarging the mean values of the flood series and weakening their dependence and reservoir regulation, decreasing the mean values and strengthening the dependence. In other words, the nonstationarities induced by these two factors may have offset each other. As a result, the nonstationarities of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ might have not been captured by the statistical method based on the Cramér–von Mises statistic. This finding highlights the significance of cause–effect analysis in judging the nonstationarities of hydrologic series (Xiong et al., 2015).
4.3 Multivariate hydrologic design characterized by average annual reliability
The multivariate hydrologic designs, characterized by AAR associated with the OR, AND and Kendall exceedance probabilities, were estimated from the predicted nonstationary multivariate distribution for $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ during the design life period from 2013 to 2100. The left columns in Figs. 7–10 show the mostlikely design events and the 90 % confidence intervals conditioned on the AAR varying from 0.01 to 0.99. The multivariate hydrologic design events associated with both the OR and Kendall exceedance probabilities exhibited the lower boundaries, whereas the design events associated with the AND exceedance probability exhibited the upper boundaries.
The design flood hydrographs were derived from the multivariate hydrologic designs against the benchmark flood hydrograph observed in 1988. Figure 11 shows the design flood hydrographs by setting AAR equal to 0.90, 0.95 and 0.99. For any given multivariate flood event, the corresponding OR exceedance probability was larger than that of AND, with the Kendall exceedance probability somewhere in between (Vandenberghe et al., 2011). These differences among the OR, AND and Kendall exceedance probabilities indicate the different design strategies. It must be noted that the choice of design strategy in engineering practice is usually a priori and is dependent on the specific design requirements and mechanisms of failure for hydraulic structures (Serinaldi, 2015; Salvadori et al., 2016).
We calculated the univariate hydrologic design events from the predicted marginal distributions to compare the design strategies under the multivariate framework with those under the univariate framework. Figures 7–10 show that the univariate hydrologic design events exactly constituted the lower boundaries of the multivariate hydrologic design events associated with the OR exceedance probability as well as the upper boundaries of the design events associated with the AND exceedance probability. Under a given AAR, the hydrologic designs under the univariate framework were generally smaller than the mostlikely design events associated with the OR exceedance probability, whereas they were larger than those associated with the AND exceedance probability; they were most approximate to those associated with the Kendall exceedance probability. The comparisons of the flood hydrographs displayed in Fig. 11 reinforced these findings.
4.4 Impacts of multivariate nonstationarity behaviours on hydrologic design values
Section 4.1 and 4.2 show the marginal distribution and dependence structure of the multivariate flood distribution of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ to be nonstationary. We estimated the multivariate hydrologic design events under an assumption of stationarity to illustrate how these nonstationarities act on the multivariate hydrologic designs, i.e. both marginal distributions and the dependence structure were treated as stationary (see the right columns in Figs. 7–10). Figure 4 suggests that both the predicted nonstationary marginal distributions for Q_{1} and V_{3} during the design life period were approximate to the stationary marginal distributions. Therefore, the nonstationary and stationary marginal distributions yielded similar design values for ${z}_{{Q}_{\mathrm{1}}}$ and ${z}_{{V}_{\mathrm{3}}}$ (see Figs. 7 and 8). The predicted nonstationary distributions for both V_{7} and V_{15} indicated smaller mean values compared to those of the stationary distributions (see Fig. 4); therefore, the corresponding hydrologic design values estimated from the nonstationary marginal distributions were generally smaller than those estimated from the stationary marginal distributions (see Figs. 9 and 10).
The nonstationary multivariate flood distribution during the design life period was also predicted to exhibit a weaker dependence structure than that of the stationary distribution (see Fig. 6). The dependence nonstationarity was expected to have a much subtler effect on the multivariate hydrologic design compared to the marginal nonstationarities (Xiong et al., 2015). To illustrate the individual effect of the dependence nonstationarity on the multivariate hydrologic design, an artificial nonstationary condition for the multivariate flood distribution was set so that only the marginal nonstationarities were considered, whereas the dependence structure was treated as stationary. The results of the multivariate hydrologic design events are shown in the middle columns in Figs. 7–10. In general, the dependence nonstationarity had less of an effect on the multivariate hydrologic designs compared the marginal nonstationarities; however, some visible differences in both the 90 % confidence intervals were still identified. In summary, the nonstationary and weaker dependence structure generally suggested wider confidence intervals for the multivariate hydrologic design values.
The statistical characteristics of both the marginal distributions and the dependence structure of multivariate flood variables can vary with time under nonstationary conditions. It is possible that the multivariate flood distribution estimated from the historical information will not reflect the statistical characteristics of flooding in the future. As a result, the stationarybased hydrologic design would not be able to deal with potential hydrologic risks of hydraulic structures. It is necessary for hydrologic designers to take into account the physical driving forces (such as human activates and climate change) responsible for the nonstationarities of multivariate flood variables.
The present study introduced possible methods for addressing multivariate hydrologic design for application in engineering practice under nonstationary conditions. A dynamic Cvine copula allowing for both timevarying marginal distributions and timevarying dependence structure was developed to capture the nonstationarities of a multivariate flood distribution. The multivariate hydrologic design under nonstationary conditions was estimated by specifying the design criterion by average annual reliability. The mostlikely design event and confidence interval were identified as the outcome of the multivariate hydrologic design. Multivariate flood series $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$ from the Xijiang River were chosen as a case study, with the main findings given below.
Urbanization and reservoir regulation were found to be the driving forces responsible for the nonstationarities of both the marginal distributions and dependence structure of the multivariate flood series $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$. The growth of the urban population generally resulted in an increased mean value of the individual flood series, whereas it weakened the dependence of $\left({Q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{15}}\right)$. The increasing reservoir index had the opposite effect on the individual flood series as well as their dependence. Under a given average annual reliability, the OR exceedance probability yielded the largest design values, followed by the Kendall and the AND exceedance probabilities. Nonstationarities in both marginal distributions and dependence structure affected the outcome of the multivariate hydrologic design. It is the marginal nonstationarities that played a dominant role in affecting the multivariate hydrologic design.
There are two remarks that can be made that are related to the practical implications of the hydrologic design methods developed in the current study that are detailed below.
The first remark relates to the length of observed flood data required for multivariate and nonstationary hydrologic designs. In theory, sufficiently long observed flood data (or other extremevalue data) are required to derive robust estimations of the distribution parameters and the correct hydrologic design values (Zheng et al., 2018). However, in reality, most data series are limited in length, thus forcing us to use what we have at hand to do research or design works without fulfilling the theoretical assumptions or requirements. Some recent studies suggested that univariate flood frequency analysis under stationary conditions usually requires flood data with a continuous period of at least 30 years (Ministry of Water Resources of People's Republic of China, 1996; Engeland et al., 2018; Kobierska et al., 2018). However, determining a definitive answer to the length of observed flood data that is required for flood frequency analysis under multivariate and/or nonstationary settings poses a challenge, since this issue has not yet been fully addressed. However, it is certain that multivariate and nonstationary hydrologic designs naturally require datasets of longer length, since they usually contain more parameters to be estimated.
The second remark related to the tradeoff between reducing estimation bias and increasing model uncertainty. Nonstationary models generally improve performance in fitting observation data by reducing estimation bias (Jiang et al., 2015b), but this is usually achieved at the expense of increasing model complexity, such as adding more model parameters and introducing more nonstationary covariates, which might induce additional sources of model uncertainty (Serinaldi and Kilsby, 2015; Read and Vogel, 2016). A careful balance between the model fitting effect and the model complexity should be maintained in practice when employing multivariate and nonstationary hydrologic design by keeping in mind the following two points: (1) the multivariate and nonstationary models should remain effective but should also be kept as simple as possible to avoid overfitting, and (2) to ensure a robust relationship between the distribution parameters and the explanatory covariates, the chosen covariates should be physically related to the flood processes and supported by a welldefined cause–effect analysis.
All the data used in this study can be requested by contacting the corresponding author Cong Jiang at jiangcong@cug.edu.cn.
A1 OR exceedance probability (formulated by Eq. 9 in the paper)
Since the cumulative distribution function $F\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\left{\mathit{\theta}}_{t}\right.\right)$ has no analytical expression, the OR exceedance probability ${p}_{t}^{\mathrm{or}}$ at time t is calculated by the Monte Carlo method as follows:

Calculate the marginal probabilities $\left({u}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15}}\right)$ of $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$.

Generate m samples $\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}i},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}i},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}i},{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}i}\right)$ ($i=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{2}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},m)$ from the Cvine copula expressed by Eq. (4).

Calculate $F\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\left{\mathit{\theta}}_{t}\right.\right)=\frac{\mathrm{1}}{m+\mathrm{1}}{\sum}_{i=\mathrm{1}}^{m}\phantom{\rule{0ex}{0ex}}\mathbf{1}\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}i}\le {u}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}i}\le {u}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}i}\le {u}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}i}\le {u}_{\mathrm{15}}\right)$.

Calculate ${p}_{t}^{\mathrm{or}}=\mathrm{1}F\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\left{\mathit{\theta}}_{t}\right.\right)$.
A2 AND exceedance probability (formulated by Eq. 10 in the paper)
The AND exceedance probability ${p}_{t}^{\mathrm{and}}$ at time t is calculated by the Monte Carlo method as follows:

Calculate the marginal probabilities $\left({u}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15}}\right)$ of $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$.

Generate m samples $\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}i},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}i},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}i},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}i}\right)$ ($i=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{2},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},m)$ from the Cvine copula expressed by Eq. (4).

Calculate ${p}_{t}^{\mathrm{and}}=\frac{\mathrm{1}}{m+\mathrm{1}}{\sum}_{i=\mathrm{1}}^{m}\phantom{\rule{0ex}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathbf{1}\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}i}\ge {u}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}i}\ge {u}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}i}\ge {u}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}i}\ge {u}_{\mathrm{15}}\right)$.
A3 The Kendall exceedance probability (formulated by Eqs. 11 and 12 in the paper)
The Kendall exceedance probability ${p}_{t}^{\mathrm{ken}}$ at time t is calculated by the Monte Carlo method as follows:

Calculate the marginal probabilities $\left({u}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15}}\right)$ of $\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\right)$.

Calculate ${\mathit{\rho}}_{t}=F\left({q}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{3}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{7}},\phantom{\rule{0.125em}{0ex}}{v}_{\mathrm{15}}\left{\mathit{\theta}}_{t}\right.\right)$ (see calculation steps 2–3 for OR exceedance probability).

Generate m samples $\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}i},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}i},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}i},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}i}\right)$ ($i=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{2},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},m)$ from the Cvine copula expressed by Eq. (4).

For $j=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{2},\phantom{\rule{0.125em}{0ex}}\mathrm{\dots}\phantom{\rule{0.125em}{0ex}},\phantom{\rule{0.125em}{0ex}}m$, calculate ${v}_{j}=\frac{\mathrm{1}}{m+\mathrm{1}}{\sum}_{i=\mathrm{1}}^{m}\phantom{\rule{0ex}{0ex}}\phantom{\rule{0.25em}{0ex}}\phantom{\rule{0.25em}{0ex}}\mathbf{1}\left({u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}i}\le {u}_{\mathrm{1},\phantom{\rule{0.125em}{0ex}}j},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}i}\le {u}_{\mathrm{3},\phantom{\rule{0.125em}{0ex}}j},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}i}\le {u}_{\mathrm{7},\phantom{\rule{0.125em}{0ex}}j},\phantom{\rule{0.125em}{0ex}}{u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}i}\le {u}_{\mathrm{15},\phantom{\rule{0.125em}{0ex}}j}\right)$.

Calculate ${K}_{C}\left({\mathit{\rho}}_{t}\right)=\frac{\mathrm{1}}{m}{\sum}_{i=\mathrm{1}}^{m}\mathbf{1}\left({v}_{i}\le {\mathit{\rho}}_{t}\right)$.

Calculate ${p}_{t}^{\mathrm{ken}}=\mathrm{1}{K}_{C}\left({\mathit{\rho}}_{t}\right)$.
To calculate the mostlikely design event and confidence interval conditioned on AAR=η, we need to generate numerous multivariate design event samples by the Monte Carlo method. Here, we give the details of generating the design event samples as follows:

Define the total number of design event samples N and the initial number of the design event sample i=0.

Generate a random integer (denoted by t_{r}) among (T_{1}, T_{1}+1, … , T_{2}).

Generate a random sample $\left({z}_{{Q}_{\mathrm{1}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}\right)$ following the multivariate distribution $F\left({z}_{{Q}_{\mathrm{1}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{3}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{7}}},\phantom{\rule{0.125em}{0ex}}{z}_{{V}_{\mathrm{15}}}\left{\mathit{\theta}}_{{t}_{r}}\right.\right)$ with the distribution parameter vector ${\mathit{\theta}}_{{t}_{r}}$.

Calculate the annual exceedance probability for each year throughout the period from T_{1} to T_{2}.

Calculate AAR during the period from T_{1} to T_{2}.

If $\left\mathrm{AAR}\mathit{\eta}\right<\mathit{\epsilon}$ (where it is a very small value, such as 0.0001), $i=i+\mathrm{1}$.

If i<N, repeat steps (2)–(6).
CJ and LX developed the main ideas. CJ and LY implemented the algorithms of the methods. CJ and JD collected the data used in the case study. CJ, LX and CYX prepared the paper.
The authors declare that they have no conflict of interest.
This research was jointly financially supported by the National Natural Science Foundation of China (NSFC grants 51809243 and 51525902), the Fundamental Research Funds for the Central Universities (grant CUG170679), the Research Council of Norway (FRINATEK Project 274310) and the “111 Project” Fund of China (B18037), all of which are greatly appreciated. The authors express their gratitude to the reviewers for their constructive comments and suggestions that have helped improve the paper.
This paper was edited by Carlo De Michele and reviewed by four anonymous referees.
Aas, K., Czado, C., Frigessi, A., and Bakken, H.: Paircopula constructions of multiple dependence, Insurance Math. Econ., 44, 182–198, https://doi.org/10.1016/j.insmatheco.2007.02.001, 2009.
Akaike, H.: A new look at the statistical model identification, IEEE Trans. Autom. Control, 19, 716–723, 1974.
Balistrocchi, M. and Bacchi, B.: Derivation of flood frequency curves through a bivariate rainfall distribution based on copula functions: application to an urban catchment in northern Italy's climate, Hydrol. Res., 48, 749–762, https://doi.org/10.2166/nh.2017.109, 2017.
Bender, J., Wahl, T., and Jensen, J.: Multivariate design in the presence of nonstationarity, J. Hydrol., 514, 123–130, https://doi.org/10.1016/j.jhydrol.2014.04.017, 2014.
Blöschl, G., Hall, J., Parajka, J., Perdigão, R. A. P., Merz, B., Arheimer, B., and Živkovic, N: Changing climate shifts timing of European floods, Science, 357, 588–590, 2017.
Bracken, C., Holman, K. D., Rajagopalan, B., and Moradkhani, H.: A Bayesian hierarchical approach to multivariate nonstationary hydrologic frequency analysis, Water Resour. Res., 54, 243–255, https://doi.org/10.1002/2017WR020403, 2018.
Bücher, A., Kojadinovic, I., Rohmer, T., and Segers, J.: Detecting changes in crosssectional dependence in multivariate time series, J. Multivariate Anal., 132, 111–128, https://doi.org/10.1016/j.jmva.2014.07.012, 2014.
Chow, V. T.: Handbook of Applied Hydrology, McGrawHill, New York, 1964.
Department of Comprehensive Statistics of National Bureau of Statistics: China Compendium of Statistics 1949–2008, China Stat. Press, Beijing, 2010 (in Chinese).
Engeland, K., Wilson, D., Borsányi, P., Roald, L., and Holmqvist, E.: Use of historical data in flood frequency analysis: A case study for four catchments in Norway, Hydrol. Res., 49, 466–486, 2018,
Favre, A. C., El Adlouni, S., Perreault, L., Thiémonge, N., and Bobée, B.: Multivariate hydrological frequency analysis using copulas, Water Resour. Res., 40, W01101, https://doi.org/10.1029/2003WR002456, 2004.
Frank, J. and Massey, J. R.: The KolmogorovSmirnov test for goodness of fit, J. Am. Stat. Assoc., 46, 68–78, 1951.
Hawkes, P. J.: Joint probability analysis for estimation of extremes, J. Hydraul. Res., 46, 246–256, https://doi.org/10.1080/00221686.2008.9521958, 2008.
He, C.: The China Modernization Report 2013, Peking University Press, Beijing, 2014 (in Chinese).
Hurvich, C. M. and Tsai, C. L.: Regression and time series model selection in small samples, Biometrika, 76, 297–307, 1989.
Jiang, C., Xiong, L., Xu, C.Y., and Guo, S.: Bivariate frequency analysis of nonstationary lowflow series based on the timevarying copula, Hydrol. Process., 29, 1521–1534, https://doi.org/10.1002/hyp.10288, 2015a.
Jiang, C., Xiong, L., Wang, D., Liu, P., Guo, S., and Xu, C.Y.: Separating the impacts of climate change and human activities on runoff using the Budykotype equations with timevarying parameters, J. Hydrol., 522, 326–338, https://doi.org/10.1016/j.jhydrol.2014.12.060, 2015b.
Kew, S. F., Selten, F. M., Lenderink, G., and Hazeleger, W.: The simultaneous occurrence of surge and discharge extremes for the Rhine delta, Nat. Hazards Earth Syst. Sci., 13, 2017–2029, https://doi.org/10.5194/nhess1320172013, 2013.
Kobierska, F., Engeland, K., and Thorarinsdottir, T.: Evaluation of design flood estimates – a case study for Norway, Hydrol. Res., 49, 450–465, 2018.
Kojadinovic, I.: npcp: Some nonparametric CUSUM tests for changepoint detection in possibly multivariate observations, R Package Version 0.19, Vienna, Austria, available at: https://cran.rproject.org/web/packages/npcp/npcp.pdf (last access: 20 March 2019), 2017.
Kundzewicz, Z. W., Pińskwar, I., and Brakenridge, G. R.: Changes in river flood hazard in Europe: a review, Hydrol. Res., 49, 294–302, 2018.
Kwon, H.H., Lall, U., and Kim, S.J.: The unusual 2013–2015 drought in South Korea in the context of a multicentury precipitation record: Inferences from a nonstationary, multivariate, Bayesian copula model, Geophys. Res. Lett., 43, 8534–8544, https://doi.org/10.1002/2016GL070270, 2016.
Kyselý, J.: A cautionary note on the use of nonparametric bootstrap for estimating uncertainties in extremevalue models, J. Appl. Meteorol. Clim., 47, 3236–3251, 2009.
Li, T., Guo, S., Liu, Z., Xiong, L., and Yin, J.: Bivariate design flood quantile selection using copulas, Hydrol. Res., 48, 997–1013, 2017.
Liang, Z., Hu, Y., Huang, H., Wang, J., and Li, B.: Study on the estimation of design value under nonstationary environment, SouthtoNorth Water Transfers, Water Sci. Technol., 14, 50–53, 2016 (in Chinese).
López, J. and Francés, F.: Nonstationary flood frequency analysis in continental Spanish rivers, using climate and reservoir indices as external covariates, Hydrol. Earth Syst. Sci., 17, 3189–3203, https://doi.org/10.5194/hess1731892013, 2013.
Loveridge, M., Rahman, A., and Hill, P.: Applicability of a physically based soil water model (SWMOD) in design flood estimation in eastern Australia, Hydrol. Res., 48, 1652–1665, 2017.
Milly, P., Betancourt, J., Falkenmark, M., Hirsch, R., Kundzewicz, Z., Lettenmaier, D., and Stouffer, R.: Climate change – Stationarity is dead: Whither water management?, Science, 319, 573–574, https://doi.org/10.1126/science.1151915, 2008.
Ministry of Water Resources of People's Republic of China: Design Criterion of Reservoir Management, Chin. Water Resour. and Hydropower Press, Beijing, 1996 (in Chinese).
Niederreiter, H.: QuasiMonte Carlo methods and pseudorandom numbers, B. Am. Math. Soc., 197, 957–1041, 1978.
Obeysekera, J. and Salas, J.: Quantifying the uncertainty of design floods under nonstationary conditions, J. Hydrol. Eng., 19, 1438–1446, https://doi.org/10.1061/(ASCE)HE.19435584.0000931, 2014.
Obeysekera, J. and Salas, J.: Frequency of recurrent extremes under nonstationarity, J. Hydrol. Eng., 21, 04016005, https://doi.org/10.1061/(ASCE)HE.19435584.0001339, 2016.
Olsen, J. R., Lambert, J. H., and Haimes, Y. Y.: Risk of extreme events under nonstationarity conditions, Risk Anal., 18, 497–510, https://doi.org/10.1111/j.15396924.1998.tb00364.x, 1998.
Parey, S., Hoang, T. T. H., and DacunhaCastelle, D.: Different ways to compute temperature return levels in the climate change context, Environmetrics, 21, 698–718, https://doi.org/10.1002/env.1060, 2010.
Qi, W. and Liu, J.: A nonstationary costbenefit based bivariate extreme flood estimation approach, J. Hydrol., 557, 589–599, https://doi.org/10.1016/j.jhydrol.2017.12.045, 2017.
Quessy, J., Saïd, M., and Favre, A. C.: Multivariate Kendall's tau for changepoint detection in copulas, Can. J. Stat., 41, 65–82, https://doi.org/10.1002/cjs.11150, 2013.
Read, L. K. and Vogel, R. M.: Reliability, return periods, and risk under nonstationarity, Water Resour. Res., 51, 6381–6398, https://doi.org/10.1002/2015WR017089, 2015.
Read, L. K. and Vogel, R. M.: Hazard function analysis for flood planning under nonstationarity, Water Resour. Res., 52, 4116–4131, https://doi.org/10.1002/2015WR018370, 2016.
Requena, A. I., Mediero, L., and Garrote, L.: A bivariate return period based on copulas for hydrologic dam design: accounting for reservoir routing in risk estimation, Hydrol. Earth Syst. Sci., 17, 3023–3038, https://doi.org/10.5194/hess1730232013, 2013.
Rootzén, H. and Katz, R. W.: Design Life Level: Quantifying risk in a changing climate, Water Resour. Res., 49, 5964–5972, https://doi.org/10.1002/wrcr.20425, 2013.
Rosner, A., Vogel, R. M., and Kirshen, P. H.: A riskbased approach to flood management decisions in a nonstationary world, Water Resour. Res., 50, 1928–1942, https://doi.org/10.1002/2013WR014561, 2014.
Salas, J. D. and Obeysekera, J.: Revisiting the concepts of return period and risk for nonstationary hydrologic extreme events, J. Hydrol. Eng., 19, 554–568, https://doi.org/10.1061/(ASCE)HE.19435584.0000820, 2014.
Salvadori, G. and De Michele, C.: Frequency analysis via copulas: theoretical aspects and applications to hydrological events, Water Resour. Res., 40, W12511, https://doi.org/10.1029/2004WR003133, 2004.
Salvadori, G. and De Michele, C.: Multivariate multiparameter extreme value models and return periods: A copula approach, Water Resour. Res., 46, W10501, https://doi.org/10.1029/2009WR009040, 2010.
Salvadori, G., De Michele, C., Kottegoda, N. T., and Rosso, R.: Extremes in Nature: An Approach Using Copulas, Springer, Dordrecht, the Netherlands, 2007.
Salvadori, G., De Michele, C., and Durante, F.: On the return period and design in a multivariate framework, Hydrol. Earth Syst. Sci., 15, 3293–3305, https://doi.org/10.5194/hess1532932011, 2011.
Salvadori, G., Durante, F., and De Michele, C.: Multivariate return period calculation via survival functions, Water Resour. Res., 49, 2308–2311, https://doi.org/10.1002/wrcr.20204, 2013.
Salvadori, G., Durante, F., Tomasicchio, G. R., and D'Alessandro, F.: Practical guidelines for the multivariate assessment of the structural risk in coastal and offshore engineering, Coastal Eng., 95, 77–83, https://doi.org/10.1016/j.coastaleng.2014.09.007, 2015.
Salvadori, G., Durante, F., De Michele, C., Bernardi, M., and Petrella, L.: A multivariate Copulabased framework for dealing with Hazard Scenarios and Failure Probabilities, Water Resour. Res., 52, 3701–3721, https://doi.org/10.1002/2015WR017225, 2016.
Salvadori, G., Durante, F., Michele, C. D., and Bernardi, M.: Hazard assessment under multivariate distributional changepoints: Guidelines and a flood case study, Water, 10, 751–765, https://doi.org/10.3390/w10060751, 2018.
Sarhadi, A., Burn, D. H., Ausín, M. C., and Wiper, M. P.: Timevarying nonstationary multivariate risk analysis using a dynamic Bayesian copula, Water Resour. Res., 52, 2327–2349, https://doi.org/10.1002/2015WR018525, 2016.
Serinaldi, F.: Dismissing return periods!, Stoch. Env. Res. Risk. A., 29, 1179–1189, https://doi.org/10.1007/s0047701409161, 2015.
Serinaldi, F. and Kilsby, C. G.: Stationarity is undead: Uncertainty dominates the distribution of extremes, Adv. Water Resour., 77, 17–36, https://doi.org/10.1016/j.advwatres.2014.12.013, 2015.
Shafaei, M., FakheriFard, A., Dinpashoh, Y., Mirabbasi, R., and De Michele, C.: Modeling flood event characteristics using Dvine structures, Theor. Appl. Climatol., 130, 713–724, https://doi.org/10.1007/s007040161911x, 2017.
Sklar, M.: Fonctions de Répartition a n Dimensions et Leurs Marges, 8 pp., Univ. Paris, Paris, 1959.
Strupczewski, W. G., Singh, V. P., and Feluch, W.: Nonstationary approach to atsite flood frequency modeling I. Maximum likelihood estimation, J. Hydrol., 248, 123–142, https://doi.org/10.1016/S00221694(01)003973, 2001.
Vandenberghe, S., Verhoest, N. E. C., Onof, C., and De Baets, B.: A comparative copula – based bivariate frequency analysis of observed and simulated storm events: A case study on Bartlett – Lewis modeled rainfall, Water Resour. Res., 47, W07529, https://doi.org/10.1029/2009WR008388, 2011.
Vezzoli, R., Salvadori, G., and De Michele, C.: A distributional multivariate approach for assessing performance of climatehydrology models, Sci. Rep., 7, 12071, https://doi.org/10.1038/s41598017123431, 2017.
Villarini, G., Serinaldi, F., Smith, J. A., and Krajewski, W. F.: On the stationarity of annual flood peaks in the Continental United States during the 20th Century, Water Resour. Res., 45, W08417, https://doi.org/10.1029/2008WR007645, 2009.
Vogel, R. M.: Reliability indices for water supply systems, J. Water Res. Pl., 113, 563–579, https://doi.org/10.1061/(ASCE)07339496(1987)113:4(563), 1987.
Vogel, R. M., Yaindl, C., and Walter, M.: Nonstationarity: Flood magnification and recurrence reduction factors in the United States, J. Am. Water Resour. As., 47, 464–474, https://doi.org/10.1111/j.17521688.2011.00541.x, 2011.
Volpi, E. and Fiori, A.: Design event selection in bivariate hydrological frequency analysis, Hydrolog. Sci. J., 57, 1506–1515, https://doi.org/10.1080/02626667.2012.726357, 2012.
Xiao, Y., Guo, S., Liu, P., Yan, B., and Chen, L.: Design flood hydrograph based on multicharacteristic synthesis index method, J. Hydrol. Eng., 14, 1359–1364, https://doi.org/10.1061/(ASCE)10840699(2009)4:12(1359), 2009.
Xiong, L. and Guo, S.: Trend test and changepoint detection for the annual discharge series of the Yangtze River at the Yichang hydrological station, Hydrolog. Sci. J., 49, 99–112, https://doi.org/10.1623/hysj.49.1.99.53998, 2004.
Xiong, L., Jiang, C., Xu, C.Y., Yu, K.X., and Guo, S.: A framework of changepoint detection for multivariate hydrological series, Water Resour. Res., 51, 8198–8217, https://doi.org/10.1002/2015WR017677, 2015.
Xu, B., Xie, P., Tan, Y., Li, X., and Liu, Y.: Analysis of flood returning to main channel influence on the flood control ability of Xijiang River, Journal of Hydroelectric Engineering, 33, 65–72, 2014 (in Chinese).
Yan, L., Xiong, L., Guo, S., Xu, C.Y., Xia, J., and Du, T.: Comparison of four nonstationary hydrologic design methods for changing environment, J. Hydrol., 551, 132–150, https://doi.org/10.1016/j.jhydrol.2017.06.001, 2017.
Yang, T., Shao, Q., Hao, Z., Chen, Xi., Zhang, Z., Xu, C.Y., and Sun, L.: Regional frequency analysis and spatiotemporal pattern characterization of rainfall extremes in the Pearl River Basin, China, J. Hydrol., 380, 386–405, https://doi.org/10.1016/j.jhydrol.2009.11.013, 2010.
Yin, J., Guo, S., Liu, Z., Chen, K., Chang, F., and Xiong, F.: Bivariate seasonal design flood estimation based on copulas, J. Hydrol. Eng., 22, 05017028, https://doi.org/10.1061/(ASCE)HE.19435584.0001594, 2017.
Zhang, L. and Singh, V. P.: Trivariate flood frequency analysis using the Gumbel–Hougaard copula, J. Hydrol. Eng., 12, 431–439, https://doi.org/10.1061/(ASCE)10840699(2007)12:4(431), 2007.
Zheng, F., Westra, S., and Sisson, S. A.: Quantifying the dependence between extreme rainfall and storm surge in the coastal zone, J. Hydrol., 505, 172–187, https://doi.org/10.1016/j.jhydrol.2013.09.054, 2013.
Zheng, F., Westra, S., Leonard, M., and Sisson, S. A.: Modeling dependence between extreme rainfall and storm surge to estimate coastal flooding risk, Water Resour. Res., 50, 2050–2071, https://doi.org/10.1002/2013WR014616, 2014.
Zheng, F., Leonard, M., and Westra, S.: Efficient joint probability analysis of flood risk, J. Hydroinform., 17, 584–597, 2015.
Zheng, F., Leonard, M., and Westra, S.: Application of the design variable method to estimate coastal flood risk, J. Flood Risk Manag., 10, 522–534, https://doi.org/10.1111/jfr3.12180, 2017.
Zheng, F., Tao, R., Maier, H. R., See, L., Savic, D., Zhang, T., Chen, O., Assumpção, T. H., Yang, P., Heidari, B., Rickermann, J., Minsker, B., Bi, W., Cai, X., Solomatine, D., and Popescu, I.: Crowdsourcing methods for data collection in geophysics: State of the art, issues, and future directions, Rev. Geophys., 56, 698–740, https://doi.org/10.1029/2018RG000616, 2018.
 Abstract
 Introduction
 Study area and data
 Methods
 Results
 Conclusion and remarks
 Data availability
 Appendix A: Calculating multivariate exceedance probabilities
 Appendix B: Generating the multivariate design event samples (formulated by Eq. 17 in the paper)
 Author contributions
 Competing interests
 Acknowledgements
 Review statement
 References
 Abstract
 Introduction
 Study area and data
 Methods
 Results
 Conclusion and remarks
 Data availability
 Appendix A: Calculating multivariate exceedance probabilities
 Appendix B: Generating the multivariate design event samples (formulated by Eq. 17 in the paper)
 Author contributions
 Competing interests
 Acknowledgements
 Review statement
 References