Multivariate hydrologic design methods under nonstationary conditions and application to engineering practice

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 inter-arrival 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 return-period concept cannot apply to engineering practice under nonstationary conditions, since by such a definition, a given multivariate flood event would correspond to a time-varying 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 C-vine copula, allowing for time-varying marginal distributions and a time-varying dependence structure. Both the most-likely 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.

1 Introduction

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 inter-arrival 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 inter-arrival 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 time-varying and non-unique return period. Consequently, the traditional return-period-based 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 risk-based or reliability-based 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 most-likely 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 most-likely 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 most-likely 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 (C-vine) 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 four-dimensional (4-D) multivariate flood series, including the annual maximum daily discharge, annual maximum 3-day flood volume, annual maximum 7-day flood volume and annual maximum 15-day 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.

2 Study area and data

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², 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.

Figure 1 Map of the Xijiang River basin (above the Dahuangjiangkou gauge).

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₁), annual maximum 3-day flood volume (V₃), annual maximum 7-day flood volume (V₇) and annual maximum 15-day flood volume (V₁₅) of the Xijiang River were defined as the multivariate flood series $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},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², 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⁹ m³ and 1.64×10⁹ m³ and catchment areas of 98 500 and 9600 km², 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.

Figure 2 Evolution of the urban population and reservoir index during both the observation and design life periods.

Table 1 Reservoir information for the Xijiang River basin.

3 Methods

The present study included the following methods: (1) the nonstationary multivariate flood distribution based on a dynamic C-vine copula, allowing for both time-varying marginal distributions and a time-varying 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₁), annual maximum 3-day flood volume (V₃), annual maximum 7-day flood volume (V₇) and annual maximum 15-day flood volume (V₁₅) 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 4-D flood series $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ at time t measured by years ($t=\mathrm{1},\mathrm{2},\mathrm{\dots },n$, and n is the length of the flood series) can be formulated through a copula C(⋅) as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }F\left(q_{\mathrm{1},t},v_{\mathrm{3},t},v_{\mathrm{7},t},v_{\mathrm{15},t}\left|{\mathit{\theta }}_t\right.\right)\\ {\displaystyle }& {\displaystyle }=C\left[F_{\mathrm{1}}\left(q_{\mathrm{1},t}|{\mathit{\theta }}_{\mathrm{1},t}\right),F_{\mathrm{3}}\left(v_{\mathrm{3},t}|{\mathit{\theta }}_{\mathrm{3},t}\right),F_{\mathrm{7}}\left(v_{\mathrm{7},t}|{\mathit{\theta }}_{\mathrm{7},t}\right)\right.,\\ {\displaystyle }& {\displaystyle }\left.F_{\mathrm{15}}\left(v_{\mathrm{15},t}|{\mathit{\theta }}_{\mathrm{15},t}\right)|{\mathit{\theta }}_{c,t}\right]\\ \text{(1)}& {\displaystyle }& {\displaystyle }=C\left(u_{\mathrm{1},t},u_{\mathrm{3},t},u_{\mathrm{7},t},u_{\mathrm{15},t}|{\mathit{\theta }}_{c,t}\right),\end{array}$$

where $F_{\mathrm{1}}\left(q_{\mathrm{1},t}|{\mathit{\theta }}_{\mathrm{1},t}\right)$, $F_{\mathrm{3}}\left(v_{\mathrm{3},t}|{\mathit{\theta }}_{\mathrm{3},t}\right)$, $F_{\mathrm{7}}\left(v_{\mathrm{7},t}|{\mathit{\theta }}_{\mathrm{7},t}\right)$ and $F_{\mathrm{15}}\left(v_{\mathrm{15},t}|{\mathit{\theta }}_{\mathrm{15},t}\right)$ denote the marginal distributions for Q₁, V₃, V₇ and V₁₅, respectively; u_1, t, u_3, t, u_7, t and u_15, t are the marginal probabilities of Q₁, V₃, V₇ and V₁₅, 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},t},{\mathit{\theta }}_{\mathrm{3},t},{\mathit{\theta }}_{\mathrm{7},t},{\mathit{\theta }}_{\mathrm{15},t},{\mathit{\theta }}_{c,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}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ defined by Eq. (1), the corresponding density function can be written as:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }f\left(q_{\mathrm{1},t},v_{\mathrm{3},t},v_{\mathrm{7},t},v_{\mathrm{15},t}\left|{\mathit{\theta }}_t\right.\right)\\ {\displaystyle }& {\displaystyle }=c\left[F_{\mathrm{1}}\left(q_{\mathrm{1},t}|{\mathit{\theta }}_{\mathrm{1},t}\right),F_{\mathrm{3}}\left(v_{\mathrm{3},t}|{\mathit{\theta }}_{\mathrm{3},t}\right),F_{\mathrm{7}}\left(v_{\mathrm{7},t}|{\mathit{\theta }}_{\mathrm{7},t}\right),\right.\\ {\displaystyle }& {\displaystyle }\left.F_{\mathrm{15}}\left(v_{\mathrm{15},t}|{\mathit{\theta }}_{\mathrm{15},t}\right)|{\mathit{\theta }}_{c,t}\right]\cdot f_{\mathrm{1}}\left(q_{\mathrm{1},t}|{\mathit{\theta }}_{\mathrm{1},t}\right)\cdot f_{\mathrm{3}}\left(v_{\mathrm{3},t}|{\mathit{\theta }}_{\mathrm{3},t}\right)\\ {\displaystyle }& {\displaystyle }\cdot f_{\mathrm{7}}\left(v_{\mathrm{7},t}|{\mathit{\theta }}_{\mathrm{7},t}\right)\cdot f_{\mathrm{15}}\left(v_{\mathrm{15},t}|{\mathit{\theta }}_{\mathrm{15},t}\right)\\ {\displaystyle }& {\displaystyle }=c\left(u_{\mathrm{1},t},u_{\mathrm{3},t},u_{\mathrm{7},t},u_{\mathrm{15},t}|{\mathit{\theta }}_{c,t}\right)\cdot f_{\mathrm{1}}\left(q_{\mathrm{1},t}|{\mathit{\theta }}_{\mathrm{1},t}\right)\\ \text{(2)}& {\displaystyle }& {\displaystyle }\cdot f_{\mathrm{3}}\left(v_{\mathrm{3},t}|{\mathit{\theta }}_{\mathrm{3},t}\right)\cdot f_{\mathrm{7}}\left(v_{\mathrm{7},t}|{\mathit{\theta }}_{\mathrm{7},t}\right)\cdot f_{\mathrm{15}}\left(v_{\mathrm{15},t}|{\mathit{\theta }}_{\mathrm{15},t}\right),\end{array}$$

