Nonstationary Design Flood Estimation in Response to Climate Change, Population Growth and Cascade Reservoir Regulation

: The hydrologic data series are nonstationary due to climate change and local anthropo ‐ genic activities. The existing nonstationary design flood estimation methods usually focus on the statistical nonstationarity of the flo w data series in the catchment, which neglect the hydraulic ap ‐ proach, such as reservoir flood regulation. In this paper, a novel approach to comprehensively con ‐ sider the driving factors of non ‐ stationarities in design flood estimation is proposed, which involves three main steps: (1) implementation of the candidate predictors with trend tests and change point detection for preliminary analysis; (2) application of the nonstationary flood frequency analysis with the principle of Equivalent Reliability (ER) for design flood volumes; (3) development of a nonsta ‐ tionary most likely regional composition (NS ‐ MLRC) method, and the estimation of a design flood hydrograph at downstream cascade reservoirs. The proposed framework is applied to the cascade reservoirs in the Han River, China. The results imply that: (1) the NS ‐ MLRC method provides a much better explanation for the nonstationary spatial correlation of the flood events in Han River basin, and the multiple nonstationary driving forces can be precisely quantified by the proposed design flood estimation framework; (2) the impacts of climate change and population growth are long ‐ lasting processes with significant risk of flood events compared with stationary distribution conditions; and (3) the swift effects of cascade reservoirs are reflected in design flood hydrographs with lower peaks and lesser volumes. This study can provide a more integrated template for down ‐ stream flood risk management under the impact of climate change and human activities.


