Multivariate Dam-Site Flood Frequency Analysis of the Three Gorges Reservoir Considering Future Reservoir Regulation and Precipitation

Abstract

Under a changing environment, the current hydrological design values derived from historical flood data for the Three Gorges Reservoir (TGR) might be no longer applicable due to the newly-built reservoirs upstream from the TGR and the changes in climatic conditions. In this study, we perform a multivariate dam-site flood frequency analysis for the TGR considering future reservoir regulation and summer precipitation. The Xinanjiang model and Muskingum routing method are used to reconstruct the dam-site flood variables during the operation period of the TGR. Then the distributions of the dam-site flood peak and flood volumes with durations of 3, 7, 15, and 30 days are built by Pearson type III (PIII) distribution with time-varying parameters, which are expressed as functions of both reservoir index and summer precipitation anomaly (SPA). The multivariate joint distribution of the dam-site flood variables is constructed by a 5-D C-vine copula. Finally, by using the criteria of annual average reliability (AAR) associated with the exceedance probabilities of OR, AND and Kendall, we derive the multivariate dam-site design floods for the TGR from the predicted flood distributions during the future operation period of the reservoir. The results indicate that the mean values of all flood variables are positively linked to SPA and negatively linked to RI. In the future, the flood mean values are predicted to present a dramatic decrease due to the regulation of the reservoirs upstream from the TGR. As the result, the design dam-site floods in the future will be smaller than those derived from historical flood distributions. This finding indicates that the TGR would have smaller flood risk in the future.

1. Introduction

As the largest hydroelectric dam in the world, the Three Gorges Dam (TGD) is located at the dividing site between the upstream and midstream of the Yangtze River, which is the largest river in China. In recent decades, the Three Gorges Reservoir (TGR) formed by the TGD has been put into operation and led to a profound impact on the hydrological regime downstream of the dam [1,2,3,4,5,6]. Since flood control is the primary goal of the TGD, the flood magnitude downstream of the dam has been found to exhibit a significant decline in recent years, especially after 2009, when the TGR began to be fully operated [7,8].

In addition to the TGD, dozens of dams with sizable reservoir capacity have been constructed or are under construction in the Upper Yangtze basin upstream of the TGD. It is indicated that the inflow flood or the dam-site flood of the TGR is unavoidably regulated by the upstream reservoirs. Due to the regulation of the reservoirs in the upstream basin of the TGR, the dam-site floods of the TGR no longer follow a natural process [8,9]. Some studies have revealed that the hydrological regime of the inflow of the TGR has presented significant variations due to the regulation of the reservoirs upstream of the TGR [9,10]. The impacts of the TGR on the downstream hydrological regime have been widely assessed by a number of studies [8,9,10,11,12]. However, very few studies have paid attention to the flood frequency of the dam-site flood of the TGR [9].

From the perspective of hydrological processes, the dam-site flood frequency of the TGR is dominated by the climatic conditions of the basin, including precipitation and temperature [7,11]. Hence, changes in some climatic conditions (such as precipitation) are also potential factors altering the flood frequency of the Upper Yangtze basin. The most relevant studies reveal that the projected temperature of this basin will present significant increase, varying from 1 °C to 4 °C, in the future [7,8,9,10,11,12]. Although different emission scenarios or climate models result in different trends for projected future precipitation, the hydrological regime, such as annual average discharge, seasonal high flow, and peak discharge of the Upper Yangtze River, will probably shift in the future [12]. In the design stage as well as the current operation stage of the TGR, the flood frequency was estimated from the historical flood records, which were not yet influenced either by reservoirs or changes of climatic conditions, and might be unable to stand for the occurrence behavior of future floods. Hence, the corresponding hydrological design of the TGR would fail to deal with future flood risk.

The current hydrological design of the TGR is based on the stationarity assumption and assessed by the design level of the return period, which is usually calculated as the inverse of exceedance probability. Under nonstationary conditions, the exceedance probability of a given flood event varies with time, and thus the return period is unsuitable as the criterion for hydrological design. Due to the failure of the traditional return-period-based method, new measurements have been introduced into hydrological design. To make the measurement of return period apply to nonstationary conditions, the concept of return period has been redefined as the expected waiting time (EWT) for an exceedance to occur [13], and the time period that results in the expected number of exceedance (ENE) equal to 1 over this period [14]. In addition to the redefinition of return period, measurement based on reliability was also employed as the hydrological design criterion for nonstationary conditions [15,16,17,18]. The reliability-based method, which concerns flood distribution during the whole design life of hydraulic structures, is more reasonable for engineering practice [15,16].

Due to a large catchment area, the flood processes of the Upper Yangtze River generally have a long duration and require multiple feature variables to describe all the characteristics of flood events, including flood peak and flood volumes with different durations. Compared to the univariate hydrological design that ignores the dependence of different flood features, multivariate flood frequency can provide a more complete understanding of flood events and should be more suitable for hydraulic structures [17,18,19]. The copula method appears to be the most popular and powerful tool in modeling the dependence of multivariate flood variables [17,18,19,20,21,22,23,24]. Li et al. [24] employed the copula method to construct the bivariate joint distribution of the dam-site flood variables of the TGR. However, using only two flood features might fail to capture the comprehensive characteristics of the long-duration flood processes.

