Assessing the Impacts of Univariate and Bivariate Flood Frequency Approaches to Flood Risk Accounting for Reservoir Operation

Flood frequency analysis plays a fundamental role in dam planning, reservoir operation, and risk assessment. However, conventional univariate flood frequency analysis carried out by flood peak inflow or volume does not account for the dependence between flood properties. In this paper, we proposed an integrated approach to estimate reservoir risk by combining the copula-based bivariate flood frequency (peak and volume) and reservoir routing. Through investigating the chain reaction of “flood frequency—reservoir operation-flood risk”, this paper demonstrated how to simulate flood hydrographs using different frequency definitions (copula “Or” and “And” scenario), and how these definitions affect flood risks. The approach was applied to the Meishan reservoir in central China. A set of flood hydrographs with 0.01 frequency under copula “Or” and “And” definitions were constructed, respectively. Upstream and downstream flood risks incorporating reservoir operation were calculated for each scenario. Comparisons between flood risks from univariate and bivariate flood frequency analysis showed that bivariate flood frequency analysis produced less diversity in the results, and thus the results are more reliable in risk assessment. More importantly, the peak-volume combinations in a bivariate approach can be adjusted according to certain prediction accuracy, providing a flexible estimation of real-time flood risk under different prediction accuracies and safety requirements.


Introduction
Floods are major natural disasters that occur across the globe.Numerous reports suggest that floods are causing massive property loss and deaths, especially in the last two decades [1,2].Flood events are complicated in nature and are primarily characterized by the flood peak, volume, duration, and hydrograph shapes [3].Most flood control measures require detailed knowledge of flood events, which have a certain frequency, also expressed as the T-year return period.However, the traditional flood frequency model is a univariate that is constructed based on a selected characteristic-peak or volume-with target flood frequency [4,5].Indeed, a single characteristic is inadequate to represent flood frequency: given a flood event with a high peak flow and low volume, its calculated frequency will be inconsistent when using the univariate method, which may lead to over-or under-estimation of its risk [6][7][8].Therefore, a multivariate method was proposed to perform a more comprehensive and objective estimation on flood frequency.In earlier work, several multivariate distributions such as lognormal bivariate distribution [9], bivariate exponential distribution [10], and bivariate gamma distribution [11] have been widely applied to flood frequency analysis.A significant limitation of these multivariate distributions is that the marginal distribution of variables must be consistent.However, it is unrealistic, since flood characteristics are random and may probably have different distribution shapes.
Recent developed copula theory provides a new means of flood frequency.Copula theory, developed by Sklar [12], is an effective approach to describing the joint probability distribution of multi-dependent variables.In copula theory, the correlated variables are not necessarily required with the same marginal distribution, which makes it widely applied in hydrological analysis.An important application of the copula theory is flood frequency analysis, considering flood peak, volume and duration, etc. [13][14][15].Based on joint distribution of flood frequency, floods at any frequency (usually at distribution tails) can be constructed [16][17][18].Besides this, the applications of copula theory have been extended to spatial and temporal relationships of floods.Previous studies [19] investigated the spatial distribution of flood events among different tributaries and sub-areas, while some constructed the joint distribution of flood events in different sub-seasons, and then optimized reservoir flood limit water level by taking the distribution as constraint to achieve better conservancy benefit [20].In the same conceptual framework, copula theory has been applied to hydrological problems such as multivariate rainfall analysis considering rainfall intensity-duration-depth [21], and rainfall peak-depth [22], a trivariate drought assessment considering drought severity-duration-frequency [23], etc. Numerous applications have proved that copula-based multivariate frequency analysis better fits the empirical observations than univariate or standard joint parametric distribution methods.
Reservoirs can effectively change flood waves to decrease flood risks.Increasingly, more natural floods are being controlled and regulated by reservoirs.In the past several decades, a number of flood risk assessment studies on reservoir operation have been proposed [24].Studies about reservoir optimal operation usually focus on reservoir overtopping risk, which is also the risk posed by the reservoir itself [25,26].In this case, the research priority is then a question of how to optimize reservoir operations for appropriate refilling and drawing down.On the other hand, flood risk assessment for a river basin is more focused on downstream flooding risk [27][28][29], because the downstream area often comprises large cities and populations to be protected.In such studies, reservoir operation procedures were often simplified, and the upstream risk was assumed constant.For instance, reservoir outflow was considered a piecewise linear function of inflow [30,31].A few studies have incorporated both upstream and downstream risks into their risk assessment, while the flood hydrographs were synthetic and simplified [32].
Flood frequency is also a key issue that determines reservoir operational response as well as flood risk.However, studies usually focus on the standard of flood frequency in dam planning phase, but seldom specify on which frequency definition the standard is based, i.e., peak or volume or entire hydrograph; in real time flood management, flood frequency on natural flood is not equal to real flood risk when accounting for reservoir operation [33].The chain effect-from flood frequency analysis to its consequences on real flood risk accounting for reservoir operations-is relatively unexplored.As part of the chain effect, Reference [34] analyzed the influence of flood hydrographs on reservoir operation and downstream damage.Four hydrographs: triangular, abrupt wave, flood pulse, and broad peak were used to compare downstream damages.Although the hydrographs were synthetic for simplification, this study proved that hydrograph shapes significantly influence downstream damage, even for hydrographs with same volumes and peak flows.Reference [35] concerned the relationship between flood frequency and reservoir operation, and constructed peak-volume bivariate flood frequency to calculate overtopping probability of dam.Reference [36] assessed the impacts of reservoir operation to downstream inundation with hydrodynamic models, but the study case is only one real flood event which happened in 2011.
The chain effect was analyzed more in depth in some other researches.Reference [37] assessed reservoir downstream flood risk using bivariate copula analysis of flood peaks in two tributaries.Three scenarios were considered and compared: 1. no reservoir built; 2. reservoir built in tributary A; 3. reservoir built in tributary B. The results showed that downstream discharge return period differs strongly among the three scenarios, which revealed that the real downstream risk is complicatedly influenced not only by whether a reservoir exists, but also by the location and storage of the reservoir.Similarly, Reference [38] highlighted the influence of reservoir operation to flood frequency, and defined a concept called "routed return period" which incorporated reservoir operation into the theoretical flood return period.This study revealed how dam characteristics such as reservoir volume and spillway length significantly influence the risk of dam overtopping.The focus of References [37,38] are close to the idea proposed in this paper.However, Reference [37] accounted for only downstream risk and reference [38] accounted for only upstream risk (the author mentioned in conclusion that downstream risk could be replicated in the same way), and both of them simplified the reservoir routing process.In summary, the impact of reservoir operation on flood frequency has attracted much attention, but overlooked points still exist.
This study aims to bridge the gap between flood frequency analysis and actual flooding risk accounting for reservoir operation.To this purpose, the following were set as two main objectives: (1) to construct multivariate flood frequency analysis model; (2) to establish the calculating chain effect of "flood frequency-reservoir operation-risk assessment".More specifically, the remainder of this paper is organized as follows: Section 2 introduces the study site-Meishan reservoir-and its historic flood characteristics.Section 3 presents (1) copula theory based multivariate flood frequency contours and generated hydrographs, (2) a real-time reservoir routing scheme to operate the constructed flood hydrographs, and (3) the derivation and integration of upstream and downstream risks based on reservoir operation results.Section 4 reports case study results with a comparative analysis and discussion.Finally, the conclusion of this study is provided in Section 5.

