Trivariate Probabilistic Assessments of the Compound Flooding Events Using the 3-D Fully Nested Archimedean (FNA) Copula in the Semiparametric Distribution Setting

Severe flooding in coastal areas can result from the joint probability of oceanographic, hydrological, and meteorological factors, resulting in compound flooding (CF) events. Recently, copula functions have been used to enable a much more flexible environment in joint modelling. Incorporating a higher dimensional copula framework via symmetric 3-D Archimedean or elliptical copulas has statistical limits, and the preservation of all lower-level dependencies would be impossible. The heterogeneous dependency in CF events can be effectively modelled via a fully nested Archimedean (FNA) copula. Incorporating FNA under parametric settings is insufficiently flexible since it is restricted by the prior distributional assumption of the function type for both marginal density and copulas in parametric fittings. If the marginal density is of a specific parametric distribution, it could be problematic if underlying assumptions are violated. This study introduces a 3-D FNA copula simulation in a semiparametric setting by introducing nonparametric marginal distributions conjoined with parametric copula density. The derived semiparametric FNA copula is applied in the trivariate model in compounding the joint impact of rainfall, storm surge, and river discharge based on 46 years of observations on the west coast of Canada. The performance of the derived model has also been compared with the FNA copula constructed with parametric marginal density. It is concluded that the performance of FNA with nonparametric marginals outperforms the FNA copula built under parametric settings. The derived model is employed in multivariate analysis of flood risks in trivariate primary joint and conditional joint return periods. The trivariate hydrologic risk is analyzed using failure probability (FP) statistics. The investigation reveals that trivariate events produce a higher FP than bivariate (or univariate) events; thus, neglecting trivariate joint analysis results in FP being underestimated. Moreover, it indicates that trivariate hydrologic risk values increase with an increase in the service time of hydraulic facilities.


