Although various techniques have been proposed to estimate the parameters of different versions of the Muskingum model, more rigorous techniques and models are still required to improve the computational precision of the calibration process. In this research, a new hybrid technique was proposed for Muskingum parameter estimation. Based on the conducted comprehensive literature review on the Muskingum flood routing models, a new improved Muskingum model with nine constant parameters was presented. Since the inflow-weighted parameter in the proposed model is a function of inflow hydrograph, it varies during the flood period and consequently can also be considered as a variable-parameter Muskingum model. The new hybrid technique was successfully applied for parameter estimation of the new version of Muskingum model for two case studies selected from the literature. Results were compared with those of other methods using several common performance evaluation criteria. The new Muskingum model significantly reduces the sum of the square of the deviations between the observed and routed outflows (SSQ) value for the double-peak case study. Finally, the obtained results indicate that not only the hybrid modified honey bee mating optimization-generalized reduced gradient algorithm somehow overcomes the shortcomings of both zero and first-order optimization techniques, but also the new Muskingum model appears to be the most reliable Muskingum version compared with the other methods considered in this study.

## INTRODUCTION

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).

*S*= total storage at time

_{t}*t*;

*I*and

_{t}*O*= inflow and outflow magnitudes at time

_{t}*t*, respectively;

*K*= storage parameter; and

*x*= a weighting parameter. The Muskingum parameters in the storage function (Equation (2)) are required to be calibrated using the historical flood events in the reach under consideration and these parameters are expected to capture the characteristics of that reach. These parameters can be utilized for flood routing for future events. Hence, parameter estimation is considered to be one of the major challenges of the Muskingum flood routing model.

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.

## METHODS

### The Muskingum model

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*), (3) the speed of the queen at the start of the mating flight (

_{Worker}*S*), (4) the speed of the queen at the end of the mating flight (

_{max}*S*), and (5) the speed reduction factor (

_{min}*χ*) (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.

#### Appraisal of new modifications for the Muskingum model

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.

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(u) is the variable exponent parameter and _{i}u_{i}=rI and _{max}L is a user-defined parameter |

7 | Easa (2013b) | 4 | |

8 | Vatankhah (2014a) | 5 | |

9 | Karahan (2014) | 10 | where βI is lateral flow and I_{ave}= [θI_{t}+ (1-θ)I_{t}_{−1}] and L = 2 (rising and falling inflow limbs) |

10 | Easa et al. (2014) | 6 | |

11 | Easa et al. (2014) | 2L + 2 | where m(u) is the variable exponent parameter and _{i}u_{i}=rI and _{max}L is a user-defined parameter |

12 | Easa (2014) | 4L + 1 | where m(u) is the variable exponent parameter and _{i}u_{i}=rI and _{max}L is a user-defined parameter |

13 | Easa (2014) | 2L + 3 | where x(u) is the variable exponent parameter and _{i}u_{i}=rI and L is a user-defined parameter _{max} |

14 | Karahan et al. (2015) | 5 | where βI is lateral flow and I_{ave}= [θI_{t}+ (1-θ)I_{t}_{−1}] |

15 | Easa (2015) | 2L + 3 | where u_{i}=I and _{i}/I_{max}m(u) _{i}=a+bln(1+cu) or _{i}m(u) _{i}=a+bexp(-cu) and _{i}L is a user-defined parameter |

16 | Haddad et al. (2015b) | 7 | |

17 | Zhang et al. (2016) | 2L + 2 | where m(u) is the variable exponent parameter and _{i}u_{i}=rI and _{max}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(u) is the variable exponent parameter and _{i}u_{i}=rI and _{max}L is a user-defined parameter |

7 | Easa (2013b) | 4 | |

8 | Vatankhah (2014a) | 5 | |

9 | Karahan (2014) | 10 | where βI is lateral flow and I_{ave}= [θI_{t}+ (1-θ)I_{t}_{−1}] and L = 2 (rising and falling inflow limbs) |

10 | Easa et al. (2014) | 6 | |

11 | Easa et al. (2014) | 2L + 2 | where m(u) is the variable exponent parameter and _{i}u_{i}=rI and _{max}L is a user-defined parameter |

12 | Easa (2014) | 4L + 1 | where m(u) is the variable exponent parameter and _{i}u_{i}=rI and _{max}L is a user-defined parameter |

13 | Easa (2014) | 2L + 3 | where x(u) is the variable exponent parameter and _{i}u_{i}=rI and L is a user-defined parameter _{max} |

14 | Karahan et al. (2015) | 5 | where βI is lateral flow and I_{ave}= [θI_{t}+ (1-θ)I_{t}_{−1}] |

15 | Easa (2015) | 2L + 3 | where u_{i}=I and _{i}/I_{max}m(u) _{i}=a+bln(1+cu) or _{i}m(u) _{i}=a+bexp(-cu) and _{i}L is a user-defined parameter |

16 | Haddad et al. (2015b) | 7 | |

17 | Zhang et al. (2016) | 2L + 2 | where m(u) is the variable exponent parameter and _{i}u_{i}=rI and _{max}L is a user-defined parameter |

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).