where $f_{\mathrm{1}}\left(q_{\mathrm{1},t}|{\mathit{\theta }}_{\mathrm{1},t}\right)$, $f_{\mathrm{3}}\left(v_{\mathrm{3},t}|{\mathit{\theta }}_{\mathrm{3},t}\right)$, $f_{\mathrm{7}}\left(v_{\mathrm{7},t}|{\mathit{\theta }}_{\mathrm{7},t}\right)$ and $f_{\mathrm{15}}\left(v_{\mathrm{15},t}|{\mathit{\theta }}_{\mathrm{15},t}\right)$ are the density functions of the marginal distributions for Q₁, V₃, V₇ and V₁₅, respectively, and c(⋅) denotes the density function of copula C(⋅). As shown by Eq. (2), the multivariate distribution of $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ can be separated into two modules, including the marginal distributions, i.e. $f_{\mathrm{1}}\left(q_{\mathrm{1},t}|{\mathit{\theta }}_{\mathrm{1},t}\right)$, $f_{\mathrm{3}}\left(v_{\mathrm{3},t}|{\mathit{\theta }}_{\mathrm{3},t}\right)$, $f_{\mathrm{7}}\left(v_{\mathrm{7},t}|{\mathit{\theta }}_{\mathrm{7},t}\right)$ and $f_{\mathrm{15}}\left(v_{\mathrm{15},t}|{\mathit{\theta }}_{\mathrm{15},t}\right)$, as well as the dependence structure expressed by the copula density function $c\left(u_{\mathrm{1},t},u_{\mathrm{3},t},u_{\mathrm{7},t},u_{\mathrm{15},t}|{\mathit{\theta }}_{c,t}\right)$. Under nonstationary conditions, both the margins and dependence structure of $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ can vary over time t.

3.1.1 Nonstationary marginal distributions based on the time-varying moment model

The time-varying 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}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ were constructed by the time-varying 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.

$$\begin{array}{}\text{(3)}& \begin{array}{ll}\mathrm{Linear}:& {\mathit{\mu }}_t={\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{1}}{\mathrm{Pop}}_t+{\mathit{\alpha }}_{\mathrm{2}}{\mathrm{RI}}_t.\\ \mathrm{Exponential}:& {\mathit{\mu }}_t=\exp \left({\mathit{\alpha }}_{\mathrm{0}}+{\mathit{\alpha }}_{\mathrm{1}}{\mathrm{Pop}}_t+{\mathit{\alpha }}_{\mathrm{2}}{\mathrm{RI}}_t\right).\end{array}\end{array}$$

Here α₀, α₁ and α₂ 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 time-varying 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 C-vine copula