Study Site
The Meishan reservoir is located in the Shi river in central China, as Figure 1 shows.The Shi river feeds into the Huai river, which is one of the largest rivers in China.The river basin is in the north-south climate transboundary zone, and thus the reservoir inflow varies greatly from year to year and is concentrated in the flood season, from June to September [39,40].The average precipitation at the Shi river basin is 1320 mm per year.About 70% of the annual precipitation occurs during the flood season, and over 70% of the flood season's precipitation is in the form of flood.As a result, flood control is the primary purpose of the Meishan reservoir.The greatest historic flood event occurred from July to August 1991: the flood peak reached 8500 m 3 /s, and the total flood volume reached 1.6 billion m 3 , which is 3 times more than the reservoir's flood storage capacity.During this extraordinary flood, the Meishan reservoir reduced 70-85% of the peak inflow using its flood storage capacity, and decreased the downstream flood depth by 4 m, demonstrating its significant flood control effects of the overall basin.
The operating rules of the Meishan reservoir are based on a series of specified reservoir water levels.During the flood season, the lower bound and upper bound of the Meishan reservoir safety water level are 125.27m (corresponding to 1.20 billion m 3 storage) and 133.00 m (corresponding to 1.70 billion m 3 storage), respectively, leaving a flood control storage capacity of 0.5 billion m 3 .A flood can be retained in the reservoir during the peak inflow period, and then released gradually to reduce downstream loss.The maximum downstream safety release is 1200 m 3 /s, but if the reservoir level exceeds 133.00 m, the dam is at risk, and the reservoir release is then permitted to exceed 1200 m 3 /s Water 2019, 11, 475 4 of 20 as a trade-off to assure dam safety.The general operating principle of the Meishan reservoir is summarized as follows, and the specific rules are introduced and simulated in Section 4.4.The operating rules of the Meishan reservoir are based on a series of specified reservoir water levels.During the flood season, the lower bound and upper bound of the Meishan reservoir safety water level are 125.27m (corresponding to 1.20 billion m³ storage) and 133.00 m (corresponding to 1.70 billion m³ storage), respectively, leaving a flood control storage capacity of 0.5 billion m³.A flood can be retained in the reservoir during the peak inflow period, and then released gradually to reduce downstream loss.The maximum downstream safety release is 1200 m³/s, but if the reservoir level exceeds 133.00 m, the dam is at risk, and the reservoir release is then permitted to exceed 1200 m³/s as a trade-off to assure dam safety.The general operating principle of the Meishan reservoir is summarized as follows, and the specific rules are introduced and simulated in Section 4.4.

Flood Event Characteristics
In this study, a total of 37 flood events during the period of 2003-2016, as well as the largest reservoir release q = inflow Q, when Q < 1200 m³/s and reservoir water level < 133 m, res water level is kept constant =1200 m³/s, when Q > 1200 m³/s and reservoir water level < 133 m, at the maximum downstream safety flow > 1200 m³/s, when reservoir water level > 133m, draw down water for dam safety, downstream may suffer as a trade-off

Flood Event Characteristics
In this study, a total of 37 flood events during the period of 2003-2016, as well as the largest historic flood in 1991, are employed as a case study.Three characteristics, flood peak inflow, volume, and duration are described in Table 1.The standard variance indicators reveal that floods with Meishan characteristics are quite diverse.Figure 2 shows scatterplots of the peak-volume relationship, peak-duration relationship, and volume-duration relationship in order to examine correlations among the three variables.The Kendall coefficient and Spearman coefficient [41] for each combination are given as well.It can be observed from the scatter trend and coefficients that the peak-volume combination has the most significant correlation compared with the other two combinations.Therefore, in this study, peak and volume are selected to model the probabilistic distribution for the Meishan reservoir.

