The Probability Density Evolution Method for Flood Frequency Analysis: a Case Study of the Nen River in China

A new approach for flood frequency analysis based on the probability density evolution method (PDEM) is proposed. It can avoid the problem of linear limitation for flood frequency analysis in a parametric method and avoid the complex process for choosing the kernel function and window width in the nonparametric method. Based on the annual maximum peak discharge (AMPD) in 54 years from the Dalai hydrologic station which is located on the downstream of Nen River in Heilongjiang Province of China, a joint probability density function (PDF) model about AMPD is built by the PDEM. Then, the numerical simulation results of the joint PDF model are given by adopting the one-sided difference scheme which has the property of direction self-adaptive. After that, according to the relationship between the marginal function and joint PDF, the PDF of AMPD can be obtained. Finally, the PDF is integrated and the frequency curve could be achieved. The results indicate that the flood frequency curve obtained by the PDEM has a better agreement with the empirical frequency than that of the parametric method widely used at present. The method based on PDEM is an effective way for hydrologic frequency analysis. 5135 Keywords: probability density evolution method (PDEM); annual maximum peak discharge (AMPD); flood frequency analysis


Introduction
The target of flood frequency analysis is to derive the probability distribution function (PDF) of annual maximum peak discharge (AMPD) according to the measured data.Then, the magnitude of a flood for a given return period can be estimated as the design basis of flood sluicing buildings, bridges, culverts, and other flood control projects [1,2].Flood frequency analysis provides not only design basis for water conservancy projects which are going to be built, but also a calculation foundation for risk analysis and systematic use of water resources for those that have been built.Hence, flood frequency analysis is of great importance in hydraulic engineering.There have been numerous studies about it and the methods for flood frequency analysis can be divided into parametric method and nonparametric method, generally.
The parametric method for flood frequency analysis mainly includes two steps: (1) select an appropriate parent distribution for the measured data and, (2) estimate the parameters for the selected distribution.In regard to parent distribution selection, the Monte Carlo Bayesian method for a comprehensive study on computing the expected probability distribution was introduced by Kuczera [3].A variety of mathematic test methods in flood parent distribution selection were adopted by Kite et al. [4][5][6][7][8].Besides, the probability plot correlation coefficient (PPCC) test which was a powerful and easy-to-use goodness of fit test in completing samples for the composite hypothesis of normality was proposed by Filliben [9].Thereafter, the PPCC test was extended for studying kinds of probability distribution types [10][11][12].With respect to the parameter estimation research for the selected parent distribution, several methods have been extensively used, such as method of moments (MOM), method of maximum likelihood (ML), and curve-fitting method.The MOM, which includes the conventional moments [13], the linear moments (L-moments) [14][15][16][17][18], and the probability-weighted moments (PWM) [19][20][21][22][23][24][25][26][27][28], is the oldest way of deriving point estimators.The method of ML can give the most likely values of the parameters for a given distribution.It was adopted in estimating parameters for various distribution types [29,30] and different situations [31].Moreover, El Adlouni, S. et al. [32] developed the generalized ML estimators for the nonstationary generalized extreme value model by incorporating the covariates into parameters.As to the curve-fitting method, it can be divided into two main categories, one is fitting the empirical frequency curve to the best by visual observation, the other is optimal curve-fitting method.Due to the characteristics of flexibility and easy manipulation, the curve-fitting method by visual observation is extensively applied in practical flood frequency analysis in China.However, the weakness of this method is prominent that different results may be obtained by different persons at different time.Obviously, it is not an objective method [33].Thus, many researchers paid attention to the optimal curve-fitting method [34][35][36].From the above studies, it can be seen that much work has been conducted on both steps of the parametric method.However, there still exists one problem in the parametric method because it is based on an assumption that the data in a given system obeys a certain distribution.It is regrettable that the limited distribution types cannot always satisfy all kinds of practice.Hence, there is the so-called linear limitation problem in parametric method.When the assumed parent distribution does not agree with the fact, the precision of corresponding calculation results cannot be guaranteed.
For the linear limitation problem in the parametric method, Tung [37] proposed the application of the nonparametric method in hydraulic frequency analysis in 1981 for the first time.The nonparametric method can solve the linear limitation problem for the reason that the process of parent distribution assumption is avoided.Therefore, it was employed and emphasized in this field gradually [38][39][40][41][42][43][44].Nevertheless, the defect of the nonparametric method is also evident that there is no specific criteria about how to select a better kernel function and window width in different situations.Because of this, it is difficult to take use of the nonparametric method in practice [45].Thereby, a method which does not need to choose kernel function and window width, meanwhile the assumption of parent distribution can be avoided, is necessarily demanded.
The probability density evolution method (PDEM) which has no need of assuming a parent distribution in flood frequency analysis takes the probability conservation principle, which means that the law of probability conservation exists during the state evolution process for a conservative stochastic system, as the theoretical foundation.The basic idea of using PDEM in flood frequency analysis is different from the nonparametric method.Therefore, the process of choosing kernel function as well as window width is avoided, and the PDEM is adopted for flood frequency analysis in this paper.Firstly, the joint PDF model about AMPD is established through PDEM.And then, by means of numerical method, the model is able to be solved and the probability density value as well as the probability distribution value of AMPD can be obtained.Finally, the frequency curve will be achieved.
This paper is organized as follows: in Section 2, the basic idea of the PDEM is introduced.In Section 3, a model about the joint PDF of AMPD is derived based on the PDEM.The solution method about the model and frequency values calculation of AMPD are suggested.Also presented in this section is the robustness study about the proposed method by statistical experiment research (Monte-Carlo method).An example and discussion will be given in Section 4, while the last section contains a brief summary of the results.