After estimating the marginal distributions, the nonstationary dependence structure of $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ as formulated by the copula density function $c\left(u_{\mathrm{1},t},u_{\mathrm{3},t},u_{\mathrm{7},t},u_{\mathrm{15},t}|{\mathit{\theta }}_{c,t}\right)$ was constructed. Given that most applied copula functions are for bivariate random variables, $c\left(u_{\mathrm{1},t},u_{\mathrm{3},t},u_{\mathrm{7},t},u_{\mathrm{15},t}|{\mathit{\theta }}_{c,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}},V_{\mathrm{3}},V_{\mathrm{7}},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₁) 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}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$, with Q₁ elected as the key variable. Thus, the density function $c\left(u_{\mathrm{1},t},u_{\mathrm{3},t},u_{\mathrm{7},t},u_{\mathrm{15},t}|{\mathit{\theta }}_{c,t}\right)$ can be decomposed into six bivariate pair copulas as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }c\left(u_{\mathrm{1},t},u_{\mathrm{3},t},u_{\mathrm{7},t},u_{\mathrm{15},t}|{\mathit{\theta }}_{c,t}\right)=c_{\mathrm{13}}\left(u_{\mathrm{1},t},u_{\mathrm{3},t}|{\mathit{\theta }}_{\mathrm{13},t}\right)\\ {\displaystyle }& {\displaystyle }\cdot c_{\mathrm{17}}\left(u_{\mathrm{1},t},u_{\mathrm{7},t}\left|{\mathit{\theta }}_{\mathrm{17},t}\right.\right)\cdot c_{\mathrm{115}}\left(u_{\mathrm{1},t},u_{\mathrm{15},t}\left|{\mathit{\theta }}_{\mathrm{115},t}\right.\right)\\ {\displaystyle }& {\displaystyle }\cdot c_{\mathrm{37}|\mathrm{1}}\left[F\left(u_{\mathrm{3},t}|u_{\mathrm{1},t}\right),F\left(u_{\mathrm{7},t}|u_{\mathrm{1},t}\right)\left|{\mathit{\theta }}_{\mathrm{37}|\mathrm{1},t}\right.\right]\\ {\displaystyle }& {\displaystyle }\cdot c_{\mathrm{315}|\mathrm{1}}\left[F\left(u_{\mathrm{3},t}|u_{\mathrm{1},t}\right),F\left(u_{\mathrm{15},t}|u_{\mathrm{1},t}\right)\left|{\mathit{\theta }}_{\mathrm{315}|\mathrm{1},t}\right.\right]\\ {\displaystyle }& {\displaystyle }\cdot c_{\mathrm{715}|\mathrm{13}}\left[F\left(u_{\mathrm{7},t}|u_{\mathrm{1},t},u_{\mathrm{3},t}\right),\right.\\ \text{(4)}& {\displaystyle }& {\displaystyle }F\left(u_{\mathrm{15},t}|u_{\mathrm{1},t},u_{\mathrm{3},t}\right)\left|{\mathit{\theta }}_{\mathrm{715}|\mathrm{13},t}\right],\end{array}$$

where ${\mathit{\theta }}_{c,t}=\left({\mathit{\theta }}_{\mathrm{13},t},{\mathit{\theta }}_{\mathrm{17},t},{\mathit{\theta }}_{\mathrm{115},t},{\mathit{\theta }}_{\mathrm{37}|\mathrm{1},t},{\mathit{\theta }}_{\mathrm{315}|\mathrm{1},t},{\mathit{\theta }}_{\mathrm{715}|\mathrm{13},t}\right)$ is the parameter vector in the C-vine copula, and

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }F\left(u_{\mathrm{3},t}|u_{\mathrm{1},t}\right)={\displaystyle \frac{\partial C_{\mathrm{13}}\left(u_{\mathrm{1},t},u_{\mathrm{3},t}|{\mathit{\theta }}_{\mathrm{13},t}\right)}{\partial u_{\mathrm{1},t}}},\\ {\displaystyle }& {\displaystyle }F\left(u_{\mathrm{7},t}|u_{\mathrm{1},t}\right)={\displaystyle \frac{\partial C_{\mathrm{17}}\left(u_{\mathrm{1},t},u_{\mathrm{7},t}\left|{\mathit{\theta }}_{\mathrm{17},t}\right.\right)}{\partial u_{\mathrm{1},t}}},\\ {\displaystyle }& {\displaystyle }F\left(u_{\mathrm{15},t}|u_{\mathrm{1},t}\right)={\displaystyle \frac{\partial C_{\mathrm{115}}\left(u_{\mathrm{1},t},u_{\mathrm{15},t}\left|{\mathit{\theta }}_{\mathrm{115},t}\right.\right)}{\partial u_{\mathrm{1},t}}},\\ {\displaystyle }& {\displaystyle }F\left(u_{\mathrm{7},t}|u_{\mathrm{1},t},u_{\mathrm{3},t}\right)=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{\partial c_{\mathrm{37}|\mathrm{1}}\left[F\left(u_{\mathrm{3},t}|u_{\mathrm{1},t}\right),F\left(u_{\mathrm{7},t}|u_{\mathrm{1},t}\right)\left|{\mathit{\theta }}_{\mathrm{37}|\mathrm{1},t}\right.\right]}{\partial F\left(u_{\mathrm{3},t}|u_{\mathrm{1},t}\right)}},\\ {\displaystyle }& {\displaystyle }F\left(u_{\mathrm{15},t}|u_{\mathrm{1},t},u_{\mathrm{3},t}\right)=\\ \text{(5)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\partial C_{\mathrm{315}|\mathrm{1}}\left[F\left(u_{\mathrm{3},t}|u_{\mathrm{1},t}\right),F\left(u_{\mathrm{15},t}|u_{\mathrm{1},t}\right)\left|{\mathit{\theta }}_{\mathrm{315}|\mathrm{1},t}\right.\right]}{\partial F\left(u_{\mathrm{3},t}|u_{\mathrm{1},t}\right)}}.\end{array}$$

