Copula-based modeling of hydraulic structures using a nonlinear reservoir model

An extreme hydrological event is considered a major cause of hydraulic structure failure, and thus, its occurrence probability provides an important basis for the evaluation of the structure return period in the design of hydraulic structures. Thus, it is important to choose appropriate variables and methods for the risk assessment of an extreme hydrological event acting on the hydraulic structures. However, the widely used univariate flood frequency analysis does not account for the interactions between hydrological events and hydraulic structures. The extreme load that can be supported by a hydraulic structure in a flood event is determined by the storage and routing processes, and changes in any characteristic variable of the flood event, such as flood peak, volume and duration, and the hydraulic structure, can have a substantial effect on the extreme load.

In univariate flood frequency analysis, the return period of a given hydrological variable is taken as the return period of structural failure. It is evident that the univariate return period of the flood peak or volume cannot provide a complete evaluation of the occurrence probability of hydrological events (Chebana & Ouarda 2011), which therefore may result in over- or underestimation of the flood risk and hydraulic design standard (Salvadori & De Michele 2004). Apparently, multivariate flood frequency analysis is preferred over univariate flood frequency analysis in hydraulic design. However, a number of drawbacks have also been identified in previous literature. For instance, a normal distribution is often assumed for non-normally distributed flood peak and volume in multivariate normal distribution and normal transformation methods. Some marginal distributions may be combined to form a joint distribution, such as bivariate lognormal distribution, bivariate exponential distribution, bivariate Gumbel distribution, and two-dimensional P-III distribution. However, it is generally assumed that the marginal distributions should come from the same family of distributions.

Recently, the Copula function has been widely used in hydrology, especially in multivariate flood frequency analysis. The copula is a function that links univariate distribution functions to a multivariate distribution function describing the dependence among correlated variables (Nelsen 2006). However, there is no universally agreed method for calculating the multivariate return period based on the copula function. Salvadori & De Michele (2004) investigated both unconditional and conditional return periods of hydrological events using copulas, and they found that it was possible to combine the marginal events in several ways using the (inclusive) OR operator and the AND operator. Salvadori et al. (2011) introduced another return period, the secondary return period (also called the Kendall return period), which was associated with the realization of dangerous events for the dam. It should be noted that the return periods mentioned above are defined based on the natural occurrence probability of floods. However, the climate changes and human activities in recent years may make the hydrological extreme value series highly non-stationary. For the frequency analysis of non-stationary hydrological extreme value series in a changing environment, the return period should be redefined and several non-stationary methods have been developed, such as expected waiting time, expected number of exceedances, equivalent reliability, and average design life level (Yan et al. 2017, 2020).

However, the return period of hydraulic structures should be defined in terms of the risk of dam overtopping or the downstream loss. In multivariate cases, structural failure is considered to be the result of the complex interactions between the structure and the hydrological loads acting on the structure, which indicates that the effects of the reservoir and dam on the inflow hydrographs should also be taken into consideration (Mediero et al. 2010). The relationships between hydrological and structure variables can be simulated using observed data. De Michele et al. (2005) generated peak-volume pairs using Gumbel copula in order to verify that the maximum water level reached at the dam by the generated hydrographs was below the crest level. Klein et al. (2010) proposed a methodology to categorize flood events for risk-based analysis and the design of flood control systems using copulas, and the maximum water level reached at the dam was estimated to graphically analyze its relation with the primary return period for each event. Requena et al. (2013) compared the multivariate return period of the hydrological loads with the structure return period by simulation.

An idealized reservoir model can also be established to simulate the relationships between hydrological and structure variables. Notably, the reservoir volume-water level and discharge-water level relationships are assumed to be linear due to the inherent difficulties in mathematical derivation. For instance, Volpi & Flori (2014) proposed an idealized reservoir model in which the reservoir volume-water level and discharge-water level relationships were assumed to be linear and the inflow process was assumed to be constant, based on which the relationships of the maximum discharge with flood peak and volume were derived. Balistrocchi et al. (2017) proposed triangular hydrographs by assuming linear reservoir volume-water level and discharge-water level relationships. However, the reservoir volume-water level and discharge-water level relationships are more likely to be nonlinear rather than linear in actual reservoir operation depending on the hydraulic structures such as sluice gates, overflow weirs, and orifices. Thus, linearization can lead to over- or under-estimation of the structure return period and flood risk.

