The Shagou reservoir is located upstream of the Shuhe River in the Huaihe River Basin in China, with a control basin area of 164 km², as shown in Figure 1. The average annual precipitation is 745 mm, and flood season (from June to September) accounts for approximately 74% of the annual precipitation. In the current flood control operation, the flood-limited water level of the Shagou reservoir is 231.5 m, and the maximum control outflow is 500 m³/s under the given return period T =20a. In this paper, the Shagou reservoir is taken as the research object to discuss the flood control risk in floodwater utilization considering the uncertainty of flood volume and peak.

The annual maximum 72 h flood volume and corresponding flood peak observed in 1964–2013 are utilized to represent the flood characterization in the Shagou reservoir. The data are provided by the Yi-Shu-Si River Basin Administration, which is responsible for the unified management of major rivers (including the Shuhe River), lakes, hubs and other projects in the Yishusi Basin. The statistical results of flood volume and flood peak are shown in Table 1.

According to the statistics of historical flood data, the highest annual maximum flood volume is 3.115 m³⨯10⁷ m³, which is approximately 3.5 times the multiyear average, and the lowest annual maximum is 70.4 m³⨯10⁴ m³. Similarly, the highest annual maximum flood peak is 1,190 m³/s, which is approximately 3.7 times the multiyear average, and the lowest annual maximum is 3.92 m³/s. The annual flood change is dramatic, and floodwater utilization is necessary for local regional water management. When the range of flood volume is [2,070, 2,204] ⨯10⁴ m³, the flood peak range is [315, 994] m³/s. Obvious uncertainty exists between flood volume and peak, which may have an adverse effect on flood control operation. Estimating the flood control risk caused by the flood volume and peak uncertainty in the Shagou reservoir is important for flood control operation decision making in floodwater utilization management.

Multivariate distribution construction using copulas has been well developed in the past years (Sklar, 1959; Shaked and Joe, 1998; Sancetta and Satchell, 2004). A bivariate joint distribution can be expressed by a copula function and its corresponding marginal distributions. A copula is a function that links two marginal distribution functions to construct a multivariate distribution function. Sklar’s theorem states that if F X Y x , y is a bivariate distribution function of 2 correlated random variables X and Y the respective marginal distributions F X x and F Y y , it is possible to write a cumulative distribution function (CDF) with two single marginal distributions as follows: F X Y x , y = C θ F X x , F Y y (1) where θ is the copula parameter. If these marginal distributions are continuous, a unique copula function C ∙ exists (Sraj et al., 2015).

There are many classes of copula functions, such as Archimedean copulas, elliptical copulas and Plackett copulas (Plackett, 1965; Fang et al., 2002). Archimedean copulas are popular because they can be easily constructed and are capable of capturing a wide range of dependence structures with several desirable properties, such as symmetry and associativity (Nelson, 2006; Hofert, 2008). The widely used bivariate Archimedean family copulas include the Clayton Copula, Frank Copula, and Gumbel-Hougaard (GH) Copula, with a parameter θ , as shown in Table 2. The copula parameter θ is usually estimated by the maximum likelihood method (Strupczewski et al., 2001). Therefore, the joint distribution function based on the copula method can be derived when the marginal distribution functions of variables are determined (Zhang and Singh, 2007a).

In current studies, root mean square error (RMSE) (Zhang and Singh, 2006), Akaike’s information criterion (AIC) (Zhang and Singh, 2007b) and Nash-Sutcliffe efficiency ( E n s ) (Xu et al., 2017) are usually employed to measure the goodness of fit of the joint distribution.

In conventional univariate analysis, the return period is usually used to represent the average time interval of a specific design flood, which is also a method used to measure the magnitude of floods. Within the copula-based framework, various definitions of the joint return period have been proposed, such as OR, AND, Kendall, dynamic, and structure-based return periods (Shiau, 2003; De Michele et al., 2005; Salvadori et al., 2011; Yin et al., 2018a). In this paper, the OR case ( T o r ) is adopted to describe the flood occurrence and can be expressed as follows: T o r = μ 1 − F q , w = μ 1 − C θ F Q q , F W w (2) where μ is the mean interarrival time between two consecutive events (in the case of annual maxima μ =1 year), and F q , w = P q ≥ Q , w ≥ W is depicted by a copula function C θ F Q q , F W w ; F Q q and F W w denote the marginal distribution functions Q and W, respectively.

