Parametric Vine Copula Construction for Flood Analysis for Kelantan River Basin in Malaysia

The multivariate approach of flood characteristics such as flood peak flow (P), volume (V), and duration (D) is much beneficial in recognizing the critical behaviour of flood episodes at a river basin scale. The incorporation of 2-dimensional copulas for establishing bivariate flood dependency frequently appears, but it could be more comprehensive if we focus all the three flood characteristic simultaneously. In such circumstances, incorporation of vine or Pair-Copula Construction (PCC) could produce a better approximation of joint probability density and much practical approach in the uncertainty analysis, in comparison with conventional trivariate copula distribution. This study demonstrated the efficacy of parametric vine copula in the modelling of trivariate flood characteristics for the Kelantan River basin in Malaysia. The D-vine tree structure is selected where the Gaussian and Frank copula is recognized for bivariate flood pairs (P-V) and (P-D) pairs in the first stage, using the maximum-pseudo-likelihood (MPL) estimation procedure. Similarly, the Gumbel copula is selected in the modelling of conditioned data obtained through the conditional distribution function of bivariate copulas selected in the previous stage based on the partial differentiation, also called h-function. Finally, the full density function of the 3-dimension structure is derived and compared with the observed flood characteristics. Furthermore, tail dependence properties and behaviour of D-vine copula are also investigated, which reveals for well capturing the general behaviour of Gaussian and Frank copula fitted to flood pair (P-V) and (V-D) and reproduces the overall flood correlation structure fairely well. Both the primary ‘OR’ and ‘AND’ joint return periods for trivariate flood events are estimated which pointing that ‘AND’ joint case produces higher return value than ‘OR’ case.