Figure 3 shows the schematic decomposition of the 4-D C-vine copula as expressed by Eq. (4). It is evident that the hierarchical structure of the 4-D C-vine 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},t}\right.\right)$, $c_{\mathrm{17}}\left(\cdot \left|{\mathit{\theta }}_{\mathrm{17},t}\right.\right)$ and $c_{\mathrm{115}}\left(\cdot \left|{\mathit{\theta }}_{\mathrm{115},t}\right.\right)$, which directly act on the marginal probabilities and describe the bivariate dependencies between the key variable Q₁ and the other three variables, i.e. V₃, V₇ and V₁₅. The second tree (T2) includes two bivariate pair copulas $c_{\mathrm{37}|\mathrm{1}}\left(\cdot \left|{\mathit{\theta }}_{\mathrm{37}|\mathrm{1},t}\right.\right)$ and $c_{\mathrm{315}|\mathrm{1}}\left(\cdot \left|{\mathit{\theta }}_{\mathrm{315}|\mathrm{1},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},t}\right.\right)$ acting on conditional distribution functions with both u_1, t and u_3, t as the conditioning variables.

Figure 3 Decomposition of the four-dimensional C-vine copula.

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 extreme-value copula widely used in hydrology, accounts for the upper tail dependence and is well-suited 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 C-vine copula formulated by Eq. (4). The bivariate Gumbel–Hougaard copula is expressed as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }C\left(u,v\right)=\exp \left\{-{\left[{\left(-\ln u\right)}^{{\mathit{\theta }}_c}+{\left(-\ln v\right)}^{{\mathit{\theta }}_c}\right]}^{\mathrm{1}/{\mathit{\theta }}_c}\right\},\\ \text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{\theta }}_c\in [\mathrm{1},\mathrm{\infty }),\end{array}$$

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}},V_{\mathrm{3}},V_{\mathrm{7}},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 time-varying copula parameters and generally formulated as follows.

$$\begin{array}{}\text{(7)}& \begin{array}{ll}\mathrm{Linear}:& {\mathit{\theta }}_{c,t}={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{\mathrm{Pop}}_t+{\mathit{\beta }}_{\mathrm{2}}{\mathrm{RI}}_t.\\ \mathrm{Exponential}:& {\mathit{\theta }}_{c,t}=\mathrm{1}+\exp \left({\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{1}}{\mathrm{Pop}}_t+{\mathit{\beta }}_{\mathrm{2}}{\mathrm{RI}}_t\right).\end{array}\end{array}$$

Here β₀, β₁ and β₂ 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},t}$, ${\mathit{\theta }}_{\mathrm{315}|\mathrm{1},t}$ and ${\mathit{\theta }}_{\mathrm{715}|\mathrm{13},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}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ following a given C-vine copula, the PIT converts the dependent flood variables into a new set of variables that are independent and uniformly distributed on [0, 1]⁴. 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 arithmetic-average 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 arithmetic-average AAR is not capable of this function. The present study employed the geometric-average method to calculate AAR, which is dominated more by the minimum than arithmetic average and is theoretically able to yield safer design values. The geometric-average 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}},v_{\mathrm{3}},v_{\mathrm{7}},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}},v_{\mathrm{3}},v_{\mathrm{7}},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}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\right)$ was calculated by the geometric-average method as follows:

$$\begin{array}{}\text{(8)}& \mathrm{AAR}\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\right)={\left[\prod _{t=T_{\mathrm{1}}}^{T_{\mathrm{2}}}\left(\mathrm{1}-p_t\right)\right]}^{\frac{\mathrm{1}}{T_{\mathrm{2}}-T_{\mathrm{1}}+\mathrm{1}}},\end{array}$$

where T₁ and T₂ 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}},v_{\mathrm{3}},v_{\mathrm{7}},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}},v_{\mathrm{3}},v_{\mathrm{7}},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}},v_{\mathrm{3}},v_{\mathrm{7}},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:

$$\begin{array}{ll}{\displaystyle }{p}_t^{\mathrm{or}}& {\displaystyle }=P\left(Q_{\mathrm{1},t}\ge q_{\mathrm{1}}\vee V_{\mathrm{3},t}\ge v_{\mathrm{3}}\vee V_{\mathrm{7},t}\ge v_{\mathrm{7}}\vee V_{\mathrm{15},t}\ge v_{\mathrm{15}}\right)\\ \text{(9)}& {\displaystyle }& {\displaystyle }=\mathrm{1}-F\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\left|{\mathit{\theta }}_t\right.\right),\end{array}$$

where “∨” stands for the OR operator and F(⋅|θ_t) is defined in Eq. (1).

The AND case for $\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{p}_t^{\mathrm{and}}=P\left(Q_{\mathrm{1},t}>q_{\mathrm{1}}\wedge V_{\mathrm{3},t}>v_{\mathrm{3}}\wedge V_{\mathrm{7},t}>v_{\mathrm{7}}\wedge V_{\mathrm{15},t}>v_{\mathrm{15}}\right)\\ {\displaystyle }& {\displaystyle }=\int \int \int {\int }_{{\mathrm{\Omega }}^{\mathrm{and}}}f\left(Q_{\mathrm{1},t},V_{\mathrm{3},t},V_{\mathrm{7},t},V_{\mathrm{15},t}\left|{\mathit{\theta }}_t\right.\right)\mathrm{d}{Q}_{\mathrm{1},t}\\ {\displaystyle }& {\displaystyle }\cdot \mathrm{d}{V}_{\mathrm{3},t}\cdot \mathrm{d}{V}_{\mathrm{7},t}\cdot \mathrm{d}{V}_{\mathrm{15},t}\\ {\displaystyle }& {\displaystyle }{\mathrm{\Omega }}^{\mathrm{and}}:q_{\mathrm{1}}<Q_{\mathrm{1},t}<\mathrm{\infty },v_{\mathrm{3}}<V_{\mathrm{3},t}<\mathrm{\infty },v_{\mathrm{7}}<V_{\mathrm{7},t}<\mathrm{\infty },\\ \text{(10)}& {\displaystyle }& {\displaystyle }{v}_{\mathrm{15}}<V_{\mathrm{15},t}<\mathrm{\infty },\end{array}$$

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}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\right)$ was first transformed into a univariate representation via Kendall's distribution function K_C(⋅) as follows:

$$\begin{array}{ll}{\displaystyle }{K}_C\left({\mathit{\rho }}_t\right)& {\displaystyle }=P\left[C\left(U_{\mathrm{1},t},U_{\mathrm{3},t},U_{\mathrm{7},t},U_{\mathrm{15},t}\left|{\mathit{\theta }}_{c,t}\right.\right)\le {\mathit{\rho }}_t\right]\\ \text{(11)}& {\displaystyle }& {\displaystyle }=P\left[F\left(Q_{\mathrm{1},t},V_{\mathrm{3},t},V_{\mathrm{7},t},V_{\mathrm{15},t}\left|{\mathit{\theta }}_t\right.\right)\le {\mathit{\rho }}_t\right],\end{array}$$

where ${\mathit{\rho }}_t=F\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\left|{\mathit{\theta }}_t\right.\right)$ is the probability level corresponding to the given flood event $\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\right)$. The exceedance probability $p_t^{\mathrm{ken}}$ in the Kendall case at time t was expressed as follows:

$$\begin{array}{}\text{(12)}& p_t^{\mathrm{ken}}=\mathrm{1}-K_C\left({\mathit{\rho }}_t\right).\end{array}$$

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 Most-likely design event and confidence interval for multivariate hydrologic design

The methods identifying both the most-likely design event, denoted by $\left(z_{Q_{\mathrm{1}}}^{*},z_{V_{\mathrm{3}}}^{*},z_{V_{\mathrm{7}}}^{*},z_{V_{\mathrm{15}}}^{*}\right)$, and the confidence interval for the multivariate hydrologic design $\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\right)$ given AAR=η are introduced below. The average annual probability density, denoted by g(⋅), of $\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\right)$ over the entire design life period from T₁ to T₂, was expressed as follows:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }g\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\right)=\\ \text{(13)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{2}}-T_{\mathrm{1}}+\mathrm{1}}}\sum _{t=T_{\mathrm{1}}}^{T_{\mathrm{2}}}f\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\left|{\mathit{\theta }}_t\right.\right).\end{array}$$

The probability distribution function for AAR≤η can be written as

$$\begin{array}{ll}{\displaystyle }\mathrm{\Phi }\left(\mathit{\eta }\right)& {\displaystyle }=\underset{\mathrm{\Omega }:\mathrm{AAR}\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\right)\le \mathit{\eta }}{\int \int \int \int }g\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\right)\\ \text{(14)}& {\displaystyle }& {\displaystyle }\mathrm{d}{q}_{\mathrm{1}}\mathrm{d}{v}_{\mathrm{3}}\mathrm{d}{v}_{\mathrm{7}}\mathrm{d}{v}_{\mathrm{15}}.\end{array}$$

By denoting the density function of Φ(η) as ϕ(η), the probability density of $\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\right)$ conditioned on AAR=η can be expressed as follows:

$$\begin{array}{ll}{\displaystyle }{g}_{|\mathrm{AAR}}& {\displaystyle }=\mathit{\eta }\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\right)\\ \text{(15)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{g\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\right)}{\mathit{\varphi }\left(\mathit{\eta }\right)}}.\end{array}$$

The most-likely design event conditioned on AAR=η was theoretically written as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left(z_{Q_{\mathrm{1}}}^{*},z_{V_{\mathrm{3}}}^{*},z_{V_{\mathrm{7}}}^{*},z_{V_{\mathrm{15}}}^{*}\right)=\\ \text{(16)}& {\displaystyle }& {\displaystyle }\arg max{g}_{|\mathrm{AAR}=\mathit{\eta }}\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\right).\end{array}$$

Unfortunately, the analytical solutions of both the most-likely design event $\left(z_{Q_{\mathrm{1}}}^{*},z_{V_{\mathrm{3}}}^{*},z_{V_{\mathrm{7}}}^{*},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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\left(z_{Q_{\mathrm{1}}}^{\mathrm{1}},z_{V_{\mathrm{3}}}^{\mathrm{1}},z_{V_{\mathrm{7}}}^{\mathrm{1}},z_{V_{\mathrm{15}}}^{\mathrm{1}}\right),\left(z_{Q_{\mathrm{1}}}^{\mathrm{2}},z_{V_{\mathrm{3}}}^{\mathrm{2}},z_{V_{\mathrm{7}}}^{\mathrm{2}},z_{V_{\mathrm{15}}}^{\mathrm{2}}\right),\mathrm{\dots },\\ \text{(17)}& {\displaystyle }& {\displaystyle }\left(z_{Q_{\mathrm{1}}}^{Nc},z_{V_{\mathrm{3}}}^{Nc},z_{V_{\mathrm{7}}}^{Nc},z_{V_{\mathrm{15}}}^{Nc}\right),\mathrm{\dots },\left(z_{Q_{\mathrm{1}}}^N,z_{V_{\mathrm{3}}}^N,z_{V_{\mathrm{7}}}^N,z_{V_{\mathrm{15}}}^N\right),\end{array}$$

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}}}^{*},z_{V_{\mathrm{3}}}^{*},z_{V_{\mathrm{7}}}^{*},z_{V_{\mathrm{15}}}^{*}\right)$ is $\left(z_{Q_{\mathrm{1}}}^{\mathrm{1}},z_{V_{\mathrm{3}}}^{\mathrm{1}},z_{V_{\mathrm{7}}}^{\mathrm{1}},z_{V_{\mathrm{15}}}^{\mathrm{1}}\right)$. The lower boundary for the confidence interval was expressed as follows:

$$\begin{array}{}\text{(18)}& \left\{\begin{array}{l}{z}_{Q_{\mathrm{1}}}^L=\min \left(z_{Q_{\mathrm{1}}}^{\mathrm{1}},z_{Q_{\mathrm{1}}}^{\mathrm{2}},\mathrm{\dots },z_{Q_{\mathrm{1}}}^{Nc}\right)\\ z_{V_{\mathrm{3}}}^L=\min \left(z_{V_{\mathrm{3}}}^{\mathrm{1}},z_{V_{\mathrm{3}}}^{\mathrm{2}},\mathrm{\dots },z_{V_{\mathrm{3}}}^{Nc}\right)\\ z_{V_{\mathrm{7}}}^L=\min \left(z_{V_{\mathrm{7}}}^{\mathrm{1}},z_{V_{\mathrm{7}}}^{\mathrm{2}},\mathrm{\dots },z_{V_{\mathrm{7}}}^{Nc}\right)\\ z_{V_{\mathrm{15}}}^L=\min \left(z_{V_{\mathrm{15}}}^{\mathrm{1}},z_{V_{\mathrm{15}}}^{\mathrm{2}},\mathrm{\dots },z_{V_{\mathrm{15}}}^{Nc}\right)\end{array}\right..\end{array}$$

The upper boundary for the confidence interval was estimated by the following:

$$\begin{array}{}\text{(19)}& \left\{\begin{array}{l}{z}_{Q_{\mathrm{1}}}^U=\max \left(z_{Q_{\mathrm{1}}}^{\mathrm{1}},z_{Q_{\mathrm{1}}}^{\mathrm{2}},\mathrm{\dots },z_{Q_{\mathrm{1}}}^{Nc}\right)\\ z_{V_{\mathrm{3}}}^U=\max \left(z_{V_{\mathrm{3}}}^{\mathrm{1}},z_{V_{\mathrm{3}}}^{\mathrm{2}},\mathrm{\dots },z_{V_{\mathrm{3}}}^{Nc}\right)\\ z_{V_{\mathrm{7}}}^U=\max \left(z_{V_{\mathrm{7}}}^{\mathrm{1}},z_{V_{\mathrm{7}}}^{\mathrm{2}},\mathrm{\dots },z_{V_{\mathrm{7}}}^{Nc}\right)\\ z_{V_{\mathrm{15}}}^U=\max \left(z_{V_{\mathrm{15}}}^{\mathrm{1}},z_{V_{\mathrm{15}}}^{\mathrm{2}},\mathrm{\dots },z_{V_{\mathrm{15}}}^{Nc}\right)\end{array}\right..\end{array}$$

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, 3-day flood volume, 7-day flood volume and 15-day 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}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\right)$ can be derived by multiplying the benchmark flood hydrograph by different amplifiers, given as described below.

The amplifier K₁ for the annual maximum daily discharge was calculated by

$$\begin{array}{}\text{(20)}& K_{\mathrm{1}}={\displaystyle \frac{z_{Q_{\mathrm{1}}}}{Q_{\mathrm{1}}^B}}.\end{array}$$

The amplifier K₃₋₁ for the 3-day flood volume except for the annual maximum daily discharge was calculated by

$$\begin{array}{}\text{(21)}& K_{\mathrm{3}-\mathrm{1}}={\displaystyle \frac{z_{V_{\mathrm{3}}}-V\left(z_{Q_{\mathrm{1}}}\right)}{V_{\mathrm{3}}^B-V\left(Q_{\mathrm{1}}^B\right)}},\end{array}$$

where V(⋅) is the operator transforming daily discharge into flood volume. The amplifier K₇₋₃ for the 7-day flood volume except for the 3-day flood volume was calculated by

$$\begin{array}{}\text{(22)}& K_{\mathrm{7}-\mathrm{3}}={\displaystyle \frac{z_{V_{\mathrm{7}}}-z_{V_{\mathrm{3}}}}{V_{\mathrm{7}}^B-V_{\mathrm{3}}^B}}.\end{array}$$

Finally, the amplifier K₁₅₋₇ for the 15-day flood volume except for the 7-day flood volume was calculated by

$$\begin{array}{}\text{(23)}& K_{\mathrm{15}-\mathrm{7}}={\displaystyle \frac{z_{V_{\mathrm{15}}}-z_{V_{\mathrm{7}}}}{V_{\mathrm{15}}^B-V_{\mathrm{7}}^B}}.\end{array}$$

