Application of Multi-Step Parameter Estimation Method Based on Optimization Algorithm in Sacramento Model

: The Sacramento model is widely utilized in hydrological forecast, of which the accuracy and performance are primarily determined by the model parameters, indicating the key role of parameter estimation. This paper presents a multi-step parameter estimation method, which divides the parameter estimation of Sacramento model into three steps and realizes optimization step by step. We ﬁrstly use the immune clonal selection algorithm (ICSA) to solve the non-liner objective function of parameter estimation, and compare the parameter calibration result of ideal artiﬁcial data with Shufﬂed Complex Evolution (SCE-UA), Parallel Genetic Algorithm (PGA), and Serial Master-slaver Swarms Shufﬂing Evolution Algorithm Based on Particle Swarms Optimization (SMSE-PSO). The comparison result shows that ICSA has the best convergence, efﬁciency and precision. Then we apply ICSA to the parameter estimation of single-step and multi-step Sacramento model and simulate 32 ﬂoods based on application examples of Dongyang and Tantou river basins for validation. It is clearly shown that the results of multi-step method based on ICSA show higher accuracy and 100% qualiﬁed rate, indicating its higher precision and reliability, which has great potential to improve Sacramento model and hydrological forecast.


Introduction
The Sacramento model is a conceptual hydrologic model, which was proposed by the Sacramento Forecast Center of National Weather Service during the period from the late 1960s to the early 1970s. It has been widely applied to the semi-humid region up to now. What is similar to the Stanford model is that the Sacramento model is a consecutive model, and it includes 17 parameters. As is known to us all, it is essential that the parameter should be correctly determined, or the model will lose its efficiency. Therefore, the parameter calibration of conceptual watershed hydrologic model has always been an important content of the flood forecasting, which belongs to the nonlinear optimization problem.
However, most of these optimization algorithms are still limited with convergence problems, and need further improvement. The Rosenbrock algorithm and SCE-UA require strict model structure and optimization guidelines, which could be significantly influenced by initial conditions and trapped in local optimum. The GA with strong robustness and versatility can generally solves the analytic expressions of objective functions and constrains, but it often runs into problems of premature convergence or slow convergence. The PSO and improved PSO algorithms also cannot avoid the weakness of slow convergence and may be easily trapped in local extremum.
As the famous conceptual watershed hydrologic model of the United States, there usually exists a large difference between the simulated result and actual measured values in the Sacramento model because of the numerous parameters and their defective calibration. Despite the fact that varied objective functions have been tried to improve the deviation of the simulated results through single-step optimization, the effect is rather limited [13]. Thus, the multi-step parameter calibration method comes into being.
The multi-step automatic parameter calibration method was developed by Terri S. Hogue et al. for the National Weather Service [14]. In the year of 2006, their team analyzed the reliability of the multi-step automatic parameter calibration method, and researched the model's availability in a variety of climatic and hydrologic conditions [15]. Meanwhile, they cooperated with the Hydrologic Laboratory of National Weather Service and then developed an advanced parameter calibration method. The "in house" National Weather Service algorithm was used to optimize the hydrologic model, and it was applied to the Sacramento Soil Moisture model and SNOW-17 model for streamflow forecasting, which had achieved desirable results. The method can be well adapted to hydrologic model, though it was only used within the forecasting system of National Weather Service before. The achievements of their studies have showed that the multi-step automatic parameter calibration method presents excellent performance in each basin, even though in the complex terrain of the western United States.
As we can see, the traditional optimization algorithm which applied to parameter calibration of the Sacramento model still has a lot of problems such as early maturing and slowly convergent speed. Besides, the traditional methods may be easy to fall into the local optimum, and its simulative fitness is poor. In view of the above limits, this paper firstly applies the immune clonal selection algorithm (ICSA) based on artificial immune system in the parameter estimation of Sacramento model, which is more accurate and has fast convergence. We pick the artificial hydrological data from the particular water basins to verify the accuracy and performance of ICSA, and compare the results with other three algorithms, i.e., SCE-UA, PGA, and Serial Master-slaver Swarms Shuffling Evolution Algorithm Based on Particle Swarms Optimization (SMSE-PSO). Then, based on the ICSA, we improved the multi-step parameter calibration method in order to increase the accuracy of the Sacramento model simulation, and give out the detailed evaluation of the simulation results which are obtained from the particular water basins.

