Nonstationary Flood Hazard Analysis in Response to Climate Change and Population Growth

The predictions of flood hazard over the design life of a hydrological project are of great importance for hydrological engineering design under the changing environment. The concept of a nonstationary flood hazard has been formulated by extending the geometric distribution to account for time-varying exceedance probabilities over the design life of a project. However, to our knowledge, only time covariate is used to estimate the nonstationary flood hazard over the lifespan of a project, which lacks physical meaning and may lead to unreasonable results. In this study, we aim to strengthen the physical meaning of nonstationary flood hazard analysis by investigating the impacts of climate change and population growth. For this purpose, two physical covariates, i.e., rainfall and population, are introduced to improve the characterization of nonstationary frequency over a given design lifespan. The annual maximum flood series of Xijiang River (increasing trend) and Weihe River (decreasing trend) are chosen as illustrations, respectively. The results indicated that: (1) the explanatory power of population and rainfall is better than time covariate in the study areas; (2) the nonstationary models with physical covariates possess more appropriate statistical parameters and thus are able to provide more reasonable estimates of a nonstationary flood hazard; and (3) the confidences intervals of nonstationary design flood can be greatly reduced by employing physical covariates. Therefore, nonstationary flood design and hazard analysis with physical covariates are recommended in changing environments.


Introduction
The conventional flood frequency analysis has served as the primary tool to estimate design flood and provide information for flood hazard management. The fundamental assumption of the conventional flood frequency analysis is that the annual maximum flood series (AMFS) is independent and identically distributed (iid), also known as the stationary assumption. However, this assumption has been challenged by climate change and human activities [1][2][3][4][5][6][7][8][9][10][11][12]. Thus, the conventional stationary flood frequency analysis should be adapted to nonstationary conditions in changing environments.
In recent decades, the nonstationary flood frequency analysis theory has been developed to address this issue. In the framework of nonstationary flood frequency analysis, the statistical parameters of the time-varying probability distribution model (TVPD) are typically modelled as a function of covariates. Therefore, the evolution of future flood distribution heavily relies on the projections of covariates. Theoretically, if the covariates are able to represent hydrological processes, the related

Time-Varying Probability Distribution Model
Researchers have documented that probability distributions in flood frequency analysis can be classified into four families: the normal family (e.g., normal, lognormal), the Pearson type III family (e.g., gamma, Pearson type III), the general extreme value (GEV) family (e.g., GEV, Gumbel), and the generalized Pareto family [37,38]. In this study, we select the lognormal (LN), gamma (GA), Gumbel (GU), and GEV distributions to represent the normal, Pearson type III and GEV families. Generalized Additive Models in Location, Scale and Shape (GAMLSS) [39] are employed to build the time-varying probability models (TVPD). In GAMLSS, for a nonstationary flood series z t (t = 1, . . . , k), the probability distribution of TVPD, denoted by G Z,t (z θ t ) , varies in different years, since the distribution type G(·) is invariant, but the statistical parameters θ t , i.e., the time-varying location parameter µ t and scale parameter σ t , are expressed as function of the covariates I t j ( j = 1, 2, . . . , n) by h(µ t ) = α 0 + n j=1 α j I t j , where h(·) represents the link function of statistical parameters, and α = (α 0 , . . . , α n ) T and β = (β 0 , . . . , β n ) T are the parameters used for describing µ t and σ t , respectively. In this study, µ t and σ t are estimated using the maximum likelihood method. It must be mentioned that, when constructing the nonstationary GEV models, the estimates of time-varying shape parameters are sensitive to the data sets and may contain a high uncertainty. Hence, the shape parameter of GEV models is kept to be constant as done by Lima et al. [24] and Du et al. [32]. Annual total rainfall (rain) and population (pop) are employed as covariates to represent the climate factors and urbanisation. In calculating future exceedance probabilities, general circulation models (GCM) are used to project rain under representative concentration pathway (RCP) 4.5 and 8.5, and the logistic growth equation is used to project the growth of future pop. See Yan et al. [19] for detailed information about the projection of future rain and pop.