#### Improving the flood routing scheme

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*in Equation (3). Although the latter modification brings about a new parameter, i.e.,

_{t}*θ*, 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.

### Proposed Muskingum model

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.

#### Simulation procedure of the proposed model

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.

*O*can be expressed as: where

_{t}*I*

_{mod}*=*

*θI*

_{t}*+*(

*1−θ*)

*I*

_{t}_{−1}and

*I*is the modified inflow at time

_{mod}*t*.

Finally, the simulation procedure particularly includes the following steps.

Step 1: Determine the nine hydrologic parameters (*K*, *x _{1}*,

*α*,

_{1}*x*,

_{2}*α*,

_{2}*x*,

_{3}*α*, and

_{3}, m*θ*) 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.

*K*,

*x*,

_{1}*α*,

_{1}*x*,

_{2}*α*,

_{2}*x*,

_{3}*α*, and

_{3}, m*θ*). This check is particularly conducted in order to avoid producing complex values for the outflow using Equation (7). According to this check, the following inequality should be satisfied to obtain acceptable and feasible outflow values.

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 new hybrid technique

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.

Parameters of the new Muskingum model | Minimum value | Maximum value |
---|---|---|

K | 0.001 | 3.3 |

x_{1} | 0.001 | 0.9 |

α_{1} | 0.001 | 2.2 |

x_{2} | 0.001 | 1.5 |

α_{2} | −0.9 | 0.9 |

x_{3} | 0.001 | 1.8 |

α_{3} | 0.001 | 2 |

m | 0.001 | 5 |

Θ | 0 | 1 |

Parameters of the new Muskingum model | Minimum value | Maximum value |
---|---|---|

K | 0.001 | 3.3 |

x_{1} | 0.001 | 0.9 |

α_{1} | 0.001 | 2.2 |

x_{2} | 0.001 | 1.5 |

α_{2} | −0.9 | 0.9 |

x_{3} | 0.001 | 1.8 |

α_{3} | 0.001 | 2 |

m | 0.001 | 5 |

Θ | 0 | 1 |

#### Performance evaluation criteria

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

## RESULTS AND DISCUSSION

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.

### First case study (smooth flood hydrograph)

*et al.*2015a) perform similar and better than the premium solver (Easa 2013b). Additionally, the premium solver (Easa

*et al.*2014) performs the best in calibration of the Muskingum model with five parameters. Although the lowest SSQ value was obtained by AGA (Zhang

*et al.*2016) for

*L*= 5 with twelve parameters, this result was achieved with the highest number of parameters. Furthermore, this model achieved larger SSQ values for four-parameter Muskingum in comparison with other four-parameter models. It should be noted that not only the optimization techniques differ in the compared scenarios for the five-parameter Muskingum models, but also these models have different formations as well. For instance, the CS algorithm (Karahan

*et al.*2015) was applied to the modified version of Muskingum which incorporates lateral flow (the 14th model in Table 1) while premium solver (Easa

*et al.*2014) was utilized for the 11th model in Table 1 for