Immune Clonal Selection Algorithm
Burnet and Talmage [16,17] proposed the famous antibody clonal selection theory in 1957, of which the basic principle is described as: the antibody is a natural product and exists on the cell surface in the form of receptor; the antigen selectively reacts with antibody. The reaction of antigens and corresponding antibody receptors leads to the clonal expansion of cells, of which some differentiate into antibody forming cells and others differentiate into immune memory cells to participate the following secondary immune response. In this process, the cells activate, differentiate and proliferate through clone to increase the number of antibodies, and finally clear the antigens by immune response. Therefore, the clonal selection is dynamic process of self-adaption to antigen stimulation in biology immune system, through which the generated biological features of study, learn and antibody diversity can be prominent references to the artificial immune system. At present, the most classic simulation method to describe the mechanism of antibody clonal selection is the immune clonal selection algorithm (ICSA) raised by De Castro in 2000 [18], and the follow-up researchers also have developed several advanced clonal selection algorithms from diverse angles [19]. The ICSA simulates the maturity process of antibody population through clone, hypermutation and selection. The principle diagram of ICSA is shown in Figure 1.
Water 2017, 9,495 3 of 21 [19]. The ICSA simulates the maturity process of antibody population through clone, hypermutation and selection. The principle diagram of ICSA is shown in Figure 1. The ICSA has splendid performance in searching and finding solutions irrelevant to problem. The essence of clone in ICSA is to produce a mutant solution group near the candidate solutions of every generation according to affinity, which enlarges the search range of solutions (i.e., it enlarges the diversity of antibodies). The simultaneous global search and local search effectively avoid being trapped in local minimum and premature convergence of evolutionary programming, while the clone selection accelerates convergence. The process of clonal selection can be recognized as: the problem in low-dimensional space is transferred into higher-dimensional space for solution, then the solutions are projected back to low-dimensional space. The detail of ICSA method can be seen in the reference [20].
The solution procedure of ICSA is showed as below.
Step 1: Population initialization: suppose the population size of antibody is N, the clonal size of antibody is c N , the mutation probability is m p , and the recombination probability is r p . Then randomly generate the initial population contains the all parameters of the Sacramento model, and set the mature condition of affinity (that is the objective function), the currently iterations is 0 k  .
Step 2: Calculate the affinity ( ( )) f A k of antibody population ( ) A k , if the mature condition of affinity is satisfied, then output the antibody with highest affinity and stop the algorithm. Otherwise, go to the step 3.
Step 3: The cloning operation of antibody population ( ) A k : Step 4: The immune gene operation of antibody population ' ( ) A k : ( ( )) C g T A k , then obtain the antibody population '' ( ) A k .
Step 5: The selection and recombination of antibody population '' ( ) A k : ( ( )) C r T A k , then obtain the antibody population '''( ) A k .
Step 7: According to the size of affinity, perform the clonal selection operation towards the antibody population '''( ) , then obtain the next generation antibody Step 8: 1 k k   , then go back to the Step 2. The ICSA has splendid performance in searching and finding solutions irrelevant to problem. The essence of clone in ICSA is to produce a mutant solution group near the candidate solutions of every generation according to affinity, which enlarges the search range of solutions (i.e., it enlarges the diversity of antibodies). The simultaneous global search and local search effectively avoid being trapped in local minimum and premature convergence of evolutionary programming, while the clone selection accelerates convergence. The process of clonal selection can be recognized as: the problem in low-dimensional space is transferred into higher-dimensional space for solution, then the solutions are projected back to low-dimensional space. The detail of ICSA method can be seen in the reference [20].
The solution procedure of ICSA is showed as below.
Step 1: Population initialization: suppose the population size of antibody is N, the clonal size of antibody is N c , the mutation probability is p m , and the recombination probability is p r . Then randomly generate the initial population A(0) = {a 1 (0), a 2 (0), · · ·, a n (0)}, in which the matrix a i (i = 1, 2, ..., n) contains the all parameters of the Sacramento model, and set the mature condition of affinity (that is the objective function), the currently iterations is k = 0.
Step 2: Calculate the affinity f (A(k)) of antibody population A(k), if the mature condition of affinity is satisfied, then output the antibody with highest affinity and stop the algorithm. Otherwise, go to the step 3.
Step 3: The cloning operation of antibody population A(k): T c C (A(k)), then obtain the antibody population A (k).
Step 4: The immune gene operation of antibody population A (k): T g C (A(k)), then obtain the antibody population A (k).
Step 5: The selection and recombination of antibody population A (k): T r C (A(k)), then obtain the antibody population A (k).
Step 6: Calculate the affinity f (A(k)) of A (k).
Step 7: According to the size of affinity, perform the clonal selection operation towards the antibody population A (k) and A(k): T s C (A(k) + A (k)), then obtain the next generation antibody population A(k + 1).
Step 8: k = k + 1, then go back to the Step 2.

Sacramento Model
The Sacramento model is a conceptual model described by a series of mathematical expressions with certain physical concepts, which is a comprehensive river rainfall-runoff model used to simulate the hydrological cycle based on the soil moisture, infiltration, drainage and evapotranspiration characteristics. Each of its variables represents an independent and understandable characteristic of the hydrological cycle, and the model parameter can be derived from the rainfall, flow data and the basin characteristic. The Sacramento model consists of three parts: evapotranspiration calculation, runoff-yield calculation and flow-concentration calculation [21]. The model structure [22][23][24] is showed in Figure 2.

Sacramento Model
The Sacramento model is a conceptual model described by a series of mathematical expressions with certain physical concepts, which is a comprehensive river rainfall-runoff model used to simulate the hydrological cycle based on the soil moisture, infiltration, drainage and evapotranspiration characteristics. Each of its variables represents an independent and understandable characteristic of the hydrological cycle, and the model parameter can be derived from the rainfall, flow data and the basin characteristic. The Sacramento model consists of three parts: evapotranspiration calculation, runoff-yield calculation and flow-concentration calculation [21]. The model structure [22][23][24] is showed in Figure 2. As we can see, the soil is divided into upper and lower zones, and the water storage of each layer can be divided into the tension water and free water [25]. The rainfall first supplements the evenly distributed tension water in the upper zone, and then supplements the free water in the upper zone. The tension water subsided for evapotranspiration, and the free water will infiltrate into the lower zone and generate the lateral interflow. Besides, the Evapotranspiration calculation mainly accounts for five kinds of evapotranspiration, including the amount of evaporation on an impervious area 1 E , tension water evapotranspiration in the upper zone 2 E , free water evapotranspiration in upper zone 3 E , evapotranspiration on the water area of the watershed 4 E , tension water evapotranspiration in the lower zone 5 E .
It is a continuum model including 17 parameters, in which PCTIM, ADIMP and SARVA are three parameters of the direct runoff, UZTWM, UZFWM and UZK are three parameters of the water capacity of upper zone, ZPERC and REXP are the parameters of percolation, PFREE, LZTWM, LZFSM, LZSK, LZFPM, LZPK, RSERV, SIDE and SSOUT are the 9 parameters of water capacity of lower zone, and it is mainly adapted to semi-arid and semi-humid regions.  As we can see, the soil is divided into upper and lower zones, and the water storage of each layer can be divided into the tension water and free water [25]. The rainfall first supplements the evenly distributed tension water in the upper zone, and then supplements the free water in the upper zone. The tension water subsided for evapotranspiration, and the free water will infiltrate into the lower zone and generate the lateral interflow. Besides, the Evapotranspiration calculation mainly accounts for five kinds of evapotranspiration, including the amount of evaporation on an impervious area E 1 , tension water evapotranspiration in the upper zone E 2 , free water evapotranspiration in upper zone E 3 , evapotranspiration on the water area of the watershed E 4 , tension water evapotranspiration in the lower zone E 5 .
It is a continuum model including 17 parameters, in which PCTIM, ADIMP and SARVA are three parameters of the direct runoff, UZTWM, UZFWM and UZK are three parameters of the water capacity of upper zone, ZPERC and REXP are the parameters of percolation, PFREE, LZTWM, LZFSM, LZSK, LZFPM, LZPK, RSERV, SIDE and SSOUT are the 9 parameters of water capacity of lower zone, and it is mainly adapted to semi-arid and semi-humid regions. Table 1 lists the value range and physical meaning of parameters in Sacramento model. Objective functions are used to evaluate the consistency of simulated flow with measured flow, and different objective functions correspond to different characteristics in hydrologic process. In this paper we choose the following objective function based on comprehensive considerations of overall water quantity balance and the consistency with flow process: In Equation (1), Q r is the parameter to control the flow magnitude (threshold value); Q 0 is the measured flow and Q 0 is the average measured flow; Q c is the calculated flow from the model and Q c is the average calculated flow; n 1 is the number of flood occurrences of which Q 0 > Q r , and n 1 is the number of flood occurrences of which Q 0 ≤ Q r ; n = n 1 + n 2 is the length of calibration data; α is the weight coefficient for different objective functions. The Equation (1) shows that different Q r and α can be used for different application requirements.