Marginal Distributions Based on the Nonparametric Method
The bivariate frequency is primarily based on the marginal distribution of each variable.The widely used parametric methods for marginal distribution estimation must assume a prior distribution of data.For hydrological statistical variables, popular marginal distributions include the Normal distribution [42], Pearson III distribution [43], Generalized Extreme Value (GEV) distribution [44], and log-logistic (log-log) distribution [45], among others.However, parametric method has its limitations.There might not be a proper distribution function for the data, and the probability might vary considerably with just a slight modification of the data.
To overcome the disadvantage of the parametric method, a less restrictive nonparametric method, kernel density estimation (KDE), was adopted to estimate the flood peak and volume marginal distribution [46,47].The KDE method can yield an empirical estimate of the probability density function without assuming any form for the distribution.For a random sample  ,  ,… , the traditional kernel probability density function (PDF) can be expressed as: where T is the sample size, h is the bandwidth that determines the smoothness of the estimate-a greater value of h will produce a smoother estimate and vice versa-and K(.) is the kernel density estimator.Commonly used kernels include the Gaussian kernel, Uniform kernel, Triangular kernel, Tricube kernel, etc. [48].For simplicity, a Gaussian kernel [49] was applied in this study.The bandwidth has a dramatic influence on the estimate.In this study, the bandwidth was estimated by enumeration, and the optimization standard is the fitting performance of the overall frequency function.

Marginal Distributions Based on the Nonparametric Method
The bivariate frequency is primarily based on the marginal distribution of each variable.The widely used parametric methods for marginal distribution estimation must assume a prior distribution of data.For hydrological statistical variables, popular marginal distributions include the Normal distribution [42], Pearson III distribution [43], Generalized Extreme Value (GEV) distribution [44], and log-logistic (log-log) distribution [45], among others.However, parametric method has its limitations.There might not be a proper distribution function for the data, and the probability might vary considerably with just a slight modification of the data.
To overcome the disadvantage of the parametric method, a less restrictive nonparametric method, kernel density estimation (KDE), was adopted to estimate the flood peak and volume marginal distribution [46,47].The KDE method can yield an empirical estimate of the probability density function without assuming any form for the distribution.For a random sample x 1 , x 2 , . . .x n , the traditional kernel probability density function (PDF) can be expressed as: where T is the sample size, h is the bandwidth that determines the smoothness of the estimate-a greater value of h will produce a smoother estimate and vice versa-and K(.) is the kernel density estimator.Commonly used kernels include the Gaussian kernel, Uniform kernel, Triangular kernel, Tricube kernel, etc. [48].For simplicity, a Gaussian kernel [49] was applied in this study.The bandwidth has a dramatic influence on the estimate.In this study, the bandwidth was estimated by enumeration, and the optimization standard is the fitting performance of the overall frequency function.

Concept of a Copula
Copulas are multivariate cumulative distribution functions (CDFs) that describe the dependence of variables marginal distributions [50].Taking bivariate copula as an example, for two random variables X and Y, their marginal CDF u and v can be represented by: Water 2019, 11, 475 6 of 20 where u, v ∈ [0, 1].Similarly, their joint CDF can be expressed as a copula function C(u, v) and a marginal distribution of variables: Apparently, C(0,0) = 0, and C(1,1) = 1. Figure 3 explains the relationships among u, v, and C(u, v).Through the copula function, the joint distribution of X and Y is simplified to the correlation between u and v.Given marginal distributions u and v, we can calculate their copula function C(u,v) (red arrowhead in Figure 3); conversely, given a certain value of a copula function C(u,v), a series of (u,v) and (x,y) with the same C(u,v) value can also be found.

Concept of a Copula
Copulas are multivariate cumulative distribution functions (CDFs) that describe the dependence of variables marginal distributions [50].Taking bivariate copula as an example, for two random variables X and Y, their marginal CDF u and v can be represented by: ( ), ( ) where , [0,1] u v∈ . Similarly, their joint CDF can be expressed as a copula function ( , ) C u v and a marginal distribution of variables: Apparently, C(0,0) = 0, and C(1,1) = 1. Figure 3 explains the relationships among u, v, and ( , ) Through the copula function, the joint distribution of X and Y is simplified to the correlation between u and v.Given marginal distributions u and v, we can calculate their copula function C(u,v) (red arrowhead in Figure 3); conversely, given a certain value of a copula function C(u,v), a series of (u,v) and (x,y) with the same C(u,v) value can also be found.Selecting proper copula function type and parameters is essential for correctly modeling a multivariate distribution.Table 2 lists several commonly used bivariate copula functions.For the parameter θ estimation, the Maximum Likelihood Estimation (MLE) method, Inference Function for margins (IFM), and Kendall correlation methods are widely utilized [51,52].In this study, the Kendall correlation method was used to calibrate parameter θ, for its fitting performance was shown to be superior to the other methods in our previous tests.Selecting proper copula function type and parameters is essential for correctly modeling a multivariate distribution.Table 2 lists several commonly used bivariate copula functions.For the parameter θ estimation, the Maximum Likelihood Estimation (MLE) method, Inference Function for margins (IFM), and Kendall correlation methods are widely utilized [51,52].In this study, the Kendall correlation method was used to calibrate parameter θ, for its fitting performance was shown to be superior to the other methods in our previous tests.
The copulas described by mathematical functions in Table 2 are called theoretical copulas.Actually, we can also construct the "true" copula by directly counting observations and calculating probabilities for each data, which is called the empirical copula.The empirical copula, denoted by Ĉ(u, v), can be expressed as the empirical distribution, as defined by Equation (4): where n denotes number of observed samples.
The empirical copula can be used to assess the fitting performance of the theoretical copula.Taking the root-mean-square error (RMSE) as an indicator [53], the fitting performance can be calculated as: 3.2.2.Derivation of Copula "Or" and "And" Scenarios In flood frequency analysis, we are interested in extreme events, so the flood frequency is mostly defined by exceedance probability, such as: where X can be the flood peak, volume, or other characteristics.Obviously, the frequency concept in hydrology is different from the non-exceedance probability in Equation (3).In order to link copula probability and flood frequency, two types of joint probabilities are considered: 1.
Copula "Or" probability, denoted by C or .C or represents P(peak ≥ x 1 or volume ≥ x 2 ).Given the concept copula C(u, v) = P(peak ≤ x 1 , volume ≤ x 2 ), C or can be described as: 2.
"And" probability, denoted by C and .C and represents P(peak ≥ x 1 , volume ≥ x 2 ).Similarly, C and can be described as: Figure 4 explains the derivation of C or and C and .The area of C or is a subset of C(u,v), as described in Equation (7).C and is the upper-right zone in Figure 4, and the area can be derived by Equation (8).Once the copula function C(u,v) is determined, C or and C and can be obtained.

