Bivariate Hydrologic Risk Assessment of Flood Episodes using the Notation of Failure Probability

Floods are becoming the most severe and challenging hydrologic issue at the Kelantan River basin in Malaysia. Flood episodes are usually thoroughly characterized by flood peak discharge flow, volume and duration series. This study incorporated the copula-based methodology in deriving the joint distribution analysis of the annual flood characteristics and the failure probability for assessing the bivariate hydrologic risk. Both the Archimedean and Gaussian copula family were introduced and tested as possible candidate functions. The copula dependence parameters are estimated using the method-of-moment estimation procedure. The Gaussian copula was recognized as the best-fitted distribution for capturing the dependence structure of the flood peak-volume and peak-duration pairs based on goodness-of-fit test statistics and was further employed to derive the joint return periods. The bivariate hydrologic risks of flood peak flow and volume pair, and flood peak flow and duration pair in different return periods (i.e., 5, 10, 20, 50 and 100 years) were estimated and revealed that the risk statistics incrementally increase in the service lifetime and, at the same instant, incrementally decrease in return periods. In addition, we found that ignoring the mutual dependency can underestimate the failure probabilities where the univariate events produced a lower failure probability than the bivariate events. Similarly, the variations in bivariate hydrologic risk with the changes of flood peak in the different synthetic flood volume and duration series (i.e., 5, 10, 20, 50 and 100 years return periods) under different service lifetimes are demonstrated. Investigation revealed that the value of bivariate hydrologic risk statistics incrementally increases over the project lifetime (i.e., 30, 50, and 100 years) service time, and at the same time, it incrementally decreases in the return period of flood volume and duration. Overall, this study could provide a basis for making an appropriate flood defence plan and long-lasting infrastructure designs.


