Parameter estimation of an improved nonlinear Muskingum model using a new hybrid method

In general, river flood routing methods can be categorized into three main distinct approaches: (1) the lumped models (e.g., linear and nonlinear Muskingum), (2) the distributed models (e.g., dynamic wave, diffusion wave, and kinematic wave models), and (3) the semi-distributed models (e.g., Muskingum–Cunge family models). The first approach is exclusively capable of routing the flood hydrograph of a river regardless of considering reach characteristics. The continuity equation accompanied with a specific storage function is utilized in this category. The true meaning of flood routing using these models is nothing but to compute the Muskingum parameters exclusively by using historical records. The distributed models utilize reach properties to route floods in any section of the reach by considering both mass conservation and transport momentum equations. Finally, the third category of models can be described as a tradeoff between the two first approaches since they predict floods along a river by using a simplified form of the Saint-Venant equations based on both physical concepts and river characteristics (Barati 2014).

Since the Muskingum model is one of the most widely used methods for flood routing, many studies have been conducted on improving both accuracy and efficiency of this model, particularly during the last two decades. However, based on the comprehensive literature review presented in the next section, more investigations seeking new improved Muskingum models are required in favor of achieving more rigorous flood routing results.

In this paper, a new Muskingum model was proposed which has nine parameters. Although all parameters of the proposed Muskingum model are constant during the flood period, this new model can be considered as a variable Muskingum model since the weighting parameter of this model varies with the inflow of the hydrograph. A new hybrid technique was also presented and applied for parameter estimation of the new Muskingum model. The new hybrid model is a combination of modified honey bee mating optimization (MHBMO) algorithm with generalized reduced gradient (GRG) algorithm. This hybrid method works in two consecutive steps in which the MHBMO algorithm performs in the first phase and the GRG algorithm acts in the second one. The parameters of the new Muskingum model were calibrated using the new hybrid model for two flood sets of data: (1) smooth and (2) multi-peak flood hydrographs. The objective function was the sum of the square of the deviations between the observed and routed outflows (SSQ). The obtained results show a reduction of 4.52% and 99.7%, respectively, with regard to the best reported results in the previous literature.

Studies on the Muskingum flood routing method can be categorized into two general groups. Those in the first category are intended to recommend different optimization techniques in favor of achieving more rigorous parameter estimation results. On the other hand, the second type attempts to somehow improve the Muskingum model or flood routing procedure. In the first category, numerous optimization algorithms have already been applied for Muskingum parameter estimation so far. These techniques generally can be categorized into four branches: (1) zero-order optimization algorithms, e.g., genetic algorithm (GA) (Mohan 1997), particle swarm optimization (PSO) (Chu & Chang 2009), and MHBMO (Niazkar & Afzali 2015a); (2) first-order optimization algorithms, e.g., Nelder–Mead simplex (NMS) (Barati 2011) and GRG (Barati 2013); (3) second-order optimization algorithms, e.g., Broyden–Fletcher–Goldfarb–Shanno (BFGS) (Geem 2006); and (4) hybrid optimization algorithms, e.g., shuffled frog leaping algorithm and Nelder–Mead simplex (SFLA-NMS) (Haddad et al. 2015a, 2015b), GA-GRG (Barati 2012).

The zero-order techniques exclusively utilize some algorithm parameters. For instance, the MHBMO algorithm has five independent control parameters: (1) size of the initial population (N_ipop), (2) the number of workers (N_Worker), (3) the speed of the queen at the start of the mating flight (S_max), (4) the speed of the queen at the end of the mating flight (S_min), and (5) the speed reduction factor (χ) (Niazkar & Afzali 2015a). This category initiates the parameter estimation process only when these controlling parameters are being evaluated. On the contrary, the first- and second-order models make the user guess one and two initial sets of values, respectively, for Muskingum parameters at the beginning of the calibration procedure, when the exact values are unknown and meant to be calibrated. The algorithm parameter evaluation in the zero-order methods and these initial-value vectors in the first- and second-order methods significantly affect the speed, and more importantly, the accuracy of Muskingum parameter estimation. The need for an initial guess for Muskingum parameters in advance of the calibration process is considered to be one of the challenges of the first- and second-order optimization algorithms. This challenge is almost addressed in the hybrid models. In these techniques, two algorithms are combined and worked one after the other in a two-phase process. In the first step, a zero-order algorithm commonly commences the parameter estimation process by only assuming delicate values for its controlling parameters. As this algorithm conducted the calibration process either after a specified number of iterations or until a stopping criterion is satisfied, the obtained results for Muskingum parameters are utilized as initial values for the corresponding parameters in a first- or second-order algorithm in the second step. Hence, not only do the hybrid techniques overcome the shortcomings of the first- and second-order algorithms, they appear to improve the accuracy of Muskingum parameter estimation in comparison with the zero-order methods as well.