Risk Assessment For Reservoir Flood Operations
Flood risk can be defined as the possibility of suffering loss or the risk probability multiplied by its consequences [54][55][56].In reservoir flood operations, the exact loss value caused by a flood may be related to massive downstream economic and social issues that are quite difficult to quantify.For

Risk Assessment For Reservoir Flood Operations
Flood risk can be defined as the possibility of suffering loss or the risk probability multiplied by its consequences [54][55][56].In reservoir flood operations, the exact loss value caused by a flood may be related to massive downstream economic and social issues that are quite difficult to quantify.For simplicity, in this study, the flood risk is denoted by "risk probability".

Integration of Upstream and Downstream Risks
There are two sources of reservoir flood operation risk: upstream and downstream.If the flood is mostly retained in the reservoir with little release, the reservoir water level will rise, so the upstream risk, which is mainly called overtopping risk [35,38], will increase, but the downstream risk will decrease.Conversely, if the reservoir releases most of the flood downstream to empty the reservoir, the upstream risk will be very low, while the downstream might suffer and is thus at risk.In this study, both upstream and downstream risks were taken into account.To describe the integrated risk R, this paper adopted the idea of flood probability integration in Reference [32], expressed as: where A and B represent flood risk in current time step and the next time step, respectively.By substituting A and B with the upstream risk (R u ) and downstream risk (R d ), the overall risk R can be calculated as: Equation ( 10) can also be written as: The implication of Equation ( 11) is that the risk probability R is the complementary set to the probability that no risk is realized.

Quantization of Upstream Risk R u
Risk upstream of the reservoir is proportional to the reservoir storage.The more water stored in the reservoir, the greater the risk when flood occurs.Therefore, the upstream risk R u can be constructed as follows: where V is the maximum reservoir storage during one round of the flood control operation, and V ub and V lb represent the upper bound and lower bound of reservoir storage, respectively.V ranges from V ub to V lb ; R u ranges from 0 to 1 and increases linearly with the increase of V.As introduced in Section 2.1, V ub and V lb are 1.70 billion m 3 and 1.20 billion m 3 in this study.

Quantization of Downstream Risk R d
Risk downstream of the reservoir is caused by excess downstream safety release q safe .If the reservoir release is less than q safe , R d is 0; when the reservoir release is greater than q safe , R d can be defined as follows: R d = max(0, (q − q safe )/(q max − q safe )) (13) where q is the maximum reservoir release during one round of the flood control operation.For the Meishan reservoir, q max and q safe are 7000 m 3 /s and 1200 m 3 /s, respectively.R d varies from 0 to 1 and increases linearly with reservoir release.Figure 5 summarized the definitions of R u , R d , and explained the derivation of R.
Water 2019, 11, 475 9 of 20 Rd = max(0, (q − qsafe)/(qmax − qsafe)) (13) where q is the maximum reservoir release during one round of the flood control operation.For the Meishan reservoir, qmax and qsafe are 7000 m³/s and 1200 m³/s, respectively.Rd varies from 0 to 1 and increases linearly with reservoir release.Figure 5 summarized the definitions of Ru, Rd, and explained the derivation of R.

Results and Discussion
Following the model procedure developed in Section 3, the flood risk of the Meishan reservoir was calculated and discussed comprehensively.The copula model computing and flood risk calculation were implemented using Matlab software, version r2014a.

Marginal Distribution and Copula Distribution of Peak and Volume
The marginal distribution of peak and volume were derived using the KDE method.Figure 6 shows the empirical distribution (dashed line) and theoretical distribution (solid line).

Results and Discussion
Following the model procedure developed in Section 3, the flood risk of the Meishan reservoir was calculated and discussed comprehensively.The copula model computing and flood risk calculation were implemented using Matlab software, version r2014a.

Marginal Distribution and Copula Distribution of Peak and Volume
The marginal distribution of peak and volume were derived using the KDE method.Figure 6 shows the empirical distribution (dashed line) and theoretical distribution (solid line).Next, copulas were used to construct the joint distribution of peak and volume.The Clayton, Gumbel, and Frank copula types were selected for this case study.RMSE was used to evaluate the goodness-of-fit of each copula, and the results are shown in Table 3.It can be observed that the Clayton copula has the lowest RMSE; thus, the Clayton copula was chosen here to fit the peakvolume joint distribution.