Introduction
Nowadays, flood events are characterized as one of the most severe and disastrous naturally occurring hydrologic consequences across the world, and the risk of their occurrence will increase in the future due to the global and regional climate-changing scenario [1][2][3]. In the operational planning, management or flood defence infrastructure design of water resources, it is often demanding to accurately estimate the flow exceedance probability or return periods for assessing the hydrologic risk. Flood frequency analysis or FFA is a statistical approach to establishing an interlink between the magnitude of flood episodes (or flood design quantiles) and their return periods (or nonexceedance probability) using the most logical and parsimonious probability distribution functions [4][5][6]. Flood is a multidimensional stochastic consequences, often defined thoroughly and comprehensively through peak discharge flow, volume and duration series [7][8][9]. Because of the multidimensional characteristics of a flood episode, the reliability of univariate return periods in the hydrologic risk assessments is questionable and thus point towards the feasibility and effectiveness of multivariate distribution modelling through integrating the multiple intercorrelated flood random vectors [10][11][12]. In actuality, the potential damage could likely be a function of several intercorrelated flood characteristics, such that ignoring the spatial dependency among them might contribute to underestimation of the uncertainty distributed over the estimated flood design quantiles, and thus, often demands more flood characteristics, i.e. based on the joint probability density function (or JPDF) and joint cumulative distribution functions (or JCDF) [13]. This is especially so from the prospect of hydraulic or flood defence infrastructure design procedures, where the accountability of multivariate design parameters could be a feasible desire [14,15]. In recent decades, the copula function has been identified as the most dynamic multivariate statistical framework or tool in comparison with traditional multivariate functions, due to its flexibility to select univariate marginal distributions and their joint dependence structure separately [16][17][18]. Numerous studies have performed the copula-based multivariate joint analysis of the flood characteristics and have estimated the design variable quantiles under different notations of return periods, i.e., based on joint distribution, conditional joint distribution or Kendall's distribution [14, 15, 19 and references therein].
The hydrologist or water experts are often interested in the evaluation of the mean inter-arrival time between the two successive occurrences of hydrologic episodes, also called return periods (or RPs). According to Brunner [20] and Latif and Firuza [21], the selection of an appropriate or justifiable RPs in solving the different structural or either nonstructural water-related queries is usually depending upon the importance of undertaken structure as well as its consequences of failure. On the otherside, the design quantiles define a higher RPs often seems a much practical approach in the hydraulic structure designs [22]. According to studies such as Salvadori et al. [23], Huang and Chen [24], Read and Vogel [25], and Moftakhari et al. [26], the traditional definition of the return periods (or RPs) would be ineffective and incapable of demonstrating the risk of potential flood hazards during the entire project design lifetime. Literature, such as Salvadori et al. [27] and Serinaldi [28] suggested the importance and necessity of the failure probability (FP) statistics, which can effectively quantify the risk of hydrologic episodes. Serinaldi [28] study pointed out that the definition of FP, which is usually defined over a given design life period, facilitates a more consistent and well-devised approach for assessing the risk of extreme flood episodes. Similarly, a study conducted by Read and Vogel [25] revealed that the definition and applicability of average return period (i.e. joint return period) did not account for the planning horizon and would be problematic. Besides this, a few other studies, such as Xu et al. [29], have performed the bivariate hydrologic risk of the multivariate flood characteristics for the Wei River basin using the concept of failure probability, which were derived using the bivariate copula distribution function. Moftakhari et al. [26] performed compound flood modelling caused by the coastal water level and fluvial flow (or river discharge) under different future sea-level scenarios, where the concept of failure probability was employed for assessing the risk. Recently, Xu et al. [30] performed copula-based joint modelling between rainfall and storm tide events, where the variation of bivariate flood risk was examined using the failure probability statistics.
Flooding has become one of the most critical hydrologic issues in the Kelantan River Basin, where the conditions are becoming more severe and intense year after year [31,32]. Examples of extreme flood scenarios due to intense and prolonged precipitation occurred in 2002, which affected a total population of about 714,287 people, or in December 2014, on the east coast of the Kelantan River basin, which affected more than 200,000 people. This river basin receives an annual rainfall of about 2500 mm, much of which occurs during the north-east monsoon (or wet season) between mid-October and mid-January. A study conducted by Hussian and Ismail [33] pointed out that the flood frequency risk is higher at the Gulliemard Bridge, Lebir and Galas gauge station. Similarly, Nashwan et al. [34] also revealed that the downstream area of the Kelantan River Basin has the highest risk of devastating flood events. Similarly, a study performed by Adnan and Atkinson [35] pointed out the existence of a significant trend in streamflow samples for both the upstream (i.e., Galas River) and downstream (i.e., Kelantan River) sub-catchments such that in the downstream area streamflow increased in the north-east monsoon circulation. Also, the expectation of the occurrence of catastrophic flooding has increased from once in every 50 to 15 years from 2004 in the Kelantan region in Malaysia [36].
The cause of frequent failure of hydrologic or flood defence infrastructure in Malaysia due to the impact of moderately severe of flood episodes might be attributed due to the lack of complete flood hydrograph, where only flood peak discharge samples often targeted in deriving flood frequency curve during the structural development. The incorporation of the copula-based methodology is still very rare in Malaysia. Recently, few copula-based multivariate analysis and modelling of flood episodes are performed for this river basin such as Latif and Firuza [37], incorporated the vine or pair-copula construction (or PPC) in the trivariate distribution modelling of the flood characteristics. In this demonstration, the D-vine tree structure was selected in establishing the trivariate dependency analysis and deriving both the primary OR and AND-joint return periods. Similarly, Latif and Firuza [38] introduced the copula-based bivariate joint analysis of the annual flood characteristics under the semiparametric framework where a distinguish varieties of nonparametric kernel density functions were introduced and tested in the modelling of the univariate marginal distribution of the flood characteristics. In another demonstration, Latif and Firuza [39] performed and tested the efficacy of the Archimedean copulas and Gaussian copula in the trivariate distribution analysis of flood characteristics for the same river basin. On the other side, Latif and Firuza [40], also performed the bivariate flood analysis using the nonparametric copula distribution framework. In such incorporation, a combination of both parametric and nonparametric marginal distribution separately conjoined by a nonparametric copula framework, which was based on the Beta kernel function. The Beta kernel copula function was incorporated to estimate bivariate copula density which further employed in deriving the joint cumulative density of flood peak-volume, volume-duration and peak-duration pairs and their associated joint as well as conditional return periods.
In this study, we incorporated and tested the Archimedean and Gaussian copulas in deriving the bivariate distribution analysis of the annual series between flood peak and volume (P-V) and flood peak flow and duration (P-D) pairs. Then the best-fitted parametric copula for each flood pair were employed in deriving the joint probability distribution and their associated return periods, i.e., 'OR' and 'AND' joint return periods, which were further employed in the estimation of failure probability or FP statistics. The FP statistics were estimated to assess the variations in the level of bivariate hydrologic risk for the flood episodes during the entire project design lifetime. In this study, the bivariate hydrologic risk analysis was carried out through two different investigation procedures: (1) demonstrating the hydrologic risk of the flood peak flow and volume pair (P-V) and flood peak flow and duration pair (P-D) in different return periods; and (2) visualizing the variation in hydrologic risk with changes of flood peak flow in different designed or synthetic flood volumes and durations under different service lifetimes. The copula-based bivariate distribution framework is applied as a case study for the daily basis streamflow discharge records, which were collected from the period 1961-2016 for the Kelantan River basin at the Gulliemard Bridge gauge station and delivered by the Department of Irrigation and Drainage (DID), Malaysia. This investigation will be useful for the hydrologist and water engineer in the prevention and management of flood episodes. Figure 1 illustrates the methodological flow diagram indicating the steps of the analysis in the present demonstration. The structure of this manuscript is organized as follows. Details of the study area and method for retrieving or delineating the trivariate flood characteristics are discussed in the next section (i.e. Section 2). A brief theoretical discussion about the copulabased bivariate joint distribution analysis of the flood characteristics also, briefing of the bivariate hydrologic risk using the concept of the failure probability (or FP) statistics are discussed in the third section (i.e. Section 3). In actuality, the FP statistics are derived for both the primary, OR-joint and AND-joint return period using the best-fitted bivariate copulas for the flood attribute pairs. The fourth section (i.e. Section 4) provides the result of the findings and discussions. The fifth section (Section 5) provides the research conclusions.