As mentioned above, the current hydrological design, which relies on historical flood records, could be unsuitable for the operation of the TGR in the future. As more and more reservoirs in the Upper Yangtze basin will be put into operation and the climatic conditions of the basin might change, it seems almost inevitable that the historical flood frequency would fail to accurately reflect the future flood risk of the TGR. In this study, we focus on the multivariate dam-site flood frequency of the TGR considering reservoir regulation and change of climatic conditions. Pearson type III (PIII) distribution [24] with time-varying parameters is employed to model the effect of reservoir regulation and climate change on the dam-site flood peak and four flood volumes with the durations of 3, 7, 15 and 30 days. Then, the multivariate joint distribution of the flood variables is constructed by a five-dimension (5-D) C-vine copula [25]. The multivariate flood distributions in the future period from 2021 to 2100 are predicted by considering the reservoirs under construction and the summer precipitation anomaly (SPA) under three typical emission scenarios. Finally, we derive the multivariate hydrological design values from the predicted multivariate flood distributions of flood peak and flood volumes using the reliability-based method.

The remainder of this paper is organized as follows. Section 2 gives a brief introduction of the study region and data set used in this study. Section 3 presents the main methodologies. Section 4 gives the results of this study. The final section displays the conclusions and discussion.

2. Study Region and Dataset

As the largest river in China, the Yangtze River has a drainage area of 1.8 million km² and the river channel length is more than 6300 km. In this study, we focus on the Upper Yangtze basin, which is actually above the Yichang hydrological gauging site and covers a catchment area of about 1 million km². Before the operation of the TGR, the dam-site flood of this reservoir can be represented by the flood observed at Yichang gauging site, given that the TGD is located approximately 40 km upstream from the Yichang gauging site and the catchment area between the TGD and the gauging site is negligible compared to the whole catchment area of the Upper Yangtze basin.

We have collected the observed daily discharge of the Yichang gauging site during the period from 1882 to 2020 as well as the historical flood records in 1560, 1613, 1788, 1796, 1860, and 1870, which were estimated by the Changjiang Water Resources Commission (CWRC) when designing the TGD [24,26]. To reconstruct the dam-site flood during the operation period of the TGR, we also collected the daily discharge observed at Cuntan and Wulong gauging sites, which control the majority of the discharge into the TGR (see Figure 1), and the daily precipitation and potential evaporation data of the meteorological stations in the catchment between Yichang gauging site and its upstream gauging sites of Cuntan and Wulong.

Table 1. Information on the reservoirs upstream the Three Gorges Reservoir.

Reservoir	Catchment Area (106 km2)	Flood control Capacity (109 m3)	Completation Year
Bikou	0.026	0.083	1976
Ertan	0.1164	0.9	1998
Baozhusi	0.0284	0.28	1998
Zipingpu	0.0227	0.167	2006
Pubugou	0.0685	0.73	2008
Silin	0.0486	0.184	2009
Pengshui	0.069	0.232	2009
Goupitan	0.0433	0.2	2011
Shatuo	0.0545	0.209	2011
Ahai	0.2354	0.215	2012
Jinanqiao	0.2374	0.158	2012
Longkaikou	0.24	0.126	2013
Ludila	0.2473	0.564	2013
Xiangjiaba	0.4588	0.903	2013
Jinping I	0.1026	1.6	2014
Xiluodu	0.4544	4.651	2014
Tingzikou	0.0611	1.44	2014
Liyuan	0.22	0.173	2015
Guanyinyan	0.2565	0.253	2015
Caojie	0.1561	0.199	2015
Changheba	0.0559	0.12	2017
Wudongde	0.4061	2.44	2021
Baihetan	0.4303	7.5	2022
Lianghekou	0.0596	0.2	2023
Shuangjiangkou	0.0393	0.663	2024

summarizes the detailed information on the reservoirs, with considerable flood control capacity in the catchment upstream the TGR. By the end of 2020, 21 reservoirs have been completed, and five reservoirs will be completed in the future. The reservoir index (RI), which is calculated as the cumulative value of the product of flood control capacity and catchment area of each reservoir [27], is used to quantify the effect of reservoirs on the dam-site flood distributions of the TGR.

The summer precipitation anomaly (SPA), which is defined as the departure of the mean value of summer (May–September) precipitation from the long-term mean value during 1961~1990, is used as a candidate covariate for the dam-site flood distributions of the TGR. To utilize the long-term historical precipitation, the SPA of the Upper Yangtze basin during the period from 1470 to 1990 is obtained from the precipitation reconstruction data set developed by Shi et al. [28]. The SPA data during the observation period from 1991 to 2020 are calculated from the precipitation observed at the meteorological stations in the Upper Yangtze basin. The SPA data during the future period from 2021 to 2100 is calculated from the mean value of the summer precipitation projected by 6 GCMs (BCC-CSM2-MR, CanESM5, CNRM-CM6-1, FGOALS-g3, IPSL-CM6A-LR, MIROC6) within the Coupled Model Intercomparison Project Phase 6 (CMIP6), which can reduce the uncertainty caused by the variability of outputs from different GCMs. For each GCM, we also consider three Shared Socioeconomic Pathways (i.e., SSP1-2.6, SSP2-4.5, and SSP5-8.5) to represent low-, medium-, and high-emission scenarios, respectively [29].

