Selection of an Optimal Distribution Curve for Non-Stationary Flood Series

: The stationarity assumption of hydrological processes has long been compromised by human disturbances in river basins. The traditional hydrological extreme-value analysis method, i


Introduction
Under the influences of climate change and human activities, significant changes are occuring to the hydrological regimes of many rivers, and the occurrence frequency of extreme hydrological events has also significantly increased [1].Schiermeier (2011) pointed out that the rainfall intensity in the northern hemisphere showed an increasing trend in the late 20th century, leading to an increase in flood risk [2].Barnett et al. (2008) found that in the second half of the 20th century, significant changes had occurred to the water cycle in the western USA, 60% of which were caused by human-induced environmental changes [3].In the context of global warming, the frequency and intensity of extreme weather events in the past 50 years showed a strengthening trend in China, mainly reflected by the extent and frequency of floods and other extreme events [4].
It was pointed out by Milly et al. (2008) that stationarity is dead because substantial anthropogenic change of the Earth's climate is altering the means and other statistics of precipitation, evapotranspiration, and rates of discharge of rivers [5].Many hydraulic structures designed by the conventional hydrologic analysis methods, which have stationarity assumption, will face the distortion risk in design frequency caused by the ever-changing environment [6,7].Compared with hydrologic frequency analysis method under the stationary assumption, which has a history of approximately 100 years, the flood frequency analysis of non-stationary hydrological series is an emerging research field, and recent research has continued to show a rapid development in this line [6][7][8][9][10].
The non-stationary flood frequency analysis methods mainly include several aspects, i.e., time-varying moments, pooled flood frequency analysis, local likelihood, quantile regression and mixture distribution model based on conditional probability [8,11].Different methods of non-stationary flood frequency analysis were suggested for eliminating the trend in the time series through one or more processing steps [8].One of the more general methods of non-stationary flood frequency analysis has been the streamflow naturalization methodology, but it needs to collect a lot of data [4].With the deepening of research on the influences of environmental change on hydrologic processes, research on hydrologic frequency analysis under non-stationary conditions has been receiving special attention in recent years [12][13][14][15][16][17].The method of time-varying moments (TVM) has become a commonly used method of non-stationary flood frequency analysis, because the TVM method can process the whole non-stationary flood series directly [18].The time-varying parameters can be estimated with a combination of the analysis on trend, periodicity and other deterministic components [13].It was found in the study of Diehl and Potter (1987) that the samples of extreme value series used for frequency analysis do not always originate from the same population, i.e., not subject to the same distribution, so that the mixture distribution model was proposed directly based on extreme sample series with non-identical distribution for frequency analysis [19].In addition, some studies have aimed to establish a correlation between environmental factors and probability distribution parameters of hydrologic extreme value, and deduce the time-varying laws based on the variation laws of these factors [20,21].
A frequency analysis model of non-stationary extreme value series proposed by Strupczewski et al. (2001) embedded the trend component into the distribution of first and second-order moments (time-varying moments), i.e., the trends of mean and variance of statistical parameters were considered [12].Recently, there have been many studies on TVM model.The study of Vogel et al. (2011) showed that the flood magnitude in many regions of the U.S. had an increasing trend, and floods with return periods of 100 years might become more common in some basins [22].The generalized addictive model for Location, Scale and Shape parameters (GAMLSS) was used to study the flood peak flow series in the Little Sugar Creek Basin with a higher degree of urbanization by Villarini et al. (2009) [23].It is shown that GAMLSS is able to describe the variability in the mean and variance of the annual maximum peak discharge by modeling the parameters of the selected parametric distribution as a smooth function of time via cubic splines.
Previous studies have shown that any violation of the basic hypotheses on independence and stationarity of observed data can significantly affect the accuracy of frequency analysis results [13,24,25].However, there are few studies that focus on flood frequency distribution curve selection for non-stationary flood series.With the aim to understand the scientific and technical issues involved in the non-stationary flood frequency analyses, the non-stationary series processing methods with TVM were compared with the traditional frequency analysis method without non-stationary series processing, to analyze the flood frequency characteristics in the changing environment of the East River basin.By the end of 2005, many reservoirs of different sizes were gradually completed in the catchment area above the control station (Boluo Station) in the East River basin.The adjustment coefficient of the large reservoirs on runoff at the Boluo Station reached 0.34 step by step.Besides, the vegetation cover continuously keeps decreasing due to urbanization process in this area.Therefore, continuous human activities and climate changes have led to hydrological change in the East River basin [26][27][28][29][30].The primary goal of this study was achieved through the following steps: (1) to reveal the mechanism of selecting the best fit distribution curves for non-stationary flood series from upstream to downstream that have been significantly affected by human activities and climate changes; and (2) to calculate and evaluate the impacts of the changing environment on the flood frequency distribution curve.Hydrological regimes in many basins in China and other developing countries are gradually changing due to climate change and human activities including the urbanization process and river regulations.The results of this study will be beneficial to rationally describe non-stationary extreme flood characteristics in the continuously changing environment of human-altered basins.