Clayton Gumbel Frank
Copula parameter θ 0.87 0.66 0.95 Next, copulas were used to construct the joint distribution of peak and volume.The Clayton, Gumbel, and Frank copula types were selected for this case study.RMSE was used to evaluate the goodness-of-fit of each copula, and the results are shown in Table 3.It can be observed that the Clayton copula has the lowest RMSE; thus, the Clayton copula was chosen here to fit the peak-volume joint distribution.

Clayton Gumbel Frank
Copula parameter θ 0.87 0.66 0.95 RMSE 0.013 0.019 0.018 Figure 7 shows the joint PDF and CDF of peak and volume based on Clayton copula.The axis in Figure 7 refers to the marginal distribution value of peak and volume.It can be observed that the PDF of peak and volume has a strong tail correlation: the probability of low-peak-low-volume and high-peak-high-volume zone combinations is extremely large.

Peak-Volume Combinations in Copula "Or" and "And" Scenarios
Since the copula function has been determined, the function of Cor and Cand can be derived using Equations ( 7) and (8). Figure 9 shows contour figures of C(u,v), Cand(u,v), and Cor(u,v) with probability contour lines from 0.05 to 0.95.For each contour line with certain copula frequency, a group of (peak, volume) samples can be generated.This is the basis for the following flood hydrograph simulation.In order to better demonstrate the peak-volume joint distribution, 10,000 groups of peak-volume samples that fit the Clayton copula distribution were generated, as shown in Figure 8.The generated 10,000 samples are highly consistent with the historic samples, and the tail correlation is also very clear: samples concentrate in the low-peak-low-volume zone and high-peak-high-volume zone, which is consistent with PDF shown in Figure 7a.These results reveal that the peak and volume are quite relevant, and the Clayton copula function can describe this relevance precisely.

Peak-Volume Combinations in Copula "Or" and "And" Scenarios
Since the copula function has been determined, the function of Cor and Cand can be derived using Equations ( 7) and (8). Figure 9 shows contour figures of C(u,v), Cand(u,v), and Cor(u,v) with probability contour lines from 0.05 to 0.95.For each contour line with certain copula frequency, a group of (peak, volume) samples can be generated.This is the basis for the following flood hydrograph simulation.

Peak-Volume Combinations in Copula "Or" and "And" Scenarios
Since the copula function has been determined, the function of C or and C and can be derived using Equations ( 7) and (8). Figure 9 shows contour figures of C(u,v), C and (u,v), and C or (u,v) with probability contour lines from 0.05 to 0.95.For each contour line with certain copula frequency, a group of (peak, volume) samples can be generated.This is the basis for the following flood hydrograph simulation.

Flood Hydrograph Generation for Peak-Volume Combinations
There are three decisive elements that characterize a flood event: peak, volume, and typical flood hydrograph [57].Previous sections discussed how to derive peak and volume under certain copula frequency.In this section, the full flood hydrograph was simulated using given peak and volume.

Peak-Volume Combinations Selection with a 0.01 Frequency
Flood risk assessment concentrates on extraordinary flood events.In this study, peak-volume combinations with a 0.01 frequency were chosen as target flood events.Clearly, there are infinite peak-volume combinations on the 0.01 frequency contour line, so it is necessary to select limited combinations for further computation.Reference [58] proposed a probability based method to narrow the line to a subset, and then select freely from the subset.In this study, we firstly exported the peak-volume combinations of the contour line using the "contour" function in Matlab software, shown in Figure 10.It is noted that the lines in Figure 10 are not as smooth as contour lines in Figure 9, for the axis here are peak and volume instead of marginal distribution u and v.

Flood Hydrograph Generation for Peak-Volume Combinations
There are three decisive elements that characterize a flood event: peak, volume, and typical flood hydrograph [57].Previous sections discussed how to derive peak and volume under certain copula frequency.In this section, the full flood hydrograph was simulated using given peak and volume.

Peak-Volume Combinations Selection with a 0.01 Frequency
Flood risk assessment concentrates on extraordinary flood events.In this study, peak-volume combinations with a 0.01 frequency were chosen as target flood events.Clearly, there are infinite peak-volume combinations on the 0.01 frequency contour line, so it is necessary to select limited combinations for further computation.Reference [58] proposed a probability based method to narrow the line to a subset, and then select freely from the subset.In this study, we firstly exported the peak-volume combinations of the contour line using the "contour" function in Matlab software, shown in Figure 10.It is noted that the lines in Figure 10 are not as smooth as contour lines in Figure 9, for the axis here are peak and volume instead of marginal distribution u and v.In Figure 10, there are 56 combinations for the "And" scenario and 4 combinations for the "Or" scenario.All the 4 combinations for the "Or" scenario were selected for better representativeness.For the "And" scenario, 10 "unrealistic" combinations, where the volumes are not proportional to peaks (the peak value is even less than average inflow calculated from volume), were deleted.For the remaining combinations, we selected 12 of them with diverse peak/volume ratios.The selected combinations are shown in Figure 11.In real time operation, the range of peak-volume combinations can be narrowed considering different flood prediction levels, which will give more definite risk assessment in the following calculation.