3. Methods

3.1. Reconstruction of the Dam-Site Floods of the TGR

During the period from 1882 to 2002 before the impoundment of the TGR, the dam-site floods of the TGR are represented by those at Yichang hydrological gauging site, which is closely located at the downstream of the TGD. After 2003, the observed floods at Yichang gauging site have been regulated by the TGR and hence cannot represent the dam-site floods of the TGR. In this study, the dam-site floods of the TGR during the post-TGD period are reconstructed using hydrological models, which are calibrated from the hydrological data during the period without TGD. As shown in Figure 1, the dam-site discharge of the TGR consists of the discharges routed from the upstream hydrological gauging sites of Cuntan and Wulong as well as the runoff generated in the TGR interval watershed, i.e., the area enclosed between these two gauges and Yichang. Hence the dam-site discharge of the TGR can be expressed as follows:

$$Q_{\mathrm{TGR}}=Q_C+Q_W+Q_{IW}$$

(1)

where Q_TGR is the daily dam-site discharge of the TGR; Q_C and Q_W are the daily discharges routed from Cuntan and Wulong gauges, respectively; and Q_IW is the daily runoff generated in the TGR interval watershed. By using Equation (1), the dam-site floods of the TGR for the period from 2003 to 2020, i.e., the period after the operation of the TGR, can be extracted from Q_TGR.

The Muskingum method [30] is used to simulate the routing processes of the discharges from both Cuntan and Wulong to Yichang, and the Xin’anjiang (XAJ) model [31] is employed to simulate the runoff generated in the TGR interval watershed from the daily precipitation and potential evaporation. The model parameters are calibrated from the observed hydro-meteorological data during the period from 1955 to 2002 by using the SCEUA algorithm with the objective function of the Nash efficiency coefficient (NSE) [32].

3.2. Time-Varying Multivariate Flood Distribution

In this study, we focus on the multivariate dam-site flood variables of the TGR, including the annual maximum daily discharge (Q₁), annual maximum 3-day flood volume ($V_3$), annual maximum 7-day flood volume ($V_7$), annual maximum 15-day flood volume ($V_{15}$), and annual maximum 30-day flood volume ($V_{30}$). According to Sklar’s Theorem [33], the probability distribution of the five-dimensional (5-D) flood series $\left(Q_1,V_3,V_7,V_{15},V_{30}\right)$ at time t can be expressed by a copula $C(\cdot )$ as follows:

$$\begin{array}{l}F\left(Q_{1,t},V_{3,t},V_{7,t},V_{15,t},V_{30,t}|{\mathsf{\theta }}_t\right)\\ =C\left[F_1\left(Q_{1,t}|{\mathsf{\theta }}_{1,t}\right),F_3\left(V_{3,t}|{\mathsf{\theta }}_{3,t}\right),F_7\left(V_{7,t}|{\mathsf{\theta }}_{7,t}\right),F_{15}\left(V_{15,t}|{\mathsf{\theta }}_{15,t}\right),F_{30}\left(V_{30,t}|{\mathsf{\theta }}_{30,t}\right)|{\mathsf{\theta }}_{c,t}\right]\\ =C\left[u_{1,t},u_{3,t},u_{7,t},u_{15,t},u_{30,t}|{\mathsf{\theta }}_{c,t}\right]\end{array}$$

(2)

where $F_1\left(Q_{1,t}|{\mathsf{\theta }}_{1,t}\right)$, $F_3\left(V_{3,t}|{\mathsf{\theta }}_{3,t}\right)$, $F_7\left(V_{7,t}|{\mathsf{\theta }}_{7,t}\right)$, $F_{15}\left(V_{15,t}|{\mathsf{\theta }}_{15,t}\right)$ and $F_{30}\left(V_{30,t}|{\mathsf{\theta }}_{30,t}\right)$ denote the marginal distributions for the univariate flood variables $Q_1$, $V_3$, $V_7$, $V_{15}$ and $V_{30}$ respectively; $u_{1,t}$, $u_{3,t}$, $u_{7,t}$, $u_{15,t}$ and $u_{30,t}$ are the corresponding marginal probabilities; ${\mathsf{\theta }}_{1,t}$, ${\mathsf{\theta }}_{3,t}$, ${\mathsf{\theta }}_{7,t}$, ${\mathsf{\theta }}_{15,t}$ and ${\mathsf{\theta }}_{30,t}$ are the corresponding marginal distribution parameter vectors; and ${\mathsf{\theta }}_{c,t}$ denotes the copula parameter vector. Hence, the parameter vector ${\mathsf{\theta }}_t$ of the multivariate flood distributions consists of both the marginal distribution parameters and the copula parameters, i.e., ${\mathsf{\theta }}_t=\left({\mathsf{\theta }}_{1,t},{\mathsf{\theta }}_{3,t},{\mathsf{\theta }}_{7,t},{\mathsf{\theta }}_{15,t},{\mathsf{\theta }}_{30,t},{\mathsf{\theta }}_{c,t}\right)$. It is theoretically possible that all parameters in ${\mathsf{\theta }}_t$ can be time-varying, but, in this study, we have to ignore the possible non-stationarity in the copula parameters ${\mathsf{\theta }}_{c,t}$ to avoid the large uncertainty in the estimation of ${\mathsf{\theta }}_{c,t}$.