Details of Study Area and Delineation of Flood Characteristics
The origin of the Kelantan River basin is the Tahan mountain range, which finally drains into the South China Sea, located in the north-eastern part of Peninsular Malaysia and occupying more than 85% of the state of Kelantan. The length of this river basin is about 248 km long and the total drain area is about 13,100 km 2 . The geographical extent of this river basin is Lat 4° 30′ N to 6° 15′ N and Long 101°E to 102° 45′ E. Annual precipitation varies between 0 mm (i.e. dry period) and 1750 mm (i.e. wet or north-eastern monsoonal period) and estimated runoff of this region is about 500 m 3 sec −1 [41]. In this study, the daily streamflow discharge observations were collected at the Gulliemard Bridge gauge station, which is located near the Kuala Kari region (i.e. downstream) of the Kelantan River and provided by the DID (see Figure 2). Agriculture activities such as paddy, rubber and palm oil are the major land use in the midstream region, and forest in the upstream region of the river basin. Due to this rapid human intervention from natural to landuse activities would be responsible for these extreme hydrologic consequences [35,36].

Figure 2. Geographical location of the study area (Gulliemard Bridge gauge station located at the downstream region of Kelantan River Basin in Malaysia)
The annual (maximum) series (AM), also called the block (annual) maxima, and peak over threshold (or POT) are the two holistic approaches of a partial data series that are widely accepted in extreme event modelling [42,43]. In actuality, the partial data series only focuses on the extreme hydrograph portion, i.e. either high flow (i.e. for flood episodes) or low flow (i.e. for drought events), instead of visualizing the entire hydrograph. In this demonstration, we prefer the AM approach in the sampling of triplet flood characteristics, i.e., the flood peak discharge flow, volume and duration series, which are retrieved from the daily streamflow discharge records at the Gulliemard Bridge gauge station between the years 1960 and 2016, and delivered by the DID (Drainage and Irrigation Department) Malaysia. The characterization of flood peak flow values is based on maximum streamflow discharge values at an annual scale, which means for each year there will be only one flood episode at the targeted site (i.e. [5,10]), using the Equation 1. Similarly, corresponding to each flood peak value of the annual basis flood events, flood volume and duration series are derived from the streamflow hydrograph using the algorithm reported by studies such as Yue et al. [10], Zhang and Singh [8], and Reddy and Ganguli [3,15], using Equations 2 and 3 as referring to Figure 3 (graphical illustration).