Introduction
The designing and planning of flood defence infrastructure system often demand the accurate estimation of the specified exceedance probabilities or design quantiles of the flow characteristics for assessing the hydrologic risk. Flood is a multivariate probabilistic consequence, usually exhibited a complex dependency shown among its triplet intercorrelated vectors, flood peak flow, volume and its durations series [1][2][3]. The frequency analysis of flood episodes is the statistical approach of establishing the interrelationship between the magnitude of flood events (or flood event quantiles) and their frequency of occurrence using the most logical and parsimonious probability distribution functions [1,4,5]. In the earlier studies, the univariate flood frequency approaches were frequently incorporated for estimating the univariate return periods or assessing flood risk [4]. However, on later, numerous studies recognized the drawback or unreliability of the single variable flood modelling approach which would be In the development of the higher dimension copula distribution framework, literature often targeted an interactive set of conventional trivariate copulas structure, but all such attempts usually surrounded with few statistical constraints and issues [19,22,23]. Few incorporations based on some standard trivariate copulas framework such as Grimaldi and Serinaldi [24] introduced the mono-parametric and asymmetric or fully nested Archimedean (or FNA) based Frank distribution function along with the Gumbel logistic distributions for flood characteristics. Similarly, Serinaldi and Grimaldi [25], also introduced and pointed out the significance and flexibility of FNA copula framework for the trivariate flood events. On another side, Genest et al. [26] pointed out the effectiveness of meta-elliptical copulas in the treatment of annual spring flood characteristics. They thus revealed that it could be flexible and effective in the handling of multidimensional observation because it exhibited a correlation matrix and which significantly preserved the pairwise dependencies exhibited among the random variables. This study also pointed out that the applicability of meta-elliptical copulas might be ineffective under the low probabilities unless the asymptotic properties of data will be justified through the strong arguments. Similarly, Reddy and Ganguli [18] study demonstrated the efficacy of the fully nested Archimedean or FNA class copula along with the Student's t copula function for the annual flood characteristics. The same study also pointed out the importance of considering the trivariate return period (derived from both joint and conditional distribution) in comparison with bivariate or univariate return periods in the hydrologic risk assessments.
Fan and Zheng [21] pointed the adequacy of the Entropy copula in conjunction with the Gibbs sampling procedure in the treatment of flood characteristics and concluded that this copulas function could be easily projected into any higher dimension distribution space. Besides the above incorporation, Kao and Govindaraju [27] demonstrated the modelling efficacy of Plackett family of copulas for the rainfall characteristics, Madadgar and Moradkhani [20] modelled the drought characteristics using the trivariate Gumble and t copula also, Song and Singh [28] performed a frequency analysis of drought events using meta-elliptical copulas.
According to Daneshkhan et al. [19], it might be quite ineffective to reflect and recognize all the possible joint concurrency exhibited among multivariate characteristics through conventional (or standard) trivariate copula distribution approach. Also, resolving the mutual dependencies of multivariate extreme via the conventional copula framework is quite complicated due to the higher degree of uncertainty as well as complex dependence structure of the hydrologic characteristics thus it often demands the estimation of tail dependence coefficient for the various tail dependency in a precise manner through adopting a comprehensive multivariate methodology [19,23].
In such circumstances, the incorporation of the vine or pair-copula construction (or PCC) in the modelling of the multidimensional complex hydrologic phenomenon could be the right choice and an effective way of uncertainty analysis in comparison with standard trivariate copulas [22,23,29]. The applicability of vine or pair-copula framework seems much popular in the field of financial data analysis and their risk management [30]. However, in the recent decades, such incorporation is also recognized and accepted in the treatment of hydrologic observations such as in the modelling of drought episodes [28,31], flood episodes [3,19], rainfall or storm distribution modelling [32,33]. The implementation of vine or pair-copula framework is entirely based on the fact that it breaks or segregates the full multivariate density through the nesting of pair-copulae, which is based on the conditional mixing procedure [29]. In other words, the multiple best-fitted bivariate (or 2-dimension) copulas are often employed for each random pair and their conditional and unconditional distribution to form a nesting framework instead of selecting a fixed copula for all the intercorrelated flood vectors and thus often reveals for higher flexibility in the construction of higher dimension hydrologic framework. Under the regular vine framework, the canonical vine, also called C-vine and D-vine are two holistic decomposition approach of the pair-copula approximation, [30,34].
The D-vine structure is highly flexible in comparison with C-vine, especially in the presence of a known variable that controls the level of mutual dependency within the samples [19,23]. The strength of dependency or degree of mutual concurrencies among the multivariate characteristics further decided in the selection of either the D-vine or Cvine structure [3,34]. According to Graler et al. [3] and Hobaek et al. [35], the manner through which the decomposition of vine copula will be performed, which will further decide their approximation capability with the multivariate observations and which further reveals that the conditioning process of the vine or PCC model is not rigid or fixed. Readers are advised to follow Czado et al. [34] for the extended details over the selection procedure of the regular vine framework also, Kurowicka and Cooke [36] and Aas and Berg [37] for the details study over the C-and D-vine copula structure.
Flood becoming a crucial and most severe natural disaster for the Kelantan River basin, especially during the wet monsoon circulation in Malaysia and perceiving for increasing in term of their frequency and magnitude [38]. Such that, in the year 2002, intense and prolonged precipitation attributed for flooding and affected 714, 287 people, similarly, in the year 2014, much heavy precipitation caused flooding in most of the part of the eastern coast of Kelantan river basin and which affected more than 200,000 people [39].
In actuality, the downstream area of the Kelantan River basin is at the highest risk of devastating flood events according to the study performed by Nashwan et al. [40]. Multivariate analysis of flood episodes could be a practical approach and much beneficial for hydrologist and water practioner in recognizing the critical behaviour of flood episodes at a river basin scale and for solving different water-related issues. In this literature, we incorporated and tested the vine or pair-copula framework to analyze the block (annual) maxima-based flood characteristics for the Kelantan river basin in Malaysia. The 3-dimensional joint probability framework is simulated using the parametric vine or pair-copula framework approach in the trivariate distribution modelling of flood characterization and estimating the multivariate (i.e., trivariate and bivariate) primary joint return periods.
The structure of this paper as follow, as a prerequisite, the brief details on the theoretical aspect in constructing the 3-dimension vine or pair-copula distribution framework is discussed in the next section (or Section-2). In the third section, the vine copula-based methodology is incorporated to a case study for the annual flood characteristics in the Kelantan River basin which further comprises, the selection of appropriate vine tree structure (i.e., D-or C-vine) and then simulation of vine structure by fitting most appropriate 2-dimensional copulas in a stage-wise hierarchical nesting manner and also, estimating the primary joint return periods for the trivariate (and bivariate) flood episodes. The fourth section provides the research conclusions. Figure 1 illustrates the methodological flow diagram, indicating the steps of analysis.