With respect to Equation (2), the probability density function of $\left(Q_1,V_3,V_7,V_{15},V_{30}\right)$ is given by:

$$\begin{array}{l}f\left(Q_{1,t},V_{3,t},V_{7,t},V_{15,t},V_{30,t}|{\mathsf{\theta }}_t\right)\\ =c\left[F_1\left(Q_{1,t}|{\mathsf{\theta }}_{1,t}\right),F_3\left(V_{3,t}|{\mathsf{\theta }}_{3,t}\right),F_7\left(V_{7,t}|{\mathsf{\theta }}_{7,t}\right),F_{15}\left(V_{15,t}|{\mathsf{\theta }}_{15,t}\right),F_{30}\left(V_{30,t}|{\mathsf{\theta }}_{30,t}\right)|{\mathsf{\theta }}_{c,t}\right]\cdot \\ f_1\left(Q_{1,t}|{\mathsf{\theta }}_{1,t}\right)\cdot f_3\left(V_{3,t}|{\mathsf{\theta }}_{3,t}\right)\cdot f_7\left(V_{7,t}|{\mathsf{\theta }}_{7,t}\right)\cdot f_{15}\left(V_{15,t}|{\mathsf{\theta }}_{15,t}\right)\cdot f_{30}\left(V_{30,t}|{\mathsf{\theta }}_{30,t}\right)\\ =c\left(u_{1,t},u_{3,t},u_{7,t},u_{15,t},u_{30,t}|{\mathsf{\theta }}_{c,t}\right)\cdot \\ f_1\left(Q_{1,t}|{\mathsf{\theta }}_{1,t}\right)\cdot f_3\left(V_{3,t}|{\mathsf{\theta }}_{3,t}\right)\cdot f_7\left(V_{7,t}|{\mathsf{\theta }}_{7,t}\right)\cdot f_{15}\left(V_{15,t}|{\mathsf{\theta }}_{15,t}\right)\cdot f_{30}\left(V_{30,t}|{\mathsf{\theta }}_{30,t}\right)\end{array}$$

(3)

where $f_1\left(Q_{1,t}|{\mathsf{\theta }}_{1,t}\right)$, $f_3\left(V_{3,t}|{\mathsf{\theta }}_{3,t}\right)$, $f_7\left(V_{7,t}|{\mathsf{\theta }}_{7,t}\right)$, $f_{15}\left(V_{15,t}|{\mathsf{\theta }}_{15,t}\right)$ and $f_{30}\left(V_{30,t}|{\mathsf{\theta }}_{30,t}\right)$ are the density functions of the marginal distributions for $Q_1$, $V_3$, $V_7$, $V_{15}$ and $V_{30}$, respectively, and $c(\cdot )$ stands for the density function of $C(\cdot )$.

The marginal distributions are characterized by the Pearson type III (PIII) distribution, which is recommended by the Design Criterion of Reservoir Management of China [24,34]. The expression of the PIII distribution is given below:

$$\begin{array}{l}f(x;\mu ,\sigma ,\nu )=\frac{1}{\sigma \left|\mu \nu \right|\Gamma \left(\frac{1}{{\nu }^2}\right)}{\left(\frac{x-\mu }{\mu \sigma \nu }+\frac{1}{{\nu }^2}\right)}^{\frac{1}{{\nu }^2}-1}{e}^{-\left(\frac{x-\mu }{\mu \sigma \nu }+\frac{1}{{\nu }^2}\right)}\\ \sigma >0,\nu \ne 0, \hspace{0.17em}\frac{x-\mu }{\mu \sigma \nu }+\frac{1}{{\nu }^2}\ge 0\end{array}$$

(4)

where $\mu$, $\sigma$ and $\nu$ are location, scale and shape parameters, respectively, referring to the mean value, coefficient of variation (Cv) and coefficient of skew (Cs) of the distribution. In this study, each of the distribution parameters is expressed as an additive function of the covariates of SPA and RI. Taking the location parameter $\mu$ as an example, it is expressed as follows:

$${\mu }_t={\alpha }_0+{\alpha }_1{g}_1\left(SP{A}_t\right)+{\alpha }_2{g}_2\left(R{I}_t\right)$$

(5)

where ${\alpha }_0$, ${\alpha }_1$ and ${\alpha }_2$ are model parameters, and are estimated using the maximum likelihood estimate (MLE) method [35]. $g_1$ and $g_2$ denote transforming functions, which are used to capture linear or nonlinear relationship between $\mu$ and the covariates. The candidates of $g_1$ and $g_2$ include identity, logarithmic and exponential functions. The selection of covariates and transforming functions are performed using the Akaike Information Criterion (AIC) [36]. The goodness-of-fit (GoF) of the PIII distribution is examined by the Kolmogorov-Smirnov test with a significance level of 0.05 [37].

3.3. Construction of Multivariate Flood Distribution Using Vine Copula

