First, a novel nonlinear Muskingum flood routing model with a variable exponent parameter and simultaneously considering the lateral flow along the river reach (named VEP-NLMM-L) was developed in this research. Then, an improved real-coded adaptive genetic algorithm (RAGA) with elite strategy was applied for precise parameter estimation of the proposed model. The problem was formulated as a mathematical optimization procedure to minimize the sum of the squared deviations (SSQ) between the observed and the estimated outflows. Finally, the VEP-NLMM-L was validated on three watersheds with different characteristics (Case 1 to 3). Comparisons of the optimal results for the three case studies by traditional Muskingum models and the VEP-NLMM-L show that the modified Muskingum model can produce the most accurate fit to outflow data. Application results in Case 3 also indicate that the VEP-NLMM-L may be suitable for solving river flood routing problems in both model calibration and prediction stages.

## INTRODUCTION

*K*and a weighting parameter

*w*are used: where is the absolute channel storage at time

*t*; and are observed rates of inflow and outflow at time

*t*, respectively;

*K*is the storage-time constant, a value usually close to the flow travel time through the routing river reach;

*w*is the dimensionless weighting factor that represents the inflow–outflow relative effects on the storage, for stream channels and for reservoir storage (Mohan 1997).

*m*has been introduced to account for the effects of nonlinearity in the following two forms of three-parameter nonlinear Muskingum model (named NLMM1 and NLMM2, respectively) (Gill 1978; Singh & Scarlatos 1987): Although the linear and nonlinear Muskingum models (Equations (1)–(4)) have been widely discussed (Tung 1985; Yoon & Padmanabhan 1993; Das 2004; Chen & Yang 2007; Geem 2011; Xu

*et al.*2012; Niazkar & Afzali 2014), they all neglect the lateral flow along the investigated river reach. However, the lateral flow really exists in many river reaches in actual flood events. Based on the assumptions that the lateral flow varied linearly along the river reach and could be expressed as a ratio of the inflow rate , O'Donnell (1985) first proposed a three-parameter linear Muskingum model which considers the lateral flow (LMM-L) as Equations (5) and (6): Apart from the necessity of considering the lateral flow, the nonlinear storage–discharge relationship and its variation (Easa 2013) are also important in Muskingum model building. Thus, a modified nonlinear Muskingum model, which assumes the exponent parameter in the NLMM2 to be varying with the inflow dividing levels and simultaneously takes the lateral flow into consideration, is presented in this paper. In addition, the precise estimation of parameters is the key point in applying the Muskingum model for real-time flood forecasting. During the past decades, various techniques have been applied to solve this problem and can be classified into two types: mathematical methods and heuristic optimization algorithms. The mathematical methods include the least-squares method (LSM) (Gill 1978), the Hooke–Jeeves (HJ) pattern search in conjunction with linear regression (HJ + LR), the conjugate gradient (HJ + CG) or Davidon–Fletcher–Powell (HJ + DFP) methods (Tung 1985), nonlinear least-squares regression (NONLR) (Yoon & Padmanabhan 1993), the nonlinear programming (NLP) algorithm (Kshirsagar

*et al.*1995), the Broyden–Fletcher–Goldfarb–Shanno (BFGS) technique (Geem 2006), the Nelder–Mead simplex (NMS) algorithm (Barati 2011), and the generalized reduced gradient (GRG) method (Hirpurkar & Ghare 2014; Easa 2015). The above-mentioned methods inevitably have the drawbacks of some special derivation conditions, a suitable guess of initial solution or being time-consuming, and they cannot guarantee convergence to the global optima in most cases (Orouji

*et al.*2013; Bozorg Haddad

*et al.*2015). Therefore, an increasing number of researchers have focused in recent years on heuristic algorithms that exhibit strong global search ability and are easy to design, such as genetic algorithms (GAs) (Mohan 1997; Chen & Yang 2007), harmony search (HS) (Kim

*et al.*2001), particle swarm optimization (PSO) algorithm (Chu & Chang 2009), differential evolution (DE) (Xu

*et al.*2012), the shuffled frog leaping algorithm (SFLA), the simulated annealing (SA) algorithm (Orouji

*et al.*2013), and the honey bee mating optimization (HBMO) (Niazkar & Afzali 2014).