Construction of 3-Dimension Vine or Pair-Copula
According to Nelsen [8], copula function interlinks the multivariate distribution to their univariate probability functions, also called the marginal distribution, and their idea was firstly given by Saklar [7]. Multivariate copula function modelled the univariate marginal density of individual intercorrelated random variables independently from their joint dependence structure. According to Sklar's theorem [7], if (X, Y) be bivariate random variables with continuous marginal distributions such as, u 1 = F X (x), and u 2 = F Y (y), then it can be characterized uniquely by using the Copula function or C, and defined on the unit square; Where, C = bivariate copula under investigation; F X (x) = F Y (y) = cumulative distribution function (or CDF) of univariate observations i.e., X and Y; H X,Y (x, y) = bivariate joint probability density function (or PDF), which can be expressed in terms of its univariate marginal functions and the associated dependence function C, as revealed from Equation 1. Also, if F X (x) and F Y (y) are continuous which indicates the copula C must be unique [11,41].
Conversely, if F X (x) , F Y (y) and the copula C[x, y] is given, then Equation (1) must defines the bivariate joint distribution functions with its marginal distributions F X (x) and F Y (y). Similarly, if f X (x) and f Y (y) are the PDF of variable X and Y, then the joint probability density can be defined through Equation (2), as given below; Where, c is the density function of bivariate copula C and which can be further expressed as The construction of vine or PCC framework is solely based on the principle of stagewise hierarchical nesting procedure, also called the conditional mixing of the bivariate (or 2-dimension) copulas where the decomposition of full multivariate density into a simple local building block are carried out through the conditional independence or pair-copula [3,23,29]. Thus, such conditional mixing approaches attributes for higher flexibility in the modelling of trivariate flood characteristics. At the beginning of vine construction, the marginal empirical distribution functions are obtained by transforming the samples of all the three flood characteristics, P, V, D, within the range [0, 1] using the rank order transformation S (see Equation 4). Figure 2 illustrates the 3-dimensional parametric vine framework. In the first tree (or Tree-1) structure, all the relevant pairwise joint association between the flood characteristics are modelled using the best-fitted bivariate copula function, i.e., between (V-P) pair and (P-D) pair (where we placed the flood variable, P between the other two variables, V and D, and which is solely based on the strength of dependencies exhibited among the flood characteristics, see section 3.2.1 for further details). In this demonstration, we introduced the different varieties of bivariate copula families in the establishing the mutual dependency between flood characteristics such as the Archimedean copulas (i.e., Clayton, Frank, Gumbel, Joe) [10,11,15] and one Elliptical copula [i.e., Gaussian or normal copula] [42,43]. In actuality, the Clayton and Gumbel family only capture the positive dependencies, but the Frank copula significantly captures the entire range of dependencies (i.e., τ θ ∈ [1, −1]).
According to Aas et al. [23] and Graler et al. [3], the bivariate copulas fitted to flood pairs (V-P) and (P-D) in the first tree (or Tree-1) can be conditioned under the flood characteristics, P, using the partial differentiation which is also called the h-function (i.e., h(v, p, θ) and h(d, p, θ)) to obtain the conditional cumulative distribution functions (CCDFs), F V|P and F D|P using the Equations 5 and 6 as given below; The copulas dependence parameter in the first tree of vine structure (see, Figure 2) i.e., (θ vp , θ pd )) are estimated from the original flood characteristics P, V and D using the maximum-pseudo-likelihood (or MPL) estimation procedure, which is based on maximizing the pseudo log-likelihood function [15,44,45]. After that, the transformed data (obtained using Equation 5 and 6), are employed to determine the unknown copula dependence parameters in the second tree (or Tree-2) and thus used to fit another bivariate copula C VD|P . The best-fitted bivariate copulas in both Tree-1 and Tree-2 are selected using the Cramer-von Mises (or CvM) test statistics which is further based on the parametric bootstrapping procedure that uses the CvM test statistics, [46,47]. Finally, the full density function c PVD (p, v, d) of the 3-dimensional copula can be derived using the Equation 7 [19,23];