The main purpose of the present study is to determine the relationships of the maximum water level at the dam with flood peak and volume by assuming nonlinear reservoir volume-water level and discharge-water level relationships and to develop a methodology for calculating the structure return period based on the copula function. The rest of this paper is organized as follows. In Section 2, we illustrate the structure-based nonlinear multivariate approach for hydraulic structure design and risk assessment. For simplicity, only the bivariate case is considered here. In Section 3, we develop an idealized nonlinear reservoir model, based on which the calculation method of the structure return period is derived. In Section 4, the structure return period based on the nonlinear reservoir model is calculated and then compared with that based on the linear reservoir model and the bivariate return period based on the Kendall distribution function. In Section 5, we discuss the effects of the discharge and storage capacity of the reservoir on the return period. Finally, the main conclusions are drawn in Section 6.

The return period is used in hydrology to estimate the recurrence interval of a hydrological event such as flood, rainstorm, and extremely high or low water level with the severity or duration greater than or equal to a given value.

In univariate analysis, X is a random variable that measures the hydrological loads such as flood peak or volume, and F_X is the probability distribution function of X. The occurrence of X is assumed to be time-independent. The return period of {X > x} is expressed as ⁠, where is the average time interval between two successive occurrences of the hydrological event X (μ = 1 year in the case of annual maxima).

Copula has emerged as a useful tool for multivariate hydrological frequency analysis, and the Archimedean and extreme value copula families are most often used for modeling flood variables. A distinct property of Archimedean copulas is that they can be constructed easily and allow for a great variety of dependence structures. Archimedean copulas such as Frank copula (Favre et al. 2004) and Clayton copula (Shiau et al. 2006) have been widely used to characterize the dependence structure between peak and volume variables. Extreme value copulas are able to connect the extreme values of the studied variables, which are very important in flood frequency analysis. The Gumbel copula is considered to best represent the relationship between peak and volume (Zhang & Singh 2006). Many bivariate return periods have been estimated by copulas in recent years, and the most representative ones are the ‘AND’ return period and the Kendall return period.

The critical level curve is the set of points corresponding to the probability level p, ⁠, and the corresponding supercritical region is given by ⁠, in which X and Y can be infinite, and the coordinate pairs of and on the critical level curve represent the probability p.

The concept of structure return period is inspired by the supercritical region of the Kendall return period. It is noted that both ‘AND’ and Kendall return periods account only for the natural occurrence of floods, while the structure return period takes into consideration the response of hydraulic structures of the reservoir to floods. In this case, the critical level depends critically on the consistent responses of the hydraulic structures, and the supercritical region is the region over the structure response isoline.

Z is defined as the design parameter characterizing a hydraulic structure (e.g., the dimension of the spillway and the elevation of the levee) or variables associated with the structure. In actual applications, Z indicates the effect of hydrological loads acting on the structure. There is a strictly monotonically increasing function between Z and hydrological load X in a univariate case. The structure function is Z = g(x), where there exists a unique inverse function g⁻¹. Note that this assumption holds for all hydraulic structures. In a univariate case, the probability distribution function of Z is the same as that of X, and thus ⁠, where ⁠. Then, the return period of Z is the same as that of X, both of which can be used in the evaluation of hydrological events.

The return period of structural failure can be easily obtained by using the standard univariate frequency analysis of the random variable Z, ⁠, which is the average time interval between two successive events in the same supercritical region ⁠. A vector of the design parameter Z can also be obtained from the return period, and the expression of may be very complex in this case. Thus, only a single design variable Z is considered in this study for the sake of simplicity.