*L*= 1. Although there is one model with seven-, eight-, nine-, eleven-, and twelve-parameter models in Table 4, these models can be relatively compared with other models that possess a similar number of parameters. For instance, AGA (Zhang

*et al.*2016) scenario for the four-, six-, eight-, ten-, and twelve-parameter models did not perform better than other similar models, except for the ten-parameter model. Therefore, this table not only compares different versions of Muskingum models for the first case study, but also demonstrates the performance of different optimization techniques applied for the related Muskingum models in the literature. According to Table 4, the hybrid MHBMO-GRG obtains the lowest value for the SSQ, SAD, MARE, and VarexQ criteria among all other scenarios while several other techniques performs as well as the recommended hybrid model based on MARE criteria. Additionally, LMM (Das 2004) performs the best in predicting the value of the maximum outflow value since it has the lowest value of DPO among all scenarios in Table 4. According to this table, the new proposed model accompanied by the new hybrid method reduced by 4.52% the SSQ value of Karahan's model (Karahan 2014) which has ten parameters. All in all, it can be concluded that the new recommended model performs the best among all other models in the literature for Wilson's data.

Algorithm | K | x_{1} | α_{1} | x_{2} | α_{2} | x_{3} | α_{3} | m | θ | SSQ |
---|---|---|---|---|---|---|---|---|---|---|

MHBMO | 2.358 | 0.027 | 0.730 | 0.236 | 0.001 | 0.763 | 0.366 | 4.366 | 0.065 | 4.951 |

Hybrid | 2.395 | 0.024 | 0.738 | 0.212 | 0.000 | 0.775 | 0.351 | 4.526 | 0.074 | 4.939 |

Change in values (%) | 1.552 | 10.764 | 1.045 | 9.916 | 100.000 | 1.684 | 4.074 | 3.674 | 13.980 | 0.234 |

Algorithm | K | x_{1} | α_{1} | x_{2} | α_{2} | x_{3} | α_{3} | m | θ | SSQ |
---|---|---|---|---|---|---|---|---|---|---|

MHBMO | 2.358 | 0.027 | 0.730 | 0.236 | 0.001 | 0.763 | 0.366 | 4.366 | 0.065 | 4.951 |

Hybrid | 2.395 | 0.024 | 0.738 | 0.212 | 0.000 | 0.775 | 0.351 | 4.526 | 0.074 | 4.939 |

Change in values (%) | 1.552 | 10.764 | 1.045 | 9.916 | 100.000 | 1.684 | 4.074 | 3.674 | 13.980 | 0.234 |

Method | No. of parameters | SSQ | SAD | DPO | MARE | VarexQ |
---|---|---|---|---|---|---|

GRG (Hirpurkar & Ghare 2014) | 3 | 178.822 | 48.96 | 1.09 | 0.053 | 98.537 |

LSM (Gill 1978) | 3 | 143.6 | 46.4 | 1.8 | 0.056 | 98.825 |

LMM (Das 2004) | 3 | 130.487 | 43.2 | 0 | 0.055 | 98.938 |

HJ + CG (Tung 1985) | 3 | 49.64 | 25.2 | 0.5 | 0.03 | 99.594 |

HJ + DEP (Tung 1985) | 3 | 45.54 | 24.8 | 0.3 | 0.03 | 99.627 |

NONLR (Yoon & Padmanabhan 1993) | 3 | 41.28 | 25.2 | 0.7 | 0.033 | 99.662 |

GA (Mohan 1997) | 3 | 38.23 | 23 | 0.7 | 0.025 | 99.702 |

PSO (Chu & Chang 2009) | 3 | 36.89 | 24.1 | 0.6 | 0.03 | 99.695 |

ICSA (Luo & Xie 2010) | 3 | 36.801 | 23.4 | 0.898 | 0.025 | 99.699 |

HS (Kim et al. 2001) | 3 | 36.78 | 23.4 | 0.92 | 0.031 | 99.627 |