According to Equation (3), the multivariate distribution of $\left(Q_1,V_3,V_7,V_{15},V_{30}\right)$ consists of five marginal distributions and a 5-D copula density function $c\left(u_{1,t},u_{3,t},u_{7,t},u_{15,t},u_{30,t}|{\mathsf{\theta }}_{c,t}\right)$. The 5-D copula density function is constructed by the pair copula method, which provides a flexible way to decompose a high dimension distribution into a series of bivariate copulas [25]. In this study, the C-vine copula with the key variable of $Q_1$ is employed to construct the dependence structure of $\left(Q_1,V_3,V_7,V_{15},V_{30}\right)$, since the flood peak $Q_1$ is usually regarded as the dominant feature for a flood event [34]. Figure 2 presents the schematic decomposition of the 5-D C-vine copula, which consists of four trees with a total 10 bivariate copulas. For more detail of the C-vine copula, readers are referred to Aas et al. [25].

The bivariate Gumbel-Hougaard (GH) copula, which inherently accounts for the upper tail dependence, is well-suited to the dependence structure of multivariate flood distribution [38]. Hence, in this study, the GH copula is employed to construct the 5-D C-vine copula. The expression of the bivariate GH copula is given as follows:

$$C\left(u,v\right)=\exp \left\{-{\left[{\left(-\ln u\right)}^{{\theta }_c}+{\left(-\ln v\right)}^{{\theta }_c}\right]}^{1/{\theta }_c}\right\}, {\theta }_c\in [1,\infty )$$

(6)

where $u$ and $v$ denote the bivariate marginal probabilities, and ${\theta }_c$ is the single parameter measuring the dependence strength and estimated by the MLE method. The GoF of the C-vine copula is examined by the Probability Integral Transform (PIT) test [25].

3.4. Multivariate Hydrological Design under Nonstationary Conditions

Under nonstationary conditions, the traditional hydrological design criteria of return period could be invalid. In this study, the criteria of average annual reliability (AAR) [18] is employed as the metric of hydrological design, and is calculated by:

$$AAR\left(q_1,v_3,v_7,v_{15},v_{30}\right)=\frac{1}{T_2-T_1+1}{\displaystyle \sum _{t=T_1}^{T_2}\left(1-p_t\right)}$$

(7)

where $\left(q_1,v_3,v_7,v_{15},v_{30}\right)$ denotes a multivariate flood event; $p_t$ is the exceedance probability at time t; and $T_1$ and $T_2$ are the beginning and end times of the design period of the hydraulic structure, which are 2021 and 2100, respectively.

For a multivariate flood event, there are multiple definitions for an exceedance probability. In this study, the exceedance probabilities in the OR case, AND case and Kendall are used to derive the multivariate hydrological designs. The OR case is that at least one of the flood features in $\left(Q_1,V_3,V_7,V_{15},V_{30}\right)$ exceeds the prescribed threshold. The AND case indicates that all flood features in $\left(Q_1,V_3,V_7,V_{15},V_{30}\right)$ exceed the prescribed threshold. The Kendall case is that the univariate representation calculated from $\left(Q_1,V_3,V_7,V_{15},V_{30}\right)$ via the Kendall’s distribution function exceeds the prescribed threshold.

The exceedance probability in the OR case at time t is denoted as $p_t^{or}$ and calculated by:

$$\begin{array}{l}{p}_t^{or}=P\left(Q_{1,t}\ge q_1\vee V_{3,t}\ge v_3\vee V_{7,t}\ge v_7\vee V_{15,t}\ge v_{15}\vee V_{30,t}\ge v_{30}\right)\\ =1-F\left(q_1,v_3,v_7,v_{15},v_{30}|{\mathsf{\theta }}_t\right)\end{array}$$

(8)

The exceedance probability in the AND case at time t is denoted as $p_t^{and}$ and calculated by:

$$p_t^{and}=P\left(Q_{1,t}>q_1\wedge V_{3,t}>v_3\wedge V_{7,t}>v_7\wedge V_{15,t}>v_{15}\wedge V_{30,t}>v_{30}\right)$$

(9)

The exceedance probability in the Kendall case at time $t$ is denoted as $p_t^{ken}$ and calculated by:

$$p_t^{ken}=1-K_C\left({\rho }_t\right)$$

(10)

where $K_C\left({\rho }_t\right)$ is the Kendall’s distribution function and calculated by

$$\begin{array}{l}{K}_C\left({\rho }_t\right)=P\left[C\left(u_{1,t},u_{3,t},u_{7,t},u_{15,t},u_{30,t}|{\mathsf{\theta }}_{c,t}\right)\le {\rho }_t\right]\\ =P\left[F\left(Q_{1,t},V_{3,t},V_{7,t},V_{15,t},V_{30,t}|{\mathsf{\theta }}_t\right)\le {\rho }_t\right]\end{array}$$

(11)

For multivariate hydrological design, a given ARR would correspond to countless groups of design events, which have different likelihood of occurrence. For the given design event $\left(q_1,v_3,v_7,v_{15},v_{30}\right)$, the multivariate joint probability density during the period from T₁ to T₂ can be calculated by:

$$g\left(q_1,v_3,v_7,v_{15},v_{30}\right)=\frac{1}{T_2-T_1+1}{\displaystyle \sum _{t=T_1}^{T_2}f\left(q_1,v_3,v_7,v_{15},v_{30}|{\mathsf{\theta }}_t\right)}$$

(12)