Model Selection and Goodness-of-Fit Test
In this study, different kinds of nonstationary models are developed considering the combination of different candidate probability distributions and variation types of statistical parameters. Thus, the Akaike Information Criterion (AIC) [40] is used to determine the optimal nonstationary model, which is given by where ζ is the maximized likelihood value and ρ is the total number of independently adjusted parameters of the model. The lower the AIC value is, the better the performance of the nonstationary model. In addition to AIC value, two kinds of diagnostic plots, i.e., the worm plot [41] and the centile curves plot, are also employed to assist the selection of optimal models and the diagnosis of the fitting quality of the optimal models. The worm plot is also known as the detrended Q-Q plot, in which the x-axis represents the standard normal quantiles, while the y-axis represents the difference value between empirical quantiles and standard normal quantiles. The centile curves plot employs five percentile curves, i.e., 5th, 25th, 50th, 75th and 95th, to assist the visual inspection of probabilistic coverage below different percentiles.

Flood Hazard Analysis
In the flood management decision-making, flood hazard over a given design life is important information to communicate event likelihood in changing environments. Under stationary conditions, for a design life starting from T 1 to T 2 , stationary flood risk R s is calculated by where p is the exceedance probability. While under nonstationary conditions, the nonstationary flood hazard that a flood event exceeds the design value z q (m) for return period m is denoted by R ns , which is given by [42,43] R ns = 1 − where p t is the time-varying exceedance probability, and G Z,t (z θ t ) is the time-varying probability distribution.

Average Design Life Level (ADLL)
The ADLL method was developed by Yan et al. [19]. Assume T 1 and T 2 are the starting year and ending year of a project, respectively, and the design life starting from T 1 to T 2 is denoted by T 1 − T 2 . Under nonstationary conditions, the annual average reliability over the design life period, denoted by RE ave , is defined as [42] RE ave In the ADLL method, it is assumed that, during the design life period T 1 − T 2 of a project, the annual average reliability for a design value z q under nonstationary conditions should be equal to the yearly reliability 1-1/m under stationary condition, i.e., RE ave In the application of the ADLL method, one should first specify a return period m and obtain the value on the right side of Equation (6), and finally the design value can be obtained by solving this equation.

Uncertainty Analysis of Design Flood
In this study, the uncertainties of design floods estimated by the ADLL method are also estimated using the nonstationary nonparametric bootstrap method. See Obeysekera and Salas [44] and Yan et al. [19] for detailed information about the nonstationary nonparametric bootstrap method.

Study Areas
The Xijiang River is located between the geographical coordinates 21 • 31 -26 • 49 N and 102 • 14 -112 • 48 E ( Figure 1) and it is the mainstream of the Pearl River, southeast China. The total length of the Xijiang River is about 2075 km and the drainage area of Xijiang River is approximately 353,100 km 2 , accounting for about 78% of the total drainage area of Pearl River basin. The average daily discharge of the Xijiang River is collected from the Dahuangjiangkou station. Being influenced by the subtropical climate, the average annual rainfall ranges from 1250 to 1750 mm, but the rainfall from April to September accounts for about 70-80% of the total rainfall. Since the 1950s, with the construction of flood control projects in the XRB, and the accelerating urbanisation process, the flood generation condition has been changed. The population of XRB has increased from about 20 million in 1954 to 50 million in 2009.
The Weihe River is the longest tributary of the Yellow River, north China, and it is located between the geographical coordinates 33 • 40 -37 • 26 N and 103 • 57 -110 • 27 E (Figure 1). The WRB has an approximate drainage area of 134,800 km 2 . The average daily discharge of the Weihe River is collected from the Huaxian station, which accounts for 79% of the total drainage area of WRB. WRB is influenced by the typical temperate continental monsoon climate, and the average annual rainfall of WRB is about 544 mm over the period 1951-2012. The WRB is the major source of water supply for the state key economic development zone, i.e., the Guanzhong Plain. Over the past decades, the total water consumption for industrial, agriculture and residential use has significantly increased. The population of WRB has doubled from 15 million in 1960 to 30 million in 2012. In addition, many reservoirs and soil and water conservation projects have been constructed in the WRB since the 1950s.
Being influenced by accelerating urbanisation and climate change, in recent decades, the nonstationarity issue of AMFS for both PRB and WRB has been reported by many researchers [6,18,37,[45][46][47]. In this study, the AMFS of Dahuangjiangkou hydrological station from XRB and Huaxian hydrological station from WRB are selected for illustration purposes ( Figure 2). The AMFS of Dahuangjiangkou station exhibits significant increasing trends while the AMFS of Huaxian station exhibits significant decreasing trends based on the Mann-Kendall test at the 0.05 significance value (Table 1).

Preliminary Data Analysis
The observed AMFS were collected from the Dahuangjiangkou and Huaxian hydrological stations, respectively. The locations and information of the two stations are shown in Figure 1 and Table 1, respectively.
Two kinds of hydrometeorological and socio-economic covariates, i.e., annual total rainfall (rain) and population (pop) are employed as covariates to model the variation of nonstationary AMFS of Dahuangjiangkou and Huaxian (Figure 3), since both climate factors and anthropogenic impacts are thought to be potential influencing factors of the nonstationary AMFS [35,48,49]. In addition, correlation tests were also conducted to examine the correlation between AMFS and pop/rain covariates from the statistical perspective, respectively. As shown from Figure 4, the Pearson correlation coefficient (PCC) between AMFS and rain is higher than that between AMFS and pop. For Huaxian, the PCCs between AMFS and rain/pop are 0.64/0.52, while, for Dahuangjiangkou, the PCCs between AMFS and rain/pop are 0.49/0.33. The observed daily total rainfall of the selected stations is obtained from the National Climate Center of China Meteorological Administration [50]. These data are aggregated to basin scale using the Thiessen polygon method for Dahuangjiangkou and Huaxian stations, respectively, and then the annual total rainfall series was calculated. The population data is employed to examine the combined impacts of urbanisation and water consumption on the nonstationarity of AMFS as done by Yan et al. [19] and Villarini et al. [36]. The population data of watersheds controlled by Huaxian and Dahuangjiangkou stations are collected from the cities located within them. The population data is provided by the Shaanxi Provincial Bureau of Statistics

Preliminary Data Analysis
The observed AMFS were collected from the Dahuangjiangkou and Huaxian hydrological stations, respectively. The locations and information of the two stations are shown in Figure 1 and Table 1, respectively.
Two kinds of hydrometeorological and socio-economic covariates, i.e., annual total rainfall (rain) and population (pop) are employed as covariates to model the variation of nonstationary AMFS of Dahuangjiangkou and Huaxian (Figure 3), since both climate factors and anthropogenic impacts are thought to be potential influencing factors of the nonstationary AMFS [35,48,49]. In addition, correlation tests were also conducted to examine the correlation between AMFS and pop/rain covariates from the statistical perspective, respectively. As shown from Figure 4, the Pearson correlation coefficient (PCC) between AMFS and rain is higher than that between AMFS and pop. For Huaxian, the PCCs between AMFS and rain/pop are 0.64/0.52, while, for Dahuangjiangkou, the PCCs between AMFS and rain/pop are 0.49/0.33. The observed daily total rainfall of the selected stations is obtained from the National Climate Center of China Meteorological Administration [50]. These data are aggregated to basin scale using the Thiessen polygon method for Dahuangjiangkou and Huaxian stations, respectively, and then the annual total rainfall series was calculated. The population data is employed to examine the combined impacts of urbanisation and water consumption on the nonstationarity of AMFS as done by Yan et al. [19] and Villarini et al. [36]. The population data of watersheds controlled by Huaxian and Dahuangjiangkou stations are collected from the cities located within them. The population data is provided by the Shaanxi Provincial Bureau of Statistics   The projected future rainfall and population of Huaxian station are provided by Yan et al. [19] and the future rainfall and population of Dahuangjiangkou station are projected using statistical downscaling method and Verhulst logistic growth equation, respectively, as done by Yan et al. [19]. It should be mentioned that the Verhulst logistic growth equation accommodates the growth restriction resulting from limited natural resources and it has been widely used to predict the growth of population. See Yan et al. [19] and Tsoularis and Wallace [53] for details. The projected future rain under RCP 4.5 (medium emission) and RCP 8.5 (high emission) scenarios and future pop are presented in Figures 5 and 6, respectively.   The projected future rainfall and population of Huaxian station are provided by Yan et al. [19] and the future rainfall and population of Dahuangjiangkou station are projected using statistical downscaling method and Verhulst logistic growth equation, respectively, as done by Yan et al. [19]. It should be mentioned that the Verhulst logistic growth equation accommodates the growth restriction resulting from limited natural resources and it has been widely used to predict the growth of population. See Yan et al. [19] and Tsoularis and Wallace [53] for details. The projected future rain under RCP 4.5 (medium emission) and RCP 8.5 (high emission) scenarios and future pop are presented in Figures 5 and 6, respectively. The projected future rainfall and population of Huaxian station are provided by Yan et al. [19] and the future rainfall and population of Dahuangjiangkou station are projected using statistical downscaling method and Verhulst logistic growth equation, respectively, as done by Yan et al. [19]. It should be mentioned that the Verhulst logistic growth equation accommodates the growth restriction resulting from limited natural resources and it has been widely used to predict the growth of population. See Yan et al. [19] and Tsoularis and Wallace [53] for details. The projected future rain under RCP 4.5 (medium emission) and RCP 8.5 (high emission) scenarios and future pop are presented in Figures 5 and 6, respectively.

Nonstationary Flood Frequency Analysis
When using time covariates, a number of nonstationary models were built considering the different combinations of candidate distributions and variation types. The optimal model was selected based on the AIC values (Table 2). When using physical covariates, i.e., rain and pop, a total of 64 models were built for Dahuangjiangkou station and Huaxian station (Figure 7). The optimal models with either time or rain and pop are summarized in Table 3. It can be seen that the AIC values of the optimal models with physical covariates are smaller than those with time covariate. Figures 8  and 9 present the goodness-of-fit of the optimal nonstationary model with time and physical covariates. The more the scatters in the worm plot located along the 0 value of the y-axis, the better is the model performance. For both stations, all scatter points in the worm plots are within the 95%

Nonstationary Flood Frequency Analysis
When using time covariates, a number of nonstationary models were built considering the different combinations of candidate distributions and variation types. The optimal model was selected based on the AIC values (Table 2). When using physical covariates, i.e., rain and pop, a total of 64 models were built for Dahuangjiangkou station and Huaxian station (Figure 7). The optimal models with either time or rain and pop are summarized in Table 3. It can be seen that the AIC values of the optimal models with physical covariates are smaller than those with time covariate. Figures 8  and 9 present the goodness-of-fit of the optimal nonstationary model with time and physical covariates. The more the scatters in the worm plot located along the 0 value of the y-axis, the better is the model performance. For both stations, all scatter points in the worm plots are within the 95%

Nonstationary Flood Frequency Analysis
When using time covariates, a number of nonstationary models were built considering the different combinations of candidate distributions and variation types. The optimal model was selected based on the AIC values (Table 2). When using physical covariates, i.e., rain and pop, a total of 64 models were built for Dahuangjiangkou station and Huaxian station (Figure 7). The optimal models with either time or rain and pop are summarized in Table 3. It can be seen that the AIC values of the optimal models with physical covariates are smaller than those with time covariate. Figures 8 and 9 present the goodness-of-fit of the optimal nonstationary model with time and physical covariates. The more the scatters in the worm plot located along the 0 value of the y-axis, the better is the model performance. For both stations, all scatter points in the worm plots are within the 95% confidence intervals (Figures 8a,b and 9a,b), but the trends of worm plots with physical covariates are gentler than those with time covariate and more scatters located around the 0 value, indicating better agreement between the nonstationary model and observations. Water 2019, 11, x FOR PEER REVIEW 9 of 20 confidence intervals (Figures 8a,b and 9a,b), but the trends of worm plots with physical covariates are gentler than those with time covariate and more scatters located around the 0 value, indicating better agreement between the nonstationary model and observations. As for centile curves, for Huaxian station, the percentages of observation points below the 5th, 25th, 50th, 75th and 95th centile curves are 3.2%, 33.9%, 45.2%, 69.4% and 95.2% using time covariate, while these percentages are 4.8%, 25.8%, 50.0%, 79.0% and 93.5% using physical covariates ( Figure  8c,d). For Dahuangjiangkou station, the percentages of observation points below the 5th, 25th, 50th, 75th and 95th centile curves are 3.7%, 27.8%, 44.4%, 72.2% and 98.1% using time covariates, while these percentages are 9.3%, 25.9%, 44.4%, 70.4% and 94.4% using physical covariates (Figure 9c,d). In addition, to assist the evaluation of the quantile interval bounds, the Root Mean Squared Error (RMSE) between the 50th centile curve and the observation data was also calculated. For Huaxian station, the RMSE is 1108.9 when using time covariate, while the RMSE is 965.8 when using physical covariates. For Dahuangjiangkou station, the RMSE for using time and physical covariates are 6983.1 and 6216.3, respectively. As for centile curves, for Huaxian station, the percentages of observation points below the 5th, 25th, 50th, 75th and 95th centile curves are 3.2%, 33.9%, 45.2%, 69.4% and 95.2% using time covariate, while these percentages are 4.8%, 25.8%, 50.0%, 79.0% and 93.5% using physical covariates (Figure 8c,d). For Dahuangjiangkou station, the percentages of observation points below the 5th, 25th, 50th, 75th and 95th centile curves are 3.7%, 27.8%, 44.4%, 72.2% and 98.1% using time covariates, while these percentages are 9.3%, 25.9%, 44.4%, 70.4% and 94.4% using physical covariates (Figure 9c,d).
In addition, to assist the evaluation of the quantile interval bounds, the Root Mean Squared Error (RMSE) between the 50th centile curve and the observation data was also calculated. For Huaxian station, the RMSE is 1108.9 when using time covariate, while the RMSE is 965.8 when using physical covariates. For Dahuangjiangkou station, the RMSE for using time and physical covariates are 6983.1 and 6216.3, respectively. Table 2. Akaike Information Criterion (AIC) values of the nonstationary models with time covariate. Note that letter "L" in the models' names represents location parameter µ and "S" represents scale parameters σ. Number "0" means the parameter is invariant while "t" means the parameter varies with time covariate. In addition, the AIC value in bold is the optimal model for each station.  Note that letter "L" in the models' names represents location parameter μ and "S" represents scale parameters σ. Number "0" means the parameter is invariant while "t" means the parameter varies with time covariate. In addition, the AIC value in bold is the optimal model for each station.   compared with the expected 95% using the cross-validation method, indicating that the selected models are well calibrated. Overall, the above results indicated that the selected optimal models are able to achieve satisfactory performance in modeling the variability of the observed flood events. In addition, the general performance of nonstationary models with physical covariates is much better than those with time covariate.  In addition to the centile curves, cross-validation procedure is also conducted for assessing the goodness-of-fit of the optimal models. See Lima et al. [24] for the standard procedure of the cross-validation method. As shown in Figure 10, the nonstationary cross-validated 50th centiles using physical covariates are able to better reproduce the observed AMFS, with the Pearson correlation coefficient being 0.53 and 0.67, respectively. In addition, the probabilistic coverage of the optimal nonstationary model of Dahuangjiangkou and Huaxian station are 90.7% and 90.3%, respectively, compared with the expected 95% using the cross-validation method, indicating that the selected models are well calibrated. Overall, the above results indicated that the selected optimal models are able to achieve satisfactory performance in modeling the variability of the observed flood events. In addition, the general performance of nonstationary models with physical covariates is much better than those with time covariate.

Flood Hazard Analysis for PRB and WRB
The stationary and nonstationary flood hazard of Dahuangjiangkou and Huaxian stations that a future flood event exceeds design flood q z corresponding to return period [2,100] m ∈ and design lives varying from 1 to 50 years were calculated using Equations (3) and (4), respectively (Figures 11 and 12). For Dahuangjiangkou station, it was found that the nonstationary flood hazard calculated using time covariates was always larger than the stationary flood hazard, since the future exceedance probabilities t p in Equation (4) were getting monotonically larger with time covariates, while the cases were complicated for the nonstationary flood hazard calculated using pop and rain covariates because t p grew in fluctuation with pop and rain covariates ( Figure 11). In addition, the nonstationary flood hazard calculated using pop and rain covariates was larger than that calculated using time covariate for [2,30]  indicating that the flood hazard in the Huaxian basin will decrease in future periods. In addition, we found that the edges of nonstationary flood hazards with physical covariates were not as smooth as that with time covariates because of the different variation types of covariates.

Flood Hazard Analysis for PRB and WRB
The stationary and nonstationary flood hazard of Dahuangjiangkou and Huaxian stations that a future flood event exceeds design flood z q corresponding to return period m ∈ [2, 100] and design lives varying from 1 to 50 years were calculated using Equations (3) and (4), respectively (Figures 11 and 12). For Dahuangjiangkou station, it was found that the nonstationary flood hazard calculated using time covariates was always larger than the stationary flood hazard, since the future exceedance probabilities p t in Equation (4) were getting monotonically larger with time covariates, while the cases were complicated for the nonstationary flood hazard calculated using pop and rain covariates because p t grew in fluctuation with pop and rain covariates ( Figure 11). In addition, the nonstationary flood hazard calculated using pop and rain covariates was larger than that calculated using time covariate for m ∈ [2,30]. For Huaxian station, the nonstationary flood hazard estimated by time or physical covariates were all smaller than stationary flood hazard because the AMFS were decreasing and the exceedance probabilities p t in Equation (4) were getting smaller (Figure 12), indicating that the flood hazard in the Huaxian basin will decrease in future periods. In addition, we found that the edges of nonstationary flood hazards with physical covariates were not as smooth as that with time covariates because of the different variation types of covariates.

Design Flood and Associated Uncertainty Using ADLL
In this study, we assume that there will be a hydrological structure to be in service for 50 years from 2015 to 2064 for illustration purposes. The ADLL method was employed to estimate the design floods for Huaxian station (AMFS with decreasing trend) and Dahuangjiangkou station (AMFS with increasing trend) based on the optimal nonstationary probability distributions with either time or physical covariates (pop and rain) ( Table 3). In addition, for the purpose of providing a comprehensive comparison between stationary and nonstationary design floods estimated using time or physical covariates, the 95% confidence intervals (CIs) of stationary and nonstationary design floods were estimated using the nonstationary nonparametric bootstrap method.
As shown in Figure 13, for Dahuangjiangkou station, the nonstationary design floods were larger than stationary design floods except for design floods with higher return periods using physical covariates. As for CIs, it was found that the CIs of nonstationary design floods were five times as large as those of stationary design floods when using time covariates, while the CIs were greatly reduced when using physical covariates. For Huaxian station, the nonstationary design floods

Design Flood and Associated Uncertainty Using ADLL
In this study, we assume that there will be a hydrological structure to be in service for 50 years from 2015 to 2064 for illustration purposes. The ADLL method was employed to estimate the design floods for Huaxian station (AMFS with decreasing trend) and Dahuangjiangkou station (AMFS with increasing trend) based on the optimal nonstationary probability distributions with either time or physical covariates (pop and rain) ( Table 3). In addition, for the purpose of providing a comprehensive comparison between stationary and nonstationary design floods estimated using time or physical covariates, the 95% confidence intervals (CIs) of stationary and nonstationary design floods were estimated using the nonstationary nonparametric bootstrap method.
As shown in Figure 13, for Dahuangjiangkou station, the nonstationary design floods were larger than stationary design floods except for design floods with higher return periods using physical covariates. As for CIs, it was found that the CIs of nonstationary design floods were five times as large as those of stationary design floods when using time covariates, while the CIs were greatly reduced when using physical covariates. For Huaxian station, the nonstationary design floods were always smaller than stationary design floods. In addition, the CIs of nonstationary design floods were dramatically reduced when using physical covariates, which were even comparable with those of stationary design floods. The much larger CIs of the design flood estimated by time covariate is mainly due to the infinite extrapolation of time covariate in far future. Therefore, physical covariates that can provide reliable future prediction products are preferred in nonstationary frequency analysis in the study areas. Table 3. Summary of the optimal nonstationary models with time and physically-based covariates. µ t and σ t are the time-varying location parameter and scale parameter, respectively. Note that letter "L" in the models' names represents location parameter and "S" represents scale parameters. Number "0" denotes the parameter is invariant, while "t", "p", "r" and "pr" denote that the parameter varies with time, pop, rain and both pop and rain, respectively.

Dahuangjiangkou Huaxian
Covariates were always smaller than stationary design floods. In addition, the CIs of nonstationary design floods were dramatically reduced when using physical covariates, which were even comparable with those of stationary design floods. The much larger CIs of the design flood estimated by time covariate is mainly due to the infinite extrapolation of time covariate in far future. Therefore, physical covariates that can provide reliable future prediction products are preferred in nonstationary frequency analysis in the study areas. Table 3. Summary of the optimal nonstationary models with time and physically-based covariates. t μ and t σ are the time-varying location parameter and scale parameter, respectively. Note that letter "L" in the models' names represents location parameter and "S" represents scale parameters. Number "0" denotes the parameter is invariant, while "t", "p", "r" and "pr" denote that the parameter varies with time, pop, rain and both pop and rain, respectively.

Discussion
Reliable prediction of future flood hazard over the design lifespan of a hydrological project is very important for local flood alert and hazard management in changing environments. It can be seen from the results of flood hazard analysis that the nonstationary flood hazard calculated using time covariate is prone to exaggerate the current trend of flood series; that is to say, it could make the flood hazard more serious for increasing flood series and make the flood hazard less serious for decreasing flood series, compared with cases using pop and rain covariates. As for the difference between the estimated stationary and nonstationary flood hazard, when using a time covariate, the nonstationary flood hazard is always larger or smaller than a stationary flood hazard for increasing and decreasing cases, respectively, while, when using pop and rain covariates, the difference between stationary and nonstationary flood hazard depends on the return period and the design lifespan of a project. Finally, when it comes to the estimated distribution parameters, the monotonous variation type of time covariate would definitely lead to monotonous variation of expectation and variance of the AMFS derived from the fitted location and scale parameters (Table 3). However, a noticeable fluctuation movement of AMFS can be observed in both stations (Figure 2a,b), which contradicts the monotonous variation patterns obtained by using time covariate, while, when using pop and rain covariates, the fitted parameters are changing in fluctuation with physical covariates, which is thought to provide better model performance and more reasonable statistical parameters. The findings of this study are consistent with other similar studies such as Condon et al. [28] and Du et al. [32], who demonstrated the evolution of future flood or low-flow hazard with only climate factors, and are different with them in terms of physical meanings, model performance and the three-dimensional flood hazard profile plot with different return periods and design life. Thus, the flood hazard analysis over the design life of hydrological projects based on the nonstationary models using physical covariates can provide more reasonable and useful information for hydrological communities and policy makers, to predict future evolution of flood hazards within the lifespan of hydrological projects.
From the perspective of design flood estimation, we also found that, for both stations, when using time covariates, the CIs of nonstationary design floods were extremely large compared with those of stationary design floods. While when using physical covariates, the CIs can be dramatically reduced for both stations, which is very important for practical application of nonstationary design strategy and thus provide an alternative for hydrological communities and policy makers when addressing the nonstationary issue.

Conclusions
The flood hazard analysis is important for flood risk reduction in changing environments. The main objective of this study is to strengthen the significance of employing physical covariates in the flood hazard analysis in changing environments. For this purpose, we compared nonstationary flood hazard calculated using time and physical covariates, i.e., pop and rain. In addition, nonstationary design floods estimated using either time covariate or pop and rain covariates were also compared. The main findings are drawn as follows: (1) The AMFS of Dahuangjiangkou and Huaxian station have exhibited significant increasing or decreasing trends due to climate change and human activities. For Dahuangjiangkou, the nonstationary lognormal model with both location and scale parameters varying with population and rainfall is the optimal model, whereas, for Huaxian station, the nonstationary lognormal model with only location parameters varying with population and rainfall is the optimal model. For both stations, the performance of the optimal nonstationary models with physical covariates is better than those with a time covariate based on the AIC values and diagnostic plots, probably indicating the better explanatory power of the physical covariates in the study areas.
(2) Different covariates employed in the construction of nonstationary models will lead to different results of nonstationary flood hazard analysis. When using time covariates, the nonstationary flood hazard is always larger or smaller than the stationary flood hazard for increasing or decreasing flood series, indicating that the current flood problem will be deteriorated even further or alleviated over time. However, the above results may exaggerate future flood hazard due to inappropriate statistical parameters of the constructed nonstationary models, which will change monotonously and indefinitely over time. In contrast, the nonstationary distributions with population and rainfall covariates have more appropriate statistical parameters and thus are able to provide more reasonable nonstationary flood risk since future exceedance probabilities are changing in fluctuation with physical covariates. Thus, nonstationary flood frequency analysis and hazard analysis with physical covariates are recommended in changing environments. (3) The use of different covariates in nonstationary frequency analysis will lead to different results of nonstationary design flood and associated CIs. The nonstationary design floods are always larger than stationary ones for increasing flood series, but the gap between stationary and nonstationary design flood can be lessened when using physical covariates. As for CIs, it is clear that the CIs of nonstationary design flood can be greatly reduced by using physical covariates regardless of whether the trend of flood series is increasing or decreasing.
Following the aforementioned conclusions, two major comments concerning the applications of nonstationary flood hazard analysis and engineering design are made as follows: Firstly, the selection of covariates is one of the most challenging issues in nonstationary frequency analysis. We should provide reliable future projection of employed covariates in nonstationary models to predict future nonstationary distributions. Two selection requirements should be satisfied for the covariates employed for nonstationary frequency analysis: (i) the covariates should own sufficient explanatory power to describe the nonstationarity of flood series; and (ii) the covariates must be projected reliably in future periods [19]. In this study, pop and rain were chosen as covariates to investigate the impacts of human activities and climatic factors, and they can also be projected using population growth models and climate models, respectively. However, it should be emphasized that pop is regarded as an indirect reflection of the intensity of human activities and can just be regarded as a simple characterization of urbanization in the developing countries where more population to some extent means more imperviousness. In the future, it is necessary to try more covariates and select the ones that are most influential on extreme floods, such as the tropical storms and land use changes, which can reflect the changing process of flood events more directly. For this purpose, efforts should also be made to improve the reliability of future projections of climate indices and land use changes in the future studies.
Secondly, only two stations showing opposite trends were selected as study areas in this study for illustration purposes. The limited number of case studies is a considerable limitation to conclude on the superiority of the nonstationary models with physical covariates. Thus, further studies are needed to improve the present work and demonstrate the superiority of physical covariates.