Compounding joint impact of rainfall, storm surge and river discharge on coastal flood risk: an approach based on 3D fully nested Archimedean copulas

Compound flooding is a multidimensional consequence of the joint impact of multiple intercorrelated drivers, such as oceanographic, hydrologic, and meteorological. These individual drivers exhibit interdependence due to common forcing mechanisms. If they occur simultaneously or successively, the probability of their joint occurrence will be higher than expected if considered separately. The copula-based multivariate joint analysis can effectively measure hydrologic risk associated with compound events. Because of the involvement of multiple drivers, it is necessary to switch from bivariate (2D) to trivariate (3D) analyses. This study presents an original trivariate probabilistic framework by incorporating multivariate hierarchal models called asymmetric or fully nested Archimedean (or FNA) copula in the joint analysis of compound flood risk. The efficacy of the derived FNA copulas model, together with symmetric Archimedean and Elliptical class copulas, are tested by compounding the joint impact of rainfall, storm surge, and river discharge observations through a case study at the west coast of Canada. The obtained copula-based joint analysis is employed in multivariate analysis of flood risks in trivariate and bivariate primary joint and conditional joint return periods. The estimated joint return periods are further employed in estimating failure probability statistics for assessing the trivariate (and bivariate) hydrologic risk associated with compound events. The statistical tests found the fully nested Frank copula outperforms symmetric 3D copulas. Our work confirms that for practical compound flood risk analysis together with bivariate or univariate return periods, it is essential to account for the trivariate joint return periods to assess the expected compound flood risk and strength of influence of different variables if they occur simultaneously or successively. The bivariate (also univariate) events produce a lower failure probability than trivariate analysis for the OR-joint cases. Thus, ignoring the compounding impacts via trivariate joint analysis can significantly underestimate failure probability and joint return period.


Introduction
Flooding is becoming the most intensive, costly and frequently occurring natural hazard worldwide (Eric et al. 2009;Lian 2016;Xu et al. 2019;Burgan and Icaga 2019). According to long-term data on natural disasters, flooding has accounted for about 40% of all loss-related natural catastrophes since 1980 (Munich Re 2020). Rapid urban expansion and climate change might further increase the severity and frequency of extreme events in many coastal areas (i.e., Milley et al. 2002). Few recent flood events on the global scale, for example, Hurricane Igor in Canada (i.e., Pirani and Najafi 2020), Cyclone Nargis in Myanmar (i.e., Fritz et al. 2009), Hurricane Katrina in the US (i.e., Jonkman et al. 2009), Hurricane Harvey in the US (i.e.,), the winter flooding in the United Kingdom in (i.e., Haigh et al. 2016 etc., all indicate severe impacts of the coastal hazards. According to Intergovernmental Panel on Climate Change (IPCC), compound events (CE) can be defined as simultaneous (i.e., co-occurrence) or successively occurrences (i.e., in close succession, a few hours to a day apart) of two or more extreme or non-extreme events (i.e., Seneviratne et al. 2012). The interaction between oceanographic, hydrological, and meteorological physical drivers can result in compound flooding (CF) events and affect low-lying coastal communities Jones 2004, 2005;Jane et al., 2020). Flooding in the coastal areas primarily arises from four main physical factors or source mechanisms, such as (1) a combination of storm surges with high predicted astronomical tides or high normal tides (also called storm tides); (2) locally or remotely generated waves (swell); (3) river discharges at the lower estuary or fluvial flooding; and (4) direct surface runoff or pluvial flooding (i.e., Hendry et al. 2019). Most existing studies individually consider these four main drivers in assessing the coastal flood risk. However, flooding in coastal areas is often characterized by multiple physical factors because they are naturally correlated or interdependent. If these concurrent events display statistical dependencies through a common forcing mechanism, the probability of their joint occurrence will be higher than that expected, considering the extremes of each variable separately, with a consequent increase in the likelihood or chance of coastal flooding.
There is a lack of a consistent mathematical description of the compound extremes phenomenon despite multiple statistical approaches suggested to study it. Such demonstrations usually account for the number of joint extremes or mutual concurrency between different flood hazard types proxy variables. The combinations of univariate analysis in demonstrating the risk of compound events are only valid if no dependency is exhibited among the targeted compound variables. However, this assumption is not reasonable in practice as it may mislead the conclusion or underestimate the hydrologic risk. The earlier studies, such as Coles (1999), Svensson and Jones (2002), and Samuels and Burt (2002), pointed out the importance of extreme dependence measures, known as chi (or ) and chi bar (or ) statistics, which helps investigate asymptotic behavior for given bivariate random observations. Besides this, Coles (2001) pointed out the existence of distinct varieties of multivariate statistical models by extending univariate extreme value theory, such as component-wise block maxima, point process methods or threshold excess model. These models are usually developed by assuming that one of the components is either completely independent or asymptotically dependent on some joint tail region of the given random pairs (Coles and Tawn 1994;Ledford and Tawn 1997). Few incorporations, such as Heffernan and Tawn (2004), proposed a semiparametric approach called the conditional method to capture the mutual concurrence between the extreme observations. Boldi and Davison (2007) pointed out the importance of the semiparametric model based on mixtures of Dirichlet distributions to model dependence between the multivariate extremes. Cooley et al. (2010) incorporated a pair-wise beta distribution model describing the asymptotic dependence between extreme observations. Zheng et al. (2013) demonstrated the adequacy of the bivariate logistic thresholdexcess model for compounding the impact of storm surge and extreme rainfall observations. Similarly, Zheng et al. (2014) examined the dependence structure between storm surge and extreme rainfall characteristics through a comparative assessment between three classes of statistical models, such as threshold excess, point process, and conditional method, in terms of their adequacy to quantify the flood risk. The investigation also revealed that if there is a strong dependency between the random variables, the performance of point process methods dominates over the threshold excess and conditional models. If lower or weak dependencies are shown between the bivariate extreme observations, none of the models performed well or provided satisfactory outcomes. Statistically, compound flooding (CF) is an extreme multidimensional consequence of multiple intercorrelated hydrologic drivers. These driving forces co-occurred or in close succession, limiting the univariate frequency distribution approach's reliability or associated return periods. Incorporating a parametric-based multivariate probability distribution framework could permit the mutual concurrencies between the contributing variables of compound events (i.e., Serinaldi 2015;Aghakouchak et al. 2014;Bevacqua et al. 2017). Compared to the empirical estimation procedure, the statistical approach minimizes the uncertainties of statistical properties that we want to examine from the given observational datasets. The outcome will depend on how well we select these parametric models to fit the given observations. The multivariate (e.g., bivariate or trivariate) dependency modeling via the conventional parametric functions (e.g., bivariate GEV or Gamma) often exhibits numerous statistical challenges (De Michele and Salvadori 2003;Nelsen 2006;Karmakar and Simonovic 2008;Madadgar and Moradkhani 2013;Saghafian and Mehdikhani 2014).
Recently, the copula-based methodology has been recognized as a highly flexible multivariate tool, widely accepted in modeling extreme hydrometeorological events (i.e., Saklar 1959;Salvadori and De Michele 2004;Karmakar and Simonovic 2009;Reddy and Ganguli 2013a, b;Kong et al. 2015;Latif and Mustafa 2020 and references therein). Recently, more attention has been given to studying compound flood modeling via a copula-based multivariate framework. For instance, Waal and Gelder (2005) modeled the bivariate dependence of extreme wave heights and periods by incorporating the Burr-Pareto-Logistic (BPL) copula. Lian et al. (2013) tested the efficacy of the Archimedean and Elliptical copulas in modeling extreme rainfall and storm surge in a coastal city. Xu et al. (2014) modeled the joint probability analysis of extreme precipitation and storm tide and its changing characteristics using the bivariate Frank, Gumbel and Clayton copula for Fuzhou City, China. Wahl et al. (2015) modeled the joint impacts of storm surge and rainfall using the 2D Extreme value (or EV) and Archimedean copulas. Masina et al. (2015) incorporated EV copulas in the bivariate joint analysis of sea levels and waves observations for the Ravenna coast of Italy. Moftakhari et al. (2017) proposed a copula-based bivariate flood hazard framework by compounding the joint impacts of the coastal sea level observations (with future SLR) and riverine flooding. Paprotny et al. (2018) incorporated distinguished varieties of the Archimedean copulas in the bivariate joint and conditional joint analysis of rainfall, storm surges, river discharge flows, and wave observation for the European region. Xu et al. (2019) introduced Elliptical (i.e., Gaussian and Student's t copula) and Archimedean copulas (Gumbel, Frank and Clayton) in building joint distribution analysis of rainfall and storm tides observations for the Haikou City's coastal region of China. Zellou and Rahali (2019) performed bivariate frequency analysis and their joint risk probability estimation by compounding the joint impact of extreme rainfall and storm surges in the urban neighborhood close to the estuary of the Bouregreg River in Morocco. Similarly, Moftakhari et al. (2019) modeled the flood hazard estimations by a joint combination of bivariate statistical analysis and hydrodynamic modeling approach for the tidal channels and estuaries. Couasnon et al. (2020) globally modeled bivariate joint distribution analysis of storm surge and daily river discharge observations to identify hot spots of the compound flooding potential from the storm surges and river discharge in the river mouths. This experiment is investigated by analyzing target variables concerning their timing, joint dependence structure, and associated joint return periods. Recently, Ghanbari et al. (2021) proposed a bivariate flood hazard distribution framework using the 2D copulas to estimate compound coastal-riverine frequency under the current and future climate scenarios to capture the impact of climate change.
All the cited studies are limited to bivariate joint distribution by compounding the mutual concurrency between two flood-driving agents via pair-wise dependency simulation. The risk of flooding in the coastal regions can be compounded by multiple driving agents and exhibit complex interdependency between them. For example, a low-pressure system like a tropical cyclone, can co-occur with storm surge, heavy rainfall, and possible river overflows that can lead to compound events (i.e., Couasnon et al. 2020;Paprotny et al. 2018). The potential damage could likely be a function of multiple relevant hydrological events where ignorance of spatial dependency among them might contribute to the underestimation of uncertainty (Renard and Lang 2007;Graler et al. 2013). A more comprehensive understanding of compound events' uncertain behavior can be achieved by simultaneously integrating all the intercorrelated random variables. Work presented in this paper focuses on trivariate joint analysis. Limited literature is available on copula-based trivariate probability analysis. For instance, Serinaldi and Grimaldi (2006) (via mono-parametric and asymmetric or fully nested Archimedean (or FNA), Grimaldi and Serinaldi (2007) (FNA); Reddy and Ganguli 2013a, b (via FNA and trivariate Student's t copula for flood characteristics), Genest et al. 2007 (via trivariate elliptical copulas); Fan and Zheng (2016) (via entropy copula). To the best of our knowledge, availability of a trivariate distribution framework is rare in the field of compound flood modeling and corresponding hydrologic risk assessments.
In this regard, the adequacy of the symmetric Archimedean copulas would be inconsistent when expanding to a higher dimension framework (i.e., n ≥ 3). Due to its monoparametric behavior, the symmetric Archimedean structure approximates the dependencies between multiple flood vector pairs by employing single dependence parameters. However, it cannot preserve all pair-wise dependencies at the lower stages (Renard and Lang 2007;Genest et al. 2007;Kao and Govindaraju 2008). In other words, all dependencies are and must be averaged to the same values. Such statistical assumption would not be feasible and practical when considering multidimensional complex hydrologic phenomena that usually exhibit different degrees of the correlation coefficients between random pairs. In these circumstances, it could be much more effective and practical to approximate each random pair individually through multiple parametric joint asymmetric functions (i.e., Serinaldi and Grimaldi 2007;Savu and Trede 2010;Reddy and Ganguli 2013a, b).
The flexibility of the multivariate hierarchical model in the higher-dimensional extreme modeling has been recognized in recent decades. These models can easily capture different mutual concurrency between and within different groups of hydrologic characteristics (Hofert and Pham 2013). The Fully Nested Archimedean (or FNA) can provide higher flexibility in building a higher-dimensional framework than the traditional symmetric Archimedean framework. In reality, the FNA is based on the joint integration of two or multiple bivariates or any dimensional Archimedean copulas structure through another Archimedean structure (Whelan 2004;Savu and Trede 2010). Few previous incorporations, such as Grimaldi and Serinaldi (2006), incorporated FNA structure in the trivariate modeling of flood peak, flood volume and duration characteristics. Reddy and Ganguli (2013a, b) introduced the FNA structure and the Elliptical copula in the trivariate analysis of flood peak flow, flood volume and duration characteristics. Based on our knowledge, the FNA copula framework has never been used in modeling compound flooding events (i.e., the joint occurrence of rainfall, storm surge, and river discharge observations) in the coastal flood risk assessments in Canada or abroad.
Estimating the risk of compound extremes is essential for policy-making, disaster risk reduction, resilience building and engineering practice. Flooding is one of Canada's most intensive hydrologic phenomena and costly natural disasters. The impact of climate change significantly increases the risk of several climate extremes across Canada, such as higher temperatures, changing precipitation patterns and storminess, rising sea levels, and diminishing sea ice (Lemmen et al. 2016;Bush and Lemmen 2019). Both, East and West Canada's coasts experience higher coastal instability due to the higher risk of coastal water levels (James et al. 2014;Atkinson et al. 2016). Moreover, on west Canada's coast of British Columbia, B.C., sea level is projected to increase by about a half-meter and one meter by the end of 2050 and 2100, respectively (B.C. Ministry Environment 2013). Canada's west (Pacific) coast is at higher risk of coastal flooding, which can be attributed to compounding the impacts of precipitation, streamflow discharge and tidal water level extremes, according to a study performed by Pirani and Najafi (2020). The risk of extreme water levels can increase the likelihood of high storm surges. When combined with other flood variables, high river discharge and rainfall events, their joint combinations can be devastating. Therefore, targeting the above variables can significantly improve the understanding of flooding impacts on the west coast of Canada. Simultaneously considering flood drivers' triplet can provide a better flood hazard assessments. The dependency level of joint probability distribution between different flood drivers can vary from coast to coast or exhibit regional differences. It is also known that the impact of climate change introduces nonstationarity or temporal variability within the statistical behavior of individual driving agents, which can further affect their joint behavior, thus increasing the risk of extreme events. The focus of the presented study is on the implementation of a trivariate distribution framework without taking into account the time-varying consequences, or dynamic statistical behavior exhibited by flood characteristics.
Based on the above-raised points, it is important for flood management practice to capture the complex interplay between the different flood-driving factors, such as rainfall, storm surge and river discharge, that can significantly impact the west coast of Canada. Simultaneous integration of these driving agents using a trivariate copula-based probability approach can improve the flood risk estimations. The motivation for this study arises from the lack of higher dimension comprehensive risk assessments for modeling compound flood extremes in the coastal areas of Canada. This study also highlights the importance of considering the failure probability (FP) method to perform trivariate hydrologic risk assessments of compounding rainfall, river discharge and storm surge events. In a few previous studies, Xu et al. (2015) used the FP method in the bivariate hydrologic risk assessments of flood peak flow, volume and duration characteristics in the Wei River in China. Moftakhari et al. (2017) used FP statistics in the bivariate compound risk analysis of coastal water level and river discharge in the United States. Xu et al. (2019) employed FP statistics in the bivariate risk analysis between storm tide and rainfall events in the coastal region of Haikou, China. These studies reveal that FP statistics can be a practical tool in performing hydrologic risk evaluation. Therefore, the present study utilized the FP statistics in the trivariate hydrologic risk assessments to investigate the compound effect of storm surge, rainfall and river discharge on the coastal flood risk. The nested form of Archimedean copula can overcome the challenges of higher-dimensional modeling and preserve all the pairwise dependencies much more effectively than the traditional symmetric 3D copula framework (i.e., Frank, Clayton, Gumbel or Gaussian etc.). Compared with symmetric Archimedean copulas, the FNA structure can use more than one dependence parameter to model the link between target variables. Our research project first incorporated a higher-dimensional copula-based framework in the modeling and hydrologic risk assessments of compound events for Canada's west (and east) coast. Each finding appears in a separate manuscript (Latif and Simonovic 2022). The present manuscript is a part of this project.
The main objective of this study is addressed by: (1) the development of trivariate distribution models using the Fully Nested Archimedean (FNA) copula and symmetric (or monoparametric) Archimedean and Elliptical class (i.e., Gaussian) copulas, and evaluation of their efficacy in compounding the joint impact of rainfall, storm surge and river discharge through a case study; and (2) evaluation of the importance of the trivariate return periods over the bivariate and univariate return periods (3) estimating the failure probability using the trivariate exceedance probability to assess the variation of trivariate (also bivariate) flood risk during the entire project design lifetime. As a prerequisite, the trivariate joint analysis includes a brief presentation of the theoretical aspects of developing 2-D copulas, 3-D symmetric copulas, and 3-D asymmetric or fully nested Archimedean copula in a trivariate joint analysis of the paper. The procedure for the estimations of copula dependence parameters, a performance measure of copulas via the goodnessof-fit test statistics, and the mathematical approach for deriving the return periods of flood characteristics are discussed in the following section manuscript. In the third section, the copula-based methodology is applied to a case study to establish a trivariate (bivariate) joint distribution relationship of multiple flood characteristics and the joint and conditional return periods. Subsequently, the trivariate return periods (and bivariate RP) and failure probability statistics are estimated to explore the trivariate (and bivariate) hydrologic risk caused by the compound rainfall, storm surge and river discharge events. The last section of the paper provides the research summary and conclusions that can be used as a reference in the hydraulic facility design and coastal flood control and mitigation practices.