Here, the objective of hydrological design is to determine the combination of x_p and y_p pairs in a given supercritical region at a given failure probability p. The supercritical regions may differ substantially in terms of their boundaries, which is closely associated with the structure under consideration and associated variables. The boundary of the supercritical region can be determined by the structure function g(x,y) described in Equation (5), which represents the hydrological load–structure interactions. Thus, it is imperative to consider the structure and the way it interacts with the hydrological load in determining the suitable structure function g(x,y). However, it is noted that the expression of g(x,y) is too complex to be described mathematically in many actual cases, and thus, there is a need to obtain the analytical solutions in a simpler and more efficient fashion.

It is clear that the structure function g(x,y) that can more accurately capture the interactions between the structure and the hydrological load acting on the structure plays a key role in calculating the structure return period. Given the complexity in mathematical derivation, an idealized reservoir model was proposed in previous studies of Volpi & Flori (2014) and Balistrocchi et al. (2017), where the reservoir volume-water level and discharge-water level relationships were assumed to be linear. However, both of them are likely to be nonlinear rather than linear in actual reservoir operation depending on the combinations of different drainage facilities such as sluice gates, overflow weirs, orifices and hydraulic turbines, as shown in Figure 1. Thus, this study focuses on the structure return period based on a nonlinear reservoir model. First, a nonlinear reservoir model is developed. Then, reservoir operation is simulated for different inflow hydrographs based on the water balance principle, and the relationships between the design variable (the maximum discharge) and two flood variables (flood peak and volume) are determined.

The nonlinear reservoir model is schematically shown in Figure 1. is the reservoir area, which is assumed to increase continuously with the increase of water level H (⁠ m at the beginning of flood discharge), and the reservoir volume V and discharge are also calculated from ⁠.

It is known that is associated with different discharge structures, which is often taken to be for orifice outflow and for weir flow in fitting the discharge capacity curves. In this study, is set to range from 0.5 to 2 due to the use of multiple discharge structures.

Let the maximum design discharge be the same as the maximum discharge in Equation (24), and thus ⁠. The function has two asymptotes: as and as ⁠. It follows from Equation (22) that the same can be derived from the limiting cases of the inflow hydrograph, which is characterized by an infinite peak flow but with a constant flood volume equal to ⁠, or by an infinite flood volume but with a constant peak flow equal to ⁠.

is the probability distribution function of the maximum design discharge ⁠; is the joint probability density function of random variables and W, where and are their marginal distribution functions, respectively. When the design variable is independent of one of the two random variables (e.g., ⁠), the boundary of can be described as a straight line ⁠, and thus ⁠.

It is important to note that the reservoir discharge-water level relationship remains approximately constant during the routing process, and both the maximum water level at the dam and the maximum reservoir volume can be achieved once the maximum discharge is achieved, which strictly corresponds to the structure return period.

The Xin'an reservoir is located on the tributary of the Qiantang River in Tongguan Xia, Zhejiang Province, China Figure 2, and the upstream river is 323 km long with a drainage area of 10,442 km². This reservoir is intended for flood control and power generation, and it has a multiyear regulation capacity with a normal pool level of 108 m and a storage capacity of 17.84 billion m³. The main discharge structure of the reservoir is the spillway with nine sluice gates, where the net width of the spillway is 117 m, and the elevation of the spillway crest is 99 m. The flood data of the Luotongbu Station, the inflow observation station of the reservoir, for the period 1961–2008 were collected, from which the annual maximum flood peak flow Q_max and the annual maximum 10-d flood volume W₁₀ were extracted.

The marginal distribution functions of Q_max and W₁₀ are derived. The generalized extreme value (GEV) distribution is used in this study as the marginal distribution function due to its wide use and good performance in flood frequency analysis. The Gumbel copula is considered to best represent the relationship between peak and volume (Zhang & Singh 2006). Three copulas (G–H, Frank, and Clayton) from the Archimedean copula family are used in this study to simulate Q_max–W₁₀ joint distribution, and the maximum likelihood method is used to estimate the parameter θ. The optimal copula function is selected using Akaike information criterion (AIC), Bayesian information criterion (BIC) and root-mean-square error (RMSE), and the smaller the values of RMSE, AIC and BIC, the better the fitting of the copula function. Table 1 shows that the smallest AIC and BIC values are obtained by G–H copula, but its RMSE value is slightly larger than that of Frank copula. Thus, the G–H copula is used to describe the Q_max–W₁₀ joint distribution in this study.