For the given AAR, the most-likely design event with the maximum probability density is used as the representative of the hydrological design values [18]. The most-likely design event has no analytical solution but can be approximately estimated through the Monte Carlo simulation. In particular, we can generate a large number of design events using the Monte Carlo simulation method and then chose the design event with the maximum probability density among all generated design events.

In China, the design flood hydrograph for a hydraulic structure is usually derived from the design flood value benchmarked against a typical flood hydrograph, which is chosen from observed flood processes [34]. In this study, the design flood hydrograph is characterized by the feature variables of flood peak and flood volumes with the durations of 3, 7, 15 and 30 days. For the given design event $\left(q_1,v_3,v_7,v_{15},v_{30}\right)$, the design flood hydrograph is derived by multiplying the typical flood hydrograph by different amplifiers, which is calculated by the ratio of the flood features of the given design event to those of the typical flood hydrograph. For details of calculating the amplifiers, readers can refer to [18].

4. Results and Analysis

4.1. Dam-Site Flood Series

Figure 3 compares the observed floods at Yichang gauging site and the simulated dam-site floods of the TGR for the period from 1955 to 2020. During the period before the impoundment of the TGR in 2003, the simulated dam-site floods suggest a good performance in fitting the observed floods. After 2003, the dam-site floods of the TGR cannot be represented by the observed floods at Yichang gauging site, which is 40 km downstream of the TGR, due to the regulation of the TGR. Thus, after 2003, the dam-site floods of the TGR will be represented by the floods simulated according to Equation (1). It is found that after 2003 the observed floods at Yichang gauging site are generally smaller than the simulated dam-site floods. This finding reveals the great role of the TGR in dampening the downstream flood magnitude, especially after 2009, when the TGR began to be fully operated. For example, the simulated dam-site flood peak of the TGR in 2020 was 71,828 m³/s, while after the regulation of the TGR, the observed flood peak at Yichang gauging site downstream the dam was 50,831 m³/s, indicating a decline of more than 20,000 m³/s.

4.2. Multivariate Dam-Site Flood Distribution of the TGR

The dam-site flood distributions of the TGR for the period from 1882 to 2020 are estimated by the PIII distribution with time-varying parameters (Equations (4) and (5)) and the results are presented in

Table 2. Results of the covariate analysis for the dam-site floods of the TGR.

Flood Variable	Distribution Parameters			p-Value of KS Test
	$\mathit{\mu }$	$\mathit{\sigma }$	$\mathit{\nu }$
Q1 (m3/s)	$50215.5+15747.3SPA-6029.2\ln \left(RI+1\right)$	0.172	0.454	0.200
V3 (109 m3)	$12.51+3.96SPA-1.42\ln \left(RI+1\right)$	0.174	0.456	0.474
V7 (109 m3)	$26.39+7.87SPA-2.66\ln \left(RI+1\right)$	0.175	0.456	0.407
V15 (109 m3)	$49.83+14.72SPA-4.84\ln \left(RI+1\right)$	0.162	0.364	0.758
V30 (109 m3)	$89.36+29.61SPA-7.34\ln \left(RI+1\right)$	0.146	0.326	0.770

. It is found that these flood variables have similar modeling results in that the location parameters are linear functions of both SPA and the logarithm of RI. The KS test also indicates that the estimated flood distributions have satisfactory fitting quality. The flood variables suggest a positive relationship with SPA as well as a negative relationship with RI. By 2020, the reservoir regulation has induced a decline of about 8900 m³/s in the mean value of the dam-site flood peak of the TGR.

In this study, we consider three different emission scenarios (SSP1-2.6, SSP2-4.5 and SSP5-8.5) to project future summer precipitation during the period from 2021 to 2100. Figure 4 displays the candidate explanatory variables of both SPA and RI for the distribution parameters of the dam-site floods of the TGR. It can be seen that the predicted SPA series under the emission scenarios of SSP1-2.6 and SSP2-4.5 exhibit similar evolutions without obvious trends, while the SPA under the emission scenario of SSP5-8.5 displays a visible increasing trend after 2050. More reservoirs with huge capacity will be completed in the next decade, hence the covariate of RI is predicted to be significantly enlarged. The dam-site flood distributions of the TGR for the period from 2021 to 2100 are obtained by using the predicted SPA and RI as the covariates of the time-varying parameters of the PIII distributions (see Table 2).

Figure 5 presents the distributions of the dam-site flood peak of the TGR during the period from 1470 to 2100. It can be seen that the estimated flood distribution is able to capture the evolution of the observed flood variables during the period from 1882 to 2020. The modeling results indicate that the mean values of the predicted flood distributions in the next decade will exhibit a sharp decline, which is mainly attributed to two huge reservoirs named Wudongde and Baihetan with a total flood control capacity of 10 billion m³ (see Table 1). Due to the increasing trend of the predicted SPA under the emission scenario SSP5-8.5, the corresponding flood distribution in the future period will present an upward trend after 2050. The distributions of the dam-site flood volumes of the TGR present the behavior similar to that of the flood peak.

After estimating the flood distribution of each individual flood variable, the 5-D C-vine copula based on the GH copula is employed to construct the multivariate joint distribution of $\left(Q_1,V_3,V_7,V_{15},V_{30}\right)$.