Probability Density Evolution Method
Due to the coupling effect between nonlinear and stochastic, the accurate prediction of nonlinear response for practical engineering structure is difficult to implement.In the light of this, Li and Chen [46,47] developed an in-depth study on the principle of probability conservation and derived the PDEM by the combination of state space description for principle of probability conservation and different physical equations, such as the classical Liouville Equation, Fokker Planck Kolmogorov (FPK) Equation, and Dostupov-Pugachev (D-P) Equation.The PDF of structure response containing the whole stochastic factors in the dynamic system can be acquired through the numerical solution of probability density evolution equation.Thus, in recent years, a systematic research of the PDEM has been made in stochastic response analysis of multi-dimensional linear and nonlinear structure systems, calculation of dynamic reliability and system reliability, as well as control aspects based on reliability.The basic idea of PDEM is briefly introduced as follows [48,49].
Suppose that an m-dimension virtual stochastic process can be described as: ( ,τ), ( , τ) where, X is the state vector, Φ represents the mapping relationship between X and Θ, τ.Θ is random vector, τ is virtual time.X l and Φl are the l-th component of X and Θ, respectively, l = 1, 2,…, m, m is an integer.According to the probability compatibility condition, when {Θ = θ}, the expression about X l (τ) can be obtained: where ( ,τ θ) P x is the conditional PDF of X l (τ) at {Θ = θ}.
According to the conditional probability formula, the joint PDF of stochastic variable (X l (τ), Θ) which is denoted as ( ,θ,τ) l P x Θ X can be derived from Equation (4): X X (5) where P Θ (θ) is the PDF of Θ.
Differentiating Equation (5) on both sides with regard to τ, and using the derivation rule of compound function will yield: Because both θ and τ are fixed values in the differential process of compound function, the partial derivative of y = x − Hl(θ,τ) is dy = dx.Therefore, substitute ∂y in Equation ( 6) with ∂x.Then, Equation ( 6) can be rewritten as follow: Equation ( 8) is the probability density evolution equation, and the corresponding initial condition can be conveniently obtained from Equation ( 5): where x0,l is the l-th component of x0.
The boundary condition of Equation ( 8) reads as: ( ,θ,τ) 0 By means of numerical difference method, the probability density evolution equation can be solved based on Equations ( 9) and (10).Thus, the joint PDF ( ,θ,τ) P x Θ l X is obtained.In accordance with the relationship between the marginal function and joint PDF, the PDF ( , τ) P x l X can be acquired: where Ωθ is the distribution domain of Θ.