Study Basin
We choose the Dongyang river basin above Balihutong hydrological station and the Tantou river basin above Yihetantou hydrological station between the Yellow River Huayuankou and Sanmenxia interval, as the study water basin, which are shown in Figure 3.
The Dongyang basin is a semi-humid and semi-arid region with the area of 537 km 2 and average annual rainfall of 550 mm, including several important rainfall stations such as Shibandao, Hengheng and Huangyuan. In this paper we take the average rainfall of the three stations as the rainfall of the entire Dongyang basin, and take the flow of Balihutong station as the basin outflow. Also, the Tantou basin is a semi-humid and semi-arid region with the area of 1695 km 2 and average annual rainfall of 580 mm, including several important rainfall stations such as Taowan, Miaozi, Baitu, and Tantou. Similarly, the average rainfall of the four stations is taken as the rainfall of the entire Tantou basin, and the flow of Tantou station is taken as the basin outflow. The Dongyang basin is a semi-humid and semi-arid region with the area of 537 km 2 and average annual rainfall of 550 mm, including several important rainfall stations such as Shibandao, Hengheng and Huangyuan. In this paper we take the average rainfall of the three stations as the rainfall of the entire Dongyang basin, and take the flow of Balihutong station as the basin outflow. Also, the Tantou basin is a semi-humid and semi-arid region with the area of 1695 km 2 and average annual rainfall of 580 mm, including several important rainfall stations such as Taowan, Miaozi, Baitu, and Tantou. Similarly, the average rainfall of the four stations is taken as the rainfall of the entire Tantou basin, and the flow of Tantou station is taken as the basin outflow.

The Generation of Artificial Data (the Ideal Data)
In order to judge whether the ICSA method is practical and available for the Sacramento model parameter calibration, we use the artificial data to calibrate the parameter for validation. The detailed procedure is showed as below.
Step 1: Suppose a set of parameter value as the real value of the model parameter, and take the observed rainfall data as the input of model, then we can obtain a set of ideal flow data, that is to say, the artificial ideal flow data.
Step 2: Take the observed rainfall data as the input, and take the corresponding ideal flow data as output of model to calibrate the parameter, then we can obtain a set of calculated parameter values and corresponding calculated flow data.
Step 3: The real parameter is compared with the calculated parameter, and the ideal flow data is compared with the calculated flow data, then conduct analysis and evaluation.
The framework diagram of artificial ideal data generation is showed in Figure 4.

Parameter Estimation
In order to compare the performances of different algorithms in parameter estimation of Sacramento model, we assume a series of artificial generated data to be the real values of the parameters for the validation of parameter estimation. The study basins are aforementioned Dongyang and Tantou river basins. Then we utilize the rainfall and evaporation data of two floods (No.820801 and No.820730) and Sacramento model to simulate these two floods, and employ SCE-UA, PGA, SMSE-PSO, ICSA to calibrate parameters. Our PC infrastructure is that the CPU is

Sacramento Model
The observed rainfall data The assumed real value of Sacramento model parameter The artificial ideal flow data

The Generation of Artificial Data (the Ideal Data)
In order to judge whether the ICSA method is practical and available for the Sacramento model parameter calibration, we use the artificial data to calibrate the parameter for validation. The detailed procedure is showed as below.
Step 1: Suppose a set of parameter value as the real value of the model parameter, and take the observed rainfall data as the input of model, then we can obtain a set of ideal flow data, that is to say, the artificial ideal flow data.
Step 2: Take the observed rainfall data as the input, and take the corresponding ideal flow data as output of model to calibrate the parameter, then we can obtain a set of calculated parameter values and corresponding calculated flow data.
Step 3: The real parameter is compared with the calculated parameter, and the ideal flow data is compared with the calculated flow data, then conduct analysis and evaluation.
The framework diagram of artificial ideal data generation is showed in Figure 4. The Dongyang basin is a semi-humid and semi-arid region with the area of 537 km 2 and average annual rainfall of 550 mm, including several important rainfall stations such as Shibandao, Hengheng and Huangyuan. In this paper we take the average rainfall of the three stations as the rainfall of the entire Dongyang basin, and take the flow of Balihutong station as the basin outflow. Also, the Tantou basin is a semi-humid and semi-arid region with the area of 1695 km 2 and average annual rainfall of 580 mm, including several important rainfall stations such as Taowan, Miaozi, Baitu, and Tantou. Similarly, the average rainfall of the four stations is taken as the rainfall of the entire Tantou basin, and the flow of Tantou station is taken as the basin outflow.

The Generation of Artificial Data (the Ideal Data)
In order to judge whether the ICSA method is practical and available for the Sacramento model parameter calibration, we use the artificial data to calibrate the parameter for validation. The detailed procedure is showed as below.
Step 1: Suppose a set of parameter value as the real value of the model parameter, and take the observed rainfall data as the input of model, then we can obtain a set of ideal flow data, that is to say, the artificial ideal flow data.
Step 2: Take the observed rainfall data as the input, and take the corresponding ideal flow data as output of model to calibrate the parameter, then we can obtain a set of calculated parameter values and corresponding calculated flow data.
Step 3: The real parameter is compared with the calculated parameter, and the ideal flow data is compared with the calculated flow data, then conduct analysis and evaluation.
The framework diagram of artificial ideal data generation is showed in Figure 4.