In this study, the maximum water level in the reservoir above the spillway crest is taken as the structure design variable. A nonlinear reservoir model is established, and then the structure return period is calculated and compared with that obtained based on a linear reservoir model.

• (1)

  Reservoir volume-water level relationship

The structure return period based on the linear reservoir model in which both reservoir volume-water level and discharge-water level relationships are assumed to be linear is also calculated following the same procedure and compared with that based on the nonlinear model (see Volpi & Flori (2014) and Balistrocchi et al. (2017) for details).

In this section, the structure return periods based on linear and nonlinear reservoir models are calculated and compared, and possible explanations of the difference are discussed.

Figure 3(a) shows the comparison of the maximum water level–return period relationships between linear and nonlinear models. It shows that when T < 360a, the of the nonlinear model is larger than that of the linear model at the same return period; while the opposite trend is observed when T > 360a. As the value increases, the T value of the nonlinear model increases more rapidly, which can be attributed to the more rapid change of and with ⁠, as shown in Figure 3(b).

Six representative return periods (T = 10, 50, 100, 200, 500, and 1,000a) are selected for further analysis, and Table 2 shows the ⁠, ⁠, and at different return periods for linear and nonlinear models. It is found that the values of the linear model are smaller than that of the nonlinear model at T = 10, 50, and 100a, while the opposite is true for the values of and ⁠. The values of ⁠, ⁠, and increase with the increase of the return period, which, however, is more pronounced for the linear model. This is because the discharge changes more rapidly with the water level for the nonlinear model, as shown in the –H curve. At a return period of 500 years, the values of ⁠, ⁠, and of the linear model are larger than those of the nonlinear model, indicating that flood risk is over-estimated by the linear model.

Figure 4 shows the comparison of the structure return period isolines at T = 10, 50, 100, 200, 500, and 1,000a between linear and nonlinear reservoir models, and gray dots are the 1,000 random combinations of flood peak and volume generated using the G–H copula. It is important to note that any flood event on the same isoline would produce the same maximum reservoir volume, maximum water level at the dam, and maximum discharge, and thus have the same impact on the reservoir. However, given their different assumptions about the reservoir volume-water level and discharge-water level relationships, the boundary of the supercritical region differs substantially between linear and nonlinear reservoir models, resulting in substantial differences in the isolines at the same return period. It is seen that at a return period of T = 10, 50, 100, and 200a, the discharge for the linear reservoir model is higher compared with that for the nonlinear reservoir model, and the nonlinear reservoir model is more sensitive to floods with high flood peak but low flood volume; whereas at a return period of T = 500a and T = 1,000a at which the maximum water level at the dam is high, the discharge and volume are approximately the same between linear and nonlinear reservoir models, and as a result the isolines are approximately coincident with each other.

It is thus concluded that reservoir discharge and volume would be increased to different extents with increasing water level for linear and nonlinear reservoir models due to their differences in the reservoir volume-water level and discharge-water level curves, which consequently have an effect on the structure return period. In addition, the reservoir volume-water level and discharge-water level curves can also have an effect on the boundary of the supercritical region, thus resulting in differences in the isolines at the same return period between linear and nonlinear reservoir models.

For the Kendall return period, the supercritical region is defined as the region enveloped by the isoline of the joint return period. Here, the structure return period which is applicable to the flood control structures and the Kendall return period which is applicable to natural occurrence of floods are compared. The isolines of the structure return period and the Kendall return period at a return period of T = 10, 50, 100, 200, 500, and 1,000a are shown in Figure 5, and gray dots in Figure 5 are the 1,000 random combination of flood peak and volume generated using the G–H copula.