4 Results

4.1 Nonstationary analysis for marginal distributions

The time-varying moment model was employed to perform nonstationary analysis for each marginal distribution of the multivariate flood series $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},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₁, whereas the Gamma distribution was chosen as the theoretical distribution for the flood volume series V₃, V₇ and V₁₅. 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.

Table 2 Results of nonstationary analysis for the marginal distributions of $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$.

The relationships between μ and covariates were built by the exponential function in Eq. (3). α₁ and α₂ 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}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$, the location parameters μ of the three short-duration flood series, i.e. Q₁, V₃ and V₇, were positively linked to Pop, whereas the RI was the driving factor reducing μ for all flood series, including Q₁, V₃, V₇ and V₁₅. Owing to the difference in covariate selections, the short-duration flood series, including Q₁, V₃ and V₇, displayed asynchronous nonstationary behaviours with the long-duration flood series V₁₅ occurring in the observation period of 1951–2012. As shown in Fig. 4, Q₁, V₃ and V₇ presented significantly increasing trends during 1951–2005, particularly since the 1980s, marking the beginning of a period of rapid urbanization in China. V₁₅ 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₁, V₃, V₇ and V₁₅, exhibited a sharp decline.

Figure 4 Nonstationary marginal distributions during both the observation and design life periods.

Table 3 Results of nonstationary analysis for the dependence structure of $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$.

The relationships between copula parameters and covariates were built by the exponential function in Eq. (7). β₁ and β₂ 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.

Figure 5 Statistical correlations between flood peak and flood volumes. Three asterisks (***) indicate the statistical correlation at the 0.01 significance level.

The predicted marginal distributions for $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ during the design life period from 2013 to 2100 were estimated using the time-varying moment model by replacing the observed covariates for μ with those predicted. Figure 4 also shows that the mean values of Q₁, V₃ and V₇ 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₁, V₃ and V₇ would be expected to reach their peaks in the mid-21th century followed by a slight declining trend because of a shrinking urban population. Since V₁₅ was only related to RI, V₁₅ would show an abrupt decline in 2023 due to the regulation of the Datengxia reservoir. In general, the predicted nonstationary marginal distributions for Q₁ and V₃ during 2013–2100 were roughly approximate to the marginal distributions under the assumption of stationarity, whereas the predicted nonstationary marginal distributions for V₇ and V₁₅ exhibited smaller mean values than those of the stationary distributions.

4.2 Nonstationary dependence structure for $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$

After estimating the nonstationary marginal distributions for $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$, the multivariate dependence structure was constructed by the dynamic C-vine copula with Q₁ elected as the key variable. Figure 5 illustrates significant correlations between the flood peak Q₁ and the flood volumes (i.e. V₃, V₇ and V₁₅). Table 3 shows the estimation results of the dynamic C-vine copula. The PIT test for the nonstationary dependence structure of $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},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 θ₁₁₅ for pair (Q₁, V₁₅) 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₁ displayed asynchronous nonstationarity behaviours with V₁₅ (see Table 2 and Fig. 4). Therefore, the dependence nonstationarity of the pair (Q₁, V₁₅) could possibly be attributed to the asynchronous marginal nonstationarities.

According to the regression function, θ₁₁₅ 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}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ but also affected the dependence of $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$. Figure 6 shows the time variations of θ₁₁₅ during the observation period of 1951–2012 as well as during the design life period of 2013–2100. Due to reservoir regulation, θ₁₁₅ presented two obvious upward change points in both 2006 and 2023. Besides this, θ₁₁₅ also exhibited an obvious decreasing trend with urban population growth from 1951 to the mid-21th century, followed by a slight increasing trend due to a shrinking urban population. During the design life period, the predicted nonstationary θ₁₁₅ suggests a weaker dependence for $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ than the dependence under the stationary assumption, since it is usually smaller than the stationary estimation.

Figure 6 Nonstationary copula parameter for pair (Q₁, V₁₅) during both the observation and design life periods.

In addition, the change-point 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}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$. Readers are referred to Bücher et al. (2014) and Kojadinovic (2017) for specific steps to implement the change-point 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}},V_{\mathrm{3}},V_{\mathrm{7}},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).

Table 4 Results of change-point detection for the marginal distributions and dependence of $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$.

Figure 7 Design values of the annual maximum daily discharge for different average annual reliability (AAR) varying from 0.01 to 0.99 under three nonstationary conditions.

Figure 8 Design values of the 3-day flood volume for different average annual reliability (AAR) varying from 0.01 to 0.99 under three nonstationary conditions.

Figure 9 Design values of the 7-day flood volume for different average annual reliability (AAR) varying from 0.01 to 0.99 under three nonstationary conditions.

Figure 10 Design values of the 15-day flood volume for different average annual reliability (AAR) varying from 0.01 to 0.99 under three nonstationary conditions.

Figure 11 Design flood hydrographs associated with OR, AND and Kendall probabilities.

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}},V_{\mathrm{3}},V_{\mathrm{7}},V_{\mathrm{15}}\right)$ during the design life period from 2013 to 2100. The left columns in Figs. 7–10 show the most-likely 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 most-likely 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}},V_{\mathrm{3}},V_{\mathrm{7}},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₁ and V₃ 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₇ and V₁₅ 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.