Parameter Estimation
In order to compare the performances of different algorithms in parameter estimation of Sacramento model, we assume a series of artificial generated data to be the real values of the parameters for the validation of parameter estimation. The study basins are aforementioned Dongyang and Tantou river basins. Then we utilize the rainfall and evaporation data of two floods (No.820801 and No.820730) and Sacramento model to simulate these two floods, and employ SCE-UA, PGA, SMSE-PSO, ICSA to calibrate parameters. Our PC infrastructure is that the CPU is

Sacramento Model
The observed rainfall data The assumed real value of Sacramento model parameter The artificial ideal flow data

Parameter Estimation
In order to compare the performances of different algorithms in parameter estimation of Sacramento model, we assume a series of artificial generated data to be the real values of the parameters for the validation of parameter estimation. The study basins are aforementioned Dongyang and Tantou river basins. Then we utilize the rainfall and evaporation data of two floods (No.820801 and No.820730) and Sacramento model to simulate these two floods, and employ SCE-UA, PGA, SMSE-PSO, ICSA to calibrate parameters. Our PC infrastructure is that the CPU is Intel Core i7-3537U CPU, the CPU Clock Speed is 2.50 GHz, the Operating System is Windows 8.1, and the Programming Language we used is Java. We set the initial population size n = 80, the largest iterations k = 300, the precision of coding ε = 0.001, Q r = 800 m 3 /s, α = 0.6, so as to ensure the comparison of the four algorithms is under the same condition. We run each algorithm for 5 times and obtain the optimal results according to objective functions. Detailed results are shown in Tables 2 and 3.
Here we use the Nash-Sutcliffe model efficiency coefficient (NSE) [26] as the evaluation criterion, and it is calculated as: (2) in which N denotes the number of flood time periods, Q sim,t denotes the simulated flow, Q obs,t denotes the measured flow, Q obs denotes the measured average flow.  Table 2 gives out the parameter estimation results of the artificial generated data. We can see that the calibration values are very close to the assumed real values with small errors, indicating that the optimal results of all four algorithms can converge to real values.  Table 3 clearly demonstrates the parameter calibration performances of the four algorithms. All four algorithms have converged to real values with NSE approaching very close to 1. The running time and iterations directly reflects the efficiency and convergence rate of the employed algorithm. For each algorithm, the running time and iterations of No.820801 Flood is smaller than that of No.820703 Flood, which is resulted from shorter last time and thus simpler calculation of the No.820801 Flood. As for the same flood, the running time and iterations of SCE-UA are the largest among four algorithms, while that of the proposed ICSA are the smallest. Clearly, the rank of efficiency and convergence rate of the four algorithms is ICSA > SMSE-PSO > PGA > SCE-UA. Besides, the ICSA owns the smallest objective function values in both two floods (2.1943 and 8.4972, respectively), implying the best consistency of the simulated results and real values.
From Tables 2 and 3 we can conclude that ICSA shows the best parameter estimation performance among all four algorithms in the experimental Sacramento model, with the highest efficiency, precision, and convergence rate. Besides, we can find that there exists some difference between the parameter calibration results obtained from each algorithm. Although the difference is quite small and the artificial ideal data were used to calibrate the parameter, it still leads to a large difference in the objective function value, and the NSE does not reach the optimal value yet. This may be caused by these two reasons, on the one hand, it may be a problem of optimal algorithm, for example, the algorithm is not good enough, or the setting termination condition is not reasonable (the maximum iterations is too small). On the other hand, the uncertainty of the model structure may be the primary cause, particularly the parameter uncertainty, in other words, although the ideal data is used for parameter calibration, if the model structure has any fault, it will still result in errors. Therefore, we will give a detailed analysis of parameter uncertainty in the next section.

The Analysis of Parameter Uncertainty
The parameter uncertainty analysis consists of two contents, the sensitivity analysis and the equifinality, and this paper will also analyze the parameter uncertainty through the two aspects.