It shows that the Kendall return period isoline is symmetric about the straight line with an angle of 45°, as it is based on the joint return period; while the symmetry axis of the structure return period isoline is more inclined to the y-axis, which is affected by the reservoir structure such that a flood event with a larger peak flow can pose a greater risk to the reservoir structure.

There are an infinite number of peak–volume pairs at the same Kendall return period, which can produce different reservoir responses and thus pose different levels of risk to the flood control structures. The Kendall return period of T = 10a is analyzed. It is noted that the isoline with a Kendall return period of T = 10a intersects that with a structure return period of 10, 50, and 100a, which are denoted as points A, B, C, and D. The occurrence probability of floods at points A, B, C, and D is 10a; the return period of the maximum water level at the dam for peak–volume pairs at points A and D is also 10a; while that at points B and C is 50 and 100a, respectively. Thus, the maximum water level in response to peak–volume pairs on the left side of point A and the right side of point D on the Kendall return period isoline of T = 10a is higher than 10a, while that between points A and D is lower than 10a. It is concluded that the reservoir would respond differently to different peak–volume pairs at the same Kendall return period, and thus flood risk is likely to be substantially underestimated in some circumstances. At a Kendall return period of T = 10a, the maximum reservoir response is 196.34a and the minimum reservoir response is 8.96a.

The Kendall return period is often used in bivariate design, and the maximum likelihood method, the most likely composition method, and the same frequency method are used for flood peak and volume estimation. The bivariate design based on the Kendall return period accounts for the natural occurrence of floods without considering the effects of hydraulic structures during the routing process, which can lead to over- or underestimation of flood risk from the perspective of hydraulic structure safety.

Although the nonlinear reservoir model has been successfully applied to the Xin'an Reservoir, other reservoirs may show quite different reservoir volume-water level and discharge-water level relationships (⁠ and ⁠). Thus, the effects of the four parameters (⁠⁠, ⁠, ⁠, and ⁠) on the results are analyzed to verify the applicability of the proposed model to other hydraulic structures.

The discharge capacity of a reservoir depends largely on the spillway length and discharge regime, which correspond to and in Equation (8), respectively. Specifically, indicates the discharge capacity of the reservoir under different discharge regimes, which can be affected not only by the structure, such as sluice gate and overflow weir, but also by flow patterns such as orifice and weir outflow through the sluice gate. In general, the larger the discharge capacity is, the larger the value will be. indicates the spillway length, and the larger the spillway length is, the larger the value will be. Here, a constant reservoir volume-water level relationship is considered ⁠, and the effects of (100, 200, and 300) and (0.5, 1, 1.5, and 2) on the discharge capacity are investigated.

Figure 6 shows the relationships between the structure return period (T) and the maximum water level (⁠⁠) for different discharge curves. It shows that at the same structure return period T, an increase in or can lead to an increase in the discharge capacity and consequently a decrease in ⁠. However, it is important to note that there is a nonlinear relationship between T and ⁠, and T increases at an increasing rate with the increase of ⁠.

As indicates the discharge capacity of the reservoir under different discharge regimes, it can have a significant effect on the return period. For instance, as is increased by 100, 50, and 33.3% from 1/m₁ = 0.5 (Q = 100H^0.5) to 1/m₂ = 1(Q = 100H), 1/m₃ = 1.5 (Q = 100H^1.5), and 1/m₄ = 2 (Q = 100H²) at K₂ = 100, the is decreased by 10.9, 23.4, and 34.3% from h₁ = 19.49 m to h₂ = 17.38 m, h₃ = 13.30 m, and h₄ = 8.74 m at T = 10a, respectively. Under the same conditions, the structure return period for the same would increase at an increasing rate with the increase of ⁠. This is because the discharge capacity of the reservoir is increased as the increases, and thus the larger the value is, the higher the increase rate of the discharge capacity of the reservoir will be.