Introduction
Traditional design flood for hydraulic structures such as reservoirs is based on a stationary hypothesis, meaning that the driving factors (e.g., climate change, urbanization and reservoir flood regulation) act to generate the hydrological variables in the same way as in the past, present and likely the future [1][2][3][4][5].However, the statistical characteristics of flood series might alter due to a changing environment [6,7].If hydrological engineers do not fully consider the nonstationarity of the hydrological series, the results of the conventional stationary flood frequency analysis would be inaccurate in practice [8].López and Francés [9] used climate and reservoir indices as external covariates in nonstationary flood frequency analysis.Yan et al. [10] considered climate change and population growth to explain the nonstationary properties of hydrological time series.Global warming, the primary factor that drives climate change, has altered the meteorological regimes in some regions [11,12].Population growth will not only lead urbanization but also lead to increased water consumption.Rapid urbanization over recent decades has significantly changed the regional hydrological characteristics of catchments [13][14][15].The river network systems have been obviously influenced by the process of urbanization, aggravating the hazards of floods and water degradation [16,17].Water consumption is broadly reflected in daily life and productions, which can also affect the hydrological regime in the catchment.Reservoirs, one of the effective facilities for flood control, hydropower generation and other social functions, have gradually formed cascade reservoir systems [18,19].The great impact brought by reservoir flood regulation has tremendously altered the hydrological characteristics in rivers [20].To sum up, it is difficult to find a river basin that is not impacted by global warming and anthropogenic activities, particularly in rapidly developing China.
Milly et al. have elucidated the challenges about how to deal with design floods and water resources management under nonstationary conditions [1].The commonly used technique to gain the changing knowledge of flood regimes is nonstationary flood frequency analysis [9,21].Strupczewski et al. [22][23][24] presented a nonstationary approach to at-site flood frequency analysis, in which the distribution parameters are represented as the functions of explanatory variables to explain the nonstationary hydrological series [25,26].Stasinopoulos and Rigby proposed "Generalized Additive Models for Location, Scale and Shape" (GAMLSS), which is a powerful implement for nonstationary frequency analysis of time series [27].The meteorological factors (such as precipitation), population growth and reservoir index (RI) are widely used as explanatory variables incorporated in GAMLSS as covariates [9,28,29].Taking RI for an example, López and Francés [9] considered the effects of RI from two aspects: the data used for flood frequency analysis are observed daily flow series' that have been affected by upstream reservoirs; some of the reservoir characteristic parameters, such as the catchment area and the reservoir total storage capacity, are integrated to form RI. However, the RI series are piecewise constants in spite of reservoir operation rules [30].Although RI has been improved (by a hand of studies) for greater performance in nonstationary model fitting, the strategies of reservoir operation are hard to be quantified with the GAMLSS-based nonstationary flood frequency analysis framework [8,31].Therefore, the covariates grounded in (or modified by) RI are unable to accurately consider the impact of cascade reservoir regulation.The flood regional composition (FRC) combined with reservoir operation rules [32][33][34] can overcome this drawback.The aim of the FRC is to study the flood generation mechanism at the investigated downstream site.The inflow of the target reservoir consists of the first upstream reservoir inflow and all interval inflows between adjacent reservoirs.Then, the outflow of the downstream reservoir can be obtained through the reservoir operation rules.Among all possible compositions based on the water balance equation, an appropriate FRC needs to be selected.Guo et al. [33] proposed the most likely regional composition (MLRC) method and derived theoretical formula for triple cascade reservoirs.The MLRC method presumes that FRCs are diverse with their occurrence probabilities, which can be quantified by the multivariate probability density function (PDF) of flood events occurring at all sub-basins, and the FRC with the largest occurrence probability should be chosen for representing the actual spatial correlation of flood events.With a rigorous statistical basis, the MLRC method has been successfully applied in the cascade reservoirs in the upper Yangtze River [35].In contrast to the nonstationary flood frequency analysis applied with RI, the natural flow data restored by the observed data is employed in the FRC framework.Although the MLRC-based approach takes into full consideration both FRC and reservoir operation rules, it is based upon a stationary assumption condition so that the other variables such as climate change and population growth are ignored.
Under nonstationary conditions, the copulas with time-varying dependence structures have been primarily applied for modeling coincidence probabilities such as the joint return periods [8,[36][37][38].In fact, the changing environments might also alter the statistical correlation of FRC, and the dependence structure of the MLRC-based method is unable to catch the nonstationary spatial correlation.Therefore, we propose a nonstationary MLRC (NS-MLRC) method and then compare it with the MLRC method under stationary distribution conditions.Equivalent Reliability (ER) [39] is employed in this study to concatenate the stationary and nonstationary design criteria grounded in a given return period [40].
In a word, this study focuses on three objectives: (1) to propose a nonstationary design flood estimation framework; (2) to develop and verify the NS-MLRC method; and (3) to estimate design flood hydrographs at downstream sites.The rest of the paper is organized as follows: Section 2 describes the methodology used in this study.Section 3 briefly introduces the study area and data acquisition.Section 4 analyzes the nonstationary design flood estimation results.Section 5 discusses the nonstationary characteristics and the worst regional flood composition.Finally, Section 6 ends with conclusions.

Methodology
The flowchart of nonstationary design flood estimation for cascade reservoirs is shown in Figure 1, of which there are three main modules.The first module preliminarily analyzes the nonstationarity of the flood volume series.The second module calculates the design flood volumes based on nonstationary flood frequency analysis and ER criteria.The last module solves the flood regional compositions using the NS-MLRC method and derives the design flood hydrograph regulated by cascade reservoirs.

Preliminary Analysis Module
This module aims to manipulate the dataset, determine the candidate univariate and multivariate distributions and covariates, and analyze the nonstationarity preliminarily.

Candidate Distribution Functions, Copulas and Covariates
For nonstationary univariate frequency analysis, the type of the probability distribution () Y f  determines which form of frequency curves will be generated.Table 1 lists five probability distributions, namely Gamma, Weibull, Log-normal, Gumbel, and Pearson type III, which are the alternative distributions for flood frequency analysis with the different link functions [41][42][43].On the other hand, the Archimedean copulas are widely used in hydrological researches [8,33].Three mono-parameter Archimedean copulas, i.e., Gumbel-Hougaard, Frank and Clayton copulas (Table 2) are selected as the candidate copulas for the time-varying dependence structure modelling.
The multiple covariates (i.e., explanatory variables, predictors or addictive terms) indicating both population growth and climate change are incorporated to explain the variation in nonstationary time series.The impact of urbanization and water consumption on the hydrologic flood series is considered by the population (Pop) [37,43], while the annual total precipitation (Prcp), one of the important hydro-climate factors, is applied to quantify the effects of climate change [10,44].

Trend Test and Change Point Detection
The Mann-Kendall trend test [45][46][47] and Pettitt's test [48] are widely used for trend test and change point detection, which can be applied to analyze the nonstationarity of time series [43,44].A nonparametric multiple change point analysis approach designed by Matteson and James [49] is used for change point analysis of multivariate time series in this study.

Hydrologic Design Module
This module is responsible for nonstationary flood frequency analysis and design flood volume estimation.

PDF Distribution Moments Link Functions
Gamma exp Three Archimedean copulas used for the time-varying dependence structure modelling in this study.

Copula
Cumulative Distribution Function Parameter Link Function Gumbel-Hougaard

Time-Varying Moments
Generally, the nonstationary flood frequency analysis can be established by the timevarying moments [50,51] based on the GAMLSS model [27].Furthermore, the time-varying moments can be improved by replacing time with other physical covariates which have physical meanings [8,10].The generalized linear model (GLM) [52] is implemented to establish the function between distribution parameters and their predictors.To facilitate understanding, a three-parameter time-varying distribution is taken as an example in the remainder of this literature.

Selection of Time-Varying Distributions
The selection of the covariates and the type of probability distribution are the two steps for seeking the best time-varying moment model.For one specific time-varying distribution, the selection of covariates contains two aspects: (1) The selection of the best GLM for a single parameter, which is conducted by the forward procedure [27].In the forward procedure, all variables not currently in the model are considered for addition at each step, while all variables currently in the model are individually considered for deletion.
(2) For a given distribution, there exist diverse strategies to select the GLMs used to model all the parameters, i.e., μ, σ, and ν.The strategy proposed by Stasinopoulos and Rigby [27] is executed for GLMs selection of all the distribution parameters.
Since the shape parameter ν is too sensitive to the data series, it is often assumed as constant as other studies have done [10,56].In practice, only the first two moments (i.e., μ and σ) are taken to consider nonstationary data series.The optimal probability distribution is selected from these five distributions mentioned in Table 1 based upon the Bayesian information criterion (BIC) [57]: where ML denotes the log-likelihood within a likelihood-based inferential procedure; df represents the total effective degrees of freedom and; n is the length of the data series.
Comparing with the widely used Akaike information criterion (AIC) [58], the BIC has a larger penalty on the over-fitting phenomenon.
After the above selection procedures is performed for each probability distribution, then the distribution corresponding to the smallest BIC value is selected as the optimal univariate distribution.

Equivalent Reliability (ER)
The hydrologic design criteria according to the definition of return period under stationary conditions might be modified to adapt nonstationarity [43,59].To ensure that the nonstationary flood frequency analysis is closely relevant to the practical condition, the flood events can be linked with the design lifespan of hydraulic projects [39].The project reliability in the no flood event exceeds the specific design value p z during the life period T T  of the project, which is defined as follows [59]: where t p is the exceedance probability at year t for design value p z ; the superscript S and NS represent stationary and nonstationary conditions, respectively.If a project is designed to withstand a flood event that occurs once in m years under stationary conditions, then the reliability within the design life period 1 2 T T  of the project is given by The ER supposes that the hydrological design values derived using stationary and nonstationary flood control design standard should have the same reliability during the life time of hydraulic structures [39].Assuming that

S N S T T T T 
 , the nonstationary design value for a specific return period m is denoted by , which is measured by the following equation: For given design life period 1 2 T T  and return period m, the stationary ER criteria is calculated by Equation (4) firstly, and then the nonstationary design value responding to m can be obtained eventually by solving Equation (5).

Design Flood Estimation Module
This module is used to establish the time-varying copula and derive the nonstationary design flood hydrographs based on the NS-MLRC method.

Time-Varying Dependence Parameter of Copulas
Traditional joint hydrological frequency analysis assumes that the parameters of both marginal distributions and copula functions are constant.Nevertheless, either the individual series or the dependence structure between the multiple series might be also nonstationary under changing environments [8,38].To consider such possibility, a general form of the joint distribution of multiple random variables 1 2 ( , ,..., ) Y Y Y at any time t , is built.According to Sklar's theorem [60], the joint PDF can be expressed in terms of its marginal distributions and the associated dependence function, i.e., the time-varying copula can be expressed as follows: where ( ) NS H  denotes the cumulative distribution function (CDF) that defines the de- pendence structure of multiple variables; ( ) C  represents the CDF of copula function; ( ) F  represents the CDF of the hydrological random variables; , , ] where are the GLM parameters.The copula-based GAMLSS with non-random sample selection is applied in this study so that the non-stationarities of marginal distributions can be considered [61,62].
The time-varying dependence parameter of copula is estimated by a penalized likelihood framework with integrated automatic multiple smoothing parameter selections [63].It is hard to calculate the distance between the fitted and the empirical frequency so that the Probability Integral Transformation (PIT) is used for the goodness-of-fit test of time-varying copula [64,65].The forward procedure is applied for the selection of t c  based on GLM, while the BIC criteria is used to select the optimal copula function.

Nonstationary Most Likely Regional Composition (NS-MLRC) Method
The sketch of flood regional composition with cascade reservoirs is shown in Figure 2, in which the random variables For the A1-A2 sub-partition, the inflow of downstream reservoir A2(x2) consists of the inflow of reservoir A1(x1) and the inflow of intermediate basin B1(y1).The inflow of downstream reservoir A2 is affected by the operation of upstream reservoir A1.Then x y is determined as the FRC of 2 x , in which inflow 1 x can be turned into discharge based on the operation strategy of reservoir A1.Thus, the discharge of reservoir A2 is relevant to the FRC of A2, i.e., ( , ) x y .For the cascade reservoir system, all the FRC  For the inherently stochastic mechanism of the flood generation, there exists multiple compositions of The output of the MLRC method is the FRC with the largest occurrence probability.Nevertheless, under changing environments, the FRC derived by the MLRC method may vary over time.According to Sklar's theorem, the time-varying copula is expressed as follows where () c  denotes the PDF of copula; ( | , , ) are the optimal time-varying marginal distributions of flood values occurring at reservoirs 1 2 A , A , ..., A n and downstream site C, respectively.
The composition is maximized by subjecting water balance constraint: is maximized when its first-order derivative equals zero so that the following equation should be satisfied: After deriving the composition x y y y can be obtained subsequently using NS-MLRC method.The Newton iteration algorithm is adopted to solve Equation ( 11) [33].
According to the flood prevention standard of cascade reservoirs, the design flood hydrographs of each sub-system are calculated by the peak and volume amplitude method [67,68] based on the results of the NS-MLRC method.River channel flood routing and reservoir flood control regulations are adopted to derive the design flood hydrograph at the downstream site C [34,69].

Study Area and Data
The Han River basin in China is located between 106-115° E and 30-35° N (see Figure 3), and has a total length of 1530 km.As one of the most important tributaries of the Yangtze River, The Han River rises in the southern of Qinling mountains, and flows from the northwest to the southeast.This mountainous region lies in the humid zone with a subtropical monsoonal climate.The annual average temperature is between 14 and 16 °C.The annual precipitation varies from 700 to 1100 mm, and about 70 to 80% of the annual precipitation occurs during the flood season (May to October), in which heavy rains in early summer and continuous rainfall in autumn often cause major floods.

Cascade Reservoirs
The Ankang (AK) and Danjiangkou (DJK) cascade reservoirs are located at the upper and middle reach of the Han River basin, respectively.The AK reservoir was built in 1992 and provides hydropower generation and flood control, while the DJK reservoir was built in 1973 and its primary functions are flood control, water supply, hydropower generation, and irrigation.These two reservoirs are selected for the case study because their inflow data is lengthy enough (over 60 years) and they have large storage capacities.
The basic information of AK and DJK cascade reservoirs is listed in Table 3.The characteristic parameter values and current flood control operation rules of the AK and DJK reservoirs are provided by the Changjiang (Yangtze River) Water Resources Commission (CWRC), Ministry of Water Resource [70].More information about AK and DJK reservoirs can be found in the references [71,72].

Dataset
Four categories of data series were collected, including restored streamflow data, observed hydro-climate data, population growth data and GCMs outputs from CMIP5.The observed data series provide information up to 2020, while a projected dataset of the future period from 2021 to 2095 is also used in this study.
(1) Restored mean daily streamflow data from both the inflow of AK reservoir and the inflow of DJK reservoir were provided by the Hydrology Bureau of the Changjiang (Yangtze River) Water Resources Commission during 1954-2020.The restoration of streamflow data can be taken as a natural flow series that eliminates the regulation impact of reservoirs.
(2) Observed daily precipitation series from 27 stations during 1951-2020 were obtained from the National Climate Center of the China Meteorological Administration (source: http://data.cma.cn/(accessed on 16 April 2021)).
(3) Given the unavailability of population data at the basin scale, the total registered population of all the prefecture-level cities amidst in the Han River basin was collected.These cities include Wuhan, Shiyan, Jingmen, Xiangyang in Hubei province, Hanzhong, Ankang and Shangluo in Shanxi province, and Nanyang in Henan province.The annual registered population data in Han River basin were obtained from the China Compendium of Statistics 1949-2008 [73], the website of the National Bureau of Statistics of China (source: http://www.stats.gov.cn/tjsj/ndsj/(accessed on 1 July 2021)) and the websites of statistical bureaus of the provinces and cities mentioned above.For future projection, the logistic growth model that adapts the growth restriction resulting from limited natural resources is applied to predict the growth of population [74].Based on the logistic growth model, the evolution of the population for the period 1950-2100 is illustrated in Figure 4.
(4) Future daily precipitation is simulated by GCM using climate change scenarios [75].Our lab research team, Tian et al. (2021) used 10 different GCMs (see Table 4) and two representative concentration pathways (RCPs) of 4.5 and 8.5 from the IPCC Fifth Assessment to project climate change for the Han River basin, which are employed in this study [76][77][78].The RCP 4.5 scenario represents a future of medium emission where climate policies limit [79].Without climate change policies, RCP 8.5 scenario presumes that high emissions of greenhouse gases continue in the future.The outputs of the GCMs not only involve the historical period before 2006 as a reference, but also cover 2021-2095 for future projection.A daily bias correction method is applied in this study for statistical downscaling [80].Six statistics containing mean, standard deviation, 85th, 90th, 95th and 99th percentiles of future precipitation series are used to test the performance of the daily bias correction method [80].Taking the BCC-CSM1-1 model as an example, the bias of raw and the corrected model for daily precipitation series during 1991-2005 are illustrated in Figure 5, where the horizontal and vertical coordinates follow the 27 meteorological stations and the bias about the six statistics, respectively.These findings indicate the great performance of the statistical downscaling method so that it can be employed for future projection.Figure 6 shows the projected future precipitation under RCP 4.5 and RCP 8.5 scenarios.The arithmetic mean values of the 10 GCMs are employed in this study for nonstationary flood frequency analysis.

Preliminary Analysis
According to the regulation characteristics of cascade reservoirs in the Han River basin, the annual maximum 15-day (denoted as W15) flood volume series of the AK and DJK reservoir are selected for flood frequency analysis.Indicated by the fitted trend lines, the overall decreasing trends of W15 series in both AK and DJK reservoirs are shown in Figure 7.The significance of trends in the flood volume series and explanatory variables (i.e., Pop and Prcp) during 1954-2020 is analyzed by the Mann-Kendall test, and the Pettitt's test is employed to detect the change points in the W15 series, which is summarized in Table 5. Findings show decreasing trends of the W15 series in both the AK and DJK reservoirs and the detected change-points are located at 1985-1986 (p-value < 0.05) for two reservoirs.The nonparametric multiple change point analysis method detects the change point of the W15 series between AK and DJK reservoirs, which takes place in 1985 (p-value < 0.05).This preliminary analysis demonstrates that both the univariate series and the dependence structure between the W15 series are all nonstationary.

Univariate Nonstationary Flood Frequency Analysis
Based on the selection procedure mentioned in Section 2.2.2, we assume that the lower bounds of forward stepwise procedure for both μ and σ are constant (i.e., μ ~1, σ ~1), while the upper bounds are Pop + Prcp).Table 6 lists the distribution parameters and BIC values for W15 at the AK and DJK reservoirs, in which the time-varying BICNS and constant BICS are also calculated.According to the minimum BICNS values, the Log-normal distribution is the best one for the W15 series in both the AK and DJK reservoirs.To demonstrate the fitting process, the detailed forward procedure applied in Log-normal distribution parameter μ is displayed in Table 7.The worm plot and centile curves for the Log-normal distribution are plotted in Figure 8.In the worm plot, all scatter points are within the 95% confidence intervals, illustrating a good fit between the optimal distribution and empirical frequency series.In terms of centile curves, the percentages of observation points within the 5th, 25th, 50th, 75th and 95th intervals are 5.97%, 26.86%, 59.70%, 70.14% and 94.02% (5.97%, 28.35%, 53.73%, 71.64% and 94.02%) for AK (DJK), respectively.The above results indicate that the selected optimal Log-normal distribution is more adequate to model the nonstationarity of the W15 series, and to express the similar type of GLMs for both AK and DJK reservoirs.As presented in Table 6, the location parameter of Log-normal distribution for modelling the W15 series is linked to Prcp and Pop while the scale parameter of W15 is constant.Furthermore, the BICNS value is always smaller than the BICS value for the same distribution.Compared to the stationary distributions, the optimal nonstationary Log-normal distribution can satisfactorily capture the nonstationarity in flood frequency analysis for both the AK and DJK reservoirs.

Design Flood Volumes in DJK Reservoir
The nonstationary design value is essential for handling changing environments in considering future flood risk management.Since the total storage capacity of the DJK reservoir is approximately 31 times larger than that of the AK reservoir (Table 3), the DJK reservoir is considered as the key flood prevention project in the Han River basin.If we assume that the design lifespan for the DJK reservoir is 100 years (from 1973 to 2072), then the years of 2021-2072 are chosen as future projections.The design flood volume under nonstationary conditions is calculated by the ER criteria.Three scenarios, i.e., stationary distribution condition (S1), nonstationary conditions based on RCP 4.5 (S2), or RCP 8.5 (S3), are compared.It is worth noting that the S1 scenario only takes cascade reservoir regulation into account while the nonstationary attributes of climate change and population growth are omitted.To produce an intuitive comparison among these scenarios, Figure 9 plots the stationary and nonstationary design flood volumes W15 for return period [2,1000] m .It should be noted that the design flood volumes in the stationary reference scenario are estimated by the Pearson type III distribution and curve-fitting estimation method recommended by the Ministry of Water Resources in China [32].It is found that the design flood volumes of S2 (S3) are lower than S1 when the return periods are less than 300 (200) years, respectively.The design flood volumes of S3 are always larger than S2.Compared with the nonstationary distribution scenarios (i.e., S2 and S3), the method used under the S1 scenario overestimates the flood volumes at a lower return period while the extreme event is underestimated at a higher return period.For the 1000year return period, the design W15 in the DJK reservoir under S1, S2 and S3 scenarios are 20.041,21.680 and 22.352 billion m 3 , respectively.These results indicate that climate change and population growth have greater impact on the features of future W15 under the S2 and S3 scenarios so that the large uncertainty and flood hazard should be considered for future hydrological design and water resources management.

Nonstationary Flood Regional Composition
After the estimation of the design flood volumes at the DJK reservoir, the flood regional composition method can be applied for analyzing the flood generation mechanism of the cascade reservoirs in the Han River basin.For nonstationary analysis of the dependence structure, the forward selection procedure is implemented to model the dependence parameter of three copulas, i.e., the Gumbel-Hougaard, Frank, and Clayton copula.The lower and upper bounds of forward stepwise procedure for θc are constant (i.e., θc ~1) and (Pop + Prcp), respectively.Table 8 summarizes the fitted dependence parameters and the goodness-of-fit results for the optimal copulas under three scenarios for nonstationary modeling of the W15 series in AK and DJK reservoirs.The p-KS(Z1) and p-KS(Z2) are pvalues of the KS test for the two Rosenblatt's probabilities integral transformations Z1 and Z2, which should be uniformly and independently distributed on [0,1].The p-Kendall is the p-value of the Kendall rank correlation test for Z1 and Z2.Table 8 shows that the p-KS(Z1), p-KS(Z2) and the p-Kendall values are larger than 0.05, which statistically supports the validity of the hypothetical models.The time-varying Gumbel-Hougaard copula is the optimal copula to model the dependence structure of the W15 series by comparing the BIC values.Obviously, the nonstationary spatial correlation of flood events can be sufficiently explained by the NS-MLRC method.Figure 3b illustrates that the flood volumes of the DJK reservoir are decomposed to the outflow of the AK reservoir and the inflow of the AK-DJK inter-basin.It is noted that the return period of the design flood is 1000 years in both the AK and DJK reservoirs, as shown in Table 3.The Log-normal distributions with explanatory variables shown in Table 6 are applied as the marginal distributions of the time-varying copulas.Following the procedures described in Section 2.3.2, the NS-MLRC results of the cascade reservoirs can be estimated.
Based on the NS-MLRC method, the nonstationary FRC of the AK and DJK reservoirs with the maximum occurrence probabilities are derived.It should be noted that there are one hundred (from 1973 to 2072) FRCs derived by the NS-MLRC method, compared with the individual result calculated by the MLRC method.Table 9 compares the FRC results under three scenarios during design lifespan with the same flood prevention standard according to ER.It is shown that the designed W15 in the AK reservoir (X1) of S2 and S3 is greater than that of nonstationary cases, while the design W15 at the (Y1) are overestimated when compared with S1.The estimation of design floods at downstream sites can be obtained using design inflow hydrographs and cascade reservoir regulations based on the derived FRC results.Different FRCs may differ in terms of their occurrence probabilities [35], and the FRCs with maximum occurrence probabilities among the lifespan of projects are analyzed.To generate the design flood at the downstream site C (Figure 3), the current flood control regulation and the commonly used Muskingum model provided by the CWRC are used.After flood control operation, 1000-year design flood hydrographs for the downstream sites under three scenarios with or without reservoir regulation are plotted in Figure 10.Table 10 lists the results of the 1000-year design floods under three scenarios with or without reservoir regulation.Figure 10 shows that the 1000-year design flood hydrographs vary considerably with less peaks and gentler flood processes at downstream sites after cascade reservoir regulation.As shown in Table 10, the design floods of the cascade reservoirs decrease significantly on account of the regulation of the AK and DJK reservoirs.Taking Qmax as an example, under the S1, S2, and S3 scenarios, the 1000-year design Qmax of AK (DJK) at downstream sites decreased by 10.97% (49.43%), 5.06% (53.25%) and 5.31% (54.66%), respectively due to cascade reservoir regulation.Furthermore, the W1 and W3 flood volumes are less than half of these without cascade reservoir regulation.

Nonstationary Characteristics
Flood events differ according to climate, reservoir regulation, water consumption, and land use, which are the potential driving forces that link to flood design.When traditional hydrologic design criteria are extended for accommodating nonstationary conditions, the future covariates should be obtained on account of the hydraulic lifespan.However, it is difficult to predict the changes of some variables, such as deforestation, so that other land use indices are not considered in this study.Based on the case study in the Han River, the influence of climate change, population growth and cascade reservoir regulation on the design flood at the downstream site are discussed as follows.
Many studies have identified how the changing pattern brings diverse control on a flow process [28,33], which can be classified as either a direct and/or indirect impact.Climate change and population growth will alter the regional hydrological characteristics of the basin, and affect the flood data series as well as the flood design (Figure 8) [43].Compared to these prolonged effects, the cascade reservoir regulation clips the flood peak and decreases the flood volumes in a transitory and swift way (Figure 10).Among the multiple driving patterns, the cascade reservoir regulation plays a dominant role in affecting the design floods at downstream sites.Furthermore, the cascade reservoir flood control strategies can be re-regulated for adapting the slow-to-change non-stationarity (such as climate change and population growth) in future work.

The Worst FRC during Reservoir Lifespan
In practical operations, the worst FRC during reservoir lifespan is more noteworthy in comparison with the most likely FRC.The reservoir water level is at the flood limited water level during flood season to provide sufficient storage for possible floods [81].When a flood event occurs, the highest water level during the period of reservoir flood regulation can quantificationally indicate the flood hazard, such as dam-break risk, downstream inundation risk and so on.In this study, the highest water level of the DJK reservoir during the whole flood regulation process is applied as an indicator to specifically represent the hazard of flood events.
Although the MLRC method has a strong statistical basis and generates a single regional composition, it is unable to evaluate the composition of the worst regional flood.The NS-MLRC method can include varying factors such as precipitation and population growth.In this study, variant FRCs derived by the NS-MLRC method are used to amplify the typical flood hydrograph, forming one hundred (from 1973 to 2072) FRC types.After cascade reservoir regulation, the comparison of the highest water levels derived from both the most likely FRCs and the worst FRCs at the DJK reservoir during its lifespan is demonstrated in Table 11.It indicates that the highest water level estimated by the most likely FRC is lower than that of the worst FRC, which means the worst FRC under the S2 and S3 scenarios is more adverse than the most likely FRC.In practice, the worst FRC in S3 is worthy of noting for future flood risk management.

Conclusions
There is an increasing need to develop an effective design flood estimation framework to deal with nonstationary data series caused by climate change and anthropogenic activities.In this study, the univariate flood frequency analysis, nonstationary hydrologic design criteria and the NS-MLRC method were adopted to derive the nonstationary design flood volumes in the Han River.The design flood hydrographs at the downstream site were estimated after cascade reservoir regulation.The main conclusions are summarized as follows: (1) The proposed NS-MLRC method can be effectively implemented as an extension of the MLRC method for explaining the nonstationary spatial correlation of the flood events.The multiple nonstationary driving forces, i.e., climate change and population growth, can be captured and precisely quantified by the proposed design flood estimation framework.
As the dominant element affecting flood design, the current cascade reservoir operation strategies can be improved to accommodate the slow-to-change nonstationarity in further research.

Figure 1 .
Figure 1.Flowchart of nonstationary design flood estimation in response to climate change, population growth and cascade reservoir regulation.
, ..., , ) n x y y y z should satisfy the principle of the water balance equation [33]:

Figure 2 .
Figure 2. Sketch diagram of the flood regional composition with cascade reservoirs.
, ..., , ) n x y y y z in compliance with the water balance equation.