The Joint PDF Model for Peak Discharge
In flood frequency analysis, the sample of AMPD, which is composed of measured data from the hydrologic stations for years, can be denoted as Z = (z1, z2, z3,…, zn), here n is the total number of data in the sample.With consideration of the stochastic characteristic of hydrologic variables, the AMPD can be expressed by a one-dimension static stochastic process Y (y).
Construct a one-dimension stochastic process: where, X(t) is the state vector of parent peak discharge.H represents the mapping relationship between X and Z, t.Z is the random vector consisting of peak discharge measurement of data.t is time.
It can be seen from Equation ( 13) that the PDF of AMPD can be obtained through the derivation of PDF for one-dimension stochastic process X(t) at t = 1.
Comparing Equation (12) with Equation (1), it is found that both forms of them are consistent.Therefore, a joint probability density evolution equation of (X, Z) corresponding to Equation (8) can be obtained as follows: (14) where P XZ (x,z,t) is the joint PDF of (X, Z).
Discretize z and bring the measured data in Equation (14).Then, the joint PDF model of peak discharge can be acquired: XZ XZ (15) Considering that the acquisition of the measured data for AMPD is an independent process; therefore; when there is no any other additional information; the probability of the measured data for AMPD which is denoted as P Z (zj) can be described as: Taking into account Equations ( 9) and ( 16), the initial condition of Equation ( 15) can be written as: where x0 is the initial value of AMPD.According to reference [50], the system state H(z,t) can be expressed as: In this way, the discrete form of velocity function in Equation ( 15) is taken as follows: ( , ) cos(2.5π ) 2.5π

Model Solution and Frequency Values Calculation of Peak Discharge
Equation ( 15) is the one-dimension variable coefficients convection equation, which can be solved by many finite difference methods.According to calculation experience, a fairly satisfactory calculation precision will be achieved by adopting the one-sided difference scheme.Considering that the sign of the velocity function may change between positive and negative with the variation of time t, in this paper, the one-sided difference scheme possessing the characteristic of direction adaptive is selected to discretize Equation (15) as follows:

P h r p h r h r p h r h r p
where Pj,m denotes the discrete form of joint PDF for XZ and is short for P XZ (xi,zj,tm), in which xi = iΔx (i = 0, ±1, ±2,…), tm = mΔt(m = 0,1, 2,…).xi and tm are the discrete values of calculation interval for time axis and space axis, respectively.Δx and Δt are the space step and time step, respectively, rL is the ratio of Δt to Δx.
In order to ensure the convergence of the calculation results, the Courant-Friedrichs-Lewy condition for the difference scheme of Equation ( 20) is needed: After establishing the joint PDF model of peak discharge and selecting a suitable difference scheme, the calculation of flood frequency based on PDEM can be described as follows: (1) Determine space step Δx and time step Δt for evaluating the velocity function and initial condition.First of all, select a suitable Δx and Δt according to Equation ( 22) by trial method.Secondly, denote the discrete points of time axis in the domain [0, 1] as tm.Thus, the velocity function can be calculated through Equation (19).Likewise, denote the discrete points of space axis in the corresponding domain as xi.After that, the discrete form of Equation ( 17) can be given as: (2) Solve Equation ( 15) by employing the difference scheme of Equation (20).Then, the discrete value of joint PDF P XZ (xi,zj,tm) at tm = 1 can be acquired.
(3) Conduct the numerical integration with respect to zj according to Equation (11).Thereby, the discrete PDF will be written as: (4) Calculate the frequency values of peak discharge through the discrete PDF.Firstly, in order to get a better calculation precision for the flood frequency, the cubic spline interpolation on the PDF of peak discharge calculated by Equation ( 24) is implemented.Secondly, the trapezoid method is adopted to calculate the area value between the adjacent points on the curve of peak discharge PDF which has already been interpolated.Then, the values of cumulative area with descending order of peak discharge are computed.Finally, the flood frequency can be obtained.
The flowchart of the flood frequency calculation based on the PDEM is depicted in Figure 1.
Figure 1.Flowchart of the flood frequency calculation based on the PDEM.