The Parameter Sensitivity Analysis
The sensitivity analysis is to study the model response caused by the model parameter variation, and focus on the consideration about the contribution of model parameter variation to the model output, in order to identify whether the parameters have a great impact on the model response or not. The sensitivity analysis method applied to hydrologic model can be generally divided into the local sensitivity analysis method and global sensitivity analysis method.
Local sensitivity analysis (LSA) method pays attention to the impact on the model output caused by single parameter, with the advantage of simple, fast and strong operability, among which the perturbation analysis method [27] is widely used because of its simple and practical character. However, the variation of one parameter value may impact the sensitivity of the other parameters, due to the poor consideration about the interaction of all parameters, the LSA method is lack of stability and there are also some limitations.
While in the global sensitivity analysis (GSA) method, the impact on the model output caused by the change of multiple parameters and the interaction of all parameters can be simultaneously taken into account, and we can obtain a comprehensive understanding about the sensitivity of each parameter, which is applicable for the hydrologic model with multiple parameters. For example, the GLUE method [28], the Sobol method [22,29], the RSA method [30] and the RSMSobol method [31], all of them are the commonly used GSA method.
In this paper, take the two floods showed in Table 2 as examples, we first apply the perturbation analysis method to the 17 parameters of Sacramento model for single-parameter sensitivity analysis, that is to say, suppose that all of parameters are at their true value, just change one parameter value according to its value range at a time, fixing the remaining 16 parameters invariable, then the corresponding NSE can be obtained and the NSE curve as shown in Figure 5.   Figure 5. The parameter sensitivity analysis results based on the perturbation analysis method.
In the above figures, red curve represents the simulated results obtained from the No.820801 flood in Doyang river basin, and the blue curve represents that of the No.820730 flood in Tandou river basin. As to these two floods, the range of corresponding NSE changes little with the variation of the parameter. For the No.820801 flood (the red curve), we can find that the parameter whose corresponding NSE with largest variation range is UZFWM, it is varying from 0.84 to 1, and the corresponding NSE of the other 16 parameters are all no less than 0.955. Thus, the parameter UZFWM is the most sensitive parameter to this flood, while the parameter PCTIM, ADIMP, SARVA, UZK, ZPERC, REXP, LZTWM, LZFPM, PFREE, RSERV, SIDE and SSOUT are almost insensitive, all of them has the NSE which is no less than 0.994. For the No.820730 flood (the blue curve), the parameter UZFWM (with NSE between 0.85 to 1) and LZSK (with NSE between 0.88 to 1) are the parameters whose NSE has largest variation range, so these two parameters are more sensitive than others to this flood. The corresponding NSE of the remaining parameters are all larger than 0.935, particularly the parameter PCTIM, ADIMP, SARVA, PFREE, RSERV, SIDE and SSOUT has the NSE no less than 0.9988, so these parameters are almost insensitive to this flood.
Besides, the Monte Carlo Sampling based GLUE method is also used for parameter sensitivity analysis from a global perspective, which means all of the 17 parameters will be changed according to their own value range, 10,000 groups of parameters will be selected randomly to conduct simulation, and the parameter value and its corresponding NSE are showed in Figure 6. insensitive, all of them has the NSE which is no less than 0.994. For the No.820730 flood (the blue curve), the parameter UZFWM (with NSE between 0.85 to 1) and LZSK (with NSE between 0.88 to 1) are the parameters whose NSE has largest variation range, so these two parameters are more sensitive than others to this flood. The corresponding NSE of the remaining parameters are all larger than 0.935, particularly the parameter PCTIM, ADIMP, SARVA, PFREE, RSERV, SIDE and SSOUT has the NSE no less than 0.9988, so these parameters are almost insensitive to this flood. Besides, the Monte Carlo Sampling based GLUE method is also used for parameter sensitivity analysis from a global perspective, which means all of the 17 parameters will be changed according to their own value range, 10,000 groups of parameters will be selected randomly to conduct simulation, and the parameter value and its corresponding NSE are showed in Figure 6.   figure), most of parameters value have corresponding NSE between 0.8 and 1, and the parameter UZFWM is still the most sensitive one to this flood, the corresponding NSE will decrease when it's greater than 25. Meanwhile, the increase of the parameter ZPERC, LZFSM, LZFPM and LZSK will result in a decrease of NSE, while the parameter REXP still holds an opposite tendency.
Note that all the simulated results were just applicable for the particular two floods. Based on the analysis results of the above two methods, while the ideal data are used for parameter calibration, we can find that the sensitivity of the same parameter is not the same when applied to the different basins or different floods. That is to say, the sensitivity of the same parameter is affected not only by flood conditions, but also by the different basins applied. In addition, the results obtained from the LSA method and the GSA method are also not the same, while multiple parameters are simultaneously taken into consideration, the single-parameter sensitivity may be changed because of the interaction or the complementarity between the parameters.

The Model Equifinality Analysis
The equifinality [32] of the model means that for the same model structure and model input, there will be multiple groups of optimal parameters which make the model output have the same accuracy. The reasons causing equifinality may be the data uncertainty, the model structure uncertainty, the complementarity between the model parameters, the randomness of the parameters, and the multi-extremum of objective function, etc.
As for the ideal data, Duan et al. [2] assumed the model is correct, the model input and output will be non-error, and there will not exist equifinality problem under this situation. While some other hydrologists like Beven et al. [33] hold the different opinion, they believe that even under the ideal condition and both the input and output data are correct, it may still results in the equifinality problem if there exists any error in the assumed model structure.
Like the opinion of Beven et al. we believe that when the ideal data are used for parameter calibration of different floods, if the assumed model structure is non-correct, it will cause the error in results. However, we can not prove that whether the model structure is correct or not. Here we assume that for a flood with the ideal data for parameter estimation, if there is no error in the input and output data and the model we chose is applicable and proven by many scholars, the , we can see that most of points fall in the area in which 0.9 ≤ NSE ≤ 1, the most sensitive of all is the parameter UZFWM, the corresponding NSE of which obviously decrease when the parameter value is greater than 30. In addition, with the increase of the parameter value of UZTWM, ZPERC, LZFSM, LZFPM, LZSK and UZFWM, the probability of their NSE falling into the area of 0.9 ≤ NSE ≤ 1 will decrease, while the parameter REXP has the opposite tendency. For the No.820730 flood (blue figure), most of parameters value have corresponding NSE between 0.8 and 1, and the parameter UZFWM is still the most sensitive one to this flood, the corresponding NSE will decrease when it's greater than 25. Meanwhile, the increase of the parameter ZPERC, LZFSM, LZFPM and LZSK will result in a decrease of NSE, while the parameter REXP still holds an opposite tendency.
Note that all the simulated results were just applicable for the particular two floods. Based on the analysis results of the above two methods, while the ideal data are used for parameter calibration, we can find that the sensitivity of the same parameter is not the same when applied to the different basins or different floods. That is to say, the sensitivity of the same parameter is affected not only by flood conditions, but also by the different basins applied. In addition, the results obtained from the LSA method and the GSA method are also not the same, while multiple parameters are simultaneously taken into consideration, the single-parameter sensitivity may be changed because of the interaction or the complementarity between the parameters.