Bivariate Joint Distribution via Copula Function
The idea of the copula function was developed by Saklar [16]. It models the multivariate distributions to their univariate marginal distribution functions, which are not necessary from the same family of probability functions, capture a broader extent of dependencies (i.e., both linear and nonlinear), and preserve their joint dependence structure [6,18,19]. Mathematically, if (M, N) are the bivariate random vectors, such that u = F M (m) = P(M ≤ m) and v = F N (n) = P(N ≤ n) is the univariate marginal distribution of individual characteristics, then it can be demonstrated through a function called copula or C, which describes the dependence structure in between them and can be defined on the unit square. It can be expressed as: Where c is the density function of bivariate copula C. It can be defined as: For the extended mathematical details and theorems about copula functions, readers are advised to refer to Saklar [16], Nelsen [18] and Genest and Favre [44], as well as the 'International Association of Hydrological Sciences (or IAHS)' for extended details and lists of their applicability in the field of hydro-climatological observations.
In this demonstration, the two-dimensional (2-D) copulas, such as the Archimedean family-based Clayton, Gumbel and Frank copula and also one Elliptical family called the Gaussian copula, are introduced and tested as candidate functions, where the best-fitted copulas are employed to model the bivariate joint dependency and estimation of joint return periods of the flood pairs (P-V) and (P-D). In actuality, different copula functions have different extents of dependency capturing capabilities. For example, both the Gumbel and Clayton copulas are only capable of modelling positively correlated random variables (i.e., 0 ≤ < ∞ & Kendall's tau τθ ≥ 0 ), while the Frank family copula accommodates the entire range of dependencies among hydrologic characteristics (i.e., τ θ ∈ [1, −1]) (i.e. [18,45]). Mathematically, we can approximate the bivariate Archimedean class copula based on the two-dimensional copula framework (i.e. [C: [0,1] 2 ⟶ [0,1]]) as given below; Where ∅ (.) and ∅ −1 represent the generator function of the specified Archimedean copulas and their inverse, such that the generator (φ: I ⟶ R+) signifies the positive, convex and decreasing function and can be approximated for (∅(1) = 0 & ∅(1) = ∞) [18]. The mathematical expression for bivariate Archimedean copulas and their associated properties (i.e., parameter range (θ), generating function (or generator) ϕ(t), and relationship of Kendall's τ and θ (τ θ )) are listed in Table 2. On other side, the Gaussian copula is an implicit copula that can easily capture both the positive as well as the negative dependency, and has almost no dependence on the tails of distribution, which are mostly distributed around the centre and can be mathematically expressed as [46,47].
where ( ) is the Debye function, for any positive integer k,