Figure 3 .
Figure 3. (a) Map of the Han River basin and (b) sketch diagram of the flood regional composition with AK and DJK cascade reservoirs.

Figure 4 .
Figure 4. Evolution of the population in Han River basin during 1950-2100.

Figure 6 .
Figure 6.Projected annual precipitation under (a) RCP 4.5 and (b) RCP 8.5 scenarios for 2021-2095.The bold red lines are the arithmetic mean value of the 10 GCMs.

Table 5 .
Results of trend test, change-point detection and multiple change point analysis for the W15 series for the years 1954-2020 (* means the maximum of the absolute value of the vector U).

Figure 7 .
Figure 7.The W15 series with the fitted linear trend lines of (a) AK and (b) DJK reservoirs.

Figure 8 .
Figure 8. Diagnostic plots (worm plot and centile curves) for evaluating the goodness-of-fit of the optimal nonstationary Log-normal distribution for the W15 series in the AK and DJK reservoirs.

Figure 9 .
Figure 9.Comparison of design W15 with different return periods and scenarios for the DJK reservoir.

Figure 10 .
Figure 10.Design flood hydrographs of 1000-year return period at downstream under three scenarios with (dash line) or without (solid line) cascade reservoir regulation.

)
The slow-to-change impacts of climate change and population growth are presented in design flood volumes according to the nonstationary flood frequency analysis method and ER criteria.The long-lasting driving factors imply the larger risks of the flood hazard.The 1000-year design W15 of the DJK reservoir under the stationary distribution scenario (S1), RCP 4.5-based (S2), and the RCP 8.5-based (S3) nonstationary scenario are 20.041,21.680 and 22.352 billion m 3 , respectively.

Table 1 .
The univariate distributions used to model the flood series in this study.

Table 3 .
Characteristic parameter values of the cascade reservoirs in the Han River basin.

Table 4 .
Description of the 10 GCMs used in this study.

Table 6 .
Distribution parameters and BIC values for W15 at the AK and DJK reservoirs.

Table 7 .
The detailed forward procedure of Log-normal distribution for W15 at the DJK reservoir.

Table 8 .
Dependence parameters and goodness-of-fit results for candidate copulas in nonstationary modeling.

Table 9 .
The FRC results derived by the MLRC (S1 scenario) and NS-MLRC (S2 and S3 scenarios) methods during design lifespan with the same flood prevention standards based on ER.X1, Y1 and X2 represent the designed flood volume W15 at the AK reservoir, inter-basin, and the DJK reservoir, respectively.

Table 10 .
Comparison of 1000-year design floods under three scenarios with or without reservoir regulation.

Table 11 .
Comparison of the highest water level derived from both the most likely FRC and the worst FRC during reservoir lifespan at the DJK reservoir for different scenarios.