Mathematical overview of the bivariate (or 2-D) copula joint simulation
This study incorporates and tests the adequacy of asymmetric or FNA copula in the trivariate analysis of storm surge, rainfall, and river discharge to investigate the compound effect on flood risk in the coastal region. The performance of the fitted FNA copulas is also compared with some frequently used Archimedean and Elliptical class copulas. This methodology also highlights the importance of trivariate return periods and failure probability statistics in the trivariate hydrologic risk assessments. Although the methodological framework outlined in this paper can be applied to different regions, we only focus on the west coast of Canada. Figure 1 shows the flow diagram of the analysis. Using the trivariate distribution framework via the 3D nested copulas framework demands a series of computational steps. After the extraction and pre-preprocessing of trivariate flood observations (i.e., rainfall, storm surges, and river discharge), the framework requires a series of analytical steps: (1) modeling of univariate marginal behavior using 1D parametric families functions, (2) establishment of bivariate joint dependence structure and estimation of joint cumulative distribution functions using the best-fitted 2D parametric copula functions; (3) modeling of trivariate joint dependence structure using the 3D FNA (also, symmetric 3D copulas structure) and estimation of joint cumulative distributions; (4) estimation of trivariate (and bivariate) joint and conditional return periods; and (5) performing multivariate hydrologic risk assessments using the failure probability statistics.
According to Nelsen (2006), the copula function connects multivariate probability distributions to univariate marginal functions. If [X, Y] represents the independent and identically distributed (i.i.d) bivariate random observations with continuous marginal distributions, i.e., u = F X (x) and v = F Y (y) , then their associated dependence function is known as copula or C, can be described as where C is any bivariate copula structure under consideration; F(x)andF(y) are the univariate marginal cumulative distribution functions (or CDF) of variables 'X' and 'Y'; H(x, y) is bivariate joint distribution expressed using the bivariate copula C (u, v) . Similarly, if f(x) and f(y) define the probability density of variables 'X' and 'Y', then joint probability density can be defined as Here, c is the density function of bivariate copula; u = F(x)andv = F(y) are the marginal distributions. Readers are advised to follow Rivest (1993) andNelsen (2006) for the extended mathematical details of copula functions.
The adequacy of the distinct varieties of 2-dimension (2-D) copulas tested in this study includes Archimedean class (i.e., Frank, Clayton), mixed version of two-parameter Archimedean functions (i.e., BB1 (Mixture of Clayton and Gumbel) (i.e., Constantino et al. 2008), BB6 (Mixture of Joe and Gumbel copula) (Manner 2010), BB7 (Mixture of Joe and Clayton copula) (Li et al. 2016), and BB8 (Mixture of Joe and Frank copula) (Tang et al. 2015)), and Extreme value class (i.e., Galambos, Husler-Reiss, Gumbel-Hougaard). Among the above classes, the Archimedean copulas are frequently used in hydrology because they cover an extensive range of joint dependency measures and have many families (Nelsen 2006;Zhang and Singh 2006;Veronika and Halmova 2014;Requena et al. 2016). (Supplementary Table (S1) lists statistical descriptions of the parametric 2-D copulas used and tested in this study.
Mathematically, the copula function (C ∶ [0, 1] → [0, 1]) can approximate the bivariate Archimedean copulas, if it follows the representation and can be expressed as where (⋅) and ϕ −1 are the generator function and the inverse of the selected Archimedean copulas, respectively. The Archimedean copulas are not derived from the standard multivariate distributions. Different families of the Archimedean copulas exhibit the different extent of dependency measure capabilities or have different ranges of Kendall's tau measures (Nelsen 2006). The Frank copula covers both positive and negative dependency measures (i.e., Kendall's tau τ θ ∈ [−1, 1]) while the Gumbel-Hougaard and Clayton families cover only the positive association (i.e., Kendall's tau τ θ ≥ 0) . On the other side, AMH copula covers both positive and negative dependency measures but are restricted for Kendall's tau τ θ ∈ [−0.181, 0.333]. The Frank copula is only a member of Archimedean families that justified radial symmetry (i.e., symmetric to x + y = 1) (Nelsen 2006). The Frank copula exhibited symmetrical behavior to secondary diagonals, which means having no tail dependency (i.e., both upper right tail 'λ UP ' and lower left tail 'λ L ' is zero). On the other hand, the GH copula is more suitable for modeling the dependency with upper-tail dependence (i.e., extreme observations). The Clayton copula exhibited a strong capability to model with lower (3) C(x, y) = ϕ −1 (ϕ(x) + ϕ(y)), for x, y ∈ [0, 1] tail dependency (Poulin et al. 2007). Accountability of the Joe copula would be justifiable when a higher positive dependency is exhibited between random variables (McNeil et al. 2015), which often plays an important role, especially in highlighting the upper tail dependence (UTD) structure in extreme value analysis (Alina 2018). The efficacy of Survival Clayton, Survival Joe, and Survival Gumbel copula is also tested in the bivariate joint dependency modeling of flood characteristics. The survival version of copulas is usually obtained by making a 180° rotation. Similarly, a rotated version of the mixed type bi-parametric Archimedean copulas is also investigated, such as Survival BB1 (rotated BB1 by 180°), Survival BB6 (rotated BB6 by 180°), Survival BB7 (rotated BB7 by 180°), and Survival BB8 (rotated BB8 by 180°).

Trivariate joint distribution analysis using symmetric 3-D copula simulation
An attempt of bivariate design quantiles and their associated return periods might be insufficient to reveal justifiable and comprehensive flood probability analysis studies. Flooding in the coastal region can result from the simultaneous or successive occurrence of multiple physical extremes (or non-extremes) events, such as rainfall, storm surge, river discharge, etc. Suppose these concurrent events display statistical dependencies (i.e., through a common forcing mechanism). In that case, the probability of their joint occurrence is higher than expected, considering each variable's extremes separately, with a consequent increase in the likelihood of coastal flooding. To visualize the probabilistic behavior of the compound flooding phenomenon in a much more comprehensive manner through accounting for more than two variables simultaneously, it is required to use a higherdimensional joint distribution framework.
The n-dimensional Archimedean copula can be formulated by extending two-dimensional form into 'n' order series, which can be as expressed by The consistency of Eq. (4) will be preserved as long as the generator function ∅( ⋅ ) is fully monotonic; otherwise, it might be inconsistent for hydrological samples due to a statistical assumption of a homogenous dependency across the variables. Also, the generator function '∅( ⋅ )' tends to be the strict generator if the pseudo-inverse function ∅ −1 ( ⋅ ) becomes an ordinary inverse function when ∅(0) tends to infinity (Grimaldi and Serinaldi 2006;Nelsen 2006;Reddy and Ganguli 2013a, b). In contrast to elliptical copulas, i.e., Gaussian copula, Archimedean copulas can capture different tail dependencies (i.e., Gumbel (upper tail) or Clayton (lower tail)).
Let us consider if the X, Y and Z series be the three intercorrelated random observations. Then the joint probability distribution, F, can be expressed as Also, the multivariate joint return period can be derived from Eq. (6) (i.e., Salvadori et al. (2010), as where F(⋅) is joint CDF or JCDF; T is return period; μ is the average inter-arrival time of sequential hydrologic or flood event = 1. In the case of trivariate distribution series, the joint distribution of random variables can be expressed as where H X,Y,Z (x, y, z) is trivariate joint distribution of random variables; F(⋅) is marginal distribution; and C is trivariate (7) H X,Y,Z (x, y, z) = C F X (x), F Y (y), F Z (z) = C u 1 , u 2 , u 3 copula function. This experiment incorporates the 3D version of monoparametric Archimedean copulas, i.e., Frank, Clayton, Gumbel-Houggard, Joe copula, and Elliptical, i.e., Gaussian (or normal) copula for establishing trivariate joint dependence structure. Supplementary Table (S2) lists a statistical description and mathematical expression of the symmetric 3-D Archimedean and Elliptical copulas incorporated in this study.

Trivariate joint analysis using 3D fully nested Archimedean (FNA) copula
The symmetric form of Archimedean copulas would be incapable of preserving all pair-wise dependency at the lower stages, often developed under the assumption of homogenous dependency among the variables of interest. This assumption could be impractical because different variables could show the different levels of dependency. In such circumstances, it must be necessary to increase modeling flexibility by allowing heterogeneous dependence via a fully nested Archimedean (FNA) copula structure (i.e., Whelan 2004;Serinaldi and Grimaldi 2007;Savu and Trede 2010). This structure is modeled or joins two or more ordinary bivariate or higher-dimensional Archimedean copulas by another Archimedean copula also called parent Archimedean copula (Hofert 2008;Reddy and Ganguli 2013a, b). Mathematically, the dependence structure of FNA for random trivariate observations can be expressed as (Savu and Trede 2010) In Eq. 8, the first derivative of φ −1 2 •φ 1 is completely monotonic where φ 1 and φ 2 are Laplace transforms, and symbol ' • ′ represents the composite function. Also, the pair ( u, v) has the bivariate marginal structure in the form of Eq. 4 with a Laplace transform φ 1 . Similarly, random pairs ( u, w) and (v, w) have the bivariate margins in the form of Eq. (4) with Laplace transformφ 2 . Figure 2 shows the schematic diagram of the nesting procedure in the 3-D copula simulation via the FNA structure. In Fig. 2 is the inner copula with copula dependence parameterθ 1 , and C 2 (C 1 , F W (w);θ 2 ) is the outer copula with dependence parameterθ 2 . Also, F U (u), F V (v), and F W (w) are the marginal cumulative distributions derived from the best-fitted univariate probability distribution functions. The formulated 3-D copula C 1 (F U (u), F V (v), F W (w)) represents the joint simulation of two bivariate structures through trivariate asymmetric Archimedean functions. Its applicability could be justified only if the dependency strength among the two variables, i.e., (u, v), will dominate the correlation structure between these variables and the third variable, i.e., (u, w) and (v,w) (i.e., Savu and Trede 2010).
This study tests the adequacy of the fully nested form of Clayton, Gumbel-Hougaard and Frank copulas in the trivariate joint distribution modeling, which are further compared with the symmetric form of Archimedean and Elliptical class of the 3D copulas (refer to "Trivariate joint distribution analysis using symmetric 3-D copula simulation") in characterizing the most justifiable joint dependence structure of compound flooding events. Table 1 lists a statistical description of the asymmetric or FNA 3D copulas structure incorporated in this study.

Estimation of copula dependence parameters
Several statistical algorithms are available for the estimation of the vector of unknown statistical or copula dependence parameters, for instance, method of moment (MOM) rank-based nonparametric measures, i.e., Kendall's tau ( ) or Spearman's rho ( ) (i.e., Reddy and Ganguli 2012), exact maximum likelihood or EML (i.e., Zhang et al. 2016), inference functions for marginal or IFM (i.e., Joe 1997), canonical maximum likelihood or CML (Genest et al. 1995), maximum pseudo-likelihood estimations (MPL) (Genest et al. 1995). This study estimates copula dependence parameters for 2-D copulas and 3-D symmetric Archimedean and Elliptical copulas using the MPL estimator. The MPL estimators use the rank-based empirical distributions for estimating copula parameters independently from their univariate marginal distributions and can be applied for both one-and multi-parameter copula functions (i.e., De Michele et al. 2005;Klein et al. 2010;Kojadinovic and Yan 2010;Reddy and Ganguli 2012). MPL estimator is the modified version of the traditional maximum likelihood method. Mathematically, it can easily estimate the copula dependence parameter by maximizing the pseudo-log-likelihood function expressed as where, is the copula dependence parameter; l( ) is the pseudo-log-likelihood function; F 1 X i,1 andF 1 X i,2 are the empirical cumulative distribution functions (CDFs). The above Eq. 9, estimates copula parameters for a bivariate joint case and can be extended for trivariate copula distribution (or any n-dimensions case) as

Goodness-of-fit (GOF) test for fitted copulas
The adequacy of the best-fitted or hypothesized copula fitted to bivariate and trivariate flood characteristics is examined Schematic diagram of constructing asymmetric or Fully Nested Archimedean (FNA) copulas in the trivariate joint distribution modeling using the Cramer-von Mises test statistics, which is often considered the most powerful model compatibility testing procedure (Genest and Rémillard 2008;Genest et al. 2009). This test is based on a parametric bootstrapping procedure and uses the Cramer-von Mises statistic S n , which can be mathematically expressed as follows: For testing the compatibility level of 2-D copulas in bivariate joint analysis For testing the fitness consistency of 3D copulas in the trivariate joint analysis where c n u 1 , u 2 , u 3 and c n u 1 , u 2 are the trivariate and bivariate empirical copulas estimated using n observational flood attribute pairs; C is the parametric copula derived under the null hypothesis H 0 , and (u 1 , u 2 , u 3 ) are the univariate marginal distributions. U 1i,n , U 2i,n or U 1i,n , U 2i,n , U 3i,n are pseudo-observations of C transformed from . N u m e r i c a l ly, t h e v a l u e o f U 1i,n , U 2i,n and U 3i,n can be estimated by using the following mathematical equation On the other side, the p values for each tested copula are estimated using the parametric bootstrapping procedure (i.e., Genest and Rémillard (2008)) which can be mathematically expressed as where N is the number of simulations. This fitness test statistics involves testing the null hypothesis H 0 against the alternative hypothesis H a as given below.
Null hypothesis ( Table 1 Mathematical expressions of the asymmetric or fully nested Archimedean (or FNA) 3-D copulas and their associated statistical properties where, O is the open subset of ℜ q for some integer value q. The acceptance or rejection of the considered copulas is based on the S n test statistics and their estimated p values. The minimum value of the S n test must indicate the minimum difference between empirical and derived parametric copulas. At the same time, their estimated p value must be larger than a significance level ( = 0.05) , which indicates the acceptance of the null hypothesis; otherwise, they are liable to rejection.
The model compatibility of the best-fitted 3-D structure is also compared analytically using the error indices statistics, for instance Mean Square Error (MSE), Root Mean Square Error (RMSE), and information criterion statistics, for instance Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). Mathematically, the RMSE test statistics can be defined as (Zhang and Singh 2007;Karmakar and Simonovic 2008).
where e(i) and t(i) are the empirical and theoretical or computed observations of the i th samples. K is the number of fitted parameters of the copula model, and N is the number of observation samples. MSE and RMSE statistics usually defined in error statistics in the units of constituents of interest. In other words, when its value is closer to zero, it can indicate optimum model performance. Minimum, the value of RMSE and MSE indicate a better fit.
Similarly, the AIC and BIC of the fitted distributions can be estimated using the MSE test value (Akaike 1974;Karmakar and Simonovic 2008) and Smaller value of AIC and BIC test statistics indicates better performance of the fitted model. Mean Absolute Error (MAE) (Willmott and Kenji 2005) test statistics are employed further to check the adequacy of fitted 3-D symmetric and asymmetric copulas as given by where C ti is the theoretical copula (or predicted) values which are developed using the 3-D FNA structure and C ei is the estimated empirical copula value, followed by Deheuvels (1979). The minimum value of MAE statistics indicates that the fitted model for joint probability distribution based on historical compound flooding observations is adequate. A study by Willmott and Matsuura (2005) confirms that MAE can be considered much more effective than RMSE. However, conversely, Chai and Draxler's (2014) study pointed out that as long as the observation exhibited unskewed (or Gaussian) distribution behavior, the performance of RMSE dominates over MAE measures.

Return periods of the flood characteristics and their risk assessments
Hydrologic and hydraulic applications are often interested in the average inter-arrival period between two design events, usually defined by the return period (Shiau 2003;Shiau 2006;Salvadori 2004). The design quantiles representing a higher return period are often used in hydraulic structure design practice (Requena et al. 2016). Concurrence probability often defines the chance that any extreme event which is either characterized through univariate or multivariate description exceeds a certain threshold level (Yue and Rassumesen 2002;Shiau 2003;Salvadori 2004). Yue and Rasmussen (2002) and Salvadori and De Michele (2004) study thoroughly discussed the basic concept of the return period.
The univariate return period can be estimated using the univariate cumulative distribution function (or CDF) value derived from the best-fitted probability distributions as where μ is the average inter-arrival times between two consecutive events, and its value must be unity (i.e., μ =1) for annual maxima-based flood modeling (Yue and Rasmussen 2002).

Derivation of primary return periods
The concept of a univariate return period might be helpful only if the concentration of single hydrologic events matches the requirements of the design process. Compound flooding resulting from multiple extreme or nonextreme factors (e.g., rainfall, storm surge or river discharge) requires multivariate exceedance probability and associated return periods. Estimating multivariate design variables quantiles under the different notations of return, i.e., based on joint probability distribution and conditional joint distribution analysis, are often an essential concern in hydrologic risk assessments (Salvadori 2004;Graler et al. 2013). Each approach to estimating return periods (i.e., joint or conditional joint) has significance and importance solely on the nature of the problem undertaken. Therefore, the selection of return periods relies on the importance of undertaken structure and its consequences of failure; their appropriate choice depends on the impact of design variables quantiles. Also, the importance of different return periods cannot be interchanged and impossible to decideon the most consistent ways (Serinaldi 2015). Literatures, such as Shiau (2003), Salvadori (2004), Salvadori and De Michele (2004) and Serinaldi (2015), pointed out a comprehensive mathematical framework for deriving different notation for return periods under copulas-based methodology.
The joint distributions for annual flood characteristics can be described through two different situations. In the first, both the hydrologic characteristic (say, X ≥ x AND Y ≥ y for the bivariate case and X ≥ x AND Y ≥ yANDZ ≥ z for the trivariate case) simultaneously exceed a certain threshold during a flood event, and their associated return period, called AND joint period, can be mathematically estimated as where F(x)andF(y) are the univariate marginal cumulative distribution functions obtained from the best-fitted probability functions, and C(F(x), F(y)) are the copula-based bivariate joint cumulative distribution functions. Similarly, the AND joint return period for trivariate joint distribution case (via 3-D copulas simulation) can be estimated by (i.e., Zhang and Singh 2007) where C(F(x), F(y), F(z)) are the copula-based trivariate joint cumulative distribution functions derived from bestfitted 3-D copulas. In the second condition of the joint return period, either of the variable (say, ( X ≥ x OR Y ≥ y ) for the bivariate case and (X ≥ x OR Y ≥ yORZ ≥ z) for trivariate case exceed given threshold. Thus, their associated return period, called OR joint return period, can be expressed as And, for trivariate joint distribution case (via 3-D copulas simulation)

Derivation of the conditional joint return period
In the design of hydraulic facilities and their risk assessments, it is essential to investigate flood events in such a manner that one can highlight the priority of one design variable over another design variable obtained using the conditional joint probability relationships (i.e., Salvadori and De Michele 2004;Shiau et al. 2006;Singh 2006, 2007;Salvadori and De Michele 2010). In this study, two different cases of conditional joint distribution and their associated return periods are considered and examined for bivariate cases and For trivariate joint distribution case (i.e., via 3-D copulas simulation), the conditional joint return period of variable X, Y given various percentile values of a third variable (Z ≤ z) in OR-joint case is expressed as Similarly, the conditional joint return period of variable X given various percentile values of variables

Evaluating multivariate hydrologic risk via failure probability statistics
According to studies made by Salvadori et al. (2011), Huang and Chen (2015), Read and Vogel (2015), and Moftakhari et al. (2017), the traditional definition of the return periods, i.e., OR-joint or AND-joint, would be ineffective and incapable of demonstrating the risk of potential flood hazards during the entire project design lifetime. Literatures, such as Salvadori et al. (2016) and Serinaldi (2014), suggested the importance and necessity of the failure probability (FP) statistic, which can quantify the risk of hydrologic events more effectively and practically. As Serinaldi (2014) pointed out, the definition of FP statistic is usually defined over a given design life period. It facilitates a more consistent and well-devised approach to assessing the risk of extreme hydrologic events. On the other side, a study conducted by Read and Vogel (2015) clearly revealed that the definition and applicability of the average return period would be problematic, which did not account for the planning horizon. Xu et al. (2015) incorporated FP statistics in the bivariate hydrologic risk evaluation for the Wei River in China. Similarly, Moftakhari et al. (2017) also utilize FPs statistics in the coastal flood risk assessments caused by fluvial flow and coastal water level. The present study incorporated FP statistics' notation to assess compound flooding risk for a trivariate joint distribution case. Let consider, (X 1 , X 2 , … … ..X T ) , (Y 1 , Y 2 , … … ..Y T ) and (Z 1 , Z 2 , … … ..Z T ) as three different series of hydrologic events where T is the arbitrary project lifetime. Again, (S 1 , S 2 , … … ..S T ) represent a sequence of bivariate hazard scenarios. According to the definition, FP statistics can be mathematically expressed as For an independent and identically distributed hydrologic occurrence, the FP can be estimated by From Eq. 25, it is clear that the FPs statistics can be estimated for both "OR" and "AND" joint cases as and (28) Similarly, the failure probability (FP) statistics for bivariate hazard scenarios (or bivariate joint cases) can be estimated as and

Case study: results and discussion
In the present study, the above-discussed methodology is used for assessing the compound effects of rainfall, storm surge and river discharge on flood risk on the west coast of Canada. The study aims to integrate or compound the joint impact of triplet flood characteristics simultaneously by comparing the efficacy of the fully nested Archimedean (FNA) copula vs symmetric or monoparametric Archimedean and Elliptical copulas in the trivariate joint distribution analysis and compound flood risk assessment.

Extraction of trivariate flood characteristics
The observed coastal water level data are collected by tidal gage station located at New Westminster (station id 7654), BC, Canada, and provided by Fisheries and Oceans Canada (https:// tides. gc. ca/ eng/ data). The coastal water level data are collected for the period 1970 to 2018. The geographic location of the tidal gage station is 49.2 • N (Latitide Decimal Degrees) and 122.91 • W (Longitude Decimal Degrees). The data are collected at local datum CD, and geolocation is Fraser River. The Fraser River is the longest river within British Columbia, flowing for 1375 km and finally draining into the Strait of Georgia (located in the South of Metro Vancouver, BC). The river's annual discharge at its mouth is about 112 cubic Kilometers (or 3550 m 3 /s). Parts of Vancouver within low-lying areas close to the ocean and Fraser River are vulnerable to coastal and river flooding. In addition, climate change will increase the region's vulnerability to flooding due to rising sea levels. The predicted astronomical high tide (or normal high tide) observations are extracted and obtained from the Canadian Hydrographic Service (CHS). Extraction of storm surge observations requires a proper time matching between normal high tide and total water level since the tidal pattern on the west coast of Canada is semidiurnal, i.e., two high peaks each day. The storm surge observations are estimated by subtracting high tide values from the coastal water level observations. The storm surge value is either positive or negative depending upon whether the predicted high tide is above or below the observed coastal water level (or total sea level).
After selecting the tidal gage station, the selection of streamflow discharge and rainfall stations is solely based on proximity within a radius of 50 km from the selected tidal gage station. In other words, the nearest rainfall gage and streamflow gage station are assigned within a specified radial distance for the same time frame. The selected stations (streamflow and precipitation) must not have more than 15% missing data. In addition to the physical distance of streamflow gages, flow routes are also tracked to ensure they are directed toward the tidal gage (Ward et al. 2018). The daily streamflow discharge records of the Fraser River are collected at the Fraser River at Hope gage station and obtained from the Environment and Climate Change Canada (https:// water office. ec. gc. ca/ search/ histo rical_e. html). The geographical location of this gage station is 49° 23ʹ 09ʹʹ N latitude and 121° 27ʹ 15ʹʹ W longitude. The rainfall data are collected at the Haney UBC RF Admin station, which is located at the geographical location 49° 15ʹ 52.100ʹʹ N latitude and 122° 34ʹ 23.400ʹʹ N longitude.
In this study, we search the dependency between storm surges, rainfall and river discharge observations to capture complex interplay for three separate cases. In Case 1 (i.e., bivariate joint dependency), we combine the annual maximum 24-h rainfall and the highest, or maximum, storm surge within a time period of ±4 days of the event. Similarly, in Case 2 (i.e., bivariate joint dependency), we combine the annual maximum 24-h rainfall series and the highest, or maximum, river discharge observation within a time period of ±4 days of the event. In case 3 (i.e., trivariate joint dependency modeling), annual maximum 24-h rainfall observation is combined with the highest value of storm surge and river discharge within a time period of ±4 days. At this stage, one can think, why the dependency of storm surge and rainfall, is observed at a time lag of ± 4 days. In the existing literature, for instance, Xu et al. (2014) and Xu et al. (2019) usually defined the compound events by selecting the second variable (e.g., highest storm surge) during the date of occurrence of the first variable (i.e., annual maximum 24-h rainfall series). In another study, Wahl et al. (2015) examined the dependency of the highest annual storm surge within a time range of ± 1 day of the highest precipitation. It can be possible that mutual concurrency between the abovediscussed flood variables is higher for time lag ± 4 days compared to time lag ± 0 day or ± 1 day or ± 2 day or ± 3 etc. or vice-versa. To do this, once we defined the highest 24-h rainfall other flood variables, maximum river discharge and maximum storm surge, are defined for the different time lag (i.e., ± 0 days, ± 1 days, ± 2 days, ± 3 days, ± 4 days and ± 5 days) from the date of annual maximum 24-h rainfall series. After that, we performed the correlation measures using Pearson, Kendall and Spearman coefficient. It is found that the maximum dependency exhibited between the annual maximum 24 h rainfall and highest rainfall (maximum storm surge) within a time lag of ± 4 days. Supplementary Table (S3) lists the descriptive statistics of trivariate flood characteristics annual maximum 24-h rainfall, maximum storm surge (time period = ±4 days), and maximum 24-h river discharge (time period = ±4 days). Different univariate statistical plots of trivariate flood characteristics are represented, such as scattergram and strip plots ( Supplementary  Fig. (S1(a-f))), histogram plots ( Supplementary Fig. (S2  (a-c))) (both in frequency and density), normal Q-Q plots ( Supplementary Fig. (S3)), and time series plots (Supplementary Fig. (S4)).

Test for stationarity (or time-trend behavior) within time series of trivariate flood characteristics
In the univariate or multivariate joint modeling or frequency analysis, individual flood characteristics must be stationary and have no serial correlation or autocorrelation behavior before introducing them into the probability distribution framework. These statistical assumptions are mandatory under stationary-based hydrological modeling and risk assessments. For this, Ljung and Box (1978)based hypothesis testing, also called Q-statistic, is implemented to investigate whether the individual time series are time-independent and have no autocorrelation (Cong and Brady 2011). Statistically, the Q-statistic usually follows a chi-square distribution with an 'h' degree of freedom under the null hypothesis H 0 (Ljung and Box 1978;Daneshkhan et al. 2016). The Q-statistics for the sample size 'n' with 't' are the total number of lags being tested with sample autocorrelation at the specific lag, i.e., ̂ t is estimated as where, Null hypothesis (H 0 ) is zero autocorrelation or independent distribution, and the alternative hypothesis (H a ) proves the existence of serial correlation (or autocorrelation). The estimated Q-statistics and their associated p value for different lag sizes (i.e., 20, 10 & 5) are shown in (Supplementary Table (S4)). Results reveal that data exhibited zero first-order autocorrelations as their estimated statistics are below their critical value for each of the flood characteristics by accepting the null hypothesis (H 0 ) at a 5% or 0.05 significance level (95% confidence interval) against their alternative hypothesis (H a ) (see Table 4). Supplementary Figure (S5) clearly illustrates the autocorrelation function (ACF) and partial ACF plot (at 95% or 0.95 confidence interval) for the individual flood characteristics. The nonparametric rank-based Mann-Kendall or M-K test is also performed (34) Q = n(n + 2) h ∑ t=1ρ to visualize the monotonic time-trend behavior within individual random characteristics under the null hypothesis H o against their alternative hypothesis H a (Mann 1945;Kendall 1975). Their estimated values are listed in (Supplementary  Table (S5)). Investigation reveals the acceptance of the null hypothesis (H 0 ), which points to zero monotonic trends at the 5% or 0.05 level of significance (or 95% confidence interval) for rainfall and river discharge. On the other side, rejection of the null hypothesis, H 0 for the variable storm surge, indicates the existence of monotonic time-trend behavior, as its calculated z value exceeds the critical value of Z-statistics (e.g., Z = ± 1.96 at 5% significance level or 95% confidence interval), refer to Table S5. In conclusion, detrending of the maximum storm surge time series observations is required to acquire a state of stationarity before proceeding with the univariate probability analyses. For the other two flood variables, no detrending steps are required. (Supplementary Fig.  (S6)) differentiated between the original sample (e.g., before detrend) of storm surge and their detrended samples after incorporating the pre-whitening procedure (e.g., Estimated Z-statistics < Critical Z-statistics at 5% significance level). Finally, the flood characteristics, e.g., rainfall, detrended version of the storm surge, and river discharge observations, are introduced into the univariate distribution framework.

Modeling of univariate marginal probability distributions of flood characteristics
Several distributions often would fit the data equally well, but each would give different estimates of a given quantile, especially in the distribution's tails. Therefore, the goodness-of-fit (GOF) procedure is used to visualize the compatibility of the fitted distributions (Karmakar and Simonovic 2008). Supplementary Table (S6) lists the probability distribution functions (or PDFs) of the candidate univariate models tested for individual flood characteristics. The Generalized Extreme Value (or GEV) distribution often exhibits significant benefits in hydrologic extreme value analysis. The GEV includes three distinct functions: Gumbel, the (Reversed) Weibull and Fréchet distributions (i.e., Coles 2001;Khaliq et al. 2006). Each function exhibits different tail behavior based on the shape parameter 'k', i.e., Gumbel is characterized by light tail behavior; Weibull distribution showed heavy tail (Graler et al. 2013). If the shape parameter 'k' is equal to 0 and corresponds to a thin upper and unbounded tail for GEV distribution, it is represented by the Gumbel function. On the other side, the normal or Gaussian distribution (e.g., Goel et al. 1998), 2-parameter (or bi-parametric) Lognormal distribution (Yue 2000), 2-parametric Logistic distribution (Bobee and Ashkar 1989), 2-parameter Gamma distribution (Daneshkhan et al. 2016) etc., are also widely and frequently used in modeling extreme hydrologic observations. The vector of unknown statistical parameters of fitted distributions is estimated via the maximum likelihood estimation (or MLE) procedure (i.e., Owen 2008). The MLE estimator provides the minor (or smallest) sampling variance of the estimated distribution parameters and the estimated quantiles (Tosunoglu and Kisi 2016) and facilitates less complexity. The estimated parameters of fitted distributions are listed in Supplementary Table (S7(a-c)).
Selecting a justifiable probability distribution in defining marginal behavior often demands higher accuracy and precision by adopting a series of quantitative and visual inspections (e.g., graphical interpretation). The empirical nonexceedance probabilities, also called observed (or empirical) cumulative frequency, are calculated using the Gringortenbased position-plotting formula (i.eGringorten 1963;Karmakar and Simonovic 2008). The estimated empirical cumulative values are compared with theoretical or modeled values (or CDF of fitted distributions) to investigate gaps or deviations between the empirical and fitted observations. Different analytical goodness-of-fit (or GOF) tests statistics are used to check the adequacy of fitted models, for instance, the Kolmogorov-Smirnov (or K-S) test (i.e., Xu et al. 2015), and Anderson-Darling (or A-D) test (i.e., Anderson and Darling 1954;Farrel and Stewart 2006), that is based on distance criteria; statistical tests that are based on information criteria statistics, e.g., Akaike Information criteria (or AIC) (i.e., Akaike 1974), Schwartz's Bayesian Information criteria (or BIC) (i.e., Schwarz 1978). The Cramer-von Mises (or CvM) criterion (i.e., Cramér 1928;von Mises 1928;Watson 1961) is also used to test the fitness level of the candidate distributions. Compared with K-S statistics, A-D statistics put extra weight on the tail of probability distribution and could be more appropriate for modeling extremes (Farrel and Stewart 2006;Alam et al. 2018). On the other side, the AIC statistics included the lack of the model's fit on the one hand and the unreliability of the model due to the number of model parameters on the other hand (Zhang and Singh 2007;Karmakar and Simonovic 2008). Lower the AIC and BIC value indicates better fit distribution. Supplementary Table (S8(a-c)) presents the quantitative evaluation or the estimated GOF test statistics of each flood characteristic's fitted univariate candidate functions. Results reveal that the performance of the GEV distribution is the best for annual maximum 24-h rainfall observations (e,g., Minimum value of K-S, A-D and CvM test value compared with other probability functions). The Normal distribution is the best fit for the approximation of storm surge series (e.g., Minimum value of A-D, CvM, AIC and BIC test statistics). The GEV distribution is the best fit for the river discharge observation (e.g., Minimum value of K-S, A-D, CvM, AIC and BIC tests). Qualitative-based model consistency checks are also conducted using the probability density function (or PDF) plots, cumulative distribution (or CDF) plots, probability-probability (or P-P) plots, and quantile-quantile (Q-Q) plots to crosscheck or validate the results obtained by the analytical-based GOF measures (refer to Supplementary  Figs. 7(a-d), 8(a-d), 9(a-d))). In summary, all the analytical and graphical testing measures point to the application of GEV distribution for rainfall series, Normal distribution for storm surge series, and GEV distribution for modeling river discharge observations.

Dependency measures of trivariate flood characteristics
Parametric-based Pearson's linear correlation (r), and two nonparametric rank-based dependence measures, such as Kendall's tau ( ) and Spearman's rho (ρ) correlation coefficient, are used to measure the strength of dependency. The Pearson correlation coefficient is not invariant to monotonic transformations of Kendall's and Spearman's correlation measures and can be strongly affected by outliers, only capturing linear dependencies (Favre 2004). Also, the Pearson correlation measure would be incompatible with modeling heavy-tailed distribution series. Kendall's tau and Spearman's rho are estimated using the ranking of random series and exhibited high resistance to outliers (Klein et al. 2011). The estimated correlation statistics are shown in Table 2, where the statistical significance of the correlation measures is checked at the significance level of 5% (or 95% confidence interval). All three random pairs (e.g., rainfall-storm surge; rainfall-river discharge, storm surge-river discharge) exhibited a positive correlation at a 5% level of significance (or 95% confidence interval), where the higher correlations are exhibited for the first two random pairs than the third pair.
Besides the analytical-based dependency measures, several graphics-based visual analyses are also carried out using different 2-D statistical plots between random pairs such as chi-plot (Fisher and Switzer 2001), Kendall's or K-plot (Genest and Bones 2003), and scatterplot, referred to Supplementary Figure (S10(a-c)). Chi-plot is a scatter plot of the pairs ( i i ) . It uses the data ranks and i value as a measure of the distance of bivariate random observations (say p i v i ) from the center of the data sets within the range of [−1, 1]. Their values will tend to be positive (in the case of positively correlated random pairs) or negative (in the case of negatively correlated random pairs). Similarly, the control limits i are another measure placed at = ±c p ∕ √ n (Fisher and Switzer 2001). Therefore, the random pairs must be outside the control limit of the chi-plot when the stronger dependency is exhibited; otherwise, when all the random pairs are inside the control limit region indicates independence between random pairs. Supplementary Figure (S10(a)) clearly shows that a large portion of random samples for the pairs rainfall-storm surge, storm surge-river discharge, and rainfall-river discharge are outside the control limit region. On the other side, Kendall's plot is analogous to the quantile-quantile (Q-Q) plot. The deviation of random pairs from the main diagonal of the K-plot indicates a dependency between random observations, otherwise revealing independence when the K-plot tends to be linear Boies 2003, Genest and. Refer to S10(b); it is clearly shown that the K-plot deviates from the main diagonal for all random pairs of the rainfall-storm surge, storm surge-river discharges, and rainfall-river discharge. Similarly, a scatterplot (refer to S10(c)) also indicates the existence of a significant level (examined at 95% confidence interval) dependency between all flood pairs the increased density of points are located away from the diagonal region, i.e., close to 45 • angle).

Bivariate joint analysis via 2-D parametric copulas
The first stage of the multivariate (i.e., trivariate) modeling includes the bivariate dependency simulation, incorporating 2-D copula functions to evaluate bivariate joint cumulative distribution functions or JCDFs. The discussion in "Theoretical research framework" of the paper provides the list of varieties of copula families (Supplementary Table (S1)) used in this work. The copula Table 2 Correlations to mutual dependencies measures among flood characteristics These correlation statistics are measured using the original extracted flood characteristics before investigating time-trend behavior and pre-whitening step, which can be required within time series of storm surge observations. Referring to "Test for stationarity (or time-trend behavior) within time series of trivariate flood characteristics" and supplementary  -c)). The Cramer-von Mises distance statistics with the parametric bootstrap procedure are estimated using Eqs. 11 and 14 and employed to analytically validate and identify the most justifiable 2-D copula in the bivariate joint distribution analysis. For this purpose, the test statistics 'S n ' and its associated p value have been computed for two different bootstrap samplings (i.e., N (number of bootstrap samples) = 1000 and 500) by the mean of the parametric bootstrap procedure, refer to S9(a-c). The smaller value of S n test statistics indicates a minor deviation or gap between an empirical and derived (or fitted) parametric copula. Investigation reveals that BB1 copula exhibited a minimum value of S n test statistics (with an estimated p value greater than 0.05 for both numbers of bootstrapping samples, N = 500 and N = 1000), see Table S9(a). Therefore, it is selected as the most justifiable copula for defining the bivariate joint dependence structure of the rainfall and storm surge. Referring to S9(b), it is found that the Husler-Reiss copula exhibited the minimum S n test statistics (with an estimated p value greater than 0.05 for both the number of bootstrap samples, N = 500 and N = 1000) is the most parsimonious copula for defining the joint structure of storm surge and river discharge. The Survival BB7 (rotated BB7 180°) copula is found the most appropriate copula for approximating the joint structure between the rainfall and river discharge; refer to S9(c). Tail dependence (TD) is of essential concern in the multivariate hydrologic frequency analysis; if not well preserved by the selected copula, it may be attributed to higher uncertainty in the estimation of extreme quantiles or wrong decisions in the design of the water-related infrastructure (Poulin et al. 2007;Reddy and Ganguli 2013a, b). The upper tail dependence coefficient (UTDC) up is often more important than the lower tail dependence coefficient (LTDC) low . To assess the utility of the selected 2-D copulas in our study (i.e., BB1, Husler-Reiss, and Survival BB7), parametric values of the UTDCs up are estimated and listed in (Supplementary Table ( S10(a-c)). The nonparametric estimates of the upper tail dependence coefficient (or empirical estimates) are evaluated by estimators suggested by Caperaa et al. (1997) (i.e., Poulin et al. 2007), CFG up , and their results are compared with the parametric coefficient of UTDC up . Results indicate that all the three selected 2-D copulas provide a reasonable estimation of the upper tail limit compared with the empirical ones, although there is some minor difference. In other words, all the selected models can effectively capture the extreme portion (or upper tail dependence behavior) of each random pair where the gap between the up (parametric) and CFG up (empirical or nonparametric) is minimum.
The graphics-based visual assessments are also carried out using the scatter plot (refer to Supplementary Figure  (S11(a-c)) to check the consistency of fitted 2-D copula, where 1000 random samples are generated from the bestfitted copula selected for each random pair. Similarly, Supplementary Fig. (S12(a-c)) illustrates the Chi-plot and K-plot of the set of 1000 random samples drawn from the best-fitted bivariate copulas derived from MPL estimators for selected pairs of flood variables. From S11(a-c), it is inferred that the BB1, Husler-Reiss, and Survival BB7 copula performed satisfactorily since the generated samples (light gray color) are adequately overlapped with the natural dependence of historical observations (red color in S11(a-c)). The Kendall's tau values are estimated from the derived 2-D copulas fitted to flood variable pairs (refer to S11(a-c)), which indicates the existence of minimum gaps between empirical and theoretical Kendall's values.
In conclusion, the performance of BB1 copula (for rainfall-storm surge pair), Husler-Reiss copula (for storm surge-river discharge pair), and Survival BB7 copula (for rainfall-river discharge pair) are satisfactory and in support of the analytical fitness test result. These bivariate functions are employed in estimating bivariate JCDFs and deriving the bivariate and trivariate joint and conditional return periods (RPs) and the failure probability (FP) statistics for assessing multivariate hydrologic risk.

Trivariate joint distribution analysis
A more comprehensive flood risk analysis can be achieved by simultaneously accounting for all the targeted random vectors instead of considering the pair-wise mutual dependency via the 2D copula joint structure. The present work's main contribution is the use of 3D FNA copula in the trivariate joint distribution modeling of the three selected flood characteristics, i.e., rainfall, storm surge and river discharge observations, to investigate the compound flooding scenario in the coastal region. Firstly, different families of 3D symmetric Archimedean and Elliptical class copulas are tested in the trivariate analysis (refer to Supplementary distribution-free Kolmogorov-Smirnov (K-S) test is also performed to validate the results obtained based on the CvM functional statistics in selecting the best-fitted 3-D copula model. The K-S test statistics are calculated for candidate copulas and used to compare empirical copula against a parametric estimate of the copula derived under the null hypothesis, where p values are estimated using the parametric bootstrap procedure. Based on the estimated statistics (shown in Table S11), it is concluded that the Gaussian (or Normal) copula performs well and is selected because it exhibited (i) the highest value of the estimated Log-likelihood (L-L) test statistic; (ii) the minimum value of the CvM functional test S n statistic with p value is higher than 0.05; and (iii) the minimum value of the K-S test statistic. Supplementary Figures (S13(a)) illustrate the scatterplot matrix using generated random samples of size N = 1000 from the 3-D symmetric Gaussian copula (fitted via MPL estimators). Also, Supplementary Fig. (S13(b)) illustrate overlapped scatterplots between observed (red color) and theoretical (light gray) flood samples extracted from the 3-D symmetric Gaussian copula.
On the other side, the Fully nested Archimedean (or FNA) copula approximates each random flood variables pair individually by incorporating multiple parametric joint asymmetric functions. It can approximate multiple copula dependence parameters between different groups of variables to capture complex inter-dependency among compound events. In this investigation, fully nested forms of Frank, Clayton and Gumbel-Houggard copulas are tested (refer to Table 1). The FNA structure is incorporated using the HAC library package of R-software (HAC 2020). In the estimation of FNA structure, all the univariate marginal CDFs are developed using the best-fitted probability distribution, i.e., (GEV, Normal, GEV for three selected variables). The dependence parameters of copulas, both inner and outer copula in the FNA structure, are calculated using the maximum likelihood estimation (or MLE) procedure, and the results are presented in Table 3. At first, the performance of the fitted nested copula families is investigated using the Loglikelihood (LL) test statistic, as shown in Table 3. A higher value of LL statistics reveals a better fit. Figure 3a shows a scatterplot using the random samples (size, N = 1000) generated from the best-fitted 3D fully nested Frank copula fitted to trivariate flood characteristics. Similarly, Fig. 3b illustrates the overlapped scatterplot between historical observations (indicated by red color) and simulated observations (indicated by light blue color) using the asymmetric Frank copula. Results show that a fully nested Frank structure recognizes and captures the mutual dependence behavior of targeted flood characteristics much more effectively than the symmetric (or traditional) 3D Gaussian copula, as illustrated in Supplementary Fig. (S13b).
The reliability and suitability of the selected best-fitted fully nested Frank copula are compared with the selected symmetric 3-D Gaussian copula by comparing the theoretical Kendall's correlation coefficient estimated from the generated random samples (sample size, N = 1000) using the selected models together with the empirical Kendall's obtained from historical flood observations (Table 4). It is concluded that the asymmetric Frank copula exhibits a minimum gap between the theoretical and empirical Kendall's value compared to the symmetric Gaussian copula. The performance of this FNA 3-D model is much more practical in these circumstances; this model regenerates the mutual dependencies of historical flood characteristics more effectively compared to the symmetric 3-D Gaussian structure. On the other side, the performance of this symmetric 3-D model is not very effective and practical in this case. It did not effectively recognize or capture all possible mutual concurrency because of the single dependence parameter associated with the fitted copula. It approximates the dependencies between multiple variable pairs by employing a single dependence parameter and is incapable of preserving all pair-wise dependencies at the lower stage.
The model performance between symmetric Gaussian and asymmetric Frank copula is compared further using the error indices statistics, such as MSE, RMSE, and information criterion, for instance, AIC and BIC criterion statistics, refer to Table 5. It is found that the selected asymmetric Frank copula performed well compared to the selected symmetric Gaussian structure. Mean absolute error (MAE) statistic between the fitted observations via asymmetric and symmetric copula and empirical copula values data is examined; refer to Eq. (19). The minimum value MAE for the asymmetric Frank copula confirms a better fit (refer to Table 5). In conclusion, the selected asymmetric Frank copula provides a much more efficient approach in trivariate joint modeling. It is thus further employed in estimating the trivariate joint and conditional return periods (RPs). The estimated joint cumulative distribution function (JCDF) estimated from the fitted asymmetric Frank copula is introduced further to estimate failure probability to assess the variation of the trivariate compound flood hazard scenario.

Estimation of primary joint return periods (JRPs)
The objective of frequency analysis of extreme hydrologic events is to establish an inter-association between the magnitude of extreme events and their frequency of occurrence through best-fitted univariate or multivariate probability distributions. The copula models can form the basis for estimating various quantities, e.g., exceedance probabilities, which can be effective and practical for risk evaluation of hydrologic extremes (i.e., estimation of joint and conditional joint return periods). The univariate return periods of the flood characteristics are derived using the CDFs of the best-fitted probability distribution functions (i.e., GEV for Annual Maximum 24-h Rainfall, Normal for Maximum Storm Surge (Time interval = ± 4 days), and GEV for Maximum River Discharge (Time interval = ± 4 days),) using Eq. 19, and represented in (Supplementary Figure (S14)).
The accountability of univariate RPs would be problematic, attributed to the overestimation or underestimation of hydrologic risk. The estimated univariate design quantiles and their associate RPs would be helpful when only one flood characteristic is significant. In practice, the considerations of multiple design variables are required for revealing effective risk practices. This is especially true from the prospect of the compound flood risk assessments where the multiple extreme or non-extreme events are occurring either in close succession (i.e., under small time lag) or simultaneously (i.e., co-incidence probability). The accountability of multivariate exceedance probabilities and their RPs is a justifiable and practical approach to better understanding flood risk using the concept of joint and conditional joint probability distribution relationships. Supplementary  Figures (S(15-17)) indicate the joint probability density function (JPDF) and joint cumulative distribution function (JCDF) of bivariate flood attributes derived from the bestfitted 2-D copula structure. The best-fitted 2-D copulas are Table 4 Investigating the performance between 3-D symmetric vs FNA copula through comparing the theoretical estimates of Kendall's coefficient (derived from best-fitted 3-D symmetric Gaussian copula and fully nested Frank copula) vs empirical estimates of Kendall's correlation coefficient The asymmetric Frank copula exhibits minimum gaps between the empirical and theoretical Kendall's tau. In other words, it can regenerate the flood dependence structure more effectively than the symmetric 3-D Gaussian structure Also, the empirical Kendall's tau estimated in the above table considered the detrended samples of storm surge together with rainfall and river discharge observations in the estimation of pair-wise dependency. In other words, Table 2 measures the correlation statistics considering the original sample (before detrend or pre-whitening procedure). Supplementary Table S5 already confirms the existence of monotonic time-trend behavior within the storm surge time series. After removing trends from the flood, it is now incorporated in developing bivariate and trivariate copula dependency

Flood attribute pairs
Kendall's estimated from historical observations (Empirical estimates) Kendall's estimated from best-fitted asymmetric or fully nested frank copula (Theoretical estimates) Kendall's estimated from the best-fitted symmetric 3D Gaussian copula structure (Theoretical estimates) Annual maximum 24-h rainfall (  employed in the estimation of bivariate primary joint RPs for both AND-joint ( T AND ) and OR-joint ( T OR ) cases for the various possible combination of flood characteristics (refer to Eqs. 20 and 22). Table 6 presents the return period estimated using univariate marginal distributions of individual flood characteristics; and joint return periods for both "OR" and "AND" joint cases for bivariate distributions. It is observed that the joint return period (JRP) in the "AND-joint" case is longer than the JRP in the "OR-joint" case when the same univariate return period is assumed. In other words, the OR-joint RPs is smaller than AND-joint RPs for every possible combination of flood characteristics, i.e., T OR < T AND . For example, at 100-year RP, the Annual Maximum 24-h Rainfall is 183.688 mm, and the Maximum Storm Surge (Time interval = ± 4 days) is 0.389 m, the AND-joint and OR-joint RP is T AND = 292. 42 years and T OR = 60.31 years. Similarly, at 100-year RP, the rainfall is 183.688 mm, and river discharge is 7847.291 m 3 /s, the AND-joint and OR-joint RP is T AND = 1039.06 years and T OR = 52.52 years, and so on. From the above discussion, it is inferred that occurrence of bivariate flood characteristics simultaneously is more frequent in the "OR" case and less frequent in the "AND" case. The 3-D asymmetric Frank copula is finally selected after making a comprehensive comparison with symmetric Gaussian copula and thus employed in estimating trivariate joint return periods for both the "AND" and "OR" joint cases, using Eqs. 21 and 23. The effect of including a third variable in the estimation of JRPs is summarized in Table 6. For all the cases considering trivariate flood characteristics, the joint return period in the "AND" case T AND RAIN,STORM,RIVER is greater than the joint return period in the "OR" case T OR RAIN,STORM,RIVER . In other words, the occurrence of trivariate flood characteristics simultaneously (i.e., co-incidence occurrence) is more frequent in the "OR" joint case and less frequent in the "AND" joint case. For example, at 100-year RP, flood events characterized with Annual Maximum 24-h Rainfall is 183.688 mm, Maximum Storm Surge (Time interval = ± 4 days) is 0.389 m, and Maximum River Discharge (Time interval = ± 4 days) is 7847.291 m 3 /s, the AND-joint and ORjoint RP is T AND = 136.34 years and T OR = 34.19 years, and so on. Based on the above results, it is concluded that for a practical flood risk analysis, it is essential for the practitioners to take accountability for both the cases of trivariate joint return periods "OR-joint" and "ANDjoint" cases, along with the bivariate and univariate return periods. Considering either OR-joint or AND-joint return periods in the flood risk analysis would be problematic and might result in underestimating or overestimating the flood risk.

Table 6
Comparison of univariate, bivariate and trivariate return periods (RPs) for flood characteristics attributing compound flooding (CF) events

Estimation of conditional joint return periods
The conditional return period relies on a conditional probability distribution relationship between the variables of interest, given that some condition is fulfilled. The bivariate conditional joint return periods of rainfall given various percentile values of storm surge and river discharge observations are estimated for two conditional cases using Eqs. 24 and 25 and presented in Supplementary Figs. (SF18(a,  b) and SF (19(a, b)). It is inferred that, for the first conditional case (i.e., T RAIN�STORM>storm ) bivariate return periods increases with an increase in the percentile values of the storm surge observations. In other words, the curves attained higher conditional return of flood events at higher conditional storm surge values as compared to the lower storm surge values for the same specified values of rainfall series (referred to Figure (S18a)). But for the second conditional case (i.e., T RAIN�STORM≤storm ), the curve attained lower conditional return of flood events at higher conditional storm surges as compared to the lower storm surges for the same specified values of rainfall observations (referred to Fig. (S18b)). It is also inferred that the higher the rainfall, the higher the return period of flood events. Figures S19a, it can be observed that for the first conditional case (i.e., T RAIN�RIVERDISCHARGE>riverdischarge ), the curve attained a higher conditional return of flood events at higher conditional river discharge observations than the lower river discharge values for the same specified values of rainfall characteristics. On the other side, for the second conditional case (i.e., T RAIN�RIVERDISCHARGE≤riverdischarge ), the curve attained a higher conditional return of flood events at lower conditional river discharge observations than the lower river discharge values for the same specified values of rainfall characteristics (refer to Fig. S19b).
The joint return periods of two flood variables, i.e., rainfall and storm surge conditional to (or for the various percentile values) river discharge observations are estimated using the Eq. 26, and they are graphically illustrated in Supplementary Figs. (S20 (a-c)). It is observed that the return periods decrease in accordance with an increment in the percentile value of the river discharge observations. In other words, smaller conditional return periods of flood events at higher conditional river discharge observations as compared to the lower conditional river discharge series for the specified values of rainfall and storm surges jointly. For example, consider a flood events with the following flood characteristics; rainfall = 145.8 mm; storm surge = 0.68 m. Using Eq. 26, it is estimated that the conditional return period is 68.4 years (when river discharge ≤ 1275 m 3 /s). The return period is 46.7 years (when river discharge ≤ 1840 m 3 /s), and the return period is 27.4 years (when river discharge ≤ 4882.5 m 3 /s), and so on. Also, at the same time, the higher the rainfall and river discharge simultaneously, the higher is the return period of the flood events. In addition, the conditional joint return period of one flood variable (i.e., rainfall) given various percentile values of the other two flood variables (i.e., storm surge and river discharge) is estimated using Eq. 27 with the results presented in Fig. 4. It is concluded that the return periods decrease for the higher percentile values of both river discharge and storm surge observation. In other words, the curve attained a higher conditional return period at lower conditional river discharge and storm surges observations than the lower river discharge and storm surge values for the same specified values of rainfall characteristics (refer to Fig. 4). Let us consider a flood event with the following flood characteristics; rainfall is 115 mm; using Eq. 27, the estimated conditional return period is 13.4 years (storm surge ≤ 0.18 m and river discharge ≤ 1275 m 3 /s). The assessed return period is 10.4 years (when storm surge ≤ 0.42 m and river discharge ≤ 2405 m 3 /s), and so on. On the other side, the conditional return period of rainfall observations given various percentile values of the river discharge and constant storm Surge (i.e., 0.55 m) are shown in Fig. 5. It is concluded that the curve attained a higher conditional return at lower conditional river discharge, with a fixed value of storm surges observations, than the lower river discharge for the same specified values of rainfall characteristics. Similarly, Fig. 6 shows the conditional return period of the rainfall observations given various percentile values of the storm surge and constant river discharge (i.e., 3395 m 3 /s). It is observed that the curve attained a higher conditional return at a higher conditional storm surge with a constant river discharge value than the lower storm surge for the same specified values of rainfall events.

Hydrologic risk assessments using failure probability
The failure probability or FP statistic usually defines the chance of potential flood events occurring at least once in a given project design lifetime (refer to "Evaluating multivariate hydrologic risk via failure probability statistics"). First, the failure probability statistics are estimated to demonstrate the bivariate hydrologic risk of flood attribute pairs, i.e., rainfall-storm surge and rainfall-river discharge, using Eqs. 32 and 32. Supplementary Figures (S21) and (S22) indicate the bivariate hazard scenario for the "OR" and "AND" joint case in the assessments of hydrologic risk for flood variables pairs, rainfall-storm surge and rainfall-river discharge estimated for different return periods (i.e., 100, 50, 20, 10, and 5 years). It is inferred that the failure probability (FP) curve (or its estimated statistics) obtained via bivariate OR-joint case are sharply rising and attain higher values at a lower return period in comparison with bivariate AND-joint case. Also, the failure probability statistics are continuously increasing in accordance with the service design lifetime of the structure (where FP statistics estimated via bivariate OR-joint case increase faster than AND-joint case (less steep slope) for the same return period). The results shown in Supplementary Figs. (S23(a-e)) and (S24(a-e) are pointing out that it is concluded that first, the univariate analysis (or events) produces lower FP statistics than the bivariate one for OR-joint cases, which further reveals that neglecting the joint impact (compounding) of the rainfall and storm surge (or rainfall and river discharge) observations would be revealed for underestimations of failure probabilities. Second, the bivariate (via ORjoint case) and univariate hydrologic risk would decrease with the increase in return periods, thus revealing that the return period is not explicitly tied to a planning period and is ineffective in characterizing the chance of events occurring during a project lifetime (Read and Vogel 2015). Also, the bivariate hydrologic risk values would increase with the increase in the service design lifetime of the hydraulic infrastructure.
Supplementary Figures (S25(a-c)) and (S26(a-c)) illustrate the variation of the bivariate hydrologic risk with changes of rainfall observations in differently designed storm surge and river discharge observations. In both Figs. S25 (a-c) and S26(a-c), the amount of storm surge and river discharge observations are considered under return periods, i.e., 100 years, 50 years, 20 years and 10 years, estimated via the best-fitted univariate marginal probability distribution, where the project design lifetime (or service time of the hydraulic facilities) is assumed to be 100 years, 50 years, and 30 years. From Fig. S25 (a-c), it is clearly revealed that the bivariate hydrologic risk (compounding joint impact of rainfall and storm surge observations) increases with an increase in the project design lifetime (or service time) and decreases with an increase in the return period of storm surge observations. From the results shown in Fig. S26(a-c), it is inferred that the bivariate hydrologic risk (compounding the joint impact of rainfall and river discharge observations) increases with an increase in the project design lifetime (or service time) and decreases with an increase in the return period of river discharge observations. It is also observed that bivariate hydrologic risk decreases instantly when the design rainfall observation is 60 mm (referred to Fig. S26 (a-c)) when observing the probability of occurrence between flood pair, rainfall-river discharge. The bivariate hydrologic risk decreases instantly when the design rainfall observation exceeds 80 mm (referred to Fig. S25 (a-c)).
Because of the bivariate distribution approach's limitation and hydrologic risk assessments, we proposed and tested the adequacy of the fully nested Archimedean copula in the trivariate distribution modeling of the rainfall, storm surge, and river discharge observations. The failure probability statistics for trivariate hydrologic risk are estimated using Eqs. 29-31 (in Supplementary Figs. (S27), Fig. 7(a-e) and 8(a-c)). From Fig. S27, it is concluded that the FP statistics estimated for the trivariate OR-joint case attained a higher   value at a lower return period than the trivariate AND-joint case. Also, the failure probability statistics are continuously increasing with the service design lifetime of the structure. They are showing higher values at lower return periods in both trivariate OR and AND-joint cases. From Fig. 7a-e, it is inferred that the bivariate (also univariate) events produce lower failure probability values than the trivariate events for OR-joint cases. These results lead to the conclusion that neglecting the joint impact through trivariate analysis of the rainfall, storm surge, and river discharge observations results in underestimation of failure probabilities. Figure 8a-c illustrates the trivariate hydrologic risk under different rainfall-storm surge-river discharge scenarios. This figure represents the variation of the probability of occurrence or trivariate hydrologic risk with the change of rainfall values for different design storm surges and river discharge observations under different return periods, i.e., 100 years, 50 years, 20 years, 10 years and 5 years. It is concluded that the trivariate hydrologic risk values would increase with an increase in in-service time for the hydraulic facilities.
In other words, increasing the designated standard for the hydraulic facility would decrease the trivariate hydrologic risk. Also, at the same time, the hydrologic risk decreases with the increase in return periods. Figure 8a-c shows that the trivariate hydrologic risk initially increases quickly when the design rainfall is above 60 mm. The value would attain the highest value at a rainfall value of 120 mm. After that, risk statistics values decrease slowly. From the above discussion, it is concluded that simultaneous accounting of all three flood characteristics, e.g., rainfall, storm surge, and river discharge, significantly impact flooding in the coastal region. It can provide a better understanding of compound flooding and more critical information for flood risk assessments. These analytical and graphical investigations are vital for flood management strategies' sustainable design and planning in coastal regions.

Summary and conclusions
Compounding through the joint distribution analysis of storm surge, rainfall, and river discharge can significantly exacerbate the impact of flooding in low-lying coastal areas like the west coast of Canada. Due to the common forcing mechanisms that can drive these flooding variables concurrently or in close succession, their ignorance can lead to devastating disasters. Trivariate joint probability analysis can increase understanding of the probabilistic behavior of compound events compared to the bivariate (or univariate) distribution case. The probability of the joint occurrence of the multiple flood drivers will be higher than expected by considering each variable independently. Therefore, the conclusions of the assessments of hydrologic risk associated with compound events would be misleading. The present study introduces asymmetric or FNA copula in the trivariate joint analysis by compounding the impact of rainfall, storm surge and river discharge on the coastal flood risk. The incorporated 3-D asymmetric copula searches the dependency of the annual maximum 24-h rainfall and the highest value of the storm surge and river discharge visualized within a time lag of ± 4 days from the date of the highest 24-h rainfall event. The west coast of Canada is used as a case study. The adequacy of some frequently used symmetric Archimedean and Elliptical copula is also tested and compared with the FNA structure in trivariate flood dependency modeling. The most parsimonious trivariate copula structure is selected and employed to assess compound events. The following are the main findings of the performed work: 1. Zero first-order autocorrelation is identified within the time series of Annual maximum 24-h rainfall, Maximum Storm surge (Time interval = ± 4 days) and Maximum River discharge (Time interval = ± 4 days). Also, no significant monotonic time trends are identified within rainfall and river discharge series. The time series of Maximum storm surge (Time interval = ± 4 days) showed a significant time-trend, which is eliminated through pre-whitening before introducing it into the probability framework. 2. The graphical and analytical-based correlation measures confirmed that the strength of dependency among the targeted variables is positive and significant at the 5% significance level (or 95% confidence interval). This investigation confirms the possibility of incorporating a 3-D copula in the joint dependence simulation. 3. The parametric-family-based marginal distributions are fitted to each flood variable using different model selection criteria statistics, such as K-S, A-D, CvM, AIC, and BIC. It has been observed that the Generalized Extreme Value (GEV) best fits the rainfall series, Normal the storm surge, and Generalized Extreme Value (GEV) distribution the river discharge observations. 4. Particular classes of 2D copulas are used in modeling bivariate joint dependence between flood pairs, where the best-fitted copulas are selected through different analytical and graphical procedures. The tail dependence assessments also examine the utility of selected copulas. BB1 copula (for rainfall-storm surge), Husler-Reiss copula (for storm surge-river discharge), and Survival BB7 copula (for rainfall-river discharge) are identified as the best-fitted functions. The selected bivariate structures are employed further in estimating bivariate joint cumulative distributions (JCDFs) and joint and conditional return periods. They are further employed in estimating failure probability to perform hydrologic risk for bivariate compound flood hazard scenarios.
5. First, three fully nested forms of Archimedean copulas (viz., Frank, Clayton, Gumbel-Hougaard (GH) copulas) are used and tested in the trivariate analysis. Also, three symmetric (or monoparametric) Archimedean copulas (viz., Frank, Clayton, Gumbel-Hougaard and Joe) and one Elliptical class, such as the Gaussian copula, are introduced, and their performance compared with asymmetric copula structure. After a comprehensive analytical model compatibility investigation and graphical visual inspections, it is confirmed that the asymmetric Frank copula performs the best. The selected fully nested Frank copula can approximate heterogenous interdependence of compound flood characteristics much better than the selected symmetric Archimedean copulas. Refer to Fig. 3b, Supplementary Table (S13b) and Table 4; it is observed that selected asymmetric Frank copula reproduced or regenerated the mutual dependence of historical flood characteristics much more comprehensively than the selected best-fitted 3-D symmetric Gaussian copula. Minimum deviations or gaps are observed between empirical and theoretical Kendall's correlation statistics. 6. The return periods for trivariate (and bivariate) OR and AND-joint cases are estimated. It is found that the AND-joint case produces a higher return period than the OR-joint case for any possible combination of flood characteristics. It is concluded that estimating trivariate return periods of flood variables is essential to examine the expected flood risks and strength of influence if these variables co-occur or are in close succession. It is also observed that the importance of different return periods solely depends upon the nature of the problem undertaken and which cannot be interchanged. Also, their appropriate choice can depend on the impact of the design quantile. Suppose the hydraulic structure design includes only considers the joint return periods in the OR case (or AND case). In that case, it is highly likely that it could be under-dimensioned (or overdimensioned). Besides this, conditional joint return periods often play a vital role in water infrastructure design; thus, their accountability cannot be ignored. First, smaller return periods for a specified value of rainfall and storm surges jointly are observed at higher conditional river discharge values. It is also found that at the lower value of conditional variables, river discharge and storm surge, return periods are higher than those obtained at a lower value of river discharge and storm surge for the same specified value of the rainfall events. In addition, the return periods of rainfall events conditioning the storm surge (or river discharge) variable with the constant value of the river discharge variable (or storm surge) is also estimated. It is found that the trivariate return periods of rainfall attained a higher return value at the lower conditional river discharge with a fixed value of storm surge event. Similarly, the trivariate return period of the rainfall attained a higher value at the higher conditional storm surge with a fixed river discharge than return periods obtained at the lower storm surge for the same specified values of rainfall events. Therefore, for an effective risk-based design of waterrelated projects, we have to take accountability for both the joint and the conditional joint probability relationship between the variable of interest. 7. The trivariate hydrologic risk of rainfall, storm surge and river discharge are examined using the failure probability statistics, which can easily reflect the flood risk level during the entire project lifetime. Investigation reveals that the bivariate (and univariate) events produce lower failure probability than trivariate events for OR-joint cases, revealing that neglecting the joint impact through trivariate analysis of the rainfall, storm surge and river discharge observations leads to underestimations of failure probabilities. It is also concluded that the trivariate hydrologic risk values would increase with an increase in in-service time for the hydraulic facilities; at the same time, hydrologic risk decrease with an increase in return periods. The variation of trivariate hydrologic risk with changes in rainfall characteristics in differently synthesize storm surge and river discharge events are also investigated. In this investigation, three different project design lifetimes of infrastructure are considered (such as 100 years, 50 years, and 30 years). The variation of the bivariate hydrologic risk with changes in rainfall events in differently designed river discharge and storm surge are also investigated. This study has a few limitations. First, the available record length of the observed data in this study was only 46 years which may be a source of uncertainty in the estimated results. The proposed copula-based methodology has a few statistical limitations, such as parametric family functions being incorporated in the modeling of the univariate marginal distribution. No one could deny that, compared with the parametric-based approach, the efficacy of nonparametric via kernel density functions could be much more effective in approximating marginal structure. On the other side, the faithful preservation of lower-level dependency presented in this study can be further improved by incorporating a vine copulabased methodology. In reality, due to the higher degree of uncertainty and complex dependency, fully describing the dependence structure of multivariate extreme is quiet. It often demands a comprehensive way of uncertainty characterization by adopting a flexible methodology such as the vine or pair-copula framework. All the above-discussed statistical issues are being considered in the current work to be prepared for future publication.