Estimating Copula Dependence Parameters and their Fitness Investigation
Different approaches are often motivated in estimating the copula dependence parameters such as the method of moment (MOM) [48], exact maximum likelihood or EML [46], inference functions for marginal (IFM) [45], canonical maximum likelihood (CML) [49], Maximum pseudo likelihood estimations (MPL) [3,49]. In this demonstration, the copula dependence parameters are derived using the method-of-moments (or MOM) estimators based on the inversion of Kendall's tau ( ) and are highly flexible and very popular, especially among the Archimedean copula families [12,49]. The MOM estimation approach is solely based on establishing an interlink between the unknown dependence parameter 'θ′ and the sample ranked correlation coefficient, i.e., Kendall's tau ( ), and if a correlation is exhibited between them, then the copula parameter can be estimated as [3,12,18]; Where ϕ(•) and ϕ ′ (t) indicate the generator function and their first derivatives for the Archimedean copula family function. In addition, Kendall's tau ( ) can be statistically defined as [50]; The mathematical range of Kendall's tau ( ) is [−1,1], such that = 1 represents concordance, = 1 represents discordance, and = 0 represents no concordance exhibited between random pairs. The adequacy of hypothesized copulas fitted to bivariate flood characteristics are identified using the Cramer-von Mises statistics, which are based on the parametric bootstrapping procedure that makes use of the Cramer-von Mises (or CvM) statistic, Sn (i.e. [3,51,52]). The simulation of copula function (i.e. estimating copula dependence parameters and their fitness test) was carried out using the R statistics platform [53]. The CvM test statistics can be mathematically expressed as; Where C n and C θ are the empirical copula (estimated using n observational flood characteristics) and parametric copula (which is derived under the null hypothesis). Under this test, the p-value is estimated using the parametric bootstrapping procedure (i.e. [51]). Overall, a minimum value of the Sn test statistic must indicate a minimum gap and deviation between an empirical and derived parametric copula.

Bivariate Risk Evaluation using the Notation of Failure Probability (FP)
In the multivariate risk framework, the return periods can be derived from the exceedance probabilities of the flood attributes pair [3,54]. In actuality, hydrology and hydraulic applications are mostly interested in the evaluation of the mean inter-arrival period between two design events, which are usually defined in a year called the return period [14,55]. The OR-joint return period of the annual flood episodes, where either of the flood vectors in a bivariate combination exceeds a specific threshold value, can be estimated [7,14] as: Where H(p, v) and H(p, d) are the joint cumulative densities estimated using the best-fitted bivariate copula C(F(p), F(v)) and C(F(p), F(d)) with marginal distribution F(p), F(v), and F (d) series. On the other side, the joint return period, where both of the flood characteristics simultaneously exceed a certain threshold value, also called primary 'AND' joint return periods, can be expressed as: , for flood pair (P − D) The notation of failure probability or FP can be practical in assessing the variation of bivariate flood risk during the entire project lifetime. The FP statistics usually define the chance of a potential flood event occurring at least once in a given project design lifetime (i.e. [26][27]30]). Mathematically, the failure probability can be expressed as: Where is the series of bivariate flood hazard scenario. In addition, the failure probability for the independent and identical distribution series can be estimated using Equation 11 [27]; However, the failure probability for both the primary 'OR' and 'AND' joint return periods can be estimated as; p T AND = 1 − (1 − P (M ≥ m AND N ≥ n)) T , for AND − joint case (20)