ES (Barati 2013) | 3 | 36.771 | 23.46 | 0.892 | 0.025 | 99.699 |

DE (Xu et al. 2011) | 3 | 36.77 | 23.46 | 0.9 | 0.026 | 99.693 |

PSF-HS (Geem 2010) | 3 | 36.768 | 23.7 | 0.9 | 0.026 | 99.693 |

BFGS (Geem 2006) | 3 | 36.768 | 23.46 | 0.9 | 0.026 | 99.693 |

HS-BFGS (Karahan et al. 2013) | 3 | 36.768 | 23.4 | 0.9 | 0.025 | 99.701 |

GRG (Barati 2013) | 3 | 36.765 | 23.47 | 0.901 | 0.025 | 99.699 |

NMS (Barati 2011) | 3 | 36.765 | 23.46 | 0.899 | 0.025 | 99.699 |

MHBMO (Niazkar & Afzali 2015a) | 3 | 36.242 | 21.88 | 0.699 | 0.028 | 99.703 |

Premium solver (Easa 2013b) | 4 | 7.67 | 10.32 | 0.31 | 0.015 | 99.937 |

GA-GRG (Barati 2012) | 4 | 7.54 | 10.2 | 0.3 | 0.015 | 99.938 |

SFLA-NMS (Haddad et al. 2015a) | 4 | 7.54 | 10.2 | 0.3 | 0.015 | 99.938 |

AGA (Zhang et al. 2016) | 4 | 10.541 | 13.05 | 0.20 | 0.017 | 99.91 |

GRG (Vatankhah 2014a) | 5 | 26.78 | 21.8 | 0.7 | 0.028 | 99.781 |

GRG (Easa 2015) | 5 | 16.57 | 14.9 | 0.3 | 0.02 | 99.864 |

CS (Karahan et al. 2015) | 5 | 9.823 | 12.6 | 0.22 | 0.016 | 99.92 |

GA-GRG (Barati 2012) | 5 | 5.53 | 6.7 | 0.1 | 0.011 | 99.955 |

SFLA-NMS (Haddad et al. 2015a) | 5 | 5.4 | 6.6 | 0.1 | 0.011 | 99.956 |

Premium solver (Easa et al. 2014) | 5 | 5.35 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

Premium solver (Easa et al. 2014) | 6 | 5.44 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

AGA (Zhang et al. 2016) | 6 | 7.475 | 10.58 | 0.31 | 0.014 | 99.94 |

SFLA-NMS (Haddad et al. 2015b) | 7 | 6.348^{b} | 7.61^{b} | 0.05 | 0.013 | 99.948 |

AGA (Zhang et al. 2016) | 8 | 5.730 | 8.58 | 0.10 | 0.012 | 99.95 |

MHBMO-GRG (this study) | 9 | 4.944 | 6.36 | 0.128 | 0.011 | 99.96 |

BFGS (Karahan 2014) | 10 | 5.178 | 9.53 | 0.35 | 0.011 | 99.958 |

AGA (Zhang et al. 2016) | 10 | 4.911 | 7.79 | 0.03 | 0.012 | 99.96 |

Premium solver (Easa 2013a) | 11 | 24.904 | 20.72 | 0.57 | 0.025 | 99.796 |

AGA (Zhang et al. 2016) | 12 | 4.535 | 7.35 | 0.26 | 0.010 | 99.96 |

Method | No. of parameters | SSQ | SAD | DPO | MARE | VarexQ |
---|---|---|---|---|---|---|

GRG (Hirpurkar & Ghare 2014) | 3 | 178.822 | 48.96 | 1.09 | 0.053 | 98.537 |

LSM (Gill 1978) | 3 | 143.6 | 46.4 | 1.8 | 0.056 | 98.825 |

LMM (Das 2004) | 3 | 130.487 | 43.2 | 0 | 0.055 | 98.938 |

HJ + CG (Tung 1985) | 3 | 49.64 | 25.2 | 0.5 | 0.03 | 99.594 |

HJ + DEP (Tung 1985) | 3 | 45.54 | 24.8 | 0.3 | 0.03 | 99.627 |