Table 3. Parameters of the 5-D vine copula for the dam-site floods of the TGR.

Bivariate Pair	Copula Parameter	Bivariate Pair	Copula Parameter	Bivariate Pair	Copula Parameter
$\left(Q_1,V_3\right)$	13.72	$F\left(Q_3|Q_1\right),F\left(Q_7|Q_1\right)$	1.77	$F\left(Q_7|Q_1,Q_3\right),F\left(Q_{30}|Q_1,Q_3\right)$	1.45
$\left(Q_1,V_7\right)$	4.62	$F\left(Q_3|Q_1\right),F\left(Q_{15}|Q_1\right)$	1.30	$F\left(Q_{15}|Q_1,Q_3,Q_7\right),F\left(Q_{30}|Q_1,Q_3,Q_7\right)$	1.60
$\left(Q_1,V_{15}\right)$	2.66	$F\left(Q_3|Q_1\right),F\left(Q_{30}|Q_1\right)$	1.19
$\left(Q_1,V_{30}\right)$	2.19	$F\left(Q_7|Q_1,Q_3\right),F\left(Q_{15}|Q_1,Q_3\right)$	1.93

summarizes the estimated parameters of the C-vine copula. The PIT test indicates that the C-vine has a satisfactory fitting effect with passing the goodness-of-fit examination at the 0.05 significance level.

4.3. Multivariate Dam-Site Design Floods for the TGR

Based on the multivariate joint distribution of the dam-site floods of the TGR, the design floods are calculated by using the criteria of AAR. In this study, we compare the dam-site design floods derived from the predicted flood distribution in the future period from 2021 to 2100 with those derived from the flood distribution in the historical period from 1470 to 1997. It is important to note that during the historical period the reservoirs in the Upper Yangtze basin have no or very little impact on the dam-site flood distribution of the TGR, while the predicted flood distribution in the future period has involved the effects of both reservoir regulation and potential climate change.

The multivariate dam-site design floods associated with the OR, AND Kendall exceedance probabilities are represented by the most-likely design values, which are given in

Table 4. Multivariate dam-site design floods associated with OR exceedance probability for the TGR.

Flood Variable	AAR = 0.9				AAR = 0.99				AAR = 0.999
	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5
Q1 (m3/s)	67,583	48,925	49,894	52,924	82,332	60,710	60,770	65,580	97,146	70,523	72,662	77,138
V3 (109 m3)	16.8	12.5	12.7	13.5	20.6	15.4	15.5	16.6	24.3	18	18.6	19.7
V7 (109 m3)	35.7	27.5	28.2	29.5	44	33.7	33.6	36.2	51.6	39.6	39.9	42.8
V15 (109 m3)	66	51	52.7	54.6	80.6	61.9	62.5	67.5	93	71	72.1	78.7
V30 (109 m3)	115.2	93.9	93.5	99.5	137.8	111.6	113.6	122.2	158.1	126.7	128.2	141.2

,

Table 5. Multivariate dam-site design floods associated with AND exceedance probability for the TGR.

Flood Variable	AAR = 0.9				AAR = 0.99				AAR = 0.999
	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5
Q1 (m3/s)	58,617	44,196	44,401	46,578	74,561	55,696	56,423	59,444	89,647	64,802	65,967	69,347
V3 (109 m3)	14.7	11.3	11.3	11.8	18.7	14.2	14.3	15.2	22.5	16.5	16.8	17.8
V7 (109 m3)	31.5	25	24.8	25.9	39.6	31.4	31.4	33.1	47.7	36.1	36.7	39.0
V15 (109 m3)	59.9	47.1	47.6	49.7	73.6	57.9	58.1	60.6	85.9	66.9	66.1	71.8
V30 (109 m3)	106.5	85.2	86	91.6	130.2	103.8	103.3	111.8	146.2	118.9	117.9	130.6

and

Table 6. Multivariate dam-site design floods associated with Kendall exceedance probability for the TGR.

Flood Variable	AAR = 0.9				AAR = 0.99				AAR = 0.999
	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5
Q1 (m3/s)	61,929	45,502	45,911	48,982	79,589	56,674	57,897	61,405	93,683	66,561	67,838	73,077
V3 (109 m3)	15.5	11.5	11.6	12.4	19.9	14.5	14.8	15.6	23.5	17	17.3	18.6
V7 (109 m3)	33	25.6	25.6	27.2	41.9	31.9	32.1	33.9	49.6	36.8	38.2	40.6
V15 (109 m3)	61.8	48.5	48.6	51.1	76.1	58.3	59.6	62.1	88.2	66.7	69.8	74.2
V30 (109 m3)	108.2	88.1	88.5	93	130.4	105.1	106.3	113.9	153.7	118.9	123.6	133.7

. In this study, we considered three typical design levels with the AARs of 0.9, 0.99 and 0.999. It is necessary to point out that, under stationary condition, these three design levels equal the return periods of 10, 100 and 1000 years, respectively. In general, the design floods associated with the OR exceedance probability are larger than those associated with other two exceedance probabilities, and the AND exceedance probability has the smallest design values. For the given design levels, the design floods derived by the predicted future flood distributions under three emission scenarios are all smaller than the design values derived by the historical flood distributions. This finding indicates that the TGR would have a smaller flooding risk in the future.