The remainder of this paper is organized as follows. First, a modified nonlinear Muskingum model called VEP-NLMM-L is described in the section below. The detailed structure and standard routing procedure of VEP-NLMM-L are also specified in this section. Second, the implementation steps of applying the improved real-coded adaptive genetic algorithm (RAGA) to estimate parameters of the VEP-NLMM-L are given in the following section. The next section defines the five evaluation criteria for assessing the quality of the calculated outflows. Then, three typical cases are employed to demonstrate the validity of VEP-NLMM-L in both model calibration and prediction stages. Finally, some conclusions about our research work are drawn.

## PROPOSED MODEL

### Model structure

As can be seen in both traditional LMM, NLMM1, NLMM2, and LMM-L, the parameters are assumed to be constant during the entire flood routing periods. It is worth noting that the exponent parameter *m* in NLMM2 has no physical meaning but only represents the average nonlinearity in the whole flood routing procedure. For the given complex, unsteady flow conditions, Easa (2013) considered *m* should be varying with the inflow levels. Hence, after integrating the later flow assumptions by O'Donnell (1985) (Equation (5)) and the variable exponent parameter suggestion by Easa (2013) into the NLMM2 (Equation (4)), a new nonlinear Muskingum model (VEP-NLMM-L) was proposed. The model structure is described below.

*L*levels by (

*L*–1) inflow dividing values , . Define the real variable , where is the maximum inflow. Thus, the lower and upper limits of are 0 and 1, respectively. For inflow level

*i*, , denoting the exponent parameter by . Variation range of are typically between 1 and 3. The constraints for real variable are given by: For a flood hydrograph that goes through time intervals

*T*, the exponent value for flood routing in time interval

*t*, , is determined by: Then, the structure of the VEP-NLMM-L can be written as:

### Routing procedure

The main steps of using the VEP-NLMM-L to route floods are as follows:

**Step 1**: Assume values of the parameters*K*,*w*, , , .**Step 2**: Calculate the storage amount using Equation (9), where the value of is determined by Equation (8) beforehand.where is the rate of change of storage volume at ; is the time interval.

where most previous studies suggested using the inflow at the previous time-point rather than at the current time-point in the calculation of the outflow (Geem 2011).

**Step 6**: Repeat Steps 2–5 for all time steps.

## IMPROVED ADAPTIVE GENETIC ALGORITHM

The genetic algorithm is a population-based, global search approach inspired by the evolutionary ideas of natural selection and survival of the fittest. GAs have been demonstrated to perform well in handling hydrological modeling problems (Wang 1997; Chang & Chen 1998; Kumar *et al.* 2012; Yang *et al.* 2013). A real-coded adaptive genetic algorithm with elite strategy for parameter estimation of the VEP-NLMM-L is described as follows.

### Individual coding and fitness assignment

*et al.*2005). Herein, each individual is coded as a vector of floating point numbers with the same length as the number of parameters in VEP-NLMM-L, i.e., individual . The total number of decision variables is , the real-value coding of variable is expressed as: where and are specified lower and upper bounds of variable ; rand(0,1) denotes a random number between 0 and 1.

*t*; is the calculated outflow at time

*t*;

*T*is the total number of time intervals.

**x**) is, the higher the fitness is. Thus, for each individual

**x**, its fitness function fit(

**x**) is defined as Equation (16):

### Adaptive probabilities of crossover and mutation

The tournament selection operator which integrates the idea of ranking is applied for selection operation. To match with the real-coded pattern, the simulated binary crossover (SBX) operator is used for crossover operation (Deb 2000). The polynomial mutation operator is adopted for maintaining diversity in the population during the evolutionary process (Deb *et al.* 2002).

*et al.*(2013) consider that the adaptive and are not only depending on the fitness value of the individual, but are also related to the dispersion degree of the population. They can be expressed as follows: where is the larger fitness value between the two selected individuals to be crossed; is the average fitness value of the population; is the maximum fitness value of the population; , , , are adaptive parameters.

### Constraint handling

**x**) to form the penalized function

*F*(

**x**), as follows: where

*J*is the number of constraints; the term is zero if and otherwise; is the penalty parameter of the

*j*th constraint; hence, if no constraint violation occurs, individual

**x**is feasible and the overall constraint violation , otherwise,

**x**is infeasible and .

As can be seen in Equation (19), the appropriate setting for values of penalty parameters is the major difficulty in using the penalty function method. This usually requires users to experiment with different values. To overcome this drawback, the following method proposed by Deb (2000) is employed.

The overall constraint violation *v*(**x**) is coupled with the tournament selection operator to implement the penalty function method according to the following three criteria when comparing two individuals or solutions:

Any feasible solution (

*v*(**x**) = 0) is preferred to any infeasible solution (*v*(**x**) > 0).Between two feasible solutions, the one having higher fit(

**x**) is preferred.Between two infeasible solutions, the one having smaller

*v*(**x**) is preferred.

Through this method, individuals are never compared in terms of both objective function value and constraint violation information. Thus, penalty parameters are not needed.

*K*,

*w*, and of some individuals are infeasible. To avoid this, they should be bounded and penalized by Equations (21) and (22): where , are the penalty constants that have very small value (for example, ) to avoid negative and ; and are the next penalized storage and outflow, respectively.

### Elite strategy

The elite strategy is the key for guaranteeing GAs converge to the global optimal solution. The idea is to directly preserve the best individuals in the current population to the next generation without crossover and mutation operations and replace the corresponding number of worst individuals of the new generated offspring generation. Here, the top 5% of individuals are preserved in each generation.

### Implementation of RAGA on VEP-NLMM-L

## EVALUATION CRITERIA FOR CALCULATED OUTFLOWS

To comprehensively assess the accuracy of the calculated (estimated or predicted) outflow hydrographs by different Muskingum models, five relevant evaluation criteria are used in this study. They are: (1) the absolute deviation of peak outflow (DPO); (2) the deviation of peak time (DPOT) (Yoon & Padmanabhan 1993); (3) the Nash–Sutcliffe criterion (Barati 2013); (4) the mean absolute error (MAE); and (5) the mean absolute relative error (MARE) (Niazkar & Afzali 2014). The criteria DPO and the DPOT check the accuracy of magnitude and occurrence time of the calculated outflow peak, respectively. The criterion (%) represents the correctness of the shape and size of the calculated outflow hydrograph, compared with the observed outflow hydrograph. The criteria MAE and MARE measure the degree of how close the calculated outflows are to the observed outflows. These five evaluation criteria are defined as follows.

### Evaluation of magnitude and occurrence time of outflow peak

### Accuracy of procedure consideration

### Mean absolute error or relative error

## CASE STUDIES