NONLR (Yoon & Padmanabhan 1993) | 3 | 41.28 | 25.2 | 0.7 | 0.033 | 99.662 |

GA (Mohan 1997) | 3 | 38.23 | 23 | 0.7 | 0.025 | 99.702 |

PSO (Chu & Chang 2009) | 3 | 36.89 | 24.1 | 0.6 | 0.03 | 99.695 |

ICSA (Luo & Xie 2010) | 3 | 36.801 | 23.4 | 0.898 | 0.025 | 99.699 |

HS (Kim et al. 2001) | 3 | 36.78 | 23.4 | 0.92 | 0.031 | 99.627 |

ES (Barati 2013) | 3 | 36.771 | 23.46 | 0.892 | 0.025 | 99.699 |

DE (Xu et al. 2011) | 3 | 36.77 | 23.46 | 0.9 | 0.026 | 99.693 |

PSF-HS (Geem 2010) | 3 | 36.768 | 23.7 | 0.9 | 0.026 | 99.693 |

BFGS (Geem 2006) | 3 | 36.768 | 23.46 | 0.9 | 0.026 | 99.693 |

HS-BFGS (Karahan et al. 2013) | 3 | 36.768 | 23.4 | 0.9 | 0.025 | 99.701 |

GRG (Barati 2013) | 3 | 36.765 | 23.47 | 0.901 | 0.025 | 99.699 |

NMS (Barati 2011) | 3 | 36.765 | 23.46 | 0.899 | 0.025 | 99.699 |

MHBMO (Niazkar & Afzali 2015a) | 3 | 36.242 | 21.88 | 0.699 | 0.028 | 99.703 |

Premium solver (Easa 2013b) | 4 | 7.67 | 10.32 | 0.31 | 0.015 | 99.937 |

GA-GRG (Barati 2012) | 4 | 7.54 | 10.2 | 0.3 | 0.015 | 99.938 |

SFLA-NMS (Haddad et al. 2015a) | 4 | 7.54 | 10.2 | 0.3 | 0.015 | 99.938 |

AGA (Zhang et al. 2016) | 4 | 10.541 | 13.05 | 0.20 | 0.017 | 99.91 |

GRG (Vatankhah 2014a) | 5 | 26.78 | 21.8 | 0.7 | 0.028 | 99.781 |

GRG (Easa 2015) | 5 | 16.57 | 14.9 | 0.3 | 0.02 | 99.864 |

CS (Karahan et al. 2015) | 5 | 9.823 | 12.6 | 0.22 | 0.016 | 99.92 |

GA-GRG (Barati 2012) | 5 | 5.53 | 6.7 | 0.1 | 0.011 | 99.955 |

SFLA-NMS (Haddad et al. 2015a) | 5 | 5.4 | 6.6 | 0.1 | 0.011 | 99.956 |

Premium solver (Easa et al. 2014) | 5 | 5.35 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

Premium solver (Easa et al. 2014) | 6 | 5.44 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

AGA (Zhang et al. 2016) | 6 | 7.475 | 10.58 | 0.31 | 0.014 | 99.94 |

SFLA-NMS (Haddad et al. 2015b) | 7 | 6.348^{b} | 7.61^{b} | 0.05 | 0.013 | 99.948 |

AGA (Zhang et al. 2016) | 8 | 5.730 | 8.58 | 0.10 | 0.012 | 99.95 |

MHBMO-GRG (this study) | 9 | 4.944 | 6.36 | 0.128 | 0.011 | 99.96 |

BFGS (Karahan 2014) | 10 | 5.178 | 9.53 | 0.35 | 0.011 | 99.958 |

AGA (Zhang et al. 2016) | 10 | 4.911 | 7.79 | 0.03 | 0.012 | 99.96 |

Premium solver (Easa 2013a) | 11 | 24.904 | 20.72 | 0.57 | 0.025 | 99.796 |

AGA (Zhang et al. 2016) | 12 | 4.535 | 7.35 | 0.26 | 0.010 | 99.96 |