The Robustness Study for PDEM
In order to analyze the robustness of the proposed method, the statistical experiment research (Monte-Carlo method) is adopted.Since the measured hydrological data of natural rivers is relatively short, usually less than 60 years, hence, the statistical experiment research with small samples has actual hydrological significance.In order to evaluate the robustness of PDEM, two commonly used parent distribution types for hydrology in this statistical test are taken to analyze and calculate two groups of samples and four parent parameters, with a total of 32 programs.The specific steps for the robustness study are as follows: 1. Give the parent parameters of statistical distribution and assume their distribution types, then extract m group of sample series with the length of n from every parent distribution.2. Calculate the design value of design frequency with PDEM.
3. Analyze various types of error of design value and study their robustness.
The specific plan is as follows: (1) The parent distribution: Pearson Type III, Log Normal.
The relative RMS error is: The relative error of the mean design value is: where, xp represents a true value of design frequency p whose distribution is known, pi x  is the design value of estimated p.The calculation results of 32 plans about statistical tests mentioned-above can be seen in Tables 2  and 3.It can be found for both distribution of Log Normal and Pearson Type III, the relative error of the mean design value ω does not exceed 15%, and the relative RMS error δ of that is not greater than 40%.In addition, the corresponding design value of same design frequency calculated by different theoretical parent distribution is relatively close.It is demonstrated that the difference of theoretical population has little influence on the calculation results.Therefore, as an estimated method of design value, the PDEM is relatively robust.

Study Area
In this study, flood frequency analysis is performed on the Nen River in Heilongjiang Province in Northeast China.The Nen River is the biggest tributary of the Songhua River, whose annual runoff and catchment are ranked third among all of the rivers in China.The area of the Nen River basin is approximately 297,000 km 2 , and the full-length of its main stream is more than 1370 km.In addition, the climate characteristic of this basin is that winter is long and cold, while summer is short and rainy.About 82% annual precipitation comes from June to September.The Dalai hydrologic station is a central control station located at the downstream of the Nen River.According to the data of AMPD in Table 4, which covers 54 years measured data from Dalai hydrologic station, the PDF and the frequency of AMPD based on the PDEM are studied.

Data Analysis Method and Parameters Selection
To show the rationality of the calculation results based on PDEM, two groups of comparisons are conducted: (1) The comparison between frequency histogram and the PDF of peak discharge calculated by PDEM as well as parametric method.
(2) The comparison between empirical frequency and the frequency values of peak discharge calculated by PDEM as well as parametric method.
All the calculations in above comparisons are carried out based on the Matlab software.The interval number of frequency histogram in this study is evaluated by the following equation as mentioned in [51]: where k is the interval number.
As to the empirical frequency, it is derived from the expectation equation wholly adopted in hydrological frequency analysis at present.Arrange the measured data in a descending order, as z1, z2,…, zm,…, zn.Then, the empirical frequency of peak discharge sample Fe(zi) can be expressed as: ( ) The related parameters of the flood frequency analysis method used in this paper are given as follows: the time step Δt and the space step Δx are selected according to Equation ( 22) by trial method.When Δt equals to 0.0001 s and Δx equals to 19.95 m 3 /s, Equation ( 22) can be satisfied.Taking into account that the value of AMPD is greater than zero, the initial sample value x0 can be selected between zero and the minimum value of measured data.From Table 4, it can be seen that the minimum value and the maximum value of AMPD are 541 m 3 /s and 16,100 m 3 /s, respectively.Hence, set x0 = 50 m 3 /s, and the calculation interval of AMPD is selected ranging from 50 to 20,000 m 3 /s.