Study Area and Data
The East River basin is one of the main sub-basins in the Pearl River basin in south China (Figure 1), with a river length of 562 km and a total basin area of 27,040 km 2 .This basin belongs to a south subtropical monsoon climate zone.The water resources in the East River basin have been highly developed for urban water use, hydroelectric power generation and farmland irrigation, where the water use efficiency of the outer channel is more than 40%.The reservoirs in the East River basin control 51% of the entire basin area, and 12 cascade hydropower Stations have been constructed in the trunk channel [26].The forest is the major land use type, whereas the urban area shows a rapid expansion, especially in the downstream areas [31].
Atmosphere 2018, 9, x FOR PEER REVIEW 3 of 16 through the following steps: (1) to reveal the mechanism of selecting the best fit distribution curves for non-stationary flood series from upstream to downstream that have been significantly affected by human activities and climate changes; and (2) to calculate and evaluate the impacts of the changing environment on the flood frequency distribution curve.Hydrological regimes in many basins in China and other developing countries are gradually changing due to climate change and human activities including the urbanization process and river regulations.The results of this study will be beneficial to rationally describe non-stationary extreme flood characteristics in the continuously changing environment of human-altered basins.

Study Area and Data
The East River basin is one of the main sub-basins in the Pearl River basin in south China (Figure 1), with a river length of 562 km and a total basin area of 27,040 km 2 .This basin belongs to a south subtropical monsoon climate zone.The water resources in the East River basin have been highly developed for urban water use, hydroelectric power generation and farmland irrigation, where the water use efficiency of the outer channel is more than 40%.The reservoirs in the East River basin control 51% of the entire basin area, and 12 cascade hydropower Stations have been constructed in the trunk channel [26].The forest is the major land use type, whereas the urban area shows a rapid expansion, especially in the downstream areas [31].The data used in the study include daily streamflow from three hydrological Stations for 1954-2009, i.e., Longchuan Station in the upper reach of East River basin, Heyuan Station in the middle reach, and Boluo Station in the lower reach, which were provided by the Hydrological Bureau of Guangdong Province.At present, three big reservoirs have been completed in the East River basin, i.e., Xinfengjiang, Fengshuba, and Baipenzhu, with a total storage capacity of 17.06 billion m 3 .The construction time of main reservoirs in the East River basin, the capacity of the reservoirs, the hydrological Stations affected by reservoirs and other information are shown in Table 1.In addition, there are thirty-six reservoirs above medium size and more than 661 small-sized reservoirs in the basin, which have been constructed continuously since 1950s to the present.The construction time of these hydraulic engineering structures spans the whole study period.Daily precipitation data for 1954-2009 were collected from 38 national standard rain gauges in the East River basin.

Methodology
In a changing environment, the flood return period estimated before the change is often unable to describe the flood frequency characteristics after environmental changes.This means that it is necessary to introduce a time reference point in describing the concept of return period.
Based on the concept of TVM, a method proposed by Strupczewski et al. (2001) to calculate the non-stationary flood frequency is used in this study [12].This method considers the trend components of first and second order moments (mean m and standard deviation σ) of a probability distribution function.For the sake of completeness, the procedure of the method is briefly introduced in this section, and more detailed information of TVM method can refer to Strupczewski et al. (2001) and Yan et al. (2017a) [12,32].

Distributions Curve Probability Density Function Relationships of Parameters and First Two Order Moments
Note: For P3, GEV and GLO, the α, ξ and k are the scale, location and shape parameters of each distribution.
There are only scale and location parameters in GMB, without shape parameter.For LN2, µ y and σ y are respectively the mean and standard deviation of Y = ln X.The m and σ are the first and second moments (mean and standard deviation), which are related to the parameters of different distributions shown above.
For the three-parameter P3 distribution, the lower bound parameter ξ in the distribution was considered to remain unchanged over time in this study [22,32,37].GEV and GLO are also three-parameter distributions, however, these methods only consider the trend components of first two moments, i.e., among three parameters of the distributions, one parameter must not change over time [6,22].The coefficient of skewness Cs is used to describe the shape of a distribution.Cs represents the extent of symmetry of a random variable with respect to the mean value.The middle part of frequency distribution shifts to the left as Cs increases with a steeper large flood part on the left and a flatter small flood part on the right.Parameter k is the shape parameter of GLO and GEV distributions, and has a certain correlation with Cs.Among GEV and GLO distributions, no lower bound parameters of random variable x exist, and the lower bound of x is determined by three parameters, while values of C s and k have a fixed function relationship [38], i.e., C s = C s (k).Therefore, the shape parameter k correlated to C s in GEV and GLO was considered to remain unchanged over time in this study.

Trend Model
In order to ensure the independence of the hydrological sequence, "extreme" refers to the annual maximum daily flow in this study, and only one sample is taken in a year.Data series used for this study are the annual maximum daily flow of 1954-2009 at the Longchuan, Heyuan and Boluo Stations.Due to the limited length of the data series, only the trends of first two moments were considered in this study with respect to the following situations: There is also a basic situation that neither of the first two moments has a trend, i.e., the distribution parameters do not change over time as a stationary option (S).
Among A, B, and C trend types, two further options exist: linear trend (L) or parabolic type (P).
For Type D, only linear trend was considered (L) in this study with the purpose of keeping the number of parameters to a minimum.
Therefore, the trend model will have seven types, of which the expressions of first two moments and the increased number of model parameters are shown in Table 3 [22,32].
Note: m is the mean of the time series, σ is the standard deviation of the time series, t is time, here we express m and σ as (m = m(t; θ (m) ) and σ = σ(t; θ (σ) )), For Type D, only linear trend was considered (L) in this study with the purpose of keeping th number of parameters to a minimum.
Therefore, the trend model will have seven types, of which the expressions of first tw moments and the increased number of model parameters are shown in Table 3 [22,32].
Note: m is the mean of the time series, σ is the standard deviation of the time series, t is time, here we express m and σ as ( ),Ө is the parameter vector, Cv is the coefficient of dispersion.m0, am, bm, σ0, aσ, bσ, Cv are elements of Ө.In different trend models, the forms of Ө are not the same.The expressions of linear trend (L) and parabolic trend (P) are given.We can see that the more complex the model, the more parameters are needed.

Parameter Estimation
The maximum likelihood estimates (MLE) method was used for parameter estimation.Mor details regarding methods on parameter estimation can be found in Yan et al. (2017a) an Strupczewski et al. (2001) [12,32].

Selection of Optimal Model
AIC (Akaike Information Criterion) not only rewards goodness of fit, but also includes penalty, which is an increasing function of the number of estimated parameters.The AIC provides means of model selection and evaluation based on the principle of maximum entropy [18].In th study, the AIC is used as a criterion of goodness of fit.The calculation equation is shown as below: Where, ML is the maximum likelihood function, and k is the number of model parameters[39 The model with the minimum AIC value should be selected as the optimum model.

Non-stationarity Detection in Flood Series
The changing trends of annual maximum peak flow of the three stations are detected b Mann-Kendall (M-K) test (Table 4) [40].It is seen that a decreasing trend at the significance level o 0.05 or less (p < 0.05) exists in all three stations.Figure 2 is a visual illusion for describing th changing details, which shows that deluge events frequently occurred for the first two decade Afterwards, the flood series of Longchuan Station showed a significant downward trend until 200 and high flood events occurred between 2005 and 2007.The flood peak flows of Heyuan and Bolu Stations show a similar but relatively gentle changing pattern.Specifically, the annual maximum flood flow had a significant downward trend from 1954 to 2009.During the same period, th precipitation amounts in the Longchuan and Heyuan Stations have displayed a slightly decreasin tendency (Table 4).In contrast, the precipitation in the Boluo Station displayed a slightly increasin For Type D, only linear trend was considered (L) in this study with the purpose of keeping the number of parameters to a minimum.
Therefore, the trend model will have seven types, of which the expressions of first two moments and the increased number of model parameters are shown in Table 3 [22,32].),Ө is the parameter vector, Cv is the coefficient of dispersion.m0, am, bm, σ0, aσ, bσ, Cv are elements of Ө.In different trend models, the forms of Ө are not the same.The expressions of linear trend (L) and parabolic trend (P) are given.We can see that the more complex the model, the more parameters are needed.

Parameter Estimation
The maximum likelihood estimates (MLE) method was used for parameter estimation.More details regarding methods on parameter estimation can be found in Yan et al. (2017a) and Strupczewski et al. (2001) [12,32].

Selection of Optimal Model
AIC (Akaike Information Criterion) not only rewards goodness of fit, but also includes a penalty, which is an increasing function of the number of estimated parameters.The AIC provides a means of model selection and evaluation based on the principle of maximum entropy [18].In this study, the AIC is used as a criterion of goodness of fit.The calculation equation is shown as below: Where, ML is the maximum likelihood function, and k is the number of model parameters [39].The model with the minimum AIC value should be selected as the optimum model.

Non-stationarity Detection in Flood Series
The changing trends of annual maximum peak flow of the three stations are detected by Mann-Kendall (M-K) test (Table 4) [40].It is seen that a decreasing trend at the significance level of 0.05 or less (p < 0.05) exists in all three stations.Figure 2 is a visual illusion for describing the changing details, which shows that deluge events frequently occurred for the first two decades.Afterwards, the flood series of Longchuan Station showed a significant downward trend until 2004, and high flood events occurred between 2005 and 2007.The flood peak flows of Heyuan and Boluo Stations show a similar but relatively gentle changing pattern.Specifically, the annual maximum flood flow had a significant downward trend from 1954 to 2009.During the same period, the precipitation amounts in the Longchuan and Heyuan Stations have displayed a slightly decreasing tendency (Table 4).In contrast, the precipitation in the Boluo Station displayed a slightly increasing .In different trend models, the forms of Atmosphere 2018, 9, x FOR PEER REVIEW For Type D, only linear trend was considered (L) in this study with th number of parameters to a minimum.
Therefore, the trend model will have seven types, of which the moments and the increased number of model parameters are shown in Tab  ),Ө is the param coefficient of dispersion.m0, am, bm, σ0, aσ, bσ, Cv are elements of Ө.In differ forms of Ө are not the same.The expressions of linear trend (L) and parabol We can see that the more complex the model, the more parameters are needed

Parameter Estimation
The maximum likelihood estimates (MLE) method was used for par details regarding methods on parameter estimation can be found in Strupczewski et al. (2001) [12,32].

Selection of Optimal Model
AIC (Akaike Information Criterion) not only rewards goodness of penalty, which is an increasing function of the number of estimated param means of model selection and evaluation based on the principle of maxim study, the AIC is used as a criterion of goodness of fit.The calculation equa 2 ln 2

= − +
Where, ML is the maximum likelihood function, and k is the number The model with the minimum AIC value should be selected as the optimum

Non-stationarity Detection in Flood Series
The changing trends of annual maximum peak flow of the three Mann-Kendall (M-K) test (Table 4) [40].It is seen that a decreasing trend a 0.05 or less (p < 0.05) exists in all three stations.Figure 2 is a visual ill changing details, which shows that deluge events frequently occurred fo Afterwards, the flood series of Longchuan Station showed a significant dow and high flood events occurred between 2005 and 2007.The flood peak flo Stations show a similar but relatively gentle changing pattern.Specificall flood flow had a significant downward trend from 1954 to 2009.Durin precipitation amounts in the Longchuan and Heyuan Stations have displa tendency (Table 4).In contrast, the precipitation in the Boluo Station displa are not the same.The expressions of linear trend (L) and parabolic trend (P) are given.We can see that the more complex the model, the more parameters are needed.

Parameter Estimation
The maximum likelihood estimates (MLE) method was used for parameter estimation.More details regarding methods on parameter estimation can be found in Yan

Selection of Optimal Model
AIC (Akaike Information Criterion) not only rewards goodness of fit, but also includes a penalty, which is an increasing function of the number of estimated parameters.The AIC provides a means of model selection and evaluation based on the principle of maximum entropy [18].In this study, the AIC is used as a criterion of goodness of fit.The calculation equation is shown as below: where, ML is the maximum likelihood function, and k is the number of model parameters [39].
The model with the minimum AIC value should be selected as the optimum model.

Non-Stationarity Detection in Flood Series
The changing trends of annual maximum peak flow of the three stations are detected by Mann-Kendall (M-K) test (Table 4) [40].It is seen that a decreasing trend at the significance level of 0.05 or less (p < 0.05) exists in all three stations.Figure 2 is a visual illusion for describing the changing details, which shows that deluge events frequently occurred for the first two decades.Afterwards, the flood series of Longchuan Station showed a significant downward trend until 2004, and high flood events occurred between 2005 and 2007.The flood peak flows of Heyuan and Boluo Stations show a similar but relatively gentle changing pattern.Specifically, the annual maximum flood flow had a significant downward trend from 1954 to 2009.During the same period, the precipitation amounts in the Longchuan and Heyuan Stations have displayed a slightly decreasing tendency (Table 4).In contrast, the precipitation in the Boluo Station displayed a slightly increasing tendency, indicating that in addition to precipitation, other factors also affect flood changes.To better understand the possible factors for the flow change, we also collected the Normalized Difference Vegetation Index (NDVI).NDVI reflects the overall information for the ground vegetation.The data are from Global Inventory Modeling and Mapping Studies (GIMMS) [41].It can be seen from Table 4 that the vegetation cover keeps decreasing due to the rapid urbanization that occurred in this area [26].The percentage of urban area was 1.94 in 1990 and this increased to 4.79 in 2010, indicating rapid urban expansion during the interim years [31].The rapid urbanization in the River basin can significantly increase the flood peak, causing more quick flow and therefore strengthening flood events, etc. (e.g., Zhang et al. 2015) [17].This is contrary to the trend of the flood sequence.Therefore, the non-stationarity of hydrological series is impacted by continuous changing environment of multi-factor synthesis.The number of reservoirs and the water storage capacity are increasing gradually.The area of urban construction land is continuously expanding, and the change of climatic factors, especially the trend of rainfall, is also a slow process.Difference Vegetation Index (NDVI).NDVI reflects the overall information for the ground vegetation.The data are from Global Inventory Modeling and Mapping Studies (GIMMS) [41].It can be seen from Table 4 that the vegetation cover keeps decreasing due to the rapid urbanization that occurred in this area [26].The percentage of urban area was 1.94 in 1990 and this increased to 4.79 in 2010, indicating rapid urban expansion during the interim years [31].The rapid urbanization in the River basin can significantly increase the flood peak, causing more quick flow and therefore strengthening flood events, etc. (e.g., Zhang et al. 2015) [17].This is contrary to the trend of the flood sequence.Therefore, the non-stationarity of hydrological series is impacted by continuous changing environment of multi-factor synthesis.The number of reservoirs and the water storage capacity are increasing gradually.The area of urban construction land is continuously expanding, and the change of climatic factors, especially the trend of rainfall, is also a slow process.

Selection of Optimal TVM Model
AIC was taken as the discrimination criterion of the optimal model (Table 5).Table 6 shows the parameter estimations for all distributions by MLE.It is seen from Table 5 that the LN2CP model at Longchuan, GMBCP model at Heyuan, and GMBCL model at Boluo Station were the best fitting models in terms of the AIC values.The performance of the optimal model selection using AIC can be analyzed from two aspects, i.e., from the characteristic of distribution curves and from the trend model, as detailed in the following sub-sections.

Selection of Optimal TVM Model
AIC was taken as the discrimination criterion of the optimal model (Table 5).Table 6 shows the parameter estimations for all distributions by MLE.It is seen from Table 5 that the LN2CP model at Longchuan, GMBCP model at Heyuan, and GMBCL model at Boluo Station were the best fitting models in terms of the AIC values.The performance of the optimal model selection using AIC can be analyzed from two aspects, i.e., from the characteristic of distribution curves and from the trend model, as detailed in the following sub-sections.

Selection of Optimal Distribution
The flood frequency curves of different distributions mainly show differences at the tails of the distribution, especially in the high flow part which generally shows big differences for different distributions.The values of coefficient of variation Cv and coefficient of skewness Cs for annual maximum flood peak flow data, gradually decreased from upstream to downstream of the catchment (Longchuan Cv = 0.775, Cs = 2.261; Heyuan Cv = 0.627, Cs = 1.211;Boluo Cv = 0.445, Cs = 1.008).In general, Cv and Cs are the important parameters for most distributions and different values of which will have significant influence on the design values.When Cv and Cs values are larger, the upper tail of frequency curve of flood series will change from "gentle" to "steep", and the heavy tailed distribution is more flexible on this part and is more suitable for fitting the data.As can be seen from the study of Adlouni et al. (2008) [42], the tail of lognormal distribution is thicker than Gumbel distribution.This conclusion coincides with the result of our model selection, which shows that LN2 distribution is suitable for Longchuan Station in upstream, and Gumbel distribution is suitable for middle and lower reach Stations of Heyuan and Boluo (Figure 3) (Table 5).
It would be important to understand the possible reasons for selecting different distributions.As rainfall is the primary source of streamflow in the study basin, the possible relationship between annual maximum streamflow and corresponding basin rainfall in the past 7 days was investigated.For Longchuan and Heyuan Stations, the flood series and the corresponding rainfall series are consistent with decreasing trend for both series, and for flood series the trend is statistically significant at the 0.05 level (Table 4), while the decreasing trend for the rainfall series is not significant at the 0.05 level (Figure 4).For Boluo Station, the flood series and the corresponding rainfall series are inconsistent with an increasing trend for rainfall series but a slightly decreasing trend for the flood series (Figure 4), meaning that the effect of rainfall on the flood is weakening.As a possible consequence, the Cs values of the peak flood series exhibit a decreasing trend from upstream to downstream, which may be caused by the variation of rainfall, urbanization activities and basin reservoir regulation storage.Due to the enhancing degree of human influence from upstream to downstream, the sensitivity of the stream flow to climate fluctuation gradually reduces from upstream to downstream in East River basin.

Selection of Optimal Trend Model
According to the trend model parameters of the optimal model (LN2CP for Longchuan, GMBCP for Heyuan, GMBCL for Boluo Station) (Figure 3) (Table 6), the change processes in means and standard deviations of flood series as showed in Figure 2 are discussed in the following paragraphs.

Selection of Optimal Trend Model
According to the trend model parameters of the optimal model (LN2CP for Longchuan, GMBCP for Heyuan, GMBCL for Boluo Station) (Figure 3) (Table 6), the change processes in means and standard deviations of flood series as showed in Figure 2 are discussed in the following paragraphs.
(1) The CP trend model for Longchuan station in the upper reaches.Referring to the seven different types of trend forms assumed in Section 3.2, in either A, B or C-type trend models, model parameters are ultimately obtained by fitting the observed series.Therefore, the potential trends of the mean and standard deviation are firstly determined by observed series intuitively.The means of observed flood flow series generally show a parabolic change trend as "decline-steady-rise" (Figure 2).Under the control of the CP model, the mean curve conforms to the change process of observed flood flow (Figure 2).It is necessary to consider the change trend of standard deviation for the hydrological extremes analysis.Standard deviation reflects the concentration degree of random variables.The smaller the standard deviation is, the more concentrated the random variables will be, while in contrast, the more easily the random variables will deviate from the average value, leading to an increase in the occurrence of extreme events.It was found that the flood peak values on both ends of the series were less concentrated; the entire series showed more intense flood peak changes in the first 20 years of the series with larger standard deviation.After 1976, the series has smaller standard deviation until 2004, followed by an increasing trend with standard deviation.In general, the change trends of the standard deviation and the mean are the same, i.e., a parabolic change trend of "decline-steady-rise".The standard deviation curve of CP model is consistent with the standard deviation of the observed annual maximum peak flow (Figure 2).As a result, the CP trend model (mean and standard deviation related by the parabolic trend model) had a better fitting effect than other trend models.period between non-stationary processing and traditional flood methods (traditional method is a frequency analysis method in which the hydrological sequence satisfies the stationarity assumption), relative difference degree E was used to quantitatively express the difference degree between the two.The equation is shown as below: where, Q 0 indicates the design flood flow of optimal curve with TVM model, and Q T indicates the design flood flow with traditional methods.
Comparing the traditional curve LN2 with the TVM model at Longchuan Station, the difference of 10-year flood flows was 28.36%, while the difference of 100-year flood flows increased to 38.65% (Table 7).From another perspective, for the same flood flow, the return period obtained by the non-stationary processing approach is longer than that of the traditional methods.This is mainly because the flood peak reduction effect of Fengshuba Reservoir led to the non-stationary flood peak flow in Longchuan Station.For the traditional curve GMB at Heyuan Station, the difference of design flood flows from TVM model and traditional methods corresponding to return periods of 10 years was 53.24%, and that with a flood peak corresponding to return periods of 100 years increased to 57.39% (Table 7).For the traditional curve GMB at Boluo Station, the difference of 10-year flood flows was 26.06%, and that of 100-year flood flows increased to 26.92% (Table 7).

Discussion
Correct selection of trend models for mean and standard deviation has a great impact on the estimated flow return period [10,17,18].Even for the same distribution curve, the selection of different trend models will cause significant differences in the calculated return period.For example, GMBBL with poorer fitting effects at the Boluo Station was non-consistent with the standard deviation for observed series that changes from large to small in the period of 50 years.The return period was significantly reduced over time, opposite with the results of the optimal model that considers the non-stationarity.Therefore, in the calculation of non-stationary flood frequency, it is required not only to consider the change trend of the mean value of flood peak, but also to consider the change trend of standard deviation.AIC and other auto-selection criteria all have deficiencies, because most of the data information is covered in the overall statistics [39].In the initial selection of the contrast model for time series, the intuitive judgments on the fitting degree and the length of time series must be considered.The first step should be to intuitively determine the potential trends of mean and standard deviation in the existing time series.The above phenomena reflect the fact that the parameter change of non-stationary flood series will cause a greater impact on the estimation of flood return period.The correct selection of change trend model for the parameters has a significant impact on flood design discharge, which should be considered in the engineering and hydrologic design.
The results of previous studies show that a changing environment in the river basins can significantly alter the flood volume and flood peak, and change the time of the peak.The impact of urbanization activities and the reservoirs on estimates of flood frequency is considerable [16,17,44].Influenced by climate change and the degree of human activities, the tail-end of the high-water frequency curve declined significantly.The occurrence probability of the same magnitude of flood was significantly reduced, as was the real design flood magnitude.Both precipitation and runoff changed smoothly during the period, and the slight difference in the changing trend of runoff is caused by factors other than precipitation.One of the most important signs of human activity is land use change.Benefitting from the significant growth of population and Gross Domestic Product (GDP), technological progress and policy support, the downstream areas such as Shenzhen and Dongguan have accelerated the pace of urbanization since 1990s [45].The percentage of urban area was 1.94 in 1990 and increased to 4.79 in 2010.The urban area downstream expanded most rapidly mainly by encroaching into the rice field, agricultural area, and forest land [31].In addition, many reservoirs of different sizes have been constructed continuously from the 1950s to the present.The adjustment coefficient of the large reservoirs on runoff at the Boluo Station reached 0.34 [26].The impact of reservoirs is similar to that of urbanization.They are not built in a short period of time, but are gradually built over the years of the whole study period.
The change factors that cause the "non-stationarity" of the flood series can be divided into two categories, i.e., continuous (e.g., urbanization, landuse change) factors and one-time (e.g., dam construction) factors.The TVM method is better suited for application on the condition of continuous change.For a situation of one-time change, TVM is of course not necessary, as we can split data into two stationary periods, i.e., pre-and post-construction periods, and fit different distributions separately.Climate change and landuse change can largely be regarded as continuous changes in the East River Basin.Multiple hydraulic engineering structures including reservoirs and locks which have been constructed year by year from the 1950s to the present in the study area also cause a continuous change tendency of flow series with some notable fluctuation.Therefore, the changed environment in the East River Basin can largely be regarded as the result of continuous changes with some notable influence of reservoir operations.If we would discuss only the impact of a single reservoir in a basin, it could not be considered as a continuous time series but instead as a one-time change time series, and it therefore would not be proper to use the TVM method.In a future study, it would be interesting to compare the influence of one-time factors vs. continuous factors on the goodness of fitting (e.g., dam volume vs. flow rate, urbanization vs. flow rate) in the East River basin or other basins influenced by multiple factors.

Conclusions
The flood series in the East River basin show non-stationarity in time impacted by continuous changing environment of multi-factor synthesis.In this study, the non-stationary processing method of time-varying moments was adopted to analyze and calculate the annual maximum daily flow series at the Longchuan, Heyuan and Boluo Stations of the East River basin.Five types of distribution curves and eight kinds of trend models, for a total of 40 models (including the case of no consideration on trend), were selected for comparison.We found that the non-stationarity of series has a significant influence on the results of frequency analysis, and this influence is station dependent, as the human activities have shown different impacts on different stations.The study has resulted in the following specific conclusions.
1. LN2 distribution with a CP model at Longchuan Station, GMB distribution with a CP model at Heyuan Station, and GMB distribution with a CL trend model at Boluo Station obtain the optimal fitting.
2. The flood series and optimal distribution curves selection in the East River basin have been significantly affected by the continuous changing environment.With the precipitation change, the continuous process of urbanization and the reservoirs storage-capacity increases, the thinner tails of

Figure 1 .
Figure 1.Location of the study area and the hydrological stations.

Figure 1 .
Figure 1.Location of the study area and the hydrological stations.

1 .
Trend in mean values-(A); 2. Trend in standard deviations-(B); 3.Both the mean value and standard deviation have a trend and the ratio of the two (i.e., Cv) is constant-(C); 4. Both the mean value and standard deviation have a trend and the ratio of the two (i.e., Cv) is not constant-(D).
is the mean of the time series, σ is the standard deviation of the time series, t is time, here we express m and σ as ( is the mean of the time series, σ is the standard deviation of the tim we express m and σ as ( et al. (2017a) and Strupczewski et al. (2001) [12,32].

Figure 2 .
Figure 2. Change process of mean value and standard deviation calculated by TVM model for annual maximum flood peak flow series in East River basin.

Figure 2 .
Figure 2. Change process of mean value and standard deviation calculated by TVM model for annual maximum flood peak flow series in East River basin.
Atmosphere 2018, 9, x FOR PEER REVIEW 9 of 16 upstream to downstream, which may be caused by the variation of rainfall, urbanization activities and basin reservoir regulation storage.Due to the enhancing degree of human influence from upstream to downstream, the sensitivity of the stream flow to climate fluctuation gradually reduces from upstream to downstream in East River basin.

Figure 3 .Figure 3 .
Figure 3.Comparison for 5 typical distributions fitted to stationary annual maximum flood peak flowseries by TVM method.

Figure 4 .
Figure 4. Change process of the sum of seven-day rainfall values before the annual maximum daily flow.

Figure 4 .
Figure 4. Change process of the sum of seven-day rainfall values before the annual maximum daily flow.

Table 1 .
Detailed information on the main reservoirs in East River basin.

Table 2 .
The probability density functions of various distributions and relationships of distribution parameters and first two order moments.

Table 3 .
The expressions of first two order moments for various types trend model.

Table 3 .
The expressions of first two order moments for various types trend model.

Table 3 .
The expressions of first two order moments for various types trend model.

Table 3 .
The expressions of first two order moments for various types

Table 4 .
Mann-Kendall trend test values of hydrological parameters and NDVI in East River basin.

Table 4 .
Mann-Kendall trend test values of hydrological parameters and NDVI in East River basin.

Table 5 .
AIC goodness-of-fit values of various distributions for flood frequency analysis in East River basin.

Table 5 .
AIC goodness-of-fit values of various distributions for flood frequency analysis in East River basin.

Table 7 .
The difference of return periods between Non-stationary processing and traditional flood frequency analysis methods.