Bivariate Dependency Modelling using Two-dimensional Copula
Approximating the most parsimonious probability density functions or PDFs for defining the univariate marginal distribution of the individual flood characteristics is often a mandatory pre-requisite before introducing the random observation into the multivariate or copula distribution framework. In this demonstration, few commonly used onedimensional parametric family functions, such as Gamma-3P, Generalized Extreme Value (GEV), Generalized Gamma-3P, Inverse Gaussian-2P, Johnson SB-4P, Lognormal-2P and Weibull-2P, are selected and tested as candidate models, and all the distribution fitting procedures are carried out using the Easyfit software (MathWave Technologies 2004-2017). In reality, no universally accepted distributions are assigned or in favour of any probability distribution functions from any literatures to model any extreme hydrologic observations (i.e. [10,56,57]). The Gringorten-based position-plotting formula is used in estimating the empirical probabilities [58]. Different goodness-of-fit test statistics, such as the Kolmogorov-Smirnov (or K-S) test and the Anderson-Darling (or A-D) test [59], based on information criteria statistics such as Akaike information criteria (or AIC) [60], Schwartz's Bayesian information criteria (or BIC) [61] and Hannan-Quinn Information criteria (HQC) [62] are used in visualizing the best-fitted marginal distribution to each individual flood characteristic. Therefore, the minimum test statistics value must indicate the best-fitted model. Table 3 illustrates the fitness test statistics of the candidate function and it is shows that flood peak flow is best modelled by the Lognormal-2P distribution, volume series by the Johnson SB-4P distribution and duration series by the Gamma-3P distribution. The strength of dependency among the triplet flood characteristics were investigated analytically using the Pearson's linear correlation (r), Kendall's tau (t) and Spearman's rho (ρ), and their estimated statistics are shown in Table 4. Based on the estimated statistics, it is clearly inferred that there is a strong positive correlation between flood pair (P-V), but a weak and negative correlation between flood pair (P-D). Three Archimedean copula families, the Clayton, Gumbel and Frank copula, and one Elliptical copula, the Gaussian, were selected. The copulas dependence parameters were derived using the method-of-moment (MOM) estimators, which are based on the inversion of Kendall's tau ( ) using Equation 9, 10 and 11, and their estimated statistics are listed in Table 5. The fitness level and adequacy of the candidate copulas for both flood pair (P-V) and (P-D) were investigated using the CvM distance statistics with the parametric bootstrap procedure, where test statistic 'Sn' and its associated p-values were estimated from 1000 and 500 simulated random samples by the mean of the parametric bootstrap procedure ( Table 5). The analytical fitness investigation reveals that the Gaussian copula exhibited the minimum 'Sn' test statistics value and also the highest p-value for both the (P-V) and (P-D) pairs, and thus was identified as the most justifiable copula for describing their joint dependence structure. A graphical visual inspection was also carried out using the scatter plot (refer to Figure 4) of the 1000 flood pairs, for both (P-V) and (P-D) observations, which were simulated from the joint distribution using the Gaussian copula function with best-fitted marginal distributions. From the comparison scatter plots of the simulated data (blue colour) with overlapped observed flood characteristics (red colour), it can be revealed that the Gaussian copula performed satisfactorily for capturing observed dependence for both flood pair (P-V) and (P-D), which means the simulated data adequately overlapped with the natural dependence of observed flood characteristics. Figure 5 and 6 illustrate the joint probability density and cumulative distribution plot of the bivariate Gaussian copula fitted to flood pair (P-V) and (P-D).  NA denotes that the Clayton and Gumbel-Houggard cannot be used for negative dependence modelling between the random variables, i.e., for flood pair (P-D), which is only applicable for positively correlated flood characteristics, i.e., for flood pair (P-V).