^{a}The routed outflow was not mentioned in the corresponding paper.

^{b}The SSQ and SAD values are reported equal to 5.44 and 6.69 in the original paper while the corresponding values in this table were computed based on the reported outflow values in the original paper.

### Second case study (multi-peak flood hydrograph)

*x*

_{1},

*α*

_{1},

*x*

_{2},

*α*

_{2}, and

*x*

_{3}are considerably changed in the second phase of the hybrid algorithm. Consequently, the routed outflow may vary because of the differences provided in the optimum values of the Muskingum model. Finally, it can be concluded that the proposed hybrid performs better than the MHBMO algorithm for both case studies considered (Tables 3 and 5). The inflow, observed, and routed outflow hydrographs using the hybrid MHBMO-GRG technique for these data are shown in Figure 4. The closeness of observed and routed outflow hydrographs in this figure obviously implies the accuracy of the utilized model.

Algorithm | K | x_{1} | α_{1} | x_{2} | α_{2} | x_{3} | α_{3} | m | θ | SSQ |
---|---|---|---|---|---|---|---|---|---|---|

MHBMO | 3.300 | 0.253 | 0.625 | 1.900 | 0.374 | 0.629 | 0.994 | 1.007 | 0.000 | 22.274 |

Hybrid | 3.501 | 0.313 | 0.410 | 1.558 | 0.582 | 0.594 | 0.994 | 1.007 | 0.000 | 22.130 |

Change in values (%) | 6.093 | 23.860 | 34.391 | 17.982 | 55.862 | 5.595 | 0.039 | 0.053 | – | 0.645 |

Algorithm | K | x_{1} | α_{1} | x_{2} | α_{2} | x_{3} | α_{3} | m | θ | SSQ |
---|---|---|---|---|---|---|---|---|---|---|

MHBMO | 3.300 | 0.253 | 0.625 | 1.900 | 0.374 | 0.629 | 0.994 | 1.007 | 0.000 | 22.274 |

Hybrid | 3.501 | 0.313 | 0.410 | 1.558 | 0.582 | 0.594 | 0.994 | 1.007 | 0.000 | 22.130 |

Change in values (%) | 6.093 | 23.860 | 34.391 | 17.982 | 55.862 | 5.595 | 0.039 | 0.053 | – | 0.645 |

*et al.*2015a), GRG (Easa 2015), GRG (Easa 2014), and GRG (Easa 2014) algorithms perform the best based on the SSQ value in four-, five-, six-, and seven-parameter Muskingum models in Table 7. Among all these scenarios considered in this table, the hybrid MHBMO-GRG with the new proposed version of Muskingum model performs the best since it significantly reduced the best SSQ available in the literature, i.e., the MHBMO algorithm (Niazkar & Afzali 2015a). Since most previous studies did not mention the routed outflow values, the comparison of scenarios in Table 7 based on other criteria may not be applicable, although the hybrid MHBMO-GRG obtains the best values for other criteria in Table 7. Finally, the application of the new proposed Muskingum model for the second case study clearly shows the superiority of this model over others for double-peak flood hydrographs.

t (hr) | I (m^{3}/s) | O (m^{3}/s) | ΔS/Δt | S | I_{mod} | O_{routed} |
---|---|---|---|---|---|---|

0 | 121 | 116 | – | – | – | 120.63 |

1 | 217 | 121 | 174.64 | 466.55 | 120.63 | 120.63 |

2 | 316 | 125 | 272.68 | 641.19 | 216.51 | 124.91 |

3 | 474 | 173 | 428.80 | 913.88 | 316.21 | 173.06 |

4 | 611 | 249 | 474.83 | 1,342.67 | 473.74 | 248.43 |

5 | 593 | 362 | 216.39 | 1,817.50 | 611.36 | 362.08 |

6 | 752 | 480 | 402.21 | 2,033.89 | 593.24 | 479.82 |

7 | 1,303 | 541 | 1,206.63 | 2,436.10 | 752.38 | 541.42 |

8 | 1,698 | 668 | 1,350.05 | 3,642.73 | 1,302.58 | 667.66 |