Results and Discussion
During the process of parent distribution selection in parametric method, three fitting methods including the chi-square test, the Kolmogorov-Smirnov test, and the Anderson-Darling test are employed to analyze six kinds of distribution types which are often used in hydrologic frequency analysis.The parameters of them are estimated by the method of ML, and the relevant results are listed in Table 5.The results of goodness of fit tests are summarized in Table 6.The calculation methods of statistic values for fit tests in Table 6 are introduced in reference [52].It can be seen from Table 6 that the order of calculation results obtained by chi-square test are somewhat different with those obtained by the other two tests.However, the top four orders of distribution types, namely Log Pearson Type III distribution, Log Normal distribution, Gamma distribution, and Pearson Type III distribution are identical.Hence, the above four distribution types are selected for subsequent investigations.According to the PDF expression of the selected distributions, the calculation results of PDF by the parametric method can be achieved as seen in Figure 2. In the same way, through the expression of PDF, the frequency values of peak discharge can also be obtained, as depicted in Figure 3.It can be seen from Figure 2 that the trends of probability density curve for AMPD calculated by parametric method with four distributions are similar, but the right tail does not agree well with the frequency histogram.For the matching result calculated by PDEM, it has a better agreement with the frequency histogram, even the mutational parts on the right tail of the probability density curve are included.Hence, compared with the parametric method, taking use of PDEM can show the preferable effect on the distribution structure of PDF for AMPD.The reason for this is that PDEM is a kind of method which can be directly estimated through the sample data and does not need to assume a parent distribution.Moreover, it can also simulate different types of measured data, including data with multi-peak distribution, which is impossible to be realized by parametric method.
Figure 3 is the frequency curve of peak discharge calculated by PDEM and the parametric method.It is shown that both methods can fit well with the empirical frequency points at the middle of the frequency interval, however, at the two ends of the frequency interval, more points will be fitted by means of the PDEM instead of the parametric method.Therefore, according to the above analysis from Figure 3, it can be concluded that employing PDEM to calculate the frequency curve of peak discharge has a better fitting performance than that with parametric method.The reason for this is analogous to the description mentioned in Figure 2.
In order to get a more clear recognition about the fitting effect of empirical frequency points and design frequency values obtained by PDEM as well as parametric method, the Equation ( 29) is selected to analyze the RMS error, which is treated as the measurement index of estimation accuracy for flood frequency.
In Equation (29), RMSE is the RMS error, Fd(zi) is the design frequency of peak discharge sample.The calculation results of Equation ( 29) are shown in Table 7. From Table 7, it can be observed that the orders of RMS error for PDEM and parametric method including four different parent distribution types are arranged as PDEM < Gamma < Log Normal < Log Pearson Type III < Pearson Type III.The results demonstrate that PDEM has a better calculation precision than that of parametric method.
The stability criteria for calculation results can be denoted by the relative error.The relative error between the 54 empirical frequency points and the corresponding design frequency obtained by the two methods is calculated and the values of relative error in certain intervals are shown in Table 8.The smaller the interval is, the higher the stability request represents.Through calculation it is found that both of the maximum relative errors for AMPD calculated by two methods appear at the maximum sample value.In reality, the larger the peak discharge is, the greater influence the projects suffer.For that reason, more attention should be paid to the maximum relative errors.Besides, the sample numbers corresponding to the relative error in certain ranges reflect the stability of the calculation results.Hence, attention should also be drawn to them.From Table 8, it can be seen that the maximum relative error of calculation results by Pearson Type III distribution is the biggest, while the sample numbers corresponding to the relative error falling in certain ranges are the least among the four distribution types in the parametric method.Therefore, the stability of the calculation results by Pearson Type III distribution is not satisfied.For Gamma and Log Normal distribution, although the general sample numbers corresponding to the relative error falling in certain ranges are more than that of Log Type III distribution, the maximum relative errors of them are also expanded.With regard to the Log Type III distribution which is contrary to Gamma and Log Normal distribution, it possesses a smaller maximum relative error, but the same fewer sample numbers.In a word, the calculation results of relative error by the parametric method are not desirable.While, it is observed that more sample numbers and a smaller maximum relative error can be obtained by PDEM than by Gamma & Log Normal distribution, and Log Type III distribution, respectively.Further analysis shows that even though there are few differences between the calculation results by PDEM and Gamma distribution as well as Log Normal distribution when the relative error is less than 3%, 5%, and 10%, the former presents advantages if the relative error is less than 1% compared with the two latter ones.The number (ratio) of the relative error less than 1% obtained by PDEM is 13 (24.074%), the corresponding results calculated by Gamma distribution as well as Log Normal distribution are 7 (12.963%)and 8 (14.815%), respectively.The differences of the calculation results between PDEM and the two latter ones are 6 (11.111%) and 5 (9.259%), respectively.That is to say the calculation results by PDEM are advantageous, especially when the stability request is high.In addition, with regard to maximum relative error, even comparing with the results of Log Pearson III distribution (61.500%), which has the smallest maximum relative error in parametric method, the performance of PDEM (47.602%) is not in the least inferior, only accounting for 77.402% of it.
From the analysis above, it can be concluded that the flood frequency in Dalai station on the Nen River calculated by PDEM not only has a remarkable advantage on calculation accuracy, but also a better stability than the results calculated by parametric method.Therefore, it is suitable to use the PDEM in flood frequency analysis due to its estimation characteristic using the sample data directly instead of assuming a parent distribution at the beginning.