Description of Study Area and Flood Characteristics
This experiment targeted the daily basis streamflow discharge records from the period (1961-2016) of the Kelantan river basin, which are collected at the Gulliemard Bridge gauge station (located at the downstream region) [see Figure  3] and delivered by the Drainage and Irrigation Department, Malaysia. The geographical location of this river basin is Lat 4° 30 ′ N to 6° 15 ′ N and Long 101°E to 102° 45 ′ E. This river is originating from Tahan mountain range to the South China Sea (a north-eastern portion of Peninsular Malaysia). The drain area of this river is about 13100 km 2 . According to DID report [38], the estimated runoff of this river basin is about 500 m 3 sec −1 where the variations of annual precipitations are in between 0 mm (dry period)-1750mm (wet or north-eastern monsoonal period). According to Hussain and Ismail [48], the frequency of flood episodes has highest in the Gulliemard Bridge, Lebir and Galas gauge stations in compare with the Nenggiri gauge station. Nashwan et al., [40] literature also revealed that the risk of devasting flood episodes in the downstream portion of this river basin is higher.
Annual (Maximum) series or block (annual) maxima and Peak over Threshold (or POT) are the two holistic technique widely accepted in the frequency analysis of hydrologic extreme [49]. In this experiment, we adopted the block (annual) maxima approach in the delineation of the triplet intercorrelated flood vectors, i.e., P, V, D series from the daily basis streamflow discharge records collected at the Gulliemard Bridge gauge station. For each observation year, there are only one flood episodes at the target site which is characterized through pointing the maximum discharge flow value, also called flood peak flow, P and their hydrograph volume, V and duration, D series are extracted corresponding to each flood peak using the algorithm as described in the literature such as Zhang and Singh [2], Yue and Wang [6]. Table 1, defines the descriptive behaviour of the derived annual flood events characteristics and also Figure 4, illustrates the histogram distribution of the derived annual flood characteristics and which indicates the existence of positively skewed distribution behaviour. In the multivariate distribution analysis, it is often an essential concern to check the existence of stationarity within individual flood characteristics before introducing into univariate or multivariate copula framework. For this, Ljung and box [50] based hypothesis testing also called Q-statistics, are undertaken [51]. Therefore, the Q-statistics for the sample size 'm' with 't' defines the total no of lag being tested with sample autocorrelations at the specific lag, i.e., ̂ are given below; Where, Null hypothesis (H 0 ) = zero autocorrelation or independent distributions; and Alternative hypothesis (H 1 ) = existence of serial correlation (or autocorrelation). The Q-statistics and their associated p-value for different lag size (i.e., 30, 20, 15, 10 and 5) are listed in Table 2 and which indicates for almost zero first-order autocorrelations because their estimated statistics are less than the critical value for each of the univariate series of flood characteristics by accepting the null hypothesis (H 0 ) at 5% or 0.05 significance level against their alternative hypothesis (H a ).

Modelling Flood Characteristics using Vine (or PCC) Model
The fitting procedure of the vine or pair-copula distribution framework to the delineated trivariate flood characteristics, i.e., flood peak (P), volume (V) and duration (D) comprise a series of computation task. The first stage of modelling comprises the selection of an appropriate vine tree structure (i.e., either D-or C-vine) that will be solely depending upon the strength of dependencies among flood characteristics pair and then after, to select the most parsimonious bivariate copulas families for pair-copulas in establishing different joint distribution pairs of flood characteristics and estimation of their dependence parameters. Different statistical tools are available that can support in the selection of most logical vine tree structure such as, graphically based, i.e., scatter plot, Kendall's plot, chi-plot and also analytically based correlation coefficients measures. Such tools can also be employed to select the most parsimonious parametric bivariate copula family for the underlying flood characteristics for describing the dependence structure.

Selection of Appropriate Vine Tree Structure
Figures 5, 6 and 7 illustrated the scatter plots, chi plots [52] and Kendall's plots [53] for analyzing the degree of interassociation and strength of dependency between flood attribute pairs. In actuality, the Chi-plot uses data ranks, also values measure the distance of bivariate random observations (say ) from the center of the data sets within the range of [-1, 1]. Thus, their values will tend to positive or negative for positively or negatively correlated random pairs, Also, the control limits are the another measuring factor in chi-plot that are placed at = ± √ ⁄ [52]. Therefore, in the case of stronger dependency exhibited, the random observations must be outside the control limit of chi-plot otherwise, inside the control limit region can be indicated for independence between random pairs. Similarly, Kendall's plot is analogous to quantile-quantile (Q-Q) plot such that, deviation of random observations from the main diagonal of K-plot reveals the existence of mutual correlation otherwise could be revealing for independence when the plot tends to be linear [14,53]. Based on the above graphical investigation, it could be revealed that flood peakvolume (P, V) pair exhibited strong and positive correlation structure while the flood pairs volume-duration (V, D) and peak-duration (P, D), both exhibited weak and negative correlation. Table 3 indicates the correlation coefficient statistics and thus, which support the similar conclusion taken from the graphical-based flood characteristics plots. Therefore, it must be concluded that the dependencies between the flood pairs (P, V) and (P, D) is higher than (V, D) that means the flood peak, P must be placed between the other two variables flood volume, V and duration, D, as graphically illustrated in Figure 2. According to Graler et al. [3], the choice of the conditioning variable (i.e. P, in Figure 2) is not unique, and different choices might lead to different results also, according to Hobaek et al. [35], different vine copula decomposition procedure differently approximates the multivariate distribution. Overall, it must be concluded that the D-vine copula could be a right choice and thus employed in this study for constructing 3dimensional flood characteristics and which are more flexible than the canonical or C-vine copula because the D-vine model construction provides higher flexibility in the selection of which flood attribute pairs to model than C-vine [23,34].