According to Eq. 2, in the bivariate case, infinite possible combinations of Q and W can be selected for a given joint return period T o r , which have the same joint probability. The different likelihood of each combination must be considered to select appropriate design scenarios (Salvadori et al., 2011; Xiong et al., 2020). Based on Sklar’s theorem, all the possible compositions of flood volume and peak differ in terms of their probability of occurrence, which can be measured by the value of the joint PDF (Salvadori et al., 2011; Gräler et al., 2013; Guo et al., 2018). The most likely flood event q * t , w * t of all possible events at a given joint return period T o r can be obtained by the following formula: q * t , w * t = arg max ⁡ f q , w = c θ F Q q , F W w ∙ f Q q ∙ f W w C θ F Q q , F W w = 1 − 1 T o r (3) where f q , w is a nonstationary joint PDF of q and w; C ∙ is the density function of the copula for nonstationary data series; and f Q q and f W w are the marginal PDFs.

Generally, the risk is simplified and defined as the probability of occurrence of a risk event (Sun et al., 2018; Li et al., 2022). In this paper, the flood control risk for the reservoir is defined as the probability of the highest level over the design flood level and can be expressed as follows: R = P H m > H d = m n (4) where H m is the highest reservoir level in the flood control operation, H d is the design flood level of the reservoir, m is the number of times that the highest water level exceeds the design water level in stochastic simulations, and n is the total number of stochastic simulations.

In this research, the framework of estimating flood control risk with the flood uncertainty of flood volume and peak is proposed based on the Monte Carlo sampling method and the most likely event selection method. It can be divided into three steps, which are shown in Figure 2.

Step 1: Establish the joint distribution function. Based on the observed flood volume and peak data with k-samples, the marginal distributions F X x and F Y y and the probability density functions f X x and f Y y are established, respectively. The copula functions can be estimated by using the historical flood series. According to the goodness of fit of the joint distribution, the optimal copula function is selected to construct the joint distribution function of flood volume and peak.

Step 2: Stochastic simulation of the flood. Based on the multilevel Monte Carlo method (Giles, 2008; Brodie, 2013; Clare et al., 2022), n sets of bivariate data with k-samples are randomly generated. Based on the generated flood series, the joint distribution function is constructed for each set. The most likely method is employed to select the appropriate q , w under the given T o r .

Step 3: Calculate the risk. According to the engineering and hydrological characteristics, various design floodwater utilization schemes have been developed. The design flood hydrograph can be obtained by n sets of q , w and typical floods. Taking n sets of flood hydrographs as input data, the simulation of flood control operation is carried out, and the flood control risk can be calculated by Eq. 4.

In this paper, the Pearson type III (P-III) distribution, which is recommended by the Chinese Ministry of Water Resources (Wang et al., 2015; Peng et al., 2017; Gao et al., 2018), was employed to obtain the single marginal distribution of flood volume and flood peak. The parameters of the P-III distribution were estimated by the moment method (Hosking, 1990) and are shown in Table 3. The K-S test ( D n ) is employed to describe how well the distributions fit the flood data. With K-S’s critical value D 0.05 = 1.36 50 = 0.1923 , the D n of flood volume and peak are both 0.0980. Therefore, flood volume data and flood peak data both failed to reject the P-III distribution.

For the flood data, Clayton, Frank and G-H Copula were used to establish the joint distribution. RMSE, AIC and E n s were obtained to measure the goodness of joint distribution as shown in Table 4. The empirical joint probabilities were estimated for the flood data by using the Gringorten formula (Gringorten, 1963; Du et al., 2019). The smaller the RMSE and AIC are, the better the joint distribution is. The larger E n s is, the better the joint distribution is. Based on the above measurement principles, the results indicated that the best-fitted joint distribution is the G-H copula function with the parameter 1.8890 in the Shagou reservoir flood volume and peak. Then, the structure of the joint distribution function of the historical flood volume and peak was determined.

According to the bivariate joint distribution of historical flood data, the most likely method was employed to obtain the bivariate flood volume and peak for a given joint return period T o r . The comparison of design flood values under the bivariate joint return period and univariate return period is shown in Table 5.