Flood Hydrograph Simulation Based on Peak-Volume Combinations
Typical flood hydrograph is another important factor in a flood hydrograph simulation, since it determines the "shape" of a hydrograph.In conventional methods, a representative hydrograph is In Figure 10, there are 56 combinations for the "And" scenario and 4 combinations for the "Or" scenario.All the 4 combinations for the "Or" scenario were selected for better representativeness.For the "And" scenario, 10 "unrealistic" combinations, where the volumes are not proportional to peaks (the peak value is even less than average inflow calculated from volume), were deleted.For the remaining combinations, we selected 12 of them with diverse peak/volume ratios.The selected combinations are shown in Figure 11.In real time operation, the range of peak-volume combinations can be narrowed considering different flood prediction levels, which will give more definite risk assessment in the following calculation.In Figure 10, there are 56 combinations for the "And" scenario and 4 combinations for the "Or" scenario.All the 4 combinations for the "Or" scenario were selected for better representativeness.For the "And" scenario, 10 "unrealistic" combinations, where the volumes are not proportional to peaks (the peak value is even less than average inflow calculated from volume), were deleted.For the remaining combinations, we selected 12 of them with diverse peak/volume ratios.The selected combinations are shown in Figure 11.In real time operation, the range of peak-volume combinations can be narrowed considering different flood prediction levels, which will give more definite risk assessment in the following calculation.

Flood Hydrograph Simulation Based on Peak-Volume Combinations
Typical flood hydrograph is another important factor in a flood hydrograph simulation, since it determines the "shape" of a hydrograph.In conventional methods, a representative hydrograph is is chosen for all scenarios, so all the simulated target hydrographs have the same "shape" [59,60].A significant drawback of this method is that it ignores the diversity of different flood processes.In order to match flood hydrograph types with the peak-volume combinations, in this study, the typical flood hydrograph was selected according to the volume/peak ratio.Given the target peak and volume, the historical hydrograph whose volume/peak ratio is the closest to target ratio is selected as typical hydrograph for this combination.
Traditional flood hydrograph simulation usually considers only one characteristic-peak or volume-and then scales the overall typical flood hydrograph by a constant ratio to fit the target characteristic.In this study, both peak and volume are fixed, so the traditional method is no longer feasible since it can meet only one target characteristic.Here, a "variable ratio amplification" rescale method [61] was adopted to simulate a flood hydrograph: where Q(t) denotes the simulated flood inflow in t time steps; Q peak and w denote the simulated target flood peak and volume, which were obtained from the 0.01 frequency contour lines; Q denotes the average inflow of the simulated hydrograph, which can be calculated by w/T; T denotes the total time steps in the typical flood hydrograph; Q D (t), Q Dmax , and Q D denote the inflow in t time steps, maximum inflow, and average inflow of the typical flood hydrograph, respectively.In this way, both peak and volume values in the simulated flood hydrograph are equal to the target value; meanwhile, the general trend of the typical flood hydrograph can still be preserved.Hydrographs using a variable ratio amplification method for C and and C or scenarios were simulated, as shown in Figures 12 and 13.As it is seen here, different simulated flood hydrographs have different "shapes", especially in the C and scenario, as the volume/peak ratios are more diverse.The diversity of simulated flood peaks, volumes, durations, and hydrographs provides full representativeness for the resulting flood risks, compared with the risk from univariate flood simulation method.
Water 2019, 11, 475 14 of 21 chosen for all scenarios, so all the simulated target hydrographs have the same "shape" [59,60].A significant drawback of this method is that it ignores the diversity of different flood processes.In order to match flood hydrograph types with the peak-volume combinations, in this study, the typical flood hydrograph was selected according to the volume/peak ratio.Given the target peak and volume, the historical hydrograph whose volume/peak ratio is the closest to target ratio is selected as typical hydrograph for this combination.
Traditional flood hydrograph simulation usually considers only one characteristic-peak or volume-and then scales the overall typical flood hydrograph by a constant ratio to fit the target characteristic.In this study, both peak and volume are fixed, so the traditional method is no longer feasible since it can meet only one target characteristic.Here, a "variable ratio amplification" rescale method [61] was adopted to simulate a flood hydrograph: where ( ) Q t denotes the simulated flood inflow in t time steps; peak Q and w denote the simulated target flood peak and volume, which were obtained from the 0.01 frequency contour lines; Q denotes the average inflow of the simulated hydrograph, which can be calculated by / w T ; T denotes the total time steps in the typical flood hydrograph; ( ) , and D Q denote the inflow in t time steps, maximum inflow, and average inflow of the typical flood hydrograph, respectively.In this way, both peak and volume values in the simulated flood hydrograph are equal to the target value; meanwhile, the general trend of the typical flood hydrograph can still be preserved.
Hydrographs using a variable ratio amplification method for Cand and Cor scenarios were simulated, as shown in Figures 12 and 13.As it is seen here, different simulated flood hydrographs have different "shapes", especially in the Cand scenario, as the volume/peak ratios are more diverse.The diversity of simulated flood peaks, volumes, durations, and hydrographs provides full representativeness for the resulting flood risks, compared with the risk from univariate flood simulation method.

Reservoir Flood Control Operations
The flood control operation of the Meishan reservoir was simulated according to its latest realtime operating rules [62].For a flood control operation, there are two key threshold parameters: flood control storage (noted as Sfc) and downstream safety storage (noted as Sds).
1. Sfc is the expected storage for the whole process during flood season, which implies that the reservoir must draw down to Sfc after a flood inflow peak to vacate storage for future flooding.2. Sds is a critical storage threshold that determines whether the downstream area should suffer: when the storage is less than Sds, the reservoir is capable of retaining more floodwater; once the reservoir's storage exceeds Sds-which means that the dam faces a serious overtopping risk-the reservoir should quickly draw down.As a trade-off, the downstream area will suffer in this case.
The Sfc and Sds for the Meishan reservoir are 1198 million m³ (corresponding to reservoir water level 125.27 m) and 1492 million m³ (corresponding to reservoir water level 133.00 m). Figure 14 illustrates reservoir routing guidance of a flood event with T time steps.