Simulation of Vine Tree by Fitting 2-Dimensional Copulas
After the selection of the appropriate vine framework, the bivariate joint dependence structure of the flood attribute pairs is defined through selecting the most appropriate 2-dimension (or bivariate) copulas function, such as C PV for (P-V) flood pairs and C PD for (P-D) pair copulas (referred to Figure 2). Similarly, the bivariate copulas, i.e., C PV and C PD are conditioned under the variable P (i.e., flood peak) [23] and thus the conditional CDFs, also called the h-functions are derived using the Equations 5 and 6, in defining the realization demands in the second tree (or Tree-2) of vine structure. Finally, a bivariate copulas C VD|P is fitted to the conditioned observations i.e., F V|P (v|p) and F D|P (d|p). Actually, the conditional cumulative distribution values are estimated for all the triplet flood characteristics in the vine tree structure.
Before initiating the copula fitting procedure for vine framework, we introduced commonly used univariate probability functions which belongs to parametric families such as Gamma-3P, Generalized Extreme Value (GEV), Generalized Gamma-3P, Inverse Gaussian-2P, Johnson SB-4P, Lognormal-2P and Weibull-2P, as a candidate function for defining the marginal structure of P, V and D series individually. All the univariate distribution fitting procedures along with their parameter estimations are carried out using the Easyfit software. The empirical observations are retrieved from commonly used Gringorten based position-plotting formula [11,54]. The best-fitted probability function is selected for each individual flood characteristics using different goodness-of-fit (GOF) test statistics such as Kolmogorov-Smirnov (or K-S) test and Anderson-Darling (or A-D) test [55], based on information criteria statistics such as Akaike Information criteria (or AIC) [56], Schwartz's Bayesian Information criteria (or BIC) [57] and Hannan-Quinn Information criteria (HQC) [58]. Table 4 illustrates the GOF test statistics, which clearly pointed out that the 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. Figure 2 illustrated the D-vine structure incorporated in the present demonstration. The pair copula exhibited three bivariate copula functions, i.e., C PV and C PD in the first tree (or Tree 1) and C VD|P in the second tree of the vine structure. In this demonstration, we introduced and investigated different copulas families which include one Gaussian copula that belongs to the Elliptical family and the Archimedean copula families such as the Gumbel, Frank, Clayton and Joe copula, as a possible candidate model. In actuality, the Gumbel, Clayton and Joe, copula cannot be used for negatively correlated flood characteristics (i.e., Kendall's tau < 0), which is only applicable in the modelling of positively correlated random variables. The copula dependence parameters fitted to the observed flood characteristics are estimated using the Maximum-Pseudo-Likelihood (MPL) estimation procedure [46,47] and their estimated values are listed in Table 5. The Cramer-von Mises (CvM) test statistics 'S n ' with the parametric bootstrapping procedure is employed to evaluate the adequacy of hypothesized copulas fitted to bivariate flood pairs (P-V) and (P-D), i.e., C PV and C PD copulas, which is often considered as one of the most potent model compatibility testing. Therefore, minimum the value of 'S n ' test statistics must indicate for minimum gap or dispensary between empirical and derived parametric copulas, in other words, this copula would be most justifiable and effective for establishing the bivariate joint relationship. Overall, from the Table 5, it must be concluded that the Gaussian copula and the Frank copula recognize as most justifiable for defining bivariate joint dependence structure of flood pairs (V, P) and (P, D) in the first tree (or Tree-1) of vine structure. Figure 8 illustrated the density and contour plot derived from the bivariate copulas fitted in the first tree between the flood pairs V-P and D-P while Figure 9 indicates the bivariate contour plot derived from joint CDFs of the Gaussian copula (C VP ) and the Frank copula (C DP ). Similarly, based on GOF test statistics, the Gumbel copula is recognized as the most parsimonious bivariate structure to model the variable (V, D) conditional on P, in the second tree (or Tree-2) (see Figure 10, which illustrates the surface plot of bivariate Gumble copula fitted with conditioned observation in the second tree of vine structure). Therefore, for the 3-dimensional case of flood characteristics, the three selected bivariate copulas such as C PV , C PD , and C VD|P are introduced and composed into the 3-dimension vine copula using the Equation 7. Simulation of the flood characteristics taken from the vine or paircopula model is based on the simulation approach as discussed by the literature such as Aas et al. [23]. Figure 11 illustrates the comparison based on scatter plot between the observed and simulated flood characteristics and which indicates a good agreement based on the fitted models, in other words, the vine model reproduces the overall flood correlation structure reasonably well.