The Model Equifinality Analysis
The equifinality [32] of the model means that for the same model structure and model input, there will be multiple groups of optimal parameters which make the model output have the same accuracy. The reasons causing equifinality may be the data uncertainty, the model structure uncertainty, the complementarity between the model parameters, the randomness of the parameters, and the multi-extremum of objective function, etc.
As for the ideal data, Duan et al. [2] assumed the model is correct, the model input and output will be non-error, and there will not exist equifinality problem under this situation. While some other hydrologists like Beven et al. [33] hold the different opinion, they believe that even under the ideal condition and both the input and output data are correct, it may still results in the equifinality problem if there exists any error in the assumed model structure.
Like the opinion of Beven et al. we believe that when the ideal data are used for parameter calibration of different floods, if the assumed model structure is non-correct, it will cause the error in results. However, we can not prove that whether the model structure is correct or not. Here we assume that for a flood with the ideal data for parameter estimation, if there is no error in the input and output data and the model we chose is applicable and proven by many scholars, the equifinality effect caused by model structure will be less than the error caused by the optimization algorithm (optimization ability and optimization condition setting).
As shown in Tables 2 and 3, in order to compare the pros and cons of the four algorithms, the initial population size and the maximum iterations we set in our algorithms were only n = 80 and k = 300. In order to validate the above assumptions, here we change them to n = 500 and k = 1000, and observe the results of simulation.
As shown above in Figure 7, when the setting condition of the optimization algorithm were changed, the value of objective function obtained by the ICSA will eventually stabilized at around 0.0035 and 0.0043, respectively.
Water 2017, 9,495 12 of 21 equifinality effect caused by model structure will be less than the error caused by the optimization algorithm (optimization ability and optimization condition setting). As shown in Tables 2 and 3, in order to compare the pros and cons of the four algorithms, the initial population size and the maximum iterations we set in our algorithms were only n = 8 0 and k = 300. In order to validate the above assumptions, here we change them to n = 500 and k = 1000, and observe the results of simulation.
As shown above in Figure 7, when the setting condition of the optimization algorithm were changed, the value of objective function obtained by the ICSA will eventually stabilized at around 0.0035 and 0.0043, respectively. Therefore, we can draw a conclusion that the main reason resulting in the difference of the parameters and the objective function values is the optimization condition setting of the optimal algorithm, which makes the objective function do not converge to the optimal value, while the equifinality is not the main reason.

Multi-Step Parameter Estimation Method
The Equation (1) gives out the traditional single-step parameter estimation method. Since each parameter of the 17 parameters in Sacramento model has its own physical meaning and combines to other ones, different types of the parameters should use corresponding estimation principles. As a result, the simulated results of single-step method often have large deviation with the actual measured results. Although using different objective functions could raise simulation accuracy, the improvement is rather limited [34]. Therefore, in this paper we present the multi-step parameter estimation method as solution. The parameter estimation of Sacramento model needs three steps as follows. Therefore, we can draw a conclusion that the main reason resulting in the difference of the parameters and the objective function values is the optimization condition setting of the optimal algorithm, which makes the objective function do not converge to the optimal value, while the equifinality is not the main reason.

Multi-Step Parameter Estimation Method
The Equation (1) gives out the traditional single-step parameter estimation method. Since each parameter of the 17 parameters in Sacramento model has its own physical meaning and combines to other ones, different types of the parameters should use corresponding estimation principles. As a result, the simulated results of single-step method often have large deviation with the actual measured results. Although using different objective functions could raise simulation accuracy, the improvement is rather limited [34]. Therefore, in this paper we present the multi-step parameter estimation method as solution.
The parameter estimation of Sacramento model needs three steps as follows.
Step 1: The following objective function is utilized for all parameters: in which Q sim,t denotes the simulated flow, and Q obs,t denotes the measured flow.
Step 1 emphasizes the underground runoff related parameters in the flood process, including LZTW M, LZFSM, LZFPM, LZSK, LZPK, PFREE. Meanwhile, estimation for all parameters benefits the relativity weakening of the surface and underground water runoff related parameters.
Step 2: Set the underground runoff related parameters obtained in Step 1 as constants, and use the following objective functions for parameter estimation: in which Q sim,t and Q obs,t are the simulated and measured flows, n is the number of simulation days or number of flood time periods.
Step 2 emphasizes the surface water runoff related parameters in the flood process, including UZTW M, UZFW M, UZK, ADI MP, ZPERC, REXP.
Step 3: Set the surface water runoff related parameters obtained in Step 2 as constants, and use the Equation (3) as objective function.