Reservoir Flood Control Operations
The flood control operation of the Meishan reservoir was simulated according to its latest real-time operating rules [62].For a flood control operation, there are two key threshold parameters: flood control storage (noted as S fc ) and downstream safety storage (noted as S ds ).

1.
S fc is the expected storage for the whole process during flood season, which implies that the reservoir must draw down to S fc after a flood inflow peak to vacate storage for future flooding.

2.
S ds is a critical storage threshold that determines whether the downstream area should suffer: when the storage is less than S ds , the reservoir is capable of retaining more floodwater; once the reservoir's storage exceeds S ds -which means that the dam faces a serious overtopping risk-the reservoir should quickly draw down.As a trade-off, the downstream area will suffer in this case.
The S fc and S ds for the Meishan reservoir are 1198 million m 3 (corresponding to reservoir water level 125.27 m) and 1492 million m 3 (corresponding to reservoir water level 133.00 m). Figure 14 illustrates reservoir routing guidance of a flood event with T time steps.

Risk Assessment and Discussion
Based on reservoir operation results, the overall risk considering both upstream and downstream can be derived from Equation (10) and written as Equation ( 16): where n denotes the number of flood events, R i is the risk of flood event i, and R u,i ,R d,i , and p i denote upstream risk, downstream risk, and the flood probability of flood event i. p i is the standardized probability for flood event i, which is the PDF of flood event i divided by the sum of PDFs of n flood events.The risk for each simulated flood event, including the 4 hydrographs for the C or scenario and 12 hydrographs of C and scenario, were calculated.Tables 4 and 5 present key specifications of the risk calculation process, including hydrograph characteristics, reservoir operation results, and risks for each flood event.
As Tables 4 and 5 shows, the overall risks in C or and C and scenarios are 0.64 and 0.31, respectively.The peak and volume values of C or scenario are much greater than that of C and scenario (also shown in Figure 11); as a result, the risk in the C or scenario is more than twice that of the C and scenario.This demonstrated that different flood frequency definitions can produce dramatically diverse flood risk, even under the same frequency and reservoir operating rules.
The impact of peak-volume combinations on flood risk can be indicated in detail.For the "Or" scenario (Table 4), the risks increase with peak inflow.In flood 3 and 4, the peak outflows exceed 1200 m 3 /s that triggered downstream risk.The Meishan reservoir performs very well in flood control: peak reduction reaches 85%, 85%, 70%, and 68% in four flood events, significantly diminishing downstream disasters.For the "And" scenario (Table 5), the highest risk is for neither the maximum peak inflow nor the maximum volume hydrographs, but for the middle hydrograph (flood 6).
Table 5 also reveals that the reservoir decreased flood peak in different cases-the greater the peak inflow, the greater the peak reduction.With reservoir flood control operation, natural flood events with the same frequency produced different flood risk consequences in real world.

Risk Assessment and Discussion
Based on reservoir operation results, the overall risk considering both upstream and downstream can be derived from Equation (10) and written as Equation ( 16): where n denotes the number of flood events, i R is the risk of flood event i, and u,i R , d,i R , and i p denote upstream risk, downstream risk, and the flood probability of flood event i. i p is the standardized probability for flood event i, which is the PDF of flood event i divided by the sum of PDFs of n flood events.The risk for each simulated flood event, including the 4 hydrographs for the Cor scenario and 12 hydrographs of Cand scenario, were calculated.Tables 4 and 5 present key specifications of the risk calculation process, including hydrograph characteristics, reservoir operation results, and risks for each flood event.In order to comprehensively investigate the effect of the copula frequency analysis method, the 0.01 frequency flood simulated with the traditional univariate method was further developed.The univariate frequency method was implemented by scaling the typical flood hydrograph to a fixed ratio to match the 0.01 frequency peak or volume value, which is called the "peak ratio" frequency and the "volume ratio" frequency, respectively.Here, 10 flood hydrographs were simulated for each frequency.The procedures of reservoir operation and risk assessment are the same as those used for the copula method.The overall risks generated from the four flood frequency definitions-the two copula definitions and two univariate definitions-are shown and compared in Figure 15.
Water 2019, 11, 475 18 of 21 ratio to match the 0.01 frequency peak or volume value, which is called the "peak ratio" frequency and the "volume ratio" frequency, respectively.Here, 10 flood hydrographs were simulated for each frequency.The procedures of reservoir operation and risk assessment are the same as those used for the copula method.The overall risks generated from the four flood frequency definitions-the two copula definitions and two univariate definitions-are shown and compared in Figure 15.All four scenarios in Figure 15 represent risks for a "0.01 flood frequency" flood, but their risks are quite diverse.The average risk for each scenario is 0.64 (copula "Or" scenario), 0.31 (copula "And" scenario), 0.65 (univariate peak-scaled scenario), and 0.77 (univariate volume-scaled scenario).Several conclusions can be drawn from the results.
First, compared with the univariate peak ratio results, the volume ratio is higher, indicating that the threaten issue for the risk of the Meishan reservoir is volume characteristic.The risk from the volume ratio is also the highest among the four scenarios, which reveals that the volume ratio method may overestimate the flood risk in this case.On the other hand, the risk resulting from the copula "And" scenario is significantly lower than the risk from the other three scenarios, for the peak and volume values are lower in the copula "And" scenario.All four scenarios in Figure 15 represent risks for a "0.01 flood frequency" flood, but their risks are quite diverse.The average risk for each scenario is 0.64 (copula "Or" scenario), 0.31 (copula "And" scenario), 0.65 (univariate peak-scaled scenario), and 0.77 (univariate volume-scaled scenario).Several conclusions can be drawn from the results.
First, compared with the univariate peak ratio results, the volume ratio is higher, indicating that the threaten issue for the risk of the Meishan reservoir is volume characteristic.The risk from the volume ratio is also the highest among the four scenarios, which reveals that the volume ratio method may overestimate the flood risk in this case.On the other hand, the risk resulting from the copula "And" scenario is significantly lower than the risk from the other three scenarios, for the peak and volume values are lower in the copula "And" scenario.
Second, the risks from the univariate methods, especially the peak ratio results, are much more diverse than copula bivariate results.This diversity, triggered by different typical flood hydrograph "shapes", is currently overlooked in dam planning and real-time flood management.Copula based bivariate methods lowered the diversity, as can be observed in Figure 15.Therefore, the copula frequency analysis can offset the risk evaluation uncertainty, although the hydrograph "shapes" are still diverse.
Third, the copula method extends flood risk assessment capability.For the univariate method, the flood characteristic can only be single (peak or volume) and deterministic, and thus the flood result is also deterministic.In the copula method, the peak-volume combinations can be set considering prediction accuracy ranges instead of uniform sampling, as shown in Figure 11.In this way, a more reliable risk assessment can be obtained by incorporating prediction accuracy.Therefore, the copula method is an effective complement to existing univariate flood frequency analysis, not only for the definition of flood frequency, but also in more versatile and adaptive applications.