It is also noted that ⁠, which is used to indicate the increase or decrease in the spillway length, can also have an effect on the ∼ T relationship. For instance, as is increased by 100 and 50% from k_2,1 = 100 (Q = 100H^0.5) to k_2,2 = 200 (Q = 200H^0.5), and k_2,3 = 300 (Q = 300H^0.5) at T = 10a and 1/m = 0.5, the is decreased by 9.1 and 8.0% from h₁ = 19.49 m to h₄ = 17.72 m and h₅ = 16.30 m, respectively. Under the same conditions, the structure return period for the same would also increase at an increasing rate with the increase of ⁠. Again, this is because the discharge capacity of the reservoir is increased as the increases, and thus the larger the value is, the higher the increase rate of the discharge capacity of the reservoir will be.

The relationship between reservoir water level and its volume is closely associated with the geographical location and characteristics of the reservoir, which can also have an effect on the structure return period. Under the conditions of ⁠, the effects of parameters (K₁ = 1, 1.5, 2, or 3) and (n/m = 1, 1.2, 1.5, or 2) in are discussed.

Figure 7 shows the ∼ T relationships for different reservoir volumes. At the same structure return period, an increase in or can lead to an increase in the storage capacity and consequently a decrease in ⁠. It is also note that the relationship between T and is nonlinear, and T increases at an increasing rate with the increase of ⁠.

Note that as the storage capacity is mainly affected by the geographical location of the reservoir, it would be much more difficult to change than reservoir discharge capacity. Nevertheless, changes in reservoir storage capacity can have a substantial impact on the safety of hydraulic structures. For instance, under the conditions of (V = H^1.1, V = 1.5H^1.1, V = 2H^1.1 and V = 3H^1.1), an increase of by 50, 100, and 200% results in a decrease in by 12.1, 22.4, and 37.4% at T = 10a, respectively. Under the same conditions, the structure return period for the same would increase at an increasing rate with the increase of ⁠. Again, this is because the storage capacity of the reservoir is increased as the increases, and thus the larger the value is, the higher the increase rate of the storage capacity of the reservoir will be.

For curves with (V = H, V = H1.1, V = H1.2, V = H1.5, and V = H2), an increase of by 10, 20, 50, and 100% results in a decrease in by 5.6, 12.6, 33.0, and 54.2% at T = 10a, respectively. Under the same conditions, the structure return period for the same would increase at an increasing rate with the increase of ⁠. Again, this is because the storage capacity of the reservoir is increased as the increases, and thus the larger the value is, the higher the increase rate of the storage capacity of the reservoir will be.

Figure 8 shows the time series of the annual maximum flood peak flow Q_max and the annual maximum 10-d flood volume W₁₀ recorded at the Luotongbu Station for the 48-year period 1961–2008. It is seen that both Q_max and W₁₀ time series show an increasing trend for the period 1961–2008.

The annual flood peak and the annual maximum 10-d flood volume records at the Luotongbu Station for the period 1961–2008 were fitted based on the GAMLSS model and the parameters were estimated. Figure 9 shows the evolution of Q_max and W₁₀ at a return period of 1,000, 500, 200, 100, 50, and 10 years for the period 1961–2008.

The distribution functions of the peak flow and volume in 1961, 1971, 1981, 1991, and 2001 under non-stationary conditions are selected, and the non-linear model is applied to the Xin'an reservoir to verify its applicability under non-stationary conditions.

Figure 10 shows the relationships between the maximum water level at the dam and the return period under non-stationary and stationary conditions, and Figure 11 shows the marginal distributions of annual maximum flood peak flow Q_max and the annual maximum 10-d flood volume W₁₀ under non-stationary and stationary conditions. It is seen from Figure 10 that under non-stationary conditions, the maximum water level at the dam H_max for the same return period also shows an increasing trend for the period 1961–2008; while the return period T for the same year increases at a higher rate with the increase of the maximum water level at the dam. However, it should also be noted that the increasing rate is higher under non-stationary conditions than that under stationary conditions as a result of the marginal distributions of annual maximum flood peak flow and annual maximum 10-d flood volume. As shown in Figure 11, higher flood peaks are more likely to occur under non-stationary conditions, which can cause more severe damage to the reservoir structure and consequently the return period.