Results
Combining the ICSA and multi-step method discussed in Chapter 2, we present the ICSA based multi-step parameter estimation method in Sacramento model. We pick the data of 10 floods of Dongyang basin    Table 4. Table 4. Results of parameters estimation of single-step (shorten as single) and multi-step (shorten as multi) method. The flood forecasting accuracy is determined by the results of parameter estimation and also reflects the effect of the parameter estimation method. Here we use three evaluation criterions, i.e., time root-mean-square of flow (TRMS), NSE, and bias (%BIAS). Lower TMRS reflects better dispersion degree, and is calculated as: The %BIAS is calculated as: Then we use the remaining floods data to validate and evaluate flood forecasting accuracy and thus the parameter estimation results by above three criterions. The flood simulation results (runoff & peak flow) and the accuracy statistics of both single-step simulation and multi-step simulation are illustrated in Figure 8. Note that the "*" represents flood for validation.
Water 2017, 9,495 14 of 21 Then we use the remaining floods data to validate and evaluate flood forecasting accuracy and thus the parameter estimation results by above three criterions. The flood simulation results (runoff & peak flow) and the accuracy statistics of both single-step simulation and multi-step simulation are illustrated in Figure 8. Note that the "*" represents flood for validation. (e) (f) (g) Figure 8. The criteria comparison between the two methods of the two river basins. (a) The runoff absolute error comparison between the two methods of the two basins; (b) The runoff relative error comparison between the two methods of the two basins; (c) The peak flow absolute error comparison between the two methods of the two basins; (d) The peak flow relative error comparison between the two methods of the two basins; (e) The TRMS comparison between the two methods of the two basins; (f) The NSE comparison between the two methods of the two basins; (g) The %BIAS comparison between the two methods of the two basins.
From the above Figure 8a-d, we can find that, for the Dongyang river basin, the absolute and relative errors of runoff, and the absolute and relative errors of peak flow obtained from the single-step method are below 1.71 mm, ±4.88%, 6.7 mm and ±3.84%, and those criteria obtained from the multi-step method are 1.1 mm, ±3.59%, 7 mm and ±4.22%. Besides, for the Tantou river basin, those criteria obtained from the single-step method are below 5.3 mm, ±14.38%, 156.39 mm and ±43.81%, and those criteria obtained from the multi-step method are below 2.83 mm, ±10.88%, 54.13 mm and ±15.08%. Therefore, we can see that the results of flood *870605 in Tantou river basin obtained from single-step method is not qualified.
Then from the Figure 8e, for the 10 floods in Dongyang river basin, the TRMS of 7 floods obtained from the multi-step method is better than the single-step. Besides, for the 32 floods in Tantou river basin, the TRMS of 20 floods obtained from the multi-step method is better.
As we can see from Figure 8f-g, for Dongyang river basin, the %BIAS from single-step method has the maximum larger than ±10% and the average value of 7.47%, the %BIAS below ±3% account From the above Figure 8a-d, we can find that, for the Dongyang river basin, the absolute and relative errors of runoff, and the absolute and relative errors of peak flow obtained from the single-step method are below 1.71 mm, ±4.88%, 6.7 mm and ±3.84%, and those criteria obtained from the multi-step method are 1.1 mm, ±3.59%, 7 mm and ±4.22%. Besides, for the Tantou river basin, those criteria obtained from the single-step method are below 5.3 mm, ±14.38%, 156.39 mm and ±43.81%, and those criteria obtained from the multi-step method are below 2.83 mm, ±10.88%, 54.13 mm and ±15.08%. Therefore, we can see that the results of flood *870605 in Tantou river basin obtained from single-step method is not qualified.
Then from the Figure 8e, for the 10 floods in Dongyang river basin, the TRMS of 7 floods obtained from the multi-step method is better than the single-step. Besides, for the 32 floods in Tantou river basin, the TRMS of 20 floods obtained from the multi-step method is better.
As we can see from Figure 8f-g, for Dongyang river basin, the %BIAS from single-step method has the maximum larger than ±10% and the average value of 7.47%, the %BIAS below ±3% account for 70%; the %BIAS from multi-step method has the average value of 3.07, and the values below ±3% account for 80%. For Tantou river basin, the %BIAS from single-step method has the maximum nearly ±15% and the average value of 2.75, the %BIAS below ±5% account for 68.75%; the %BIAS from multi-step method has the average value of 2.05, and the values below ±5% account for 71.88%.
We compare the error results of single-step and multi-step flood simulations with the permissible error of "Forecasting norm for hydrology intelligence of China" (SL250-2000). For the Dongyang river basin, the evaluation criterions of single-step and multi-step simulations are all within the accepted error range (±20%), thus the model parameters can be recognized as basically reasonable. Note that, all evaluation criterions of multi-step simulations are better than that of single-step simulations. Through the comparison of NSE, it is found that among the ten simulated floods for validation, seven floods reach the first level and three reach the second level (according to SL250-2000) for both single-step and multi-step methods. That is to say, both single-step and multi-step method for Dongyang river basin reach qualified level.
However, for the Tantou river basin, the relative error of flood flow peak of No.870605 flood in single-step simulation exceeds the permissible error, indicating that the flood forecasting of No.870605 by single-step method is unqualified. Also, the other evaluation criterions of single-step method are all worse than that of multi-step method. Through the comparison of NSE, we can see that among the 32 simulated floods, for single-step method, 15 floods reach the first level, 16 reach the second level, and 1 is unqualified; for multi-step method, 17 floods reach the first level, 15 reach the second level, and none is unqualified.
The simulation results comparison of flood forecasting between single-step method and multi-step method is showed in Figure 9. Pre1 means the simulation results which obtained from the multi-step method, while Pre2 means the simulation results which obtained from the single-step method.
for 70%; the %BIAS from multi-step method has the average value of 3.07, and the values below ±3% account for 80%. For Tantou river basin, the %BIAS from single-step method has the maximum nearly ±15% and the average value of 2.75, the %BIAS below ±5% account for 68.75%; the %BIAS from multi-step method has the average value of 2.05, and the values below ±5% account for 71.88%.
We compare the error results of single-step and multi-step flood simulations with the permissible error of "Forecasting norm for hydrology intelligence of China" (SL250-2000). For the Dongyang river basin, the evaluation criterions of single-step and multi-step simulations are all within the accepted error range (±20%), thus the model parameters can be recognized as basically reasonable. Note that, all evaluation criterions of multi-step simulations are better than that of single-step simulations. Through the comparison of NSE, it is found that among the ten simulated floods for validation, seven floods reach the first level and three reach the second level (according to SL250-2000) for both single-step and multi-step methods. That is to say, both single-step and multi-step method for Dongyang river basin reach qualified level.
However, for the Tantou river basin, the relative error of flood flow peak of No.870605 flood in single-step simulation exceeds the permissible error, indicating that the flood forecasting of No.870605 by single-step method is unqualified. Also, the other evaluation criterions of single-step method are all worse than that of multi-step method. Through the comparison of NSE, we can see that among the 32 simulated floods, for single-step method, 15 floods reach the first level, 16 reach the second level, and 1 is unqualified; for multi-step method, 17 floods reach the first level, 15 reach the second level, and none is unqualified.
The simulation results comparison of flood forecasting between single-step method and multi-step method is showed in Figure 9. Pre1 means the simulation results which obtained from the multi-step method, while Pre2 means the simulation results which obtained from the single-step method. (e) (f)

Discussion
In this paper, both the artificial generated data and the actual flood data are used to calibrate the parameter of Sacramento model, although we have improved the flood prediction accuracy by analyzing the parameter calibration process under these two conditions, there are still some difficult issues required further discussion.
(1) We first use the ICSA and the other three optimization algorithms for Sacramento model parameter calibration based on the artificial generated data, the initial aim of which is to validate the practicability and better performance of ICSA by comparing the simulated results of the four algorithms. During this process, although the artificial generated ideal data are used for parameter calibration, there still exists some bias in the parameter value, objective functions and the corresponding NSE, which results in the search for the origin of the bias. Therefore, the analysis of the bias sources under the ideal condition should be an aspect that needs to be further studied. Many scholars have validated their own opinions from the view of their specific example, but there