9 | 1,635 | 988 | 595.80 | 4,992.78 | 1,697.88 | 988.21 |

10 | 1,356 | 1,323 | −193.30 | 5,588.58 | 1,635.02 | 1,322.39 |

11 | 976 | 1,457 | −792.09 | 5,395.27 | 1,356.09 | 1,457.42 |

12 | 613 | 1,391 | −1,074.30 | 4,603.19 | 975.80 | 1,390.66 |

13 | 982 | 1,176 | 107.83 | 3,528.89 | 613.34 | 1,175.57 |

14 | 1,279 | 925 | 595.70 | 3,636.73 | 982.03 | 925.51 |

15 | 1,391 | 966 | 515.84 | 4,232.43 | 1,279.36 | 966.78 |

16 | 1,169 | 1,120 | −132.67 | 4,748.27 | 1,391.49 | 1,120.86 |

17 | 958 | 1,238 | −453.19 | 4,615.60 | 1,169.20 | 1,238.74 |

18 | 581 | 1,195 | −923.50 | 4,162.41 | 957.96 | 1,195.36 |

19 | 417 | 1,064 | −782.03 | 3,238.90 | 580.78 | 1,064.08 |

20 | 324 | 826 | −578.67 | 2,456.88 | 416.82 | 825.91 |

21 | 263 | 627 | −413.30 | 1,878.20 | 323.83 | 626.46 |

22 | 222 | 479 | −291.73 | 1,464.90 | 263.18 | 479.28 |

23 | 176 | 374 | −235.43 | 1,173.17 | 221.75 | 374.26 |

24 | 172 | 299 | −130.55 | 937.74 | 176.36 | 299.41 |

t (hr) | I (m^{3}/s) | O (m^{3}/s) | ΔS/Δt | S | I_{mod} | O_{routed} |
---|---|---|---|---|---|---|

0 | 121 | 116 | – | – | – | 120.63 |

1 | 217 | 121 | 174.64 | 466.55 | 120.63 | 120.63 |

2 | 316 | 125 | 272.68 | 641.19 | 216.51 | 124.91 |

3 | 474 | 173 | 428.80 | 913.88 | 316.21 | 173.06 |

4 | 611 | 249 | 474.83 | 1,342.67 | 473.74 | 248.43 |

5 | 593 | 362 | 216.39 | 1,817.50 | 611.36 | 362.08 |

6 | 752 | 480 | 402.21 | 2,033.89 | 593.24 | 479.82 |

7 | 1,303 | 541 | 1,206.63 | 2,436.10 | 752.38 | 541.42 |

8 | 1,698 | 668 | 1,350.05 | 3,642.73 | 1,302.58 | 667.66 |

9 | 1,635 | 988 | 595.80 | 4,992.78 | 1,697.88 | 988.21 |

10 | 1,356 | 1,323 | −193.30 | 5,588.58 | 1,635.02 | 1,322.39 |

11 | 976 | 1,457 | −792.09 | 5,395.27 | 1,356.09 | 1,457.42 |

12 | 613 | 1,391 | −1,074.30 | 4,603.19 | 975.80 | 1,390.66 |

13 | 982 | 1,176 | 107.83 | 3,528.89 | 613.34 | 1,175.57 |

14 | 1,279 | 925 | 595.70 | 3,636.73 | 982.03 | 925.51 |

15 | 1,391 | 966 | 515.84 | 4,232.43 | 1,279.36 | 966.78 |

16 | 1,169 | 1,120 | −132.67 | 4,748.27 | 1,391.49 | 1,120.86 |

17 | 958 | 1,238 | −453.19 | 4,615.60 | 1,169.20 | 1,238.74 |

18 | 581 | 1,195 | −923.50 | 4,162.41 | 957.96 | 1,195.36 |

19 | 417 | 1,064 | −782.03 | 3,238.90 | 580.78 | 1,064.08 |

20 | 324 | 826 | −578.67 | 2,456.88 | 416.82 | 825.91 |