Conclusions
Estimation of flood risk is essential for flexible reservoir operation and basin management.Several studies have already been conducted to analyze the flood frequency simulation and reservoir operation risks, but the impacts of these issues on each other were relatively unexplored.This study combined bivariate flood frequency analysis to explicitly simulate its impact on upstream-downstream risk accounting for the effect of reservoir operation.By applying the model to the Meishan reservoir, two definitions of flood frequency, known as copula "Or" and "And" scenarios, were proposed and simulated.Integrated flood risks were carried out in four scenarios with the same frequency: univariate peak-scaled, univariate volume-scaled, bivariate copula "Or", and bivariate copula "And" scenarios.Their differences were comprehensively compared and discussed.
It is concluded that the overall risk from the univariate volume-scaled method is the highest among the four results, which demonstrates that volume is the key factor to the flood risk of Meishan reservoir.Risks from the univariate peak-scaled methods are quite diverse and unstable due to different typical hydrographs, while risks from copula-based results are much more reliable.It is worth remarking in this paper that different flood frequency analysis approaches do have a significant "chain reaction" influence on overall risks to the reservoir basin.Since dam construction standard are determined by hydrological frequency in most countries, this paper provides useful references for recognizing different consequences caused by different frequency definitions.In real-time operation process, the copula-based flood frequency analysis provides a more convincing and adaptive reference for flood risk assessment than the conventional univariate method, especially in the case of imperfect flood prediction.

= 21 Figure 1 .
Figure 1.Location of the Meishan reservoir and its river basin.

Figure 1 .
Figure 1.Location of the Meishan reservoir and its river basin.

Figure 2 .
Figure 2. Bivariate scatter figures and coefficients of historic flood peak, volume, and duration.

Figure 2 .
Figure 2. Bivariate scatter figures and coefficients of historic flood peak, volume, and duration.

Figure 3 .
Figure 3. Schematic diagram of the concepts of copula.

Figure 3 .
Figure 3. Schematic diagram of the concepts of copula.

Figure 4 .
Figure 4. Sketch of the relationships between C(u,v), C or (u,v), and C and (u,v).

Figure 5 .
Figure 5. Illustration of upstream risk Ru, downstream risk Rd, and their integration R.

Figure 5 .
Figure 5. Illustration of upstream risk R u , downstream risk R d, and their integration R.

Figure 6 .
Figure 6.Empirical and kernel density estimation (KDE)-based theoretical distribution for peak and volume.

Figure 6 .
Figure 6.Empirical and kernel density estimation (KDE)-based theoretical distribution for peak and volume.

Figure 9 .
Figure 9. Contour map of C(u,v), C and (u,v), and C or (u,v).

Figure 10 .
Figure 10.Peak-volume combinations on the 0.01 frequency contour line.

Figure 11 .
Figure 11.Peak-volume combinations with a 0.01 frequency in Cand and Cor scenarios.

Figure 10 .
Figure 10.Peak-volume combinations on the 0.01 frequency contour line.

Figure 11 .
Figure 11.Peak-volume combinations with a 0.01 frequency in Cand and Cor scenarios.

Figure 11 .
Figure 11.Peak-volume combinations with a 0.01 frequency in C and and C or scenarios.4.3.2.Flood Hydrograph Simulation Based on Peak-Volume CombinationsTypical flood hydrograph is another important factor in a flood hydrograph simulation, since it determines the "shape" of a hydrograph.In conventional methods, a representative hydrograph

Figure 12 .
Figure 12.Simulated hydrographs of C or scenario.

Figure 13 .
Figure 13.Simulated hydrographs of C and scenario.

Figure 14 .
Figure 14.Flowchart of Meishan reservoir flood control operation routing.
s and reservoir water level < 133 m, reservoir water level is kept constant

Table 1 .
Statistical indicators of historic flood events.

Table 1 .
Statistical indicators of historic flood events.

Table 2 .
Some popular and important copulas.

Table 3 .
Copula parameters and root mean square error (RMSE) for Clayton, Gumbel, and Frank copula types.

Table 3 .
Copula parameters and root mean square error (RMSE) for Clayton, Gumbel, and Frank copula types.

Table 4 .
Risk results of the Meishan reservoir for the "Or" scenario.

Table 5 .
Risk results of the Meishan reservoir for the "And" scenario.