The results show that the flood volume and peak obtained by the bivariate joint return for a given joint return period T o r are larger than the values obtained by the univariate return period. As shown in Table 5, the design flood volume for the univariate return period is 7.59% more than that for the bivariate joint return period, and the design flood peak for the univariate return period is 8.22% more than that for the bivariate joint return period. The value of flood volume and peak for joint return period T o r = 100 a will make a larger return period in the respective univariate frequency analysis, which is nearly 144 a. Considering the correlation between flood volume and peak, multivariable flood events can be described more reasonably and may demonstrate a new theoretical basis for flood control operations.

Based on the joint distribution function of flood volume and peak, the multilevel Monte Carlo method was employed to generate bivariate design data. Considering the number of historical flood data, 50 samples are randomly generated in each set. In the same way, the joint distribution was constructed, and the most likely design values were obtained. The number of stochastic simulation sets is 10,000 in this paper and the design flood values based on the stochastically simulated flood for given T o r are shown in Figure 3. The statistics of bivariate design data under different bivariate joint return periods are shown in Figure 4.

The results reveal that the generated design flood values are scattered around the most likely design value calculated based on historical flood data. As the bivariate joint return period increases, the generated design values also increase. For the T o r = 20 a , the range of flood volume and flood peak is [1,412, 4,476] ⨯10⁴ m³ and [509, 1,692] m³/s, respectively. And for the T o r = 200 a the range of flood volume and flood peak is [1990, 7,036] ⨯10⁴ m³ and [721, 2,788] m³/s, respectively. The Shagou reservoir is designed with the univariate return period T=100a, and the joint return period of the simulated joint design value spans from 10a to 1000a compared to the historical data under the same bivariate joint return period T o r = 100 a . The joint design value is related to the economy of engineering construction and the reliability of flood control operation for reservoirs. If a small joint design value is adopted, the scale of the reservoir will be small, which is not conducive to providing its own benefits and may lead to a huge loss of life and property in flood control. The uncertainty of flood volume and peak should be given more attention in the field of floodwater utilization.

As mentioned in Section 2.1, the Shagou reservoir was designed with the standard T=100a, and the corresponding design flood water level is 233.13 m, which was obtained by univariate flood volume frequency analysis. In the current scheduling rules, the water level is limited to 231.5 m in the flood season. The essence of floodwater utilization is transferring more floodwater into ordinary water resources in the flood season for use in the non-flood season, where the flood-limited water level plays an important role. In recent years, numerous studies have been carried out to scientifically raise the flood-limited water level in the flood season without decreasing flood control standards or damaging the ecological environment of rivers (Li et al., 2010; Liu et al., 2015; Chang et al., 2017; Xie et al., 2018; Ye et al., 2019; Wang et al., 2020). In this paper, according to the current operation, the flood-limited water level is designed to be from 231.5 m to 232.0 m with an interval of 0.1 m to indicate various floodwater utilization schemes.

The flood in August 1974 was selected as a typical flood, and various flood hydrographs were obtained by using the design bivariate design value in Section 4.4. By simulating the current flood control operation, the maximum water level ( H s ), the count of the maximum water level exceeding the current design water level ( m s ), and the flood control risk ( R s ) were calculated under different flood-limited water levels ( H c ), which are shown in Table 6. The statistical results of the maximum water level are shown in Figure 5 under different flood-limited water levels.

In the flood control operation of a reservoir, the flood-limited water level is considered to be the initial water level in a flooding process. With the initial water level rise, the maximum water level in the flood control operation process will also increase accordingly. According to the statistical results, the averages of the maximum water levels are all less than the design water level of 233.13 m under flood utilization with flood-limited water levels of 231.5, 231.6, 231.7 and 231.8 m. With the increase in flood-limited water level, the counts of the maximum water level exceeding the current design water level similarly increased. When the flood-limited water level increased from 231.5 m to 232, H c and R s nearly doubled. The flood control risk increased with the rising flood-limited water level, and the benefit of floodwater utilization grew correspondingly. Assuming that the water level of the reservoir can be kept at the flood-limited water level when the flood season ends, the benefit growth rate and risk growth rate can be calculated quantitatively compared to the current operation and are shown in Figure 6.

The results reveal that different increased rates occurred between benefit and risk with the same increased flood-limited water level. In particular, the growth rates of benefit and risk are 0.08% and 91% with H c =232.0 compared to the current operation, respectively. Achieving a lower benefit increment brings a greater risk increment to water management, which shows that the marginal benefit may be contained in floodwater utilization.