Discussion
In this paper, both the artificial generated data and the actual flood data are used to calibrate the parameter of Sacramento model, although we have improved the flood prediction accuracy by analyzing the parameter calibration process under these two conditions, there are still some difficult issues required further discussion.
(1) We first use the ICSA and the other three optimization algorithms for Sacramento model parameter calibration based on the artificial generated data, the initial aim of which is to validate the practicability and better performance of ICSA by comparing the simulated results of the four algorithms. During this process, although the artificial generated ideal data are used for parameter calibration, there still exists some bias in the parameter value, objective functions and the corresponding NSE, which results in the search for the origin of the bias. Therefore, the analysis of the bias sources under the ideal condition should be an aspect that needs to be further studied. Many scholars have validated their own opinions from the view of their specific example, but there is no confirmation of others' results by using their own view. This indicates that different models, basins and input conditions may have certain impacts on the bias, in one case, the bias caused by the model structure is larger than that caused by the optimization algorithm, whereas in another case it can be the opposite. Thus, the causes of the bias should be analyzed according to the specific application conditions.
(2) While in the actual flood simulation and prediction, the ICSA is combined with the multi-step parameter calibration method, and multiple criteria are used for evaluation, which have improved the prediction accuracy of Sacramento model. However, we haven't made the analysis of the prediction uncertainty, the main reason is that we have made fewer studies in this field, and there are only some immature views. First of all, the current parameter uncertainty analysis are mainly focused on the simulation process based on daily streamflow data [22][23][24]35] and then give the calibrated parameter set, while the probability forecasting method are used for the flood event simulation process [36], in which the parameter uncertainty analysis are avoided. The flood event model parameter calibration is to select a group of parameters from multiple floods by using a kind of evaluation criterion, although the group of parameters is receivable to the multiple floods, it may not be credible for one of the floods. From the view of equifinality, this kind of parameter should also be parameter set, as for the multiple different floods, the parameter variation range will be larger than that of the single flood, thus, it's difficult to simultaneously analyze the parameter uncertainty of multiple floods. Because of the model equifinality, the model parameter become the parameter set, in order to reduce the selection range of the parameter and improve the parameter selection efficiency, Guo et al. [37][38][39] proposed to adopt the multi-objective algorithms for parameter calibration and uncertainty analysis. This method provides an idea for us to solve the model equifinality problem of multiple floods, we can utilize multiple objective functions and criteria to evaluate the simulation results of each flood, and simultaneously adopt multiple criteria to conduct comprehensive evaluation of the simulated results of all floods. Then the double-evaluation will be formed to reduce the space of parameter set so that we can improve the parameter selection accuracy.

Conclusions
According to the analysis and studies about the commonly used algorithms and methods for parameter calibration of Sacramento model, we can find that there are still some problems in the practical application of the optimization algorithms, for example, the premature convergence, easy to fall into local optimum, or slow convergence. The traditional single-step parameter calibration method also has disadvantages of lower simulation and prediction accuracy, and poor performance. In view of all these problems, this paper presented a Multi-step parameter calibration method based on the ICSA for Sacramento model in order to improve the accuracy and stability of parameter calibration.
First of all, we give a brief introduction about the ICSA and Sacramento model basic theory, and the 17 parameters of the Sacramento model as well as their value range are also showed in the paper. Besides, according to the different flow magnitude, different objective functions of optimization algorithm for Sacramento model parameter calibration are set up. Meantime, combined with the specific characteristics of the Sacramento model, the ICSA is used to solve the parameter calibration problem for Sacramento model, and the detailed solution steps are proposed.
Then, take the Dongyang river basin and Tantou river basin as typical study basin, and a flood is selected from each of basins for parameter calibration based on the artificial generated data. The detailed procedure of artificial data generation is also introduced in this paper, we adopt the corresponding rainfall, evapotranspiration and the initial soil moisture data of the flood as the input of Sacramento model, and assumed a group of parameter value as the real value of parameter to generate a group of flow data. After that, we still use the previous input data, and take the artificial generated flow data as the output to calibrate the parameters. Under the same condition, the ICSA, SMSE-PSO, PGA and SCE-UA are all used to solve the parameter calibration problem of the selected two floods based on the artificial generated data, and the objective functions and the NSE are taken as the evaluation criteria. The results showed that to these two floods, whether from the view of convergence speed, the objective functions or the NSE, the ICSA is the most excellent algorithm. Thus, we can validate that the ICSA is applicable and practical for Sacramento model parameter calibration and provide a better solution algorithm for the multi-step parameter calibration method.
Subsequently, aimed at the deviation of parameter value, NSE and objective functions showed in Tables 2 and 3, we still take the two floods in the table as examples to make further analysis about Sacramento model parameter uncertainty (from the view of parameter sensitivity and model equifinality). Based on the analysis result of perturbation analysis method and the GLUE method, while using the artificial generated ideal data for parameter calibration, the parameter sensitivity will be influenced not only by the specific conditions of different floods and the applied basins, but also by the interactions and complementarity between the parameters. In addition, the optimal conditions setting of ICSA algorithm for parameter calibration are improved, the objective functions of the two floods are found to be stabilized at around 0.0035 and 0.0043 respectively, they have been obviously improved compared to the results in Table 3. Therefore, as for the bias in Tables 2 and 3, our assumption can be validated that the error caused by the model equifinality are less than that caused by the optimization algorithm, so the model equifinality is not the main reason for this bias.
Finally, based on the above studies, we introduce the multi-step parameter calibration method and combine with the ICSA for Sacramento model parameter calibration. We selected 7 floods from Dongyang river basin and 25 floods from Tantou river basin for simulation, took 3 floods in Dongyang river basin and 7 floods in Tantou river basin for validation, and the TRMS, NSE, %BIAS are used as evaluation criteria. The results showed that the simulated results of multi-step method can fully satisfied the requirement, while the result of one flood calibrated by the single-step method is unqualified, and the simulation performance of multi-step method is obviously better than the single-step method. Besides, based on the 10 floods selected from the two basins, these two methods are used for flood forecasting simulation. It is easy to find that the simulated results obtained from the multi-step method are closer to the corresponding observed value, and the prediction accuracy is apparently higher compared to the single-step method. In a word, the combination of the ICSA and multi-step parameter calibration method for Sacramento model can not only improve the optimization ability of algorithm but also improve the accuracy and stability of the parameter calibration.