The second-type efforts have tried to search for sufficient ways on how to obtain more precise flood routing results specifically as shown in the two sub-sections below.

These modifications were basically conducted by: (1) introducing new constant-parameter Muskingum models, e.g., four-parameter Muskingum model (Easa 2013b) and seven-parameter Muskingum model (Haddad et al. 2015b); and (2) utilizing variable-parameter Muskingum models, e.g., m-variable model (Easa 2013a) and VPNLMM-L (Karahan 2014).

Table 1 shows chronological improvement of the Muskingum flood routing model. In this table, α is a flow parameter, β is an exponent parameter, r is a dimensionless parameter for dividing flood period into sub-periods, and θ is a weighted parameter in weighted-average inflow modification.

Table 1

Chronological review of Muskingum flood routing models

No. .	Researchers .	Number of parameters .	Storage equation .
1	McCarthy (1938)	2
2	Chow (1959)	3
3	Gill (1978)	3	where m is the constant exponent parameter
4	Gavilan & Houck (1985)	4
5	O'Donnel (1985)	3	where βI is lateral flow and β is the lateral inflow factor
6	Easa (2013a)	2L + 1	where m(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter
7	Easa (2013b)	4
8	Vatankhah (2014a)	5
9	Karahan (2014)	10	where βI is lateral flow and Iave= [θIt+ (1-θ)It−1] and L = 2 (rising and falling inflow limbs)
10	Easa et al. (2014)	6
11	Easa et al. (2014)	2L + 2	where m(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter
12	Easa (2014)	4L + 1	where m(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter
13	Easa (2014)	2L + 3	where x(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter
14	Karahan et al. (2015)	5	where βI is lateral flow and Iave= [θIt+ (1-θ)It−1]
15	Easa (2015)	2L + 3	where ui=Ii/Imax and m(ui) =a+bln(1+cui) or m(ui) =a+bexp(-cui) and L is a user-defined parameter
16	Haddad et al. (2015b)	7
17	Zhang et al. (2016)	2L + 2	where m(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter

No. .	Researchers .	Number of parameters .	Storage equation .
1	McCarthy (1938)	2
2	Chow (1959)	3
3	Gill (1978)	3	where m is the constant exponent parameter
4	Gavilan & Houck (1985)	4
5	O'Donnel (1985)	3	where βI is lateral flow and β is the lateral inflow factor
6	Easa (2013a)	2L + 1	where m(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter
7	Easa (2013b)	4
8	Vatankhah (2014a)	5
9	Karahan (2014)	10	where βI is lateral flow and Iave= [θIt+ (1-θ)It−1] and L = 2 (rising and falling inflow limbs)
10	Easa et al. (2014)	6
11	Easa et al. (2014)	2L + 2	where m(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter
12	Easa (2014)	4L + 1	where m(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter
13	Easa (2014)	2L + 3	where x(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter
14	Karahan et al. (2015)	5	where βI is lateral flow and Iave= [θIt+ (1-θ)It−1]
15	Easa (2015)	2L + 3	where ui=Ii/Imax and m(ui) =a+bln(1+cui) or m(ui) =a+bexp(-cui) and L is a user-defined parameter
16	Haddad et al. (2015b)	7
17	Zhang et al. (2016)	2L + 2	where m(ui) is the variable exponent parameter and ui=rImax and L is a user-defined parameter

View Large

In Table 1, the simplest version of the flood routing Muskingum model was proposed by McCarthy (1938), which has two constant parameters, i.e., K and x. In this model, storage variation with weighted inflow and outflows is presumed to be linear. Since the corresponding relation is practically found to be nonlinear most times, Chow (1959), Gill (1978), and Gavilan & Houck (1985) separately introduced three different nonlinear Muskingum models. Among these nonlinear models, Gill's nonlinear Muskingum model is considered to be the most widely used variant of Muskingum models and much research has been specifically devoted to improve the accuracy of this model in the literature. In all the aforementioned models, it is assumed that no lateral flow is conveyed throughout the reach under investigation.

In contrast, O'Donnel (1985) ceased no-lateral-flow assumption and presented a new three-parameter Muskingum model (Table 1). In the recommended model, lateral flow is hypothesized to be not only a linear multiplication of inflow at each time step but also constantly distributed along the reach. After a distinct gap in the time table of Muskingum model improvement, a variable-exponent Muskingum model was presented by Easa (2013a). This model has the same formation as Gill's model while the exponent parameter in this model, i.e., m, can possess the desirable number of discrete values in the flood period. In this variable-exponent model, the flood period is divided into a user-defined number of sub-periods using a new non-dimensional parameter, i.e., r. The exponent parameter can have a different value in each sub-period. In other words, unlike previous models in which all Muskingum parameters are kept constant during the flood period, the exponent parameter is a discontinuously variable in this model. The lack of a standard approach for selecting an appropriate number of sub-periods is one of its critical shortcomings (Karahan 2014). The next modified version of the Muskingum flood routing model was presented by combining Chow's and Gill's nonlinear Muskingum models, which has four constant parameters (Easa 2013b). Vatankhah (2014a) proposed a new five-parameter Muskingum model using Gill's and Galvin and Houck's nonlinear models. The idea of considering a linear lateral flow with respect to inflow in each time interval in the linear Muskingum model of O'Donnell was extended for Gill's nonlinear model (Karahan et al. 2015). This new model, which has five parameters, was also proposed with a modification for flood routing scheme, which will be covered later in this section. Additionally, Karahan (2014) divided the flood period in respect to peak time inflow hydrograph into two sub-periods, and applied this model with different parameter values for the rising and falling limbs. Since he also utilized the modification for flood routing scheme proposed in Karahan et al. (2015), the model can be considered as a ten-parameter model in which Muskingum parameters can discontinuously change after inflow peak-time. A new constant six-parameter Muskingum model was suggested by Easa et al. (2014). The exponent-variable model (Easa 2013a) was combined with the four-parameter model (Easa 2013b) and a new exponent-variable model was introduced (Easa et al. 2014). Moreover, the obtained model was modified by Easa (2014) by extending the idea of variable-exponent model. The latter has the same formation as the four-parameter Muskingum model (Easa 2013b), in which each of the four parameters can discontinuously vary in the flood period. The technique utilized for dividing the flood period into several sub-periods is the same as the one suggested in Easa (2013a). Easa (2014) proposed a new variable-weighting parameter Muskingum model, which conceptually considers lateral flow the same as the one in Karahan et al.’s model (2015). The variable exponent-parameter model in Table 1 was modified, where the exponent parameter can continuously vary by logarithmic or exponential functions during flood period (Easa 2015). Each of these functions has three new parameters and the number of parameters for this model consequently is equal to five. Additionally, a new seven-parameter Muskingum model was proposed by Haddad et al. (2015b). More recently, Zhang et al. (2016) proposed a new Muskingum model which has a discontinuous variable exponent. This exponent varies in different sub-periods of the flood routing period in the same way as the one suggested by Easa (2013a). Although the recommended model considers the linear lateral flow, it has the same shortcomings as Easa's model since it is a modified form of Easa's model. Finally, it should be noted that Muskingum models with a greater number of parameters achieve better results since they efficiently provide more superior fitting capacity in comparison with Gill's nonlinear model (Table 1).

In favor of achieving better routing results, some researchers utilized I_t in this equation while others attempted to use I_t−1 instead. Moreover, an arithmetic average, 0.5(I_t+I_t−1) (Niazkar & Afzali 2015a) and weighted average, i.e., θI_t+ (1−θ)I_t−1 (Karahan et al. 2015) were suggested for substitution of I_t in Equation (3). Although the latter modification brings about a new parameter, i.e., θ, to the parameter calibration list, its assessment clearly shows that this modification performs sufficiently well for all utilized data sets.

Another modification in flood routing procedure was conducted in calculating the storage value. Traditionally, Euler's method was utilized for computation of new storage value (Tung 1985), while other explicit time integration schemes, e.g., improved Euler's method (Vatankhah 2014b), fourth-order Runge–Kutta method (Vatankhah 2014a, 2014b), Runge–Kutta–Fehlberg method (Vatankhah 2014b), were also applied for this purpose. Although the application of these higher order explicit schemes generally improves the accuracy of Muskingum flood routing results, they did not perform better than conventional Euler's method for all the utilized data sets. According to the state-of-the-art of the Muskingum model, more investigations on Muskingum models and new techniques for more accurate parameter estimation are still required.

Hence, the new Muskingum model has a variable inflow-weighted parameter.

This new model, unlike the 6th, 12th, 13th, and 14th models in Table 1, considers variable parameter, i.e., inflow-weighted parameter, without dividing the flood period into sub-periods. These modifications divide flood period into a user defined number of sub-periods, i.e., L, while selection of an appropriate number of sub-divisions is not standardly addressed in application of these models. This definitely produces confusion for users in practice. Not only the way in which these models divide flood period bringing additional parameter(s) to the list of parameters required to be calibrated, but also the technique utilized for this division is quite questionable. Moreover, although these modifications used a greater number of parameters, they did not perform as efficiently as more complex constant-parameter models. To be more specific, the 6th model in Table 1, a variable-exponent Muskingum method, for L = 10 has 21 parameters. It achieves the SSQ value equal to 24.203 by the Premium Solver for Wilson's data (Easa 2013a), whereas the 7th model in Table 1, a four-constant-parameter Muskingum model, obtains 7.67 for SSQ for the corresponding data by any algorithm (Easa 2013b). It also may be implied that the variation of the exponent parameter in the Muskingum model does not have a significant effect on the accuracy of routing results. In contrast, the proposed model not only considers variable inflow-weighted parameter, which is a function of inflow, but it also predicts floods with no requirement for flood period division. The applicability and performance of the proposed model in the simulation process is presented in the following sections.

In the routing procedure for the nine-parameter Muskingum model, Euler's method and the weighted average technique proposed by Karahan et al. (2015) are employed. In order to estimate the parameters of the proposed Muskingum model, the routing procedure can be conducted as follows.

Finally, the simulation procedure particularly includes the following steps.

Step 1: Determine the nine hydrologic parameters (K, x₁, α₁, x₂, α₂, x₃, α₃, m, and θ) by implementing the new hybrid optimization algorithm. The selection of values for these parameters can be obtained differently in the hybrid techniques as follows:

• In the first phase: The MHBMO algorithm randomly chooses initial values for the parameters in the defined bounds.

• In the second phase: The GRG algorithm utilizes the initial values of these parameters from the obtained result in the first phase.

Step 4: Compute the time rate of change of the channel storage using Equation (8).

It should be noted that negative values for reach storage is not acceptable and should be penalized in the routing procedure.

Step 6: Calculate the magnitude of the outflow at the next time using Equation (7). The third step to the sixth step should be repeated for all time intervals.

The HBMO algorithm was inspired from the mating process of honey bees, which makes this algorithm a phenomenon-mimicking algorithm. Not only is this algorithm a well-documented optimization technique, but it has also already been applied for several optimization problems in water resources engineering and especially Muskingum parameter estimation (Esmi Jahromi & Afzali 2014; Niazkar & Afzali 2015a, 2015b, 2016b; Afzali et al. 2016). The GRG algorithm, as a mathematical technique embedded in Excel spreadsheets, was also utilized for calibration of Gill's nonlinear Muskingum parameters (Barati 2013, Niazkar & Afzali 2016a). The new proposed hybrid model in this paper utilizes the combination of these two algorithms as a more powerful tool for efficient Muskingum parameter estimation. This new hybrid method not only overcomes the shortcomings of the zero-order algorithm, i.e., the MHBMO algorithm, but also provides a useful initial guess for the first-order method, i.e., the GRG algorithm. Therefore, the combination of these methods, i.e., the new hybrid model, performs more promisingly in finding optimum Muskingum parameter values.

This new hybrid approach works in two successive steps: (1) the MHBMO algorithm initiates the calibration process by assuming random values for Muskingum parameters and precedes the simulation process for several iterations; and (2) the GRG algorithm continues the routing procedure using the final values of Muskingum parameters obtained in the first phase.

In the first step of the new hybrid technique, a written code of the MHBMO algorithm in the MATLAB program was utilized. In this regard, the controlling parameters of the MHBMO algorithm were evaluated exactly as the ones utilized in Niazkar & Afzali's study (2015a). The utilized range for Muskingum parameters is shown in Table 2. Finally, the obtained results for Muskingum parameters in the first step were used as the initial vector for the second phase. In the second phase, the routing procedure was simulated in Excel spreadsheets. The built-in optimization tool in Excel, i.e., the GRG algorithm, was utilized for the second phase of calibration of Muskingum parameters. Since the GRG algorithm is a deterministic technique, the calibration process will terminate based on the inserted initial guess. The final values of Muskingum parameters will be achieved at the end of this step.

The utilized performance evaluation criteria adopted for verifying the efficiency of the hybrid MHBMO-GRG algorithm are as follows (Niazkar & Afzali 2015a; Afzali 2016).

The performance of the new hybrid method, i.e., MHBMO-GRG, in the Muskingum parameter estimation of the new proposed model was demonstrated using two case studies: (1) experimental (Wilson 1974) and (2) real (Viessman & Lewis 2003) data sets. Since these examples have already been applied in many studies on parameter estimation of Muskingum models, the applicability of the new Muskingum model potentially could be compared with most of the previous models in the literature.

In order to conduct sensitivity analysis of new Muskingum parameters in terms of SSQ, the ranges presented in Table 2 were considered using 0.001 intervals. In the achieved results, the minimum value for the selected parameter is obtained equal to the reported optimum values in Tables 3 and 5 for the first and second case studies, respectively. SSQ is relatively less sensitive to θ and more sensitive to x₁ and α₁ in the sensitivity analysis for the first case study. Furthermore, the result of sensitivity analysis for the second case study shows that SSQ is relatively less sensitive to α₂ and more sensitive to α₃. According to the results of sensitivity analysis for the two case studies, the greatest and the least parameter to which SSQ is sensitive differs for the selected case studies. Finally, it is necessary to note that the obtained results from this analysis clearly demonstrate that the parameters to which the final results are more sensitive depends upon the case study under consideration, and no specific parameter can be reported as the one to which the final results are mostly sensitive.

The Muskingum model is one of the most widely utilized hydrologic flood routing approaches although calibration of its parameters is considered to be challenging. In this study, a new hybrid optimization technique is introduced as a new tool to address this challenge. This technique combines the MHBMO and the GRG algorithms in two successive steps. In the first step, the MHBMO algorithm, as a zero-order method, commences the estimation process of Muskingum parameters for several specific iterations. Afterwards, the GRG algorithm, as a first-order approach, continues the calibration procedure using the obtained result from the first step as an improved initial guess for Muskingum parameters. This hybrid technique not only overcomes the requirement for appropriate initial guess for deterministic Excel-embedded optimization technique, but it also brings about more efficiency in the whole calibration process while it prevents sticking in probable local optimums by stochastic optimization techniques like the MHBMO algorithm. Additionally, a new improved version of Muskingum model, which comprises nine constant parameters, is proposed. Since the inflow-weighted parameter in this new Muskingum model varies with inflow of flood hydrograph, this model can be allocated as a variable-parameter Muskingum model, too. The new hybrid model was applied for parameter estimation of the new nonlinear Muskingum model using a simulation process based on Euler's method. The proposed hybrid algorithm was utilized to minimize the most common objective function for Muskingum calibration procedure for two case studies including smooth and double-peak flood hydrographs. The performance of the recommended algorithm for the new improved Muskingum model was compared with various utilized techniques and different Muskingum versions in the literature by applying several common criteria in flood routing studies. The hybrid MHBMO-GRG algorithm reached the best values for SSQ, SAD, MARE, and VarexQ among all other techniques for both case studies. For instance, this new hybrid method particularly reduced the best SSQ value available in the literature by 99.7% for the utilized case study. The obtained results demonstrate that the new hybrid method can be efficiently utilized for rigorous nonlinear Muskingum parameter estimation. Furthermore, the proposed modified Muskingum flood routing model appeared to be the most reliable version of Muskingum model among all versions which were considered in this study.