The SPA series under the emission scenarios SSP1-2.6 and SSP2-4.5 are similar to the historical series, hence the declines in the design floods in the future period are mainly attributed to the regulation of the reservoirs upstream the TGR. In next years the Baihetan and Wudongde reservoirs, which both have huge capacities of flood control, are anticipated to play an important role of declining the inflow flood of the TGR. Compared to the hydrological design value derived from the historical flood distribution, the design flood peak with the AAR of 0.999 generally exhibits an absolute decrease of about 25,000 m³/s. For the emission scenarios of SSP5-8.5, though the precipitation in flooding season is predicted to increase, the design flood peak with the AAR of 0.999 would still present a decrease by about 20,000 m³/s. Hence in the future, the regulation of reservoirs will be of great benefit in reducing the flood risk of the TGR.

The dam-site flood hydrographs with different design levels are derived from the corresponding design values benchmarked against the typical flood hydrograph in 1954, and displayed in Figure 6, Figure 7 and Figure 8. It can be seen that the flood hydrographs derived from the predicted future flood distributions tend to become flatter than those derived from the historical flood distribution. This finding also indicates that the reservoirs should have a more profound impact on flood peak.

In the current engineering practice, the dam-site design floods are usually estimated based on the univariate hydrological design method, which is recommend by the Design Criterion of Reservoir Management of China [34]. The univariate hydrological design method ignores the dependence of multivariate flood variables and assumes that the design values of all flood variables have a same occurrence probability. As shown in

Table 7. Univariate dam-site design floods for the TGR.

Flood Variable	AAR = 0.9				AAR = 0.99				AAR = 0.999
	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5	Historical	SSP1-2.6	SSP2-4.5	SSP5-8.5
Q1 (m3/s)	64,221	47,119	47,568	50,338	80,236	58,766	59,420	63,533	94,651	69,258	70,128	75,508
V3 (109 m3)	16.0	12.0	12.1	12.8	20.1	15.0	15.1	16.2	23.7	17.7	17.9	19.2
V7 (109 m3)	33.8	26.1	26.3	27.7	42.2	32.6	32.9	34.9	49.8	38.4	38.8	41.4
V15 (109 m3)	62.9	49.1	49.5	52.1	76.9	59.9	60.5	64.2	89.3	69.3	70.1	74.9
V30 (109 m3)	111.8	90.5	91.3	96.6	135.1	108.7	109.9	117.6	155.2	124.4	126	135.9

, the dam-site design flood peak with the AAR of 0.999 (which is equal to the return period of 1000 years under stationary condition) is 94,651 m³/s, which is close to the design value (98,000 m³/s) used in the engineering practice of the TGR. It can be seen that the univariate design flood values are generally smaller than the multivariate design values associated with the OR and Kendall exceedance probabilities, and larger than those associated with the AND exceedance probability.

5. Conclusions and Discussion

Under a changing environment, both the reservoirs in the Upper Yangtze basin and the changes of climatic conditions can lead to non-stationarity in the dam-site flood frequency of the TGR. The current hydrological design based on the historical flood records might fail to deal with the operation of the TGR. In this study, we perform an analysis of the multivariate dam-site flood frequency of the TGR by considering future reservoir regulation and precipitation, and then derive a multivariate hydrologic design using the design criterion of AAR. The main conclusions of this study are as follow:

After the regulation of the TGR, the dam-site flood variables are found to present obvious departure from the observed floods downstream the TGD. This finding indicates that the TGR has played a profound role in dampening the downstream flood magnitude, especially in the recent decade.

For all at-site flood variables of the TGR, the location parameters (i.e., mean values) of the PIII distributions are linked to the covariates of RI and SPA. In particular, the flood mean values present a negative relationship with RI and a positive relationship with the SPA. In the future period, the flood mean values are predicted to present a significant decline due to the regulation of the newly-built reservoirs in the Upper Yangtze basin.

The dam-site design floods of the TGR in the future are smaller than the design values derived from the historical flood records. Different emission scenarios will result in different design floods. The design flood peaks under the high emission scenario (SSP5-8.5) are larger than those under the lower emission scenarios (SSP1-2.6 and SSP2-4.5).

Compared to previous studies, which usually only concerned the impact of the TGR on downstream floods, this paper focuses on the dam-site flood frequency of the TGR. This study reveals that the reservoirs upstream of the TGR will be of great benefit in reducing the flood risk of the TGR. As more and more reservoirs will be put into operation, the joint regulation of the TGR and these reservoirs upstream should further reduce the flood risk of both the TGR and the region downstream of the TGD. For a better economic benefit, it is necessary to make an adjustment for the current operation strategy of the TGR. For example, making a proper lift of the water level of the TGR in the flooding season will not lead to exceedance in the design risk level of the dam, but can increase electricity generation.

In this study, the effect of the upstream reservoirs on the dam-site flood of the TGR is identified by a statistical method. The sample size of the observed series would be likely to induce a large uncertainty in the modeling results since only in recent years have the majority of the reservoirs in the Upper Yangtze basin been put into operation. It requires an extension of the observed flood records or a process-based simulation to build a more robust relationship between the dam-site flood distributions of the TGR and the reservoirs in the Upper Yangtze basin.