21 | 263 | 627 | −413.30 | 1,878.20 | 323.83 | 626.46 |

22 | 222 | 479 | −291.73 | 1,464.90 | 263.18 | 479.28 |

23 | 176 | 374 | −235.43 | 1,173.17 | 221.75 | 374.26 |

24 | 172 | 299 | −130.55 | 937.74 | 176.36 | 299.41 |

Method | No. of parameters | SSQ | SAD | DPO | MARE | VarexQ |
---|---|---|---|---|---|---|

Premium solver (Easa 2013a) | 3 | 71,114 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

MHBMO (Niazkar & Afzali 2015a) | 3 | 7,751.662 | 304.675 | 8.717 | 0.037 | 99.837 |

GA-GRG (Barati 2012) | 4 | 76,785 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

SFLA-NMS (Haddad et al. 2015a) | 4 | 73,379 | 1,034 | NM^{a} | NM^{a} | NM^{a} |

GA-GRG (Barati 2012) | 5 | 70,738 | 1,018 | NM^{a} | NM^{a} | NM^{a} |

SFLA-NMS (Haddad et al. 2015a) | 5 | 69,860 | 993 | NM^{a} | NM^{a} | NM^{a} |

Premium solver (Easa et al. 2014) | 5 | 69,726.75 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

GRG (Easa 2015) | 5 | 57,391 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

Premium solver (Easa et al. 2014) | 6 | 68,015 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

GRG (Easa 2014) | 6 | 52,748 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

SFLA-NMS (Haddad et al. 2015b) | 7 | 73,399 | 1,034 | 50 | NM^{a} | NM^{a} |

GRG (Easa 2014) | 7 | 38,823 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

MHBMO-GRG (this study) | 9 | 22.13 | 9.293 | 0.03 | 0.002 | 100 |

BFGS (Karahan 2014) | 10 | 1,277,047.73^{b} | 4,664.23 | 24.23 | 0.37 | 73.16 |

Premium solver (Easa 2013a) | 11 | 44,894 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

Method | No. of parameters | SSQ | SAD | DPO | MARE | VarexQ |
---|---|---|---|---|---|---|

Premium solver (Easa 2013a) | 3 | 71,114 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

MHBMO (Niazkar & Afzali 2015a) | 3 | 7,751.662 | 304.675 | 8.717 | 0.037 | 99.837 |

GA-GRG (Barati 2012) | 4 | 76,785 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

SFLA-NMS (Haddad et al. 2015a) | 4 | 73,379 | 1,034 | NM^{a} | NM^{a} | NM^{a} |

GA-GRG (Barati 2012) | 5 | 70,738 | 1,018 | NM^{a} | NM^{a} | NM^{a} |

SFLA-NMS (Haddad et al. 2015a) | 5 | 69,860 | 993 | NM^{a} | NM^{a} | NM^{a} |

Premium solver (Easa et al. 2014) | 5 | 69,726.75 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

GRG (Easa 2015) | 5 | 57,391 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

Premium solver (Easa et al. 2014) | 6 | 68,015 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

GRG (Easa 2014) | 6 | 52,748 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

SFLA-NMS (Haddad et al. 2015b) | 7 | 73,399 | 1,034 | 50 | NM^{a} | NM^{a} |

GRG (Easa 2014) | 7 | 38,823 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

MHBMO-GRG (this study) | 9 | 22.13 | 9.293 | 0.03 | 0.002 | 100 |

BFGS (Karahan 2014) | 10 | 1,277,047.73^{b} | 4,664.23 | 24.23 | 0.37 | 73.16 |

Premium solver (Easa 2013a) | 11 | 44,894 | NM^{a} | NM^{a} | NM^{a} | NM^{a} |

^{a}The routed outflow was not mentioned in the corresponding paper.

^{b}The SSQ value is reported equal to 26,185 in the original paper while the corresponding value in this table was computed based on the reported outflow values in the original paper.

### Sensitivity analysis

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_{1} and α_{1} 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 *α*_{2} and more sensitive to *α*_{3}. 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.

## CONCLUSION

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.