Bivariate Hydrologic Risk Estimation
The selection of an appropriate or justifiable return period solely depends upon the importance of undertaken structure as well as its consequences of failure [20]. The best-fitted Gaussian copulas were employed in deriving primary 'OR' and 'AND' joint return periods or RPs for flood pairs (P-V) and (P-D) using Equations 13 to 16 for the given flood characteristics (see Table 6). It can be noticed that for the same values of the peak flow and volume pair (P-V) or the flood peak flow and duration pair (P-D), the joint return period in the 'AND' case is much greater than 'OR' joint cases, i.e., TOR (OR-joint case) < TAND (AND-joint case. = 563.126478 years and T OR PD = 9.25889919 years for P-D pairs. In addition, from the same results it was also concluded that the univariate RP is greater than the OR-joint case but smaller than the AND-joint case. We can further conclude that the univariate approach in the hydrologic risk would be inappropriate and ineffective. The failure probability statistics (see section 3.2) were estimated to demonstrate the bivariate hydrologic risk of flood characteristics for flood pair (P-V) and (P-D), which usually measure the chance of flood episodes potentially occurring once in a given project design lifetime. Figures 7 and 8 illustrate the risk statistics of the flood episodes for different return periods (i.e., 100, 50, 20, 10, 5 years) where the failure probabilities are estimated based on the univariate approach (indicated by blue and green colour), and bivariate events based on the OR-joint scenario for the flood pair (P-V) (Figure 7 a, b, c, d, and e) and flood pair (P-D) (Figure 8 a, b, c, d, and e). The results revealed that the risk statistics value would increase by an increment in the service lifetime. For example, for the 100-year univariate return period events, the failure probability attained the value, FP = 0.179654 when the service lifetime is 10 years, which is lower in comparison with 50-year RP events, i.e., FP = 0.3086, at the same service lifetime. Similarly, for the 5-year RP events, the FP is 0.95831, which is higher than the FP statistics estimated for 10-and 20year RP events at the same service lifetime. On the other side, the univariate events produce lower FP values than the bivariate events (see both Figure 7 and 8). Thus, it could be revealed that neglecting the joint impact or mutual dependencies between multiple flood characteristics would be problematic and underestimate the FPs in the assessments of flood risk. From the same figure, it is also inferred that with the increase in the value of return period statistics, both the bivariate and univariate hydrologic risk would be decreased. Similarly, Figures 9 and 10 illustrates bivariate hydrologic risk statistics for the flood peak flow and volume pair (P-V) and the flood peak flow and duration pair (P-D) under the different scenarios where we considered the 5-year, 10-year, 20-year, 50-year, and 100-year RP of the flood volume and duration series. In this analysis, variation of the bivariate hydrologic risk with the change of flood peak flow in different synthetic or designed flood durations was studied for different project lifetimes (i.e. 30, 50, and 100 years) service time. The investigation revealed that the value of the bivariate hydrologic risk statistics is increasing with the increment in the project lifetime and at the same time, it is decreasing with the increment in the RP of flood volume and duration.

Conclusion
Water resources operational planning, managements or either flood defence infrastructure designs often demand the accurate estimations of flow exceedance probability for visualizing the risk of flood episodes. The Kelantan River basin in Malaysia is often affected by the most intensive and severe monsoonal flooding episodes. In this study, we incorporated the copula function in deriving the joint distribution of annual flood characteristics and failure probabilities for assessing the bivariate hydrologic risk. The Gaussian copula is recognized as the most parsimonious function for the flood pairs (P-V) and (P-D). The selected copula was then employed in deriving the joint return periods (RPs) for both the 'OR' and 'AND' case for a different possible combination of flood characteristics. It is revealed that the univariate RPs are higher than the 'OR' joint case but less than the 'AND' joint case, which further reveals that the RPs derived from univariate flood characteristics do not facilitate a sufficient and comprehensive approach to flood risk assessments. The bivariate hydrologic risk of flood episodes based on the joint relationship between flood pairs (P-V) and (P-D) was estimated using the definition of the failure probability (FP), which could be (b) (c) useful in determining the risk level of extreme floods during the entire project design lifetime. The hydrologic risk analysis was carried out through two different investigation procedures. From the first investigation, it is inferred that the bivariate hydrologic risk statistics would increase by an increment in the service lifetime for the hydraulic facilities and at the same instant, both the bivariate and univariate hydrologic risk would incrementally decrease in return periods value. However, it could be revealed that the definition of return period is not capable of recognizing the probability or chance of the occurrence of extreme hydrologic events during the project lifetime.
In addition, univariate events produce lower FP statistics values in comparison with the bivariate events. Thus, it could be revealed that neglecting the mutual concurrency between multiple flood characteristics would be problematic and might underestimate the FP in the assessments of flood risk. Similarly, the variation of the bivariate hydrologic risk with the change of flood peak flow in different synthetic or designed flood volumes and flood durations for different service times are also demonstrated. It is concluded that hydrologic risk statistics incrementally increase by the project lifetime and at the same instant, incrementally decrease in the RP of the designed flood volume and duration series. Overall, this study provides a basis for making an appropriate flood defence plan and long-lasting flood defence infrastructure designs. According to Zhang [6], Bender et al. [63], Sarhadi et al. [64] the validity or accuracy in the estimated flood design quantiles under the time-invariant or stationary risk framework might be questionable, due to the ignorance of the accountability of the changing environment (Climate change or LULC change) in the flood distribution analysis. The existence of the non-stationarity usually tries to interrupt the hydrological behaviour within the catchment region thus might alter the expectation of such extreme happening under the stationary hydrologic risk assessment [65]. In other words, the actual associated flood risk either greater or smaller than the hazards statistics accounted under stationary risk concept and might reveal for under-dimensioned or overdimensioned in the designing of hydraulic structure [66]. All the above raised issues can be considered as a future research purpose for the same river basin in order to achieve much practical and justifiable flood hazard assessments.