Introduction
Extreme coastal flooding episodes are among the most significant hazards threatening infrastructure and human life (Bates et al. 2005;Padgett et al. 2008). Rapid urban expansion and conditions of climate change increase their severity or risk in many coastal cities worldwide (Milly et al. 2002;Zheng et al. 2014). Compound flooding (CF) events have a multidimensional stochastic character that can be comprehensively characterized through consideration of the joint probability relationship between multiple physical drivers. According to the Intergovernmental Panel on Climate Change (IPCC), a compound event (CE) can be defined as the occurrence of two or more extreme or non-extreme events in close succession or simultaneously (Seneviratne et al. 2012). For example, storm surges, high predicted astronomical tides, locally or remotely generated waves (swell), fluvial flooding, and pluvial flooding are four commonly occurring primary mechanisms that, when combined, can increase the risk of coastal flooding (Hendry et al. 2019). Most existing studies individually consider these four main drivers in coastal flood risk assessments. However, flooding in coastal regions is often characterized by multiple physical factors because they are naturally intercorrelated. If these concurrent events display statistical dependencies through a common forcing mechanism, the probability of their joint occurrence will be higher than that expected from separately considering the variables of each extreme.
Sea-level rise (SLR) can significantly impact coastal areas through four different phenomena: coastal flooding, coastal inundation, coastal erosion, and saltwater intrusion (Church et al. 2013). The SLR will cause changes in the intensity and frequency of extreme events due to the combined effects of high spring tides, storm surges, extreme precipitations, surface waves, and increased river discharge. Storm surges are the primary driver of coastal flooding (Resio and Westerink 2008). Combination with other hydrometeorological factors, such as high river discharge and/or precipitation, can result in devastation of coastal regions. More extreme coastal flood events are likely to result when rainfall and storm surges co-or sequentially occur (Zheng et al. 2013). Interdependency of driving factors (e.g., storm surge, river discharge, etc.) may arise as a result of having common meteorological forcing mechanisms. For example, rainfall and storm surges or high river discharge can be affected or driven by the same meteorological systems (tropical or extra-tropical cyclones) and result in compound events (CEs) (Svensson and Jones 2002;Paprotny et al. 2018). Ignoring their dependencies can significantly impact flood risk estimation. In previous studies, more attention has been given to investigating CEs through multivariate statistical methods, either accounting for the number of joint extremes or the existence of mutual concurrence between the proxy variables of different flood hazard types (e.g., Coles et al. 1999, Coles 2001, Svensson and Jones 2002Zheng et al. 2013 and references therein).
Different conventional parametric distributions, e.g., bivariate normal, gamma, etc., have been extensively studied in the modelling of extreme events and the analysis of their frequency (Krstanovic and Singh 1987;Yue 1999Yue , 2001 despite the existence of numerous statistical constraints and the lack of modelling flexibility. In current practice, the copulabased approach is widely accepted and has gained broader attention in modelling CEs and their hydrologic risk (Lian et al. 2013;Xu et al. 2014;Wahl et al. 2015;Moftakhari et al. 2017;Bevacqua et al. 2017;Paprotny et al. 2018;Xu et al. 2019;Moftakhari et al. 2019;Ghanbari et al. 2021, and references therein). In the existing relevant literature, mostly different classes of 2-D copulas have been incorporated; for instance, Lian et al. (2013) (between extreme rainfall and storm surge), Xu et al. (2014) (between extreme precipitation and storm tide), Wahl et al. (2015) (between rainfall and storm surge), Masina et al. (2015) (between sea levels and waves), Moftakhari et al. (2017) (between coastal sea level (with future SLR) and riverine flooding), Paprotny et al. (2018) (between rainfall, river discharge, storm surges and wave observations), Xu et al. (2019) (between rainfall and storm tides), and Zellou and Rahali (2019) (between storm surges and extreme rainfall). In a recent demonstration, Couasnon et al. (2020) modelled the risk of compound flooding potential arising from the joint probability of storm surges and river discharge in river mouths. On the other side, Ghanbari et al. (2021) introduced 2-D copulas in estimating compound coastal riverine frequency under current and future climate scenarios to capture the impact of climate change.
Bivariate joint analysis would be insufficient and less practical in hydrological risk assessments. CEs are usually more strongly correlated with multiple flooding drivers. For instance, a low-pressure system, such as a tropical cyclone, can co-occur with storm surges, heavy rainfall, and possible river overflows, leading to flooding (Couasnon et al. 2020). A better and more comprehensive understanding of uncertain behaviour of the CE can be achieved by simultaneously integrating all the intercorrelated drivers instead of visualizing their bivariate joint probability occurrence. The potential damage could likely be a function of multiple relevant vectors of specified hydrological events, where ignorance of spatial dependency might be attributed to underestimating the uncertainty. A limited number of studies have performed trivariate distribution via 3-D copula density and pointed to the importance of a trivariate return period over bivariate or univariate return periods. For instance, Serinaldi and Grimaldi (2007) (via parametric-based fully nested Archimedean (FNA) for flood characteristics), Genest et al. (2007) (via trivariate meta-elliptical copulas for annual spring flood analysis), Kao and Govindaraju (2008) (incorporated Plackett family for rainfall events), Wong et al. (2010) (3-D Gumbel and Student's t copulas for drought characteristics), Reddy and Ganguli (2013) (FNA and trivariate Student's t copula for flood characteristics), and Fan and Zheng (2016) (entropy copula for flood characteristics).
The adequacy of symmetric 3-D Archimedean copulas in higher-dimension modelling is questionable because of the involvement of single copula dependence parameters. It would be incapable of preserving of all pairwise mutual concurrency at the lower stages, where all dependencies are averaged (Genest et al. 2007). CF is a mixture of multiple extreme or non-extreme events in which complex dependency between the variables of interest is exhibited. The asymmetric FNA structure is highly flexible compared with symmetric Archimedean copula density in constructing higher dimension dependence modelling, and each random pair is individually approximated through multiple parametric joint asymmetric functions (Serinaldi and Grimaldi 2007;Savu and Trede 2010;Reddy and Ganguli 2013). FNA structure comprises the joint integration of two or more 2-D Archimedean copula structures through another Archimedean structure. It can thus easily capture different mutual concurrency between and within different groups of random variables (Whelan 2004;Savu and Trede 2010;Hofert and Pham 2013).
Before pointing out the existing limitation of FNA copula distribution modelling, which lacks flexibility, there is a need to highlight some critical facts about the copula-based methodology. Copula distribution can be developed using parametric, semiparametric, and nonparametric fitting procedures, depending on how univariate marginals and multivariate join dependence structures are defined. Most of the above-referenced approaches are often implemented under parametric settings, with prior subjective assumption of the function type for marginal PDFs and 2-D copula density in parametric fitting. Very few studies have incorporated 2-D parametric copulas with nonparametric marginal density, i.e., semiparametric distribution settings (Karmakar and Simonovic 2009;Reddy and Ganguli 2012;Latif and Mustafa 2021). A third approach, called nonparametric copula framework, comprises the modelling of univariate marginals via nonparametric distribution and their joint dependence structures using the nonparametric copula density, such as beta kernel copula density (Charpentier et al. 2006;Rauf and Zeephongsekul 2014;Latif and Mustafa 2020).
Nonparametric kernel density estimation (KDE) is a data-driven model with higher flexibility than parametric distributions and can be directly retrieved from a given distribution series. It represents a much more stable data smoothing approach for characterizing the uncertainty in hydrologic characteristics without making any prior hypothesis about the form of the probability density function PDF (unimodal or multimodal) (Adamowski 1985;Moon and Lall 1994). It has also been shown that, in practice, it is difficult to approximate the distribution tail beyond the most significant values under the parametric distribution framework (Bardsley 1988). More significantly, in the case of skewed or multimodal distributions, the parametric models might lead to inconsistencies in the estimated extreme quantiles. In the existing literature, for instance, storm surge observation modelling is rarely incorporated using a nonparametric distribution approach. Besides this, a few examples in the literature, such as Reddy and Ganguli (2012) and Latif and Mustafa (2020), have already demonstrated how effectively the nonparametric KDE approach captures univariate marginal probability distribution of flood peak flow, flood volume and flood duration, outperforming parametric models.
In the previous studies, the FNA copula framework has been constructed using the standard parametric distribution approach (Grimaldi and Serinaldi 2006;Reddy and Ganguli 2013). In other words, the parametric distribution fitting estimates the univariate marginals PDFs and composition of simple parametric bivariate copulas (the inner and outer copula in the FNA structure) must belong to the Archimedean class. They hold a prior subjective assumption for marginal PDFs; thus, it would be difficult and questionable to fully acknowledge the flexibility of FNA copulas in the trivariate joint density. The copulas already eliminate the restriction of selecting flood marginals from the same parametric family and allow different types of marginals to coexist and join together to form a multivariate joint distribution. Therefore, in our study, we make a methodological attempt to construct FNA copulas by allowing nonparametric KDE to be used in modelling marginal density. Introducing the distribution free-based KDE in the FNA copula simulation, also called the semiparametric fitting approach, can facilitate higher flexibility and more accurately approximate the multivariate joint density of hydrologic events.
In Canada, flooding events are becoming one of the costliest natural extremes (Public Safety Canada 2022). For instance, Hurricane Igor in the year 2010 (due to the joint impact of excessive rainfall, swollen rivers, storm surge, and wave heights, with peaks up to 2.5 m) impacted 97 communities and caused at least CAD 100 million in damage to public infrastructure in Newfoundland, Canada (Pirani and Najafi 2020). According to the studies by James et al. (2014) and Atkinson et al. (2016), most of the east and west coast regions and the Beaufort Sea coast of the north coast region in Canada are experiencing a sea-level rise (SLR) that will, throughout this century, have significant impacts through coastal changes. On the west coast of Canada, most of the residents of the province of British Columbia (BC) reside nearer to the coastline. The problem of SLR can increase the risk of tidal or extreme water level events, which could pose a high risk to low-lying coastal communities.
There is a need to develop and investigate a trivariate hazard framework that can effectively capture and outline the likelihood of compound events due to the joint occurrence, either simultaneously or successively, of multiple drivers of the risk of CF at Canada's coasts. The primary methodological importance provided in this manuscript is the construction of an asymmetric or FNA copula framework through semiparametric distribution fitting to investigate the compound effects of rainfall, storm surges, and river discharge series on coastal flood risk. This proposed modelling approach uses the nonparametric marginal PDFs of flood-contributing variables. In eliminating prior subjective hypotheses about the flood marginal distribution type, co-joined by parametric Archimedean class copulas, it can provide better flexibility than traditional FNA copula simulation. By integrating the efficacy of nonparametric kernel approximation with a lack of distribution assumption, the resulting FNA copula can better model the uncertainty characteristics of compound extreme.
The main objective of this study is addressed through (1) the development of the 3-D FNA copula framework in the semiparametric fitting procedure and comparison of its performance with that of FNA copula with the parametric fitting approach in compounding the joint impact of rainfall, storm surge, and river discharge observations based on 46 years of observations of the west coast of Canada and (2) evaluation of the trivariate return period (RP) and failure probability (FP) in assessing the trivariate hydrologic risk. These estimated trivariate joint and conditional RPs can provide vital information for flood risk assessments in low-lying coastal regions. The statistical procedure for the nonparametric KDE in the univariate marginal distribution and fitting a 3-D FNA copula framework under semiparametric settings is discussed in the following section of the manuscript. In the third section, the copula-based methodology is applied in a case study to establish trivariate joint dependency modelling of multiple flood characteristics and the joint and conditional return periods. The performance of the proposed semiparametric FNA model is compared with that of the FNA copulas with parametric distribution settings. Subsequently, the trivariate return periods and failure probability statistics are estimated. The last section of the paper provides the research summary and conclusions, which can serve as a reference for the design of hydraulic facilities and establishing coastal flood control and mitigation practices.

Nonparametric Estimation of Marginal Flood Behaviour
Estimating the risk of CF events in a coastal region is necessary for developing policy decisions, disaster risk reduction, resilience building, and engineering practice. In the present study, we developed and tested the adequacy of the semiparametric approach in 3-D FNA copula simulation. The proposed 3-D model is applied in investigating the compound effects of rainfall, river discharge, and storm surge events on trivariate flood risk in the coastal zone. The fitted nonparametric kernel marginal PDFs have been employed in constructing fully nested 3-D copulas by combining two parametric Archimedean copulas that were obtained by parametric fitting. The performance of the developed semiparametric FNA model is compared with that of the FNA copulas constructed using the parametric fitting approach (both marginal PDFs and copulas are in parametric settings). Figure 1 illustrates the methodological workflow details of this study.
In parametric functions, it has always been assumed that the random characteristics come from known populations whose probability density functions PDFs are predefined. However, the best-fitted marginal distributions may not be from the same probability distribution family (Adamowski 1985;Kim and Heo 2002). Nonparametric distribution approaches do not require prior assumptions and have higher flexibility than parametric density estimators (Adamowski 1996;Moon and Lall 1994).

Univariate Kernel Density Estimation or KDE
The univariate KDE is a nonparametric approach to approximate the PDFs, f (x), of given random observations of ′ x ′ or compound flooding drivers such as rainfall, storm surges, or river discharge observations, where inferences about the flood populations are made based on the finite samples having the following statistical property where K(x) represents the univariate kernel density function and can be used as a PDF. The kernel function can be approximated through general equations (Hardle 1991) as where 'h' is the bandwidth of the kernel function, which helps to regulate the smoothness and roughness of the shape of the estimated density functions or PDF (Moon and Lall 1994). Mathematically, if X 1 , X 2 , X 3 , … , X n are the n th set of independent and identically distributed (or i.i.d) random series having the density function f (x) , then the univariate KDE of f (x) are obtained by averaging Eq. (2) of the given random observations (Kim and Heo 2002) as where 'n' is the number of random observations; X i is the i th observation and f h (x) is the kernel density estimate.
The efficiency of the KDE often depends upon two factors: (1) the appropriate choice of the kernel bandwidth; and (2) the selection of a proper kernel function. Table 1 lists the different types of commonly used 1-D kernel functions (along with their PDFs) used in this study.

Estimation of Kernel Bandwidth
The appropriate choice of the kernel smoothing parameter or bandwidth 'h' is crucial because the shape of kernel density estimates are dependent on this (Moon and Lall 1993;Efromovich 1999). Insufficient smoothing can result in a rough density, while oversmoothing can bypass or smooth over essential features Santhosh & Srinivas 2013). Different articles have pointed out different statistical approaches for kernel bandwidth estimation; for instance, Tarboton et al. (1998) (least square cross-validation LSCV), Scott and Terrell (1987) (biased cross-validation), Wand and Jones (1995) (Plug-in estimate), and Silverman (1986) (rule of thumb, ROT), among others.
The plug-in bandwidth estimators are promising, conceptually simple, theoretically justifiable, and perform well (Wand and Jones 1995). In this study, the direct plug-in (DPI) method is used to select the bandwidth N, where unknown functionals that appear in expressions for the asymptotically optimal bandwidths are replaced by kernel estimates (Sheather and Jones 1991;Chen 2015). Mathematically, a natural estimate of R(γ(y)f (y)) is given as In addition, ∅ = ∫ γ(y)f �� (y)f (y)dy can be estimated by Hence, the direct plug-in bandwidth selector for U = ∫ (y)f 2 (x)dx can be obtained by replacing R( (y)f (y) with R , and ∫ γ(y)f �� (y)dy with ∅ , which is estimated by It must be noted that the choice of g, also called the pilot bandwidth, is decided by the kernel density estimation obtained in Eqs. (4) and (5). Chen's (2015) study reveals that it could be justifiable to calculate g using a normal scale bandwidth selector.

3-D FNA Copula Simulation in Semiparametric Settings
The idea of the copula function was developed by Saklar (1959). It models the multivariate distributions to their univariate marginal distribution functions, which are not necessarily from the same family of probability functions. It captures a broader extent of dependencies (both linear and nonlinear) and preserves their joint dependence structure (Salvadori and De Michele 2004;Nelsen 2006). If H X,Y (x, y) represents the joint distribution function of bivariate random observations (X, Y) along with their univariate marginal distributions u = F X (x) and v = F Y (y) , then there is a copula function C(⋅) between these marginals, which can be determined as In Eq. 7, if univariate marginals are continuous, copula C(⋅) is unique for the particular bivariate joint distribution. The probability density functions (PDFs) can be obtained by integrating Eq. (7) as Symmetric 2-D Archimedean copulas are frequently applicable in the case of bivariates. Different copulas are available in this class and have a more comprehensive range of dependencies measuring capabilities (Nelsen 2006). The performance monoparametric Archimedean class copulas is inconsistent and restricted when projecting into highdimensional distribution. For instance, the symmetric 3-D Archimedean copulas would be ineffective in preserving all pairwise mutual concurrency at the lower stages, since it typically depends on a single dependence parameter of the copula generator function. Therefore, it could be more practical to approximate the complex joint structure of compound events by incorporating multiple parametric joint asymmetric functions called FNA copulas (Whelan 2004). In this approach, each random pair is approximated (4) R = � γ 2 (y)f(y)dF n (y) = 1 n(n − 1) individually instead of modelling the dependency of all the intercorrelated variables using the same dependence parameters.
The traditional approach to building the asymmetric or FNA copula framework is usually based on parametric distribution settings, where parametric family distributions are used to estimate the univariate marginals PDFs, and the bivariate parametric copulas in nesting form, which must belong to the Archimedean class, are employed in establishing trivariate dependency. Using this methodological concept, it would be challenging to fully acknowledge the flexibility of FNA copulas without eliminating the restriction of approximating their marginal behaviour only via parametric family functions with prior distributional assumption. The main contribution of our study is in developing a 3-D FNA framework under the semiparametric distribution fitting procedure. The proposed model allows the nonparametric KDE to capture the flood marginals instead of parametric PDFs, which are further co-joined by the composition of two Archimedean class 2-D copulas in the parametric fitting approach.
The heterogeneous dependency between the CF variables can be effectively modelled using the nested form for Archimedean FNA structure, which would be impractical via a symmetric 3-D copula structure. The FNA structure is usually modelled by two or more ordinary bivariates, or any higher dimensional Archimedean copulas, and by another Archimedean copula, also called parent Archimedean copula (Reddy and Ganguli 2013). Mathematically, the dependence structure of FNA for random trivariate observations can be expressed as (Savu and Trede 2010) where variables S, R, and RD represent storm surge, rainfall, and river discharge, respectively. From Eq. 9, the first derivative of −1 2 • 1 is completely monotonic where 1 and 2 are Laplace transforms, and the symbol ' • ′ represents the composite function. Moreover, the pair ( s, r) has the bivariate marginal structure in the form of Eq. 7 with Laplace transform 1 . Similarly, random pairs ( s, rd) and (r, rd) have the bivariate margins in the form of Eq. (7) with Laplace transform 2 . Figure 2 shows the schematic diagram of the nesting procedure in constructing the 3-D FNA copula in the semiparametric distribution setting. In Fig. 2, C 1 (F R (r), F S (s);θ 1 ) is the inner copula with copula dependence parameter 1 . The marginal behaviour of both flood characteristics is modelled using the nonparametric kernel density estimation KDE, and the 2-D copula is approximated under standard parametric settings. Similarly, C 2 (C 1 , F R.D. (rd);θ 2 ) is the outer copula with copula dependence parameterθ 2 , where nonparametric kernel density is used to approximate the marginal behaviour of third variables RD. The formulated 3-D copula C 1 (F R (r), F S (s), F R.D. (rd)) represents the joint simulation of two bivariate parametric copulas with nonparametric marginal structures through trivariate asymmetric Archimedean functions. The applicability of the FNA copula is only valid if the dependency strength among the two variables, (s, r), dominates the correlation structure between these variables and the third variable, (s, r) and (r, rd) (Savu and Trede 2010).
The performance of the parametric-based FNA copula is analytically and graphically compared with the derived semiparametric FNA model in building a trivariate flood dependence structure. Besides the methodological importance in this study, the estimation of failure probability statistics in the trivariate compound flood hazard can provide a practical approach to coastal hazard assessments for hydrologists and water experts.

Selecting Compound Flooding Variables
The intention of this research is to improve the flexibility of the multivariate FNA copula framework by introducing the nonparametric KDE in compounding the joint probability occurrence of storm surge, rainfall, and river discharge observations. It is assumed that the performance of FNA copulas incorporated in the previous studies (under parametric distribution settings) can be improved by removing the restriction of prior distributional assumptions in approximating marginal PDFs. The proposed framework is applied to a case study on the west coast of Canada. First, the observed sea-level (or coastal water level) data are collected by a tidal gauge station located at New Westminster (station ID 7654), BC, Canada, with geographical location coordinates of 49.2 • N latitude and 122.91 • W longitude , and provided by Fisheries and Oceans Canada (https:// tides. gc. ca/ eng/ data) from 1970 to 2018. Second, the predicted coastal water level or tidal data are obtained from the Canadian Hydrographic Service (CHS). The storm surge data are obtained by differencing high tide values from the observed CWL. The calculated storm surge values are either positive or negative depending on whether the predicted high tide is above or below the observed CWL. The calculation of storm surge often demands a proper time synchronization between observed CWL and high tide values. The CWL data are collected at the local datum, called chart datum (CD), with the geolocation at Fraser River. The streamflow discharge and rainfall gauge station selection are solely based on having a proximity within the radial distance of 50 km from the targeted tidal gauge station, where the time frame must be the same. It should be noted that the selected station must not have less than 15% of missing observations. Daily streamflow discharge observations are collected at the Fraser River at Hope gauge station and delivered by the Environment and Climate Change Canada (https:// water office. ec. gc. ca/ search/ histo rical_e. html). The geographical location of this gauge 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.
The annual maximum 24 h rainfall events are defined using the daily-basis 24 h rainfall observations. The highest (or maximum) storm surge and river discharge records are observed within a period (or time lag) of ±1 day from the date of the annual maximum 24 h rainfall event for each calendar year and then used in establishing trivariate joint dependency modelling. Supplementary Table ST

Modelling Univariate Flood Marginals Using Nonparametric and Parametric Distribution
Before fitting the univariate probability density, conducting stationarity and serial correlation tests is a prerequisite (Daneshkhan et al. 2016). For this, Ljung and Box (1978)-based hypothesis testing, also known as Q-statistics, is used for investigating the autocorrelation behaviour within individual series. Test results revealed that the data exhibited zero firstorder autocorrelation (refer to Supplementary Table ST 2). No pre-whitening is required for all three random observations. The nonparametric Mann-Kendall (M-K) test is used to examine monotonic timetrend behaviour within individual time series (Mann 1945;Kendall 1975). Investigation reveals the presence of monotonic time trend behaviour only within storm surge observations, as its calculated z value exceeds the critical value of Z-statistics ( Z estimated = 2.85(p − value = 0.004) > Z critical = ±1.96 at 5% significance level or 95% confidence interval-refer to Supplementary Table ST 3). In Supplementary Fig. SF 3, a positive trend in storm surge observations is clearly seen. The storm surge observations require a differencing procedure to detrend before being introduced into the univariate probability distribution framework.
Several distributions often fit the data equally well, but each would give different estimates of a given quantile, especially in the distribution tails. Table 1 lists some standard 1-D kernel functions employed and tested in modelling the marginal distribution. The bandwidth of the candidate kernel functions is estimated using the direct-plug-in (DPI) method (refer to Eqs. (4-6)) (Sheather and Jones 1991;Wand and Jones 1995). The nonparametric kernel density approximations do not facilitate a closed form of the probability density function PDF and cumulative distribution function CDF; thus, CDFs are estimated through an empirical procedure that is based on numerical integration, as pointed out by Kim and Heo (2002) and Kim et al. (2006). The fitted univariate kernel density functions are listed in Table 2. The data-reproducing capabilities of the fitted KDE with observational samples are tested by comparing the theoretical cumulative density of each flood series with the empirical non-exceedance probabilities. The empirical non-exceedance probability is estimated using the Gringorten plotting position formulae (refer to Gringorten (1963)). Error index statistics such as mean square error (MSE) and root mean square error (RMSE) (Moriasi et al. 2007) are estimated to check the consistency of the fitted nonparametric models (refer to Table 2). In addition, information criterion statistics, such as Akaike information criterion AIC (i.e., Akaike 1974), Schwartz's Bayesian information criteria (or BIC) (i.e., Schwarz 1978) and Hannan-Quinn information criteria (HQC) (Hannan and Quinn 1979) are also used to check the fitness level of the fitted nonparametric models (refer to Table 2). It is observed that the normal kernel functions outperform others (exhibited minimum value of MSE, RMSE, AIC, BIC, and HQC statistics) and are thus selected for fitting the univariate marginal distribution of all three flood variables. A qualitative assessment by visual inspection is also carried out for each flood characteristic using cumulative distribution plots and probability-probability (P-P) plots (refer to Supplementary Figs. SF 4a-c and SF 5a-c). The findings of the analytical-based fitness test results are confirmed through visual inspection.
The adequacy of some frequently used parametric family distributions is also tested and compared to the nonparametric KDE. The vector of unknown statistical or model parameters is calculated using the maximum likelihood estimation (MLE) procedure, and their estimated parameters are listed in Table 3. The MLE estimators usually facilitate minor sampling variance in their estimated parameters of fitted distributions (Owen 2008). Different analytical goodness-of-fit (GOF) test statistics are used to check the consistency of the fitted candidate parametric models, such as the Kolmogorov-Smirnov (K-S) test (Xu et al. 2015) and Anderson-Darling (A-D) test (Anderson and Darling 1954;Farrel and Stewart 2006) and Cramer-von Mises (CvM) criterion (Cramér 1928;von Mises 1928). The results are summarized in Table 4. GEV distribution is the best-fitted distribution for defining the marginal distribution of annual maximum 24 h rainfall, normal distribution for storm surge, and GEV distribution for river discharge observations (the minimum value of K-S, A-D, and CvM test statistics). Supplementary Table ST 4 shows a comparison of the fitness level between the best-fitted nonparametric KDE and parametric functions. It is observed that, for all three flood variables, nonparametric kernel functions dominate parametric functions in defining their marginals (the minimum value of MSE, RMSE, AIC, BIC and HQC statistics). This investigation concludes that the KDE approximation of univariate marginals is much more effective without imposing any predefined density or distributional assumption than parametric distribution settings.

Bivariate Joint Exceedance Probability via Parametric 2-D Copulas with Nonparametric Marginals
In estimating trivariate joint return periods (JRPs), one must first evaluate bivariate joint cumulative distribution functions (CDFs) between each flood attribute pair. At the initial step, the strength of dependency between target flood variables is measured using the parametric-based Pearson's linear correlation (r) and two nonparametric rank-based dependence measures, such as Kendall's tau ( ) and Spearman's rho (ρ) correlation coefficients (Klein et al. 2011). Kendall's correlation measures can describe a broader class of mutual concurrency and is more suitable for describing the dependence structure of extreme hydrological events (Mirabbasi et al. 2012). Table 5 lists dependence measures between flood characteristics (estimated at a significance level of 5% (0.05) or 95% (0.95) confidence interval), concluding that all three flood variables exhibit positive dependency. In this study, 15 different bivariate copulas such as monoparametric Archimedean copulas (Frank, Clayton, Joe, and Gumbel); rotated (by 180 degrees) versions of Archimedean copulas (survival Clayton, survival Joe, and survival Gumbel); mixed versions of Archimedean copulas (BB1, BB6, BB7, and BB8); rotated (by 180 degrees) versions of the biparametric copulas (survival BB1, survival BB6, survival BB7, and survival BB8) are implemented and tested in modelling of bivariate joint dependence structure, which is further employed in deriving bivariate joint CDFs. The copulas' dependence parameters are calculated using the maximum pseudo-likelihood (MPL) estimation procedure (Kojadinovic and Yan 2010;Reddy and Ganguli 2012), and their values are listed in Table 6(a-c). Mathematically, MPL estimates the copula dependence parameter by maximizing the pseudo-log-likelihood function, which can be 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 (CDF). The fitness level of 2-D parametric copulas is investigated using the Cramer-von Mises (CvM) distance statistics via the parametric bootstrap procedure with two different bootstrap samplings (number of bootstrap samples N = 1000 and 500) (Genest and Rémillard 2008;Genest et al. 2009). The CvM test is based on a parametric bootstrapping approach and uses the Cramer-von Mises functional statistic S n , which can be mathematically estimated by where the p-values for each candidate copula are estimated using the parametric bootstrapping procedure (Genest and Rémillard 2008), which is mathematically expressed by The acceptance or rejection of the targeted parametric 2-D copulas is solely based on the S n test statistics and their estimated p-values. The minimum value of the S n test must indicate the minimum gap between empirical and derived parametric copulas. In addition, the estimated p-value must be larger than a significance level ( = 0.05) . From the results presented in Table 6(a-c), it is observed that BB7, Gumbel-Hougaard, and survival BB7 (BB7 rotated 180 degrees) copulas are selected as the best-fitted copulas in defining the bivariate joint dependence structure between rainfall and storm surge, storm surge and river discharge, and rainfall and river discharge observations. All three selected copulas exhibited minimum S n test values, where the estimated p-value is greater than 0.05 for both numbers of bootstrapping samples, N = 500 and 1000.
The adequacy of the best-fitted 2-D copulas is further examined via tail dependence coefficient (TDC) estimations (refer to Table 6(a-c)). Tail dependence assessment is mandatory in multivariate extreme event modelling or frequency analysis. If not well preserved by the selected 2-D copula, it may be attributed to higher uncertainty in estimated quantiles or wrong decisions in the design of the water-related infrastructure (Poulin et al. 2007). In the hydrologic frequency analysis, the upper tail dependence coefficient (UTDC) assessments are often of more interest than the lower tail dependence coefficient (LTDC). The nonparametric estimates of the UTDC, also called empirical estimates, are evaluated by estimators suggested by Capéraà et al. (1997) (Poulin et al. 2007), CFG up , and their results are compared with the parametric coefficient of UTDC up . Referring to Table 6(a-c), it is observed that all the selected 2-D copulas can effectively capture the extreme portion   (upper tail dependence behaviour) of each flood attribute pair, where the gap between the up and CFG up is at minimum. The selected copulas better estimate the upper tail limit than empirical estimates, although some minor gaps still exist between them. Finally, the selected 2-D copulas with nonparametric normal KDE marginals are employed to estimate the bivariate joint cumulative distribution function (JCDFs) and their associated joint return periods.

Trivariate Flood Dependence Using FNA Copulas with Nonparametric Marginals
The higher dimensional expressions in joint modelling are only available for a few copula families and do not possess sufficient flexibility in modelling multidimensional complex hydrological phenomena. The FNA copula can individually approximate each pair of random flood variables by incorporating multiple parametric joint asymmetric functions. It would be challenging to better assess the flexibility of the FNA copula framework without exempting the restriction of approximating their marginal behaviour only via the parametric distribution approach. From Section 3.2, the performance of nonparametric KDE is found to be much better than that of the parametric models. This study tests five asymmetric or fully nested forms of Frank, Clayton, Joe, AMH (Ali-Mikhail-Haq), and Gumbel-Hougaard (G-H) copulas to model trivariate flood dependence structure using the semiparametric fitting approach. For the extended mathematical details regarding the densities of nested Archimedean copulas incorporated in this study, readers are referred to Hofert and Pham (2013).
The normal KDE are selected as the most appropriate density for all three flood variables (refer to Table 2 and Section 3.2). The dependence parameters of copulas, both inner and outer, in the FNA structure are calculated using the full maximum likelihood estimation procedure (Okhrin and Ristig 2014;Okhrin et al. 2015), where all the computation is carried out using the HAC library package of R software (Okhrin 2020). The performance of the fitted 3-D FNA copulas is tested using the Akaike information criterion (AIC) (i.e., Akaike 1974), Schwartz's Bayesian information criteria (or BIC) (i.e., Schwarz 1978) and Hannan-Quinn information criteria (HQC) (Hannan and Quinn 1979). Table 7 shows that the fully nested Frank copula has a minimum value of AIC, BIC, and HQC test statistics (also, higher value of log-likelihood (L-L) statistics), which indicates a better fit. A visual comparison is also carried out using a scatter plot between the simulated (N = 1000) observations derived from the best-fitted 3-D nested Frank copula (with nonparametric marginals PDF), which confirms the finding (refer to Supplementary Fig. SF 6b). It is clear that the selected 3-D copula performed satisfactorily, since the generated samples (light blue colour) adequately overlap with the natural dependence of historical observations (red colour in Supplementary Fig. 6b). Supplementary Fig. (SF 6a) shows a 3-D scatterplot using random samples (size, N = 1000) generated from the best-fitted FNA Frank copula with nonparametric marginals.
The performance of parametric FNA copulas is also tested and compared with the above-incorporated semiparametric FNA copula density. The marginal distribution of individual flood characteristics is modelled under the parametric fitting procedure, i.e., the MLE estimator (Tables 3 and 4). The dependence parameters of fitted FNA parametric copulas are estimated using the full maximum likelihood estimation (or full MLE) procedure. The obtained values are listed in Table 8. It is concluded that 3-D fully nested via Gumbel-Hougaard copulas (with parametric marginals such as GEV, normal, and GEV distributions) results in the minimum value of AIC, BIC, and HQC statistics (with Table 6 Estimating dependence parameters of fitted 2-D parametric copulas via Maximum pseudo-likelihood (MPL) estimation procedure, their goodness-of-fit (GOF) test statistics and tail dependence assessments for flood pairs (a) rainfall-storm surges (b) storm surges-river discharge (c) rainfall-river discharge   ).
Also, BB7 copula exhibited a minimum gap between parametric and nonparametric estimates of upper tail dependence coefficient (where nonparametric or empirical estimates of upper tail dependence λ CFG up estimators are suggested by (Capéraà et al. 1997 andFrahm et al. 2005) (b) Gumbel-Hougaard copula (indicated by bold letter with an asterisk) exhibit a minimum value of S n Goodness-of-fit test statistics (p-value greater than 0.05) is identified as the most parsimonious copula in deriving the bivariate joint probability simulation between the Maximum Storm surges (time interval = ±1days ) and Maximum River dis-  (-12.757). In other words, nested Frank copula with nonparametric marginal PDF exhibits the minimum value of AIC, BIC, and HQC statistics. Similarly, graphical-based visual inspection is used to compare the empirical copula and the theoretical copula developed under the parametric and semiparametric approaches. Mathematically, the d-dimension empirical copula can be formulated as In Eq. (13), F j (⋅) indicate the estimated univariate marginal distribution of the random variable Y i .
The empirical probabilities are estimated using Eq. (13) and then compared with the fitted probabilities obtained from the FNA copulas derived in both semiparametric and parametric settings (refer to Fig. 3a, b). The semiparametric FNA model fitted to flood characteristics exhibits a minimum gap between empirical and theoretical probabilities; most observations lie near the diagonal inclined at 45 degrees to the axis. The performance of both developed models is also analytically compared using the error index statistics mean square error (MSE) and root mean square error (RMSE) (refer to Singh et al. 2004;Moriasi et al. 2007;and Table 9). Besides this, mean absolute error (Singh et al. 2004) and Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe 1970) are also estimated to measure the model performance (refer to Table 9). The MAE statistic is solely based on the absolute value to minimize bias towards large events relative to the RMSE (Bennett et al. 2013). The study of Willmott and Matsuura (2005) revealed that use of MAE represents a more effective approach than that based on RMSE. Chai and Draxler (2014) pointed out that RMSE performance can only dominate over MAE statistics in samples characterized by Gaussian distribution. Minimum, MSE, RMSE, and MAE values must indicate a better fit. Apart from this, NSE-based model evaluation criteria are used for comparison of the data and residual variance structure numerically defined within a range of -∞ (i.e., inferior performance) to 1 (i.e., ideal fitness level) (Nash and Sutcliffe 1970). Examples in the literature, such as Sevat and Dezetter (1991) and Legates and McCabe (1999), demonstrate the fitness measuring strength of NSE statistics in model performance evaluation, especially during hydroclimatic model simulations. Table 9 shows that use of the semiparametric approach in the FNA copula results in a minimum value of MSE, RMSE, and MAE and a high value of NSE test statistics.  The reliability of the developed nested Frank copula is further investigated by comparing the theoretical Kendall's correlation coefficient obtained from the best-fitted copula with the empirical Kendall's coefficient obtained from the historical flood characteristics (refer to Table 10). It is observed that the nested Frank copula exhibits a minimum gap between empirical and theoretical Kendall's value in comparison with the 3-D nested Gumbel copula developed under parametric fitting. In other words, compared with the other multivariate copula joint approaches, the selected 3-D Frank FNA copula accurately regenerates the historical flood characteristics' mutual dependence (or correlation) structure with nonparametric marginals. Figures 4 illustrates the derived nested Frank with KDE marginal density in approximating the trivariate dependence. In conclusion, the developed FNA copula with KDE marginal PDFs provides a much more efficient modelling approach by eliminating the restriction of distributional assumptions (or predefined PDFs). Therefore, it is recommended for modelling in cases of a trivariate flood dependence structure. The derived 3-D copula is further employed to estimate joint and conditional return periods and failure probability for assessing the hydrologic risk associated with compound flooding events.

Estimating Joint Return Periods
In the frequency analysis of extreme hydrologic consequences, the interassociation between extreme event quantiles and their non-exceedance probability is measures by fitting the most justifiable distributions (Yue and Rasmussen 2002;Salvadori 2004). The univariate return period (RP) is estimated using the CDFs of the best-fitted univariate distributions obtained via the normal KDE as The use of univariate exceedance probabilities and their return periods would be problematic and impractical in the design of hydraulic facilities. In the design of hydraulic structures, considering more than one variable is desirable. Compound flooding (CF) events are a multidimensional consequence of multiple extreme or non-extreme factors. This could demand multivariate return periods, such as joint RP and conditional joint RP. In hydrologic risk assessments, each approach for estimating return periods (joint and conditional) has its importance, which cannot be interchanged is and solely dependent upon the nature of the problem under consideration (Serinaldi 2015;Brunner et al. 2016). Their selection relies on the significance of the problem structure and its consequences of failure.
This study considers two significant classes of return periods: joint and conditional joint RP (Shiau 2003(Shiau , 2006Salvadori 2004;Salvadori and De Michele 2004;Zhang and Singh 2007). The joint RP can be described in two forms: OR-joint and AND-joint. In the ANDjoint case, two flood characteristics ( X ≥ x AND Y ≥ y for the bivariate case and X ≥ x AND Y ≥ y AND Z ≥ z for the trivariate case) simultaneously exceed a certain threshold, and their associated return period is thus estimated for the bivariate case (via 2-D copula) by T AND x,y = 1 P(X ≥ xANDY ≥ y) = 1 (1 − F(x) − F(y) + H(x, y) = 1 (1 − F(x) − F(y) + C(F(x), F(y))

Table 10
Comparing theoretical and empirical Kendall's coefficient derived from the best-fitted nested Frank copula in semiparametric approach (with nonparametric marginal density) and nested Gumbel copula in parametric approach (with parametric marginal density) Flood attribute pairs T AND X,Y,Z (x, y, z) = where F(x), F(y), and F(z) are the univariate flood marginals derived using best-fitted nonparametric kernel density estimation KDE; C(F(x), F(y)), C(F(y), F(z)), C(F(z), andF(x)) are the copula-based bivariate CDFs; and C(F(x), F(y), andF(z)) is the trivariate joint CDF derived from 3-D FNA copula under semiparametric settings (nonparametric marginals). Similarly, in the OR-joint case, either of the flood characteristics ( X ≥ x OR Y ≥ y ) for the bivariate case and (X ≥ x OR Y ≥ yORZ ≥ z) for the trivariate case exceed the given threshold for bivariate case (via 2-D copula) as and for trivariate case (via 3-D copula) as Using Eqs. (15) and (17), bivariate AND-joint and OR-joint RPs are estimated using the nonparametric marginal CDFs, and bivariate joint CDFs derived from the best-fitted parametric 2-D copulas for flood pairs rainfall-storm surge and rainfall-river discharge (refer to Table 11). Similarly, the trivariate AND-joint and OR-joint RPs are estimated from the trivariate joint CDFs obtained via FNA Frank copula constructed under semiparametric settings using Eqs. (16) and (18) (refer to Table 11). The results show that the joint return period in the "AND" case, T AND RAIN,STORM,RIVER is higher than the joint return period in the "OR" case T OR RAIN,STORM,RIVER . For example, at 100-year RP, flood events characterized by rainfall of 147.54 mm, storm surge of 0.65 m, and river discharge of 5427.17 m 3 /s, the AND-joint and OR-joint RP are T AND = 161.16 years and T OR = 34.16 years, respectively. Similarly, 200-year RP has flood characteristics of rainfall of 149.80 mm, storm surge of 0.67 m, and river discharge of 5484.14 m 3 /s; the AND-joint and OR-joint RP are T AND = 313.88 years and T OR = 67.50 years, and so on. It is concluded that the simultaneously occurrence of trivariate flood characteristics is more frequent in the "OR-joint" case and less frequent in the "AND-joint" case.

Estimating the Conditional Joint Return Period of Flood Events
This study employs the derived 3-D FNA copula joint framework (under semiparametric settings) to estimate conditional joint return periods. The conditional RP must rely on the conditional distribution relationship between the variables of interest (Salvadori andDe Michele 2004, 2010;Zhang and Singh 2007;Reddy and Ganguli 2013). Two cases of conditional joint distribution are considered, and their associated return periods are derived for the case of trivariate flood distribution.
Case 1: The conditional joint return period of flood variables X (rainfall) and Y (storm surge) given various percentile values of a third variable (Z ≤ z) (river discharge) in ORjoint case is estimated as where C(x, y, z) is the trivariate joint CDF derived from the selected FNA copula, and F(z) is the univariate marginal CDF derived from the best-fitted normal KDE. Figure 5a-c illustrate the joint return periods of the rainfall and storm surge observations given various percentiles (25, 50, 75, 90, 95, and 99) values of river discharge observations. It is clearly observed that the return periods decrease with an increase in the percentile values of the conditional variable (river discharge). For example, referring to Fig. 5a-c and Eq. (19), a flood event has the following flood characteristics: rainfall = 145.8 mm, storm surge = 0.52 m, the conditional return period is 55.18 years (when river discharge ≤ 1085 m 3 /s at 25 th percentile), 34.74 years (when river discharge ≤ 1615 m3/sec at 50 th percentile), and the return period is 23.47 years (when river discharge ≤ 2162.5 m3/sec at 75 th percentile), and so on. Case 2: The conditional joint return period of variable X (rainfall) given various percentile values of flood variables (Y ≤ y)(storm surge) and (Z ≤ z) (river discharge) is estimated as where C(F(y), F(z)) is the bivariate joint CDFs estimated from the best-fitted 2-D Gumbel copula (refer to Table 6b) with nonparametric marginals estimated via the normal KDE (for both storm surge and river discharge series) (refer to Table 2). From Fig. 6, it is concluded that a higher return period results from lower conditional storm surges and river discharge observations than the lower river discharge and storm surge values for the same specified values of rainfall characteristics. For example, referring to Fig. 6 and Eq. (20), flood events with flood characteristics of rainfall = 145.8 mm, conditional return period of 30 years (when storm surge ≤ 0.36 m and river discharge ≤ 2162.5 m3/sec), and return period of 26.97 years (when storm surge ≤ 0.44 m and river discharge ≤ 3100 m3/sec).In other words, the return period decreases when the percentile values of both conditional variables (storm surge and river discharge) increase. Moreover, the conditional return period of rainfall given various percentile values of storm surge with constant river discharge (river discharge = 5377 m 3 /s, 99 th percentile value) is also estimated using Eq. (20) and illustrated in Fig. 7. Similarly, Fig. 8 illustrates the joint return period of the rainfall conditional to river discharge series with a constant value of storm surge observations (storm surge = 0.62605 m, 99 th percentile value). From Fig. 7, it is observed that the return period is higher at lower storm surges, with a fixed value of river discharge observations than the lower storm surge for the same specified values of rainfall characteristics. Similarly, referring to Fig. 8, it is found that the RP is higher at a higher conditional river discharge with a constant storm surge value than the lower river discharge for the same specified values of rainfall events. (20)

Evaluating Multivariate Hydrologic Risk via Failure Probability
Using the traditional definition of joint return periods, it would be impossible to describe the risk of potential flood hazards during the entire project design lifetime, which has already been pointed out in the existing literature, such as by Salvadori et al. (2011), Huang and Chen (2015), and Xu et al. (2015). The traditional return period approach does not account for the planning horizon (Read and Vogel (2015). In such circumstances, the failure probability FP statistics can facilitate a more consistent and well-devised approach that quantifies the risk of hydrologic events. Previous studies, such as Xu et al. (2015), estimated FP statistics for assessing the hydrologic risk of bivariate flood events. Similarly, Moftakhari et al. (2017) introduced FP statistics in the bivariate flood risk assessments by compounding the joint impact of sea-level and fluvial flooding events. Our study estimates FP statistics for a trivariate joint distribution case to assess the risk of compound flooding events. Mathematically, if(X 1 , X 2 , … , X T ) , (Y 1 , Y 2 , … , Y T ), and (Z 1 , Z 2 , … , Z T ) are three different series of hydrologic events, T is the arbitrary project design lifetime and (S 1 , S 2 , … , S T ) represents a sequence of trivariate hazard scenarios. Then, the FP statistics can be expressed as and Using Eqs. (22) and (23), the FP statistics are estimated to demonstrate the hydrologic risk of trivariate flood characteristics (refer to Figs. 9a-e and 10a-b). The trivariate events produce higher failure probability values than the bivariate events for OR-joint cases (or univariate events). Therefore, neglecting the trivariate analysis by compounding the joint impact of the rainfall, storm surge, and river discharge observations results in underestimation of the failure probabilities. In addition, it is observed from the same Fig. 9a-e that the trivariate and bivariate (or univariate) hydrologic risk decreases with the increase in return periods, which further clearly points out 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. It is also concluded that trivariate hydrologic risk (also bivariate risk) values would increase with the increase in the service design lifetime of the hydraulic infrastructure. Figure 10a-b illustrates the variation in the occurrence probability or trivariate hydrologic risk with changes in rainfall events in differently designed storm surges and river discharge observations estimated under 100-year, 50-year, 20-year, and 10-year RPs for two different project design lifetimes (or service time of the hydraulic facilities) of 100 and 50 years. It is observed that the trivariate hydrologic risk increases with increased    in-service time of the hydraulic facilities; also, the hydrologic risk decreases with the increase in return periods. In other words, increasing the designated standard for the hydraulic facility would reduce the trivariate hydrologic risk. It is concluded from the above discussion that simultaneously accounting for all three flood characteristics, e.g., rainfall, storm surge, and river discharge, results in a better understanding of compound flooding phenomena and more critical information for flood risk assessments. The above analytical and graphical analyses are essential for the design and planning of sustainable flood management strategies in coastal regions.

Research Summary and Conclusions
The complex interplay between rainfall, storm surge, and river discharge, called compound flooding (CF) events, can exacerbate the effects of flooding events in coastal areas. A more comprehensive understanding of extreme hydrologic behaviour can be achieved by simultaneously integrating the targeted intercorrelated random variables instead of assessing their bivariate joint dependence structures. Only a few copula families exhibited higher dimensional mathematical structures in the joint modelling of the multidimensional hydrologic consequences. For instance, the symmetric 3-D Archimedean copulas would be ineffective in preserving all pairwise mutual concurrency at the lower stages. Hydrological phenomena can exhibit complex interdependency between targeted random pairs, rendering the application of symmetric 3-D copulas impractical in trivariate analysis. The FNA copula framework can be used to model each pair individually by incorporating multiple parametric joint asymmetric functions or joint integration of two or multiple bivariates or any dimensional Archimedean copula through another Archimedean structure. The FNA copula framework is often constructed under parametric distribution settings. Univariate marginals and their mutual dependence structures are modelled by employing parametric functions and their fitting procedure. Therefore, it would be quite complex to acknowledge the flexibility of copula in the parametric-based FNA joint distribution framework to a desirable degree without eliminating the restriction of subjective assumption of the function type for univariate marginal distribution. In this study, the multivariate FNA copula framework was constructed using the semiparametric fitting procedure in the trivariate joint modelling of rainfall, storm surge, and river discharge observations to assess compound flooding risk. The proposed trivariate distribution was applied to a case study on the west coast of Canada. The following are the main findings of the presented work: 1. The distribution-free normal KDE provides a much more accurate approximation of marginal distribution for rainfall, storm surge, and river discharge series than the parametric probability distributions. 2. BB7 copula is selected as best-fitted for rainfall-storm surge pair, Gumbel-Hougaard for storm surge-river discharge pair, and survival BB7 for rainfall-river discharge pair based on GOF test statistics. The adequacy of these selected copulas is examined via upper tail dependence coefficient (UTDC) assessments. 3. The trivariate joint analysis introduces five fully nested forms of Archimedean copulas with selected normal KDE marginals, called the semiparametric framework. The dependence parameters of the fitted FNA structure are calculated using the full maximum likelihood estimation procedure. Similarly, the above-selected FNA copulas are also incorporated under the parametric settings with parametric marginal distributions (GEV for rainfall series, normal for storm surge, and GEV for river discharge; refer to Table 4). Different analytical and graphical GOF test measures show that the FNA copula with nonparametric marginals PDF dominates other incorporated multivariate distribution frameworks (3-D asymmetric or FNA Gumbel copula with parametric marginal PDFs). The investigation reveals that using the selected nested Frank copula with KDE marginals, the correlation structure of the historical flood characteristics is satisfactorily regenerated compared with other multivariate copula joint approaches. In conclusion, by eliminating the restriction of distributional assumption in defining their univariate marginal behaviour, the nested Archimedean copulas under semiparametric settings exhibited improved performance compared with the traditional approach of constructing multivariate FNA copulas (under parametric settings). 4. In engineering-based hydraulic or flood defence infrastructure design, the accountability of both joint return periods (OR-joint and AND-joint cases) and conditional joint return periods is essential. The trivariate joint and conditional return periods are estimated (along with bivariate joint return periods) using the best-fitted derived semiparametric FNA Frank copula framework. On the other side, this study used the concept of FP statistics in the trivariate hydrologic risk assessments of rainfall, storm surge, and river discharge observations. It is concluded that neglecting the trivariate distribution analysis by compounding the mutual concurrency between rainfall, storm surge, and river discharge observations results in FP being underestimated because trivariate events produce higher FP than bivariate or univariate events. Moreover, it is concluded that both trivariate and bivariate (or univariate) hydrologic risk decrease with the increase in return periods and increase with the increase in the service lifetime of the hydraulic infrastructure. The variation in the trivariate hydrologic risk with changes in rainfall events in differently designed storm surges and river discharge observations was estimated under 100-year, 50-year, 20-year, and 10-year RPs for two different project design lifetimes (or service time of the hydraulic facilities) of 100 and 50 years. It is also observed that trivariate hydrologic risk increases with increased in-service time of the hydraulic facilities, and the hydrologic risk decreases with the increase in return periods.
The asymmetric copulas exhibit a few statistical challenges. For instance, the applicability of this framework can be justified when two correlation structures are near or equal to and smaller than the third correlation structure. It might also be quite ineffective in preserving all the lower-stage dependencies between multiple intercorrelated random observations, only limited to the positive range. No one could deny that justifiable conservation of all the lower-level dependencies is often challenging in copula-based multivariate joint distribution analysis, especially when complex dependencies are exhibited within multidimensional observations such as compound flooding events (Joe 1997;Kurowicka and Cooke 2006;Aas and Berg 2009). Besides this, the observed data length used in this study is 46 years (annual maximum sampling procedure). The level of uncertainty in the estimated quantiles may be affected by the length of data; results would be more consistent and reliable if a very long term was used (i.e., Haylock and Nicholls 2015;Bisht et al. 2017). Furthermore, this study used storm surge, rainfall, and streamflow data collected from one station because of limited data availability. A more comprehensive investigation would be possible if more stations were considered using the above-proposed methodology. In actuality, the assessments of coastal flood events using multiple station data will provide spatial and temporal variation, or in other words, nonuniformity due to variation in data that needs to be addressed. All such limitations will be addressed in our subsequent work for some other coastal regions in Canada.