To reasonably evaluate the performance of the VEP-NLMM-L, the data set from Wilson (1974), a flood in the River Wyre in October 1982 (O'Donnell 1985), and two different flood events that occurred in the same river reach of the River Wye (Natural Environment Research Council (NERC) 1975) were selected as three illustrative case studies (Cases 1–3). Search spaces for the 2*L* + 2 parameters of the VEP-NLMM-L are specified in Table 1. The parameter represents the contribution of lateral flow to the outflow and is related to the area under the inflow–outflow hydrograph. As is shown in Table 1, the lateral flow is very small in Case 1 and Case 3, but it is large and makes a significant contribution to the outflow in Case 2.

. | K
. | α
. | w
. | m (_{i}i = 1,…, L)
. | r (_{i}i = 1,…, L-1)
. |
---|---|---|---|---|---|

Case 1 | [0.01, 1.00] | [−0.10, 0.10] | [0.01, 0.50] | (1.00, 3.00) | (0.00, 1.00) |

Case 2 | [0.01, 6.00] | [−3.00, 3.00] | |||

Case 3 | [0.01, 1.00] | [−0.10, 0.10] |

. | K
. | α
. | w
. | m (_{i}i = 1,…, L)
. | r (_{i}i = 1,…, L-1)
. |
---|---|---|---|---|---|

Case 1 | [0.01, 1.00] | [−0.10, 0.10] | [0.01, 0.50] | (1.00, 3.00) | (0.00, 1.00) |

Case 2 | [0.01, 6.00] | [−3.00, 3.00] | |||

Case 3 | [0.01, 1.00] | [−0.10, 0.10] |

In all cases, the VEP-NLMM-L was solved over 50 runs using the RAGA to search for the optimal solution and the parameter settings of the RAGA were as follows: population size pop = 200; distribution indexes for the SBX and the polynomial mutation are 1.0 and 100.0, respectively; adaptive controlling parameters , , , ; the maximum number of iterations MaxIter is 1,000.

### Case 1: Data set of Wilson (1974)

The data set from Wilson (1974), which has been demonstrated to have an obvious nonlinear relationship between the storage and the weighted-flow (Yoon & Padmanabhan 1993; Mohan 1997), was selected as the first case. These flood data are a single peak hydrograph and have been previously studied by many researchers (O'Donnell 1985; Yoon & Padmanabhan 1993; Al-Humoud & Esen 2006; Niazkar & Afzali 2014), where and *T* = 21.

For this case, optimal solutions for the VEP-NLMM-L with different numbers of inflow dividing levels (*L* = 1–5) are listed in Table 2, and corresponding results by the LMM-L (O'Donnell 1985) and the NLMM2 (Niazkar & Afzali 2014) are also given for comparisons. As is shown in Table 2, the SSQ values by LMM-L and NLMM2 are 815.680 and 36.242, respectively. The value 36.242 is the best existing value for the NLMM2 so far in the literature. The SSQs by estimating the VEP-NLMM-L of different *L* have obvious reductions, at least (*L* = 1) 98.71% and 70.91% lower than those obtained by the LMM-L and NLMM2, respectively. Comparing the VEP-NLMM-L with the classical NLMM2, the variation range of the exponent parameter is small and the optimal values of other parameters are close, but the improvement on accuracy of flood routing is excellent. For example, the optimal inflow dividing values are , , , when *L* = 5.

Model . | K
. | w
. | α
. | m/m (_{i}i = 1,…,L)
. | SSQ . | |
---|---|---|---|---|---|---|

LMM-L (O'Donnell 1985) | 5.4300 | 0.1393 | 0.0235 | – | 815.680 | |

NLMM2 (Niazkar & Afzali 2014) | 0.6589 | 0.3399 | – | 1.8456 | 36.242 | |

VEP-NLMM-L (this study) | L = 1 | 0.5362 | 0.3005 | −0.0215 | 1.8634 | 10.541 |

L = 2 | 0.4921 | 0.2978 | −0.0195 | 1.8996, 1.8833 | 7.475 | |

L = 3 | 0.5875 | 0.3048 | −0.0199 | 1.8583, 1.8446, 1.8361 | 5.730 | |

L = 4 | 0.6176 | 0.3115 | −0.0199 | 1.8474, 1.8344, 1.8299, 1.8231 | 4.911 | |

L = 5 | 0.5565 | 0.2938 | −0.0192 | 1.8889, 1.8697, 1.8540, 1.8573, 1.8505 | 4.535 |

Model . | K
. | w
. | α
. | m/m (_{i}i = 1,…,L)
. | SSQ . | |
---|---|---|---|---|---|---|

LMM-L (O'Donnell 1985) | 5.4300 | 0.1393 | 0.0235 | – | 815.680 | |

NLMM2 (Niazkar & Afzali 2014) | 0.6589 | 0.3399 | – | 1.8456 | 36.242 | |

VEP-NLMM-L (this study) | L = 1 | 0.5362 | 0.3005 | −0.0215 | 1.8634 | 10.541 |

L = 2 | 0.4921 | 0.2978 | −0.0195 | 1.8996, 1.8833 | 7.475 | |

L = 3 | 0.5875 | 0.3048 | −0.0199 | 1.8583, 1.8446, 1.8361 | 5.730 | |

L = 4 | 0.6176 | 0.3115 | −0.0199 | 1.8474, 1.8344, 1.8299, 1.8231 | 4.911 | |

L = 5 | 0.5565 | 0.2938 | −0.0192 | 1.8889, 1.8697, 1.8540, 1.8573, 1.8505 | 4.535 |

The five prescribed evaluation criteria in the section ‘Evaluation criteria for calculated outflows' are calculated according to the estimated outflows and shown in Table 3. In general, the VEP-NLMM-L produced the most accurate routing outflows in terms of all criteria and the forecasting accuracy increased with the number of inflow dividing levels. All the models except the LMM-L estimate the peak outflow on the correct time interval (DPOT = 0). The VEP-NLMM-L with *L* = 4 estimates the closest peak outflow to the observed value among all models.

Model . | DPO . | DPOT . | η
. | MARE . | MAE . | |
---|---|---|---|---|---|---|

LMM-L (O'Donnell 1985) | 5.80 | −1 | 93.33 | 0.137 | 4.918 | |

NLMM2 (Niazkar & Afzali 2014) | 0.70 | 0 | 99.70 | 0.028 | 0.994 | |

VEP-NLMM-L (this study) | L = 1 | 0.20 | 0 | 99.91 | 0.017 | 0.593 |

L = 2 | 0.31 | 0 | 99.94 | 0.014 | 0.481 | |

L = 3 | 0.10 | 0 | 99.95 | 0.012 | 0.390 | |

L = 4 | 0.03 | 0 | 99.96 | 0.012 | 0.354 | |

L = 5 | 0.26 | 0 | 99.96 | 0.010 | 0.334 |

Model . | DPO . | DPOT . | η
. | MARE . | MAE . | |
---|---|---|---|---|---|---|

LMM-L (O'Donnell 1985) | 5.80 | −1 | 93.33 | 0.137 | 4.918 | |

NLMM2 (Niazkar & Afzali 2014) | 0.70 | 0 | 99.70 | 0.028 | 0.994 | |

VEP-NLMM-L (this study) | L = 1 | 0.20 | 0 | 99.91 | 0.017 | 0.593 |

L = 2 | 0.31 | 0 | 99.94 | 0.014 | 0.481 | |

L = 3 | 0.10 | 0 | 99.95 | 0.012 | 0.390 | |

L = 4 | 0.03 | 0 | 99.96 | 0.012 | 0.354 | |

L = 5 | 0.26 | 0 | 99.96 | 0.010 | 0.334 |

*L*= 4 for Wilson's data is shown in Figure 2.

### Case 2: River Wyre October 1982 flood

This case study employed the River Wyre flood event in October 1982 (O'Donnell 1985), which exhibits a considerable increase in the flood volume (lateral flow) between the inflow and outflow sections (some 25 km apart). The flood data have a multi-peaked inflow, and *T* = 31.

Table 4 lists the optimal solution vectors for the VEP-NLMM-L of different inflow dividing levels (*L* = 1–5). The SSQs calculated by the VEP-NLMM-L of different *L* have obvious reductions, at least 88.23% (*L* = 1) lower than the value obtained by the LMM-L. For VEP-NLMM-L with *L* = 5, the optimal inflow dividing values are , , , .

Model . | K
. | w
. | α
. | m/m (_{i}i = 1,…,L)
. | SSQ . | |
---|---|---|---|---|---|---|

LMM-L (O'Donnell 1985) | 5.5400 | 0.1370 | 2.6200 | – | 468.840 | |

VEP-NLMM-L (this study) | L = 1 | 5.1485 | 0.2282 | 2.5285 | 1.0000 | 55.193 |

L = 2 | 4.9473 | 0.2524 | 2.5212 | 1.0120, 1.0059 | 28.462 | |

L = 3 | 4.7858 | 0.2518 | 2.5225 | 1.0188, 1.0102, 1.0134 | 26.239 | |

L = 4 | 4.2020 | 0.2730 | 2.5142 | 1.0135, 1.0498, 1.0458, 1.0398 | 20.490 | |

L = 5 | 3.7303 | 0.2856 | 2.5166 | 1.0477, 1.0841, 1.0772, 1.0722, 1.0658 | 17.574 |

Model . | K
. | w
. | α
. | m/m (_{i}i = 1,…,L)
. | SSQ . | |
---|---|---|---|---|---|---|

LMM-L (O'Donnell 1985) | 5.5400 | 0.1370 | 2.6200 | – | 468.840 | |

VEP-NLMM-L (this study) | L = 1 | 5.1485 | 0.2282 | 2.5285 | 1.0000 | 55.193 |

L = 2 | 4.9473 | 0.2524 | 2.5212 | 1.0120, 1.0059 | 28.462 | |

L = 3 | 4.7858 | 0.2518 | 2.5225 | 1.0188, 1.0102, 1.0134 | 26.239 | |

L = 4 | 4.2020 | 0.2730 | 2.5142 | 1.0135, 1.0498, 1.0458, 1.0398 | 20.490 | |

L = 5 | 3.7303 | 0.2856 | 2.5166 | 1.0477, 1.0841, 1.0772, 1.0722, 1.0658 | 17.574 |

Table 5 shows the evaluation criteria of the estimated outflow hydrographs by the LMM-L and VEP-NLMM-L. Similarly to Case 1, the VEP-NLMM-L still produces more accurate routing outflows in terms of all criteria than the LMM-L and the accuracy increases with the number of inflow dividing levels *L*. All models can estimate the peak outflow on the correct time interval (DPOT = 0). The VEP-NLMM-L with *L* = 4 also gives the closest peak outflow to the observed value.

Model . | DPO . | DPOT . | η
. | MARE . | MAE . | |
---|---|---|---|---|---|---|

LMM-L (O'Donnell 1985) | 1.10 | 0 | 98.39 | 0.101 | 3.231 | |

VEP-NLMM-L (this study) | L = 1 | 1.88 | 0 | 99.81 | 0.030 | 0.950 |

L = 2 | 0.29 | 0 | 99.90 | 0.024 | 0.665 | |

L = 3 | 0.17 | 0 | 99.91 | 0.025 | 0.689 | |

L = 4 | 0.07 | 0 | 99.93 | 0.023 | 0.662 | |

L = 5 | 0.16 | 0 | 99.94 | 0.021 | 0.624 |

Model . | DPO . | DPOT . | η
. | MARE . | MAE . | |
---|---|---|---|---|---|---|

LMM-L (O'Donnell 1985) | 1.10 | 0 | 98.39 | 0.101 | 3.231 | |

VEP-NLMM-L (this study) | L = 1 | 1.88 | 0 | 99.81 | 0.030 | 0.950 |

L = 2 | 0.29 | 0 | 99.90 | 0.024 | 0.665 | |

L = 3 | 0.17 | 0 | 99.91 | 0.025 | 0.689 | |

L = 4 | 0.07 | 0 | 99.93 | 0.023 | 0.662 | |

L = 5 | 0.16 | 0 | 99.94 | 0.021 | 0.624 |

*L*= 4 for the River Wyre flood, October 1982, is shown in Figure 3.

### Case 3: River Wye December 1960 and January 1969 floods

*L*= 4) for Scenario I were estimated. Then, the obtained parameters were applied to predict the outflow hydrograph for Scenario II. The optimal parameter vectors by optimizing the LMM and VEP-NLMM-L (

*L*= 4) for Scenario I are listed in Table 6. The corresponding estimated outflows of Scenario I are shown in Figure 4.

Model . | K
. | w
. | α
. | m/m (_{i}i = 1,…,L)
. |
---|---|---|---|---|

LMM-L | 5.0500 | 0.1730 | 0.0720 | – |

VEP-NLMM-L (L = 4) | 0.0981 | 0.3168 | 0.0617 | 1.8190, 1.7742, 1.8131, 1.8169 |

Model . | K
. | w
. | α
. | m/m (_{i}i = 1,…,L)
. |
---|---|---|---|---|

LMM-L | 5.0500 | 0.1730 | 0.0720 | – |

VEP-NLMM-L (L = 4) | 0.0981 | 0.3168 | 0.0617 | 1.8190, 1.7742, 1.8131, 1.8169 |

*L*= 4) is much lower than the value obtained by LMM-L (SSQ = 26,113). The peak of predicted outflows by the proposed model could be controlled with a 2.3% error (DPO = 7.46 and ) and the peak occurrence time was correct (DPOT = 0), and also performed much better than the LMM-L.

## CONCLUSION

This paper proposed a new multi-parameter nonlinear Muskingum model (VEP-NLMM-L). The VEP-NLMM-L includes a variable exponent parameter accounting for the nonlinearity of flood wave and a flow parameter that represents the lateral flow along the investigated river reach. The variable exponent parameter varies depending on the inflow level. An improved real-coded adaptive genetic algorithm (RAGA) with elite strategy was designed and applied for the parameter calibration of the VEP-NLMM-L. From the comparative results of different Muskingum models in three typical case studies, several conclusions can be drawn as below:

The VEP-NLMM-L was successfully modeled and applied in three watersheds (Cases 1–3) with different characteristics. The lateral flow contribution in both Case 1 and 3, i.e., the data set of Wilson (1974) and the two floods in the River Wye, are minor. However, an obvious nonlinear storage–discharge relationship exists in Case 1. The River Wyre, October 1982 flood (Case 3) has a major lateral flow contribution.

For Case 1 and Case 2, optimal results by estimating the VEP-NLMM-L indicate that the SSQ values have obvious reductions, at least 70.91% and 88.23%, compared to the best SSQs by traditional models (NLMM2 and LMM) reported in the literature, respectively. Although the improvement on accuracy of flood routing is excellent, the variation range of the optimal variable exponent parameter is small and the optimal values of all parameters in the VEP-NLMM-L are not substantially different from the NLMM2 constant exponent parameter. Comparisons of the five evaluation criteria based on the estimated outflows also show that the VEP-NLMM-L can produce the most accurate flood routing hydrograph than the other methods in the model calibration stage and its accuracy will increase with the number of inflow dividing levels.

Furthermore, in Case 3, the River Wye, December 1960 flood was used for estimating model parameters (model calibration stage) and the River Wye, January 1969 flood for model forecasting validation (model prediction stage); the results show that the proposed VEP-NLMM-L has excellent applicability and reliability in solving flood routing problems.

## ACKNOWLEDGEMENTS

This study was financially supported by the Hubei Support Plan of Science and Technology (No. 2015BCA291), the Wuhan Planning Project of Science and Technology (No. 2014060101010062) and the Natural Science Foundation of China (Grant No. 51509099).