Estimating Return Period of Flood Characteristics
Accounting of multiple interactive flood characteristics, i.e., P, V, D through the multivariate analysis of extreme episodes is much beneficial and useful in recognizing the critical behaviour of flood episodes at a river basin scale and for making a justifiable risk-based decision making in both structural and non-structural water-related issues. According to Salvadori and De Michele [9] and Salvadori [59], the primary joint return periods are estimated using the inclusive probability such as, "OR" and "AND" case. The primary joint return period for "OR" case for the trivariate and bivariate events can be estimated as; Where x 1 = F V (v), x 2 = F P (p) and x 3 = F D (d) is the univariate marginal distributions of the flood characteristics and C(x 1 , x 2 , x 3 ) is the trivariate copula function.
Similarly, the "OR-joint" case for bivariate events; T (V,P) OR (v, p) = , for V-D pair (12) Similarly, the "AND-joint" return periods for trivariate and bivariate flood characteristics are estimated as; T (P,D) At first, the univariate return periods are obtained from the best fitted marginal distribution of the flood characteristics, i.e., lognormal-2P for peak(P), Johnson SB-4P for volume (V) and Gamma-3P for duration (V) series (see Table 6). The joint return periods for 'OR' and 'AND' case for the trivariate flood distributions using the PCC model as well as bivariate distribution are derived from the Equation (9)(10)(11)(12)(13)(14)(15)(16), and their value are listed in the same Table 6. It is clearly revealed that the primary joint return periods in "AND" cases produce a higher return than the joint return in 'OR' cases for both the bivariate and trivariate distribution (based on PCC model). Investigating the tail dependence properties and their coefficient is an essential aspect in the extreme value analysis otherwise might be attributed for underestimation of hydrologic risk [22]. In this study, we also investigated the tail dependence behaviour of the derived D-vine copula for the flood characteristics extracted from the fitted copulas through both the graphical and analytical observation. In the graphical approach, we visualized and compared the scatter plot, Chi-plot and Kendall or K-plot of the random observations derived from the bivariate Normal (or Gaussian) copula, C VP (as fitted to (V, P) flood pair) and the Frank copula C PD (as fitted to (P, D) flood pair) with the D-vine copula fitted to (V, P) and (P, D) pairs (refer to Figures 12 and 13). From the Figures 12 and 13, it is revealed that vine copula is much consistent and well capturing the general behaviour of both the Normal copula fitted for (V, P) pair and Frank copula fitted for (P, D) pair, based on the scatter plot (i.e., the first row of Figures 12 and 13). Also, based on Chi-plot and K-plot it indicates the quite similar tail dependencies behaviour of both the Normal copula and vine copula for (V, P) and the Frank copula and vine copula for (P, D) flood characteristics (see 2 nd and 3 rd row of Figures 12 and 13). The non-parametric estimation of tail dependence coefficient, based on the Caperaa-Fougeres-Genest estimator λ U CFG [60] is employed to estimate the tail dependence of the fitted vine model, as given below; Where ( 1 , 1 ), ( 2 , 2 ), … … . . , ( , ) are the random observations simulated from the PCC model fitted to the flood characteristics. Thus, using the Equation (16), the tail dependence coefficients (see Table 7) are estimated for each flood attributes pairs i.e., (V-P), (P-D), (V-D) based on the sample generated from fitted vine or pair-copula model.