Six representative return periods (T = 10, 50, 100, 200, 500, and 1,000a) are considered for comparison Table 3. Different Q_max–W₁₀ joint distributions are used to obtain the nonlinear return periods under non-stationary conditions, and the increasing trend of Q_max and W₁₀ for the period 1961–2008 implies that there is a higher probability of the joint occurrence of large Q_max and W₁₀. Under conditions of constant reservoir volume-water level and discharge-water level relationships, a higher maximum water level at the dam is more likely to occur, which in turn may cause more severe damage to the reservoir. At a return period of T = 10a, the probability of the Q_max–W₁₀ joint distribution under non-stationary conditions is similar to that under stationary conditions. Thus, under conditions of constant reservoir volume-water level and discharge-water level relationships, there would be no substantial difference in the return period between non-stationary and stationary conditions. However, the probability of the joint occurrence of large Q_max and W₁₀ would increase as the return period increases under non-stationary conditions because of the effects of the marginal distributions of annual maximum flood peak flow and annual maximum 10-d flood volume. In this case, the return periods obtained under stationary conditions are more conservative.

In this study, a general nonlinear reservoir model was proposed for the design and risk assessment of hydraulic structures in a bivariate context. A case study was then performed with the Xin'an reservoir, and the nonlinear structure return period was calculated and compared with the linear structure return period and the Kendall return period. The main conclusions are summarized as follows.

1. The supercritical region D_Z plays a determinant role in the design and risk assessment of hydraulic structures, and its boundary is affected by the structure and hydrological load acting on the structure. Thus, it is imperative to consider the structure and the way it interacts with the hydrological load in determining the structure function. In this study, a nonlinear reservoir model is established, and the relationship between structure and hydrological load is determined using the approximate polynomial integration method and then the supercritical region D_Z is obtained. Unlike the linear reservoir model, nonlinear reservoir volume-water level and discharge-water level relationships are assumed in this study, thus making it possible to capture the routing processes of the reservoir more accurately and flexibly.

2. The maximum water level at the dam is taken as the structure design variable, and the structure return period based on the nonlinear reservoir model is calculated and compared with that based on the linear reservoir model. Given their different assumptions about the reservoir volume-water level and discharge-water level relationships, the boundary of the supercritical region differs substantially between linear and nonlinear reservoir models, resulting in substantial differences in the structure return period. The change rates of discharge and reservoir volume with reservoir water level may also differ between linear and nonlinear reservoir models, and thus they show different levels of sensitivity to floods with different peak–volume combinations and thus lead to differences in the structure return period.

3. The structure-based multivariate return period is distinctively different from the hydrological event-based multivariate return period. The structure return period is more applicable than the Kendall return period to flood control structures as it takes into consideration the interactions of the structure and hydrological load acting on the structure. However, it is important to note that the Kendall return period can lead to over- or underestimation for hydraulic structures in multivariate flood frequency analysis.

4. For the nonlinear reservoir model, the effects of ⁠, ⁠, ⁠, and involved in the reservoir volume-water level and discharge-water level relationships (⁠ and ⁠) are investigated to verify the applicability of the proposed model. It is found that an increase in or leads to an increase in reservoir storage capacity; while an increase in or leads to an increase in reservoir discharge capacity. Different reservoirs would respond differently to the same flood event due to their different storage and discharge capacity, which can have an effect on the structure return period.

However, as it is difficult to derive the general solution of the nonlinear differential equation, the inflow hydrograph is simplified to be a constant process in this study, which may make the results less accurate. Efforts should also be made to better describe the fluctuation of the inflow hydrograph and to investigate the effects of the operation rules on the routing processes of the reservoir.

This research is funded by the National Natural Science Foundation of China (52079040), the Fundamental Research Funds for the Central Universities (B210201023), and Jiangsu Water Science and Technology Project (2019027 and 2020065).

All relevant data are included in the paper or its Supplementary Information.