5 Conclusion and remarks

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 stationary-based 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 C-vine copula allowing for both time-varying marginal distributions and time-varying 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 most-likely design event and confidence interval were identified as the outcome of the multivariate hydrologic design. Multivariate flood series $\left(Q_{\mathrm{1}},V_{\mathrm{3}},V_{\mathrm{7}},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}},V_{\mathrm{3}},V_{\mathrm{7}},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}},V_{\mathrm{3}},V_{\mathrm{7}},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 extreme-value 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 trade-off 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 well-defined cause–effect analysis.

Appendix A: Calculating multivariate exceedance probabilities

A1 OR exceedance probability (formulated by Eq. 9 in the paper)

Since the cumulative distribution function $F\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},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:

1. Calculate the marginal probabilities $\left(u_{\mathrm{1}},u_{\mathrm{3}},u_{\mathrm{7}},u_{\mathrm{15}}\right)$ of $\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\right)$.

2. Generate m samples $\left(u_{\mathrm{1},i},u_{\mathrm{3},i},u_{\mathrm{7},i},u_{\mathrm{15},i}\right)$ ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },m)$ from the C-vine copula expressed by Eq. (4).

3. Calculate $$\begin{gathered} F\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\left|{\mathit{\theta }}_t\right.\right)=\frac{\mathrm{1}}{m+\mathrm{1}}{\sum }_{i=\mathrm{1}}^m \\ \mathbf{1}\left(u_{\mathrm{1},i}\le u_{\mathrm{1}},u_{\mathrm{3},i}\le u_{\mathrm{3}},u_{\mathrm{7},i}\le u_{\mathrm{7}},u_{\mathrm{15},i}\le u_{\mathrm{15}}\right) \end{gathered}$$.

4. Calculate $p_t^{\mathrm{or}}=\mathrm{1}-F\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},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:

1. Calculate the marginal probabilities $\left(u_{\mathrm{1}},u_{\mathrm{3}},u_{\mathrm{7}},u_{\mathrm{15}}\right)$ of $\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\right)$.

2. Generate m samples $\left(u_{\mathrm{1},i},u_{\mathrm{3},i},u_{\mathrm{7},i},u_{\mathrm{15},i}\right)$ ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },m)$ from the C-vine copula expressed by Eq. (4).

3. Calculate $$\begin{gathered} p_t^{\mathrm{and}}=\frac{\mathrm{1}}{m+\mathrm{1}}{\sum }_{i=\mathrm{1}}^m \\ \mathbf{1}\left(u_{\mathrm{1},i}\ge u_{\mathrm{1}},u_{\mathrm{3},i}\ge u_{\mathrm{3}},u_{\mathrm{7},i}\ge u_{\mathrm{7}},u_{\mathrm{15},i}\ge u_{\mathrm{15}}\right) \end{gathered}$$.

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:

1. Calculate the marginal probabilities $\left(u_{\mathrm{1}},u_{\mathrm{3}},u_{\mathrm{7}},u_{\mathrm{15}}\right)$ of $\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\right)$.

2. Calculate ${\mathit{\rho }}_t=F\left(q_{\mathrm{1}},v_{\mathrm{3}},v_{\mathrm{7}},v_{\mathrm{15}}\left|{\mathit{\theta }}_t\right.\right)$ (see calculation steps 2–3 for OR exceedance probability).

3. Generate m samples $\left(u_{\mathrm{1},i},u_{\mathrm{3},i},u_{\mathrm{7},i},u_{\mathrm{15},i}\right)$ ($i=\mathrm{1},\mathrm{2},\mathrm{\dots },m)$ from the C-vine copula expressed by Eq. (4).

4. For $j=\mathrm{1},\mathrm{2},\mathrm{\dots },m$, calculate $$\begin{gathered} v_j=\frac{\mathrm{1}}{m+\mathrm{1}}{\sum }_{i=\mathrm{1}}^m \\ \mathbf{1}\left(u_{\mathrm{1},i}\le u_{\mathrm{1},j},u_{\mathrm{3},i}\le u_{\mathrm{3},j},u_{\mathrm{7},i}\le u_{\mathrm{7},j},u_{\mathrm{15},i}\le u_{\mathrm{15},j}\right) \end{gathered}$$.

5. 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)$.

6. Calculate $p_t^{\mathrm{ken}}=\mathrm{1}-K_C\left({\mathit{\rho }}_t\right)$.

Appendix B: Generating the multivariate design event samples (formulated by Eq. 17 in the paper)

To calculate the most-likely 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:

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

2. Generate a random integer (denoted by t_r) among (T₁, T₁+1, … , T₂).

3. Generate a random sample $\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\right)$ following the multivariate distribution $F\left(z_{Q_{\mathrm{1}}},z_{V_{\mathrm{3}}},z_{V_{\mathrm{7}}},z_{V_{\mathrm{15}}}\left|{\mathit{\theta }}_{t_r}\right.\right)$ with the distribution parameter vector ${\mathit{\theta }}_{t_r}$.

4. Calculate the annual exceedance probability for each year throughout the period from T₁ to T₂.

5. Calculate AAR during the period from T₁ to T₂.

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

7. If i<N, repeat steps (2)–(6).