Conclusions
The necessity of multivariate flood distribution analysis is often an essential concern among the hydrologist, which usually reveals an effective outcome during the assessments of hydrologic risk. The higher dimension flood distribution modelling in the light of the standard trivariate copula function usually surrounded with few statistical constraints and issues. In such circumstances, the incorporation of the vine or pair-copula construction (or PCC) could be the logical choice for revealing a comprehensive and practical approach in the uncertainty analysis and which is solely based on stagewise hierarchical nesting procedure. In this demonstration, the parametric vine copula distribution framework is incorporated for the block (annual) maxima flood characteristics collected at the Gulliemard Bridge gauge station in the Kelantan River basin in Malaysia. Tests for stationarity or existence of serial correlation (or autocorrelation) within the individual time series of flood characteristics is often a pre-requisite before introducing the random samples into a univariate or a multivariate framework. It was found that the time series of the flood characteristics do not exhibit any significant trend. Different dependency measures graphical approaches such as scatter plot, Kendall's plot, chi-plot and also analytically based correlation coefficients statistics are estimated to select the most appropriate vine tree structure. Investigation revealed that the D-vine copula could be a right choice and thus employed in constructing the hierarchical nesting structure of 3-dimensional flood characteristics where the flood variable, P is placed in between variable, V and D because the dependency between (P-V) and (P-D) is higher than (V-D). In actuality, the D-vine copula usually exhibited higher flexibility in the selection of which flood attribute pairs to model, in comparison with the canonical or C-vine copula [23,34]. Different copulas families which include the Gaussian copula that belongs to the Elliptical family and the Archimedean copula families such as the Gumbel, Frank, Clayton and Joe copula, are introduced and test as a possible candidate model in demonstrating the pair-wise (or bivariate) joint dependency modelling of flood characteristics, i.e., (P-V) and (P-D) in Tree-1 and also, the bivariate copulas fitted to the conditioned observations in the Tree-2 of the selected D-vine structure (see Figure 2). The vector of unknown copula dependence parameters, in both Tree-1 & Tree-2 are estimated using the Maximum-Pseudo-Likelihood (MPL) estimation procedure where the best-fitted bivariate copulas are selected using the Cramer-von Mises (CvM) test statistics 'S n ' with the parametric bootstrapping procedure.
The Lognormal-2P distribution, Johnson SB-4P distribution and Gamma-3P distribution are identified as possible univariate marginal distribution for the flood peak flow, volume and duration series based on the GOF test statistics. Firstly, the Gaussian copula and Frank copula are selected for defining the bivariate joint dependence structure of flood pairs (V, P) (i.e., 2-D copula fitted to peak-volume pair, C PV ) and (P, D) (i.e., 2-D copula fitted to peak-duration pair, C PD ) in the first tree of vine structure (see Figure 2). Similarly, the Gumbel copula is recognized as the most justifiable bivariate structure in the modelling of the variable (V, D) conditional on P (i.e., the conditioned observations, F V|P (v|p) and F D|P (d|p)), in the second tree of vine framework. The scatterplot between observed and simulated flood characteristics are investigated and which indicating on an excellent agreement or in words the fitted D-vine distribution reproduces overall flood correlation structure reasonably well. The tail dependence behaviour of the fitted D-vine copula is also investigated for flood characteristics extracted from the fitted copulas through both the graphical and analytical observation and which reveals that vine copula is well capturing the general behaviour of both the fitted Normal and Frank copula for flood pairs (V, P) and (P, D). The primary joint return periods for trivariate (also, bivariate) flood events are examined and revealed that "AND" return period produces higher return value than 'OR' joint case. Unless the modelling efficacy of vine copula model is investigated in the above demonstration and from few existed kinds of literature, but it is still surrounded with the few modelling issues and constraints such that selection of best-fitted parametric copulas in the vine framework is often a tough challenge [19]. In such circumstances, the non-informative vine methodological approach, also called the minimum information-based vine construction could facilitate higher flexibility and accentuate the potential of traditional parametric vine approach, especially when the length of the observation samples under consideration are limited [13]. Daneshkhan et al. [19] study already demonstrated the performance of minimum information pair-copula and thus revealed that it could be possible to approximate multivariate density to any level of desired precision. All the above-raised issues can be considered as a future research purpose for the same river basin in order to achieve much comprehensive approximation of joint probability density and practical flood hazard assessments.