Conclusions
To avoid the disadvantage of the parametric method and nonparametric method in flood frequency analysis, the calculation of frequency curve for peak discharge based on the PDEM is presented in this paper.Take the example of flood frequency analysis on the Nen River, the numerical analysis results indicate that the frequency curve of peak discharge calculated by PDEM fit better with the empirical frequency points than that by parametric method.The detail conclusions are as follows: (1) The PDEM can avoid the problem of linear limitation in flood frequency analysis using parametric method.It has a more rational theoretical background due to the fact that there is no need to give a parent distribution when calculating the flood frequency.Furthermore, the PDEM in flood frequency analysis can also avoid the complex problem of choosing the kernel function and window width in nonparametric method.The joint PDF model about AMPD built based on the PDEM is solved through numerical method with sample data directly.Thus, compared with the nonparametric method, the adoption of PDEM method in frequency analysis is easier to operate.
(2) It is obvious that the frequency curve of peak discharge calculated by PDEM fit better with the empirical frequency points than by parametric method with four distribution types.When the RMS error is treated as the accuracy index of calculation results, the corresponding value for PDEM is no more than 90.476% of that for parametric method.While the stability of calculation results is judged by relative error less than 1%, the former stability value is more than 1.625 times of the latter ones.In addition, when the value of AMPD is huge, in another word, which has great influence in the engineering design, the relative error calculated by PDEM is only 77.402% of the best calculation results in parametric method.Overall, the results of estimation accuracy and stability calculated with PDEM are better than those obtained with the parametric method.Therefore, it can be concluded that adopting PDEM to calculate the frequency curve of peak discharge is feasible and effective.And it is suggested to consider PDEM in more frequency analysis in hydraulic engineering.

Table 5 .
Estimated parameter values of different distribution types.

Figure 2 .
Figure 2. Probability density curve of AMPD from Dalai hydrologic station.

Figure 3 .
Figure 3. Frequency curve of AMPD from Dalai hydrologic station.

Table 2 .
The results calculated by PDEM with theoretical population of Log Normal distribution.

Table 3 .
The results calculated by PDEM with theoretical population of Pearson Type III distribution.

Table 4 .
The measured data of annual maximum peak discharge (AMPD) from Dalai hydrologic station.

Table 6 .
The results of goodness of fit tests.

Table 7 .
Calculation results of RMS error.

Table 8 .
Calculation results of relative error.