Flood Routing in River Reaches Using a Three-Parameter Muskingum Model Coupled with an Improved Bat Algorithm

Design of hydraulic structures, flood warning systems, evacuation measures, and traffic management require river flood routing. A common hydrologic method of flood routing is the Muskingum method. The present study attempted to develop a three-parameter Muskingum model considering lateral flow for flood routing, coupling with a new optimization algorithm namely, Improved Bat Algorithm (IBA). The major function of the IBA is to optimize the estimated value of the three-parameters associated with the Muskingum model. The IBA acts based on the chaos search tool, which mainly enhances the uniformity and erogidicty of the population. In addition, the current research, unlike the other existing models which consider flood routing, is based on dividing one reach to a few intervals to increase the accuracy of flood routing models. Three case studies with lateral flow were considered for this study, including the Wilson flood, Karahan flood, and Myanmar flood. Seven performance indexes were examined to evaluate the performance of the proposed Muskingum model integrated with IBA, with other models that were also based on the Muskingum Model with three-parameters but utilized different optimization algorithms. The results for the Wilson flood showed that the proposed model could reduce the Sum of Squared Deviations (SSD) value by 89%, 51%, 93%, 69%, and 88%, compared to the Genetic Algorithm (GA), Particle Swarm Optimization (PSO) algorithm, Pattern Search (PS) algorithm, Harmony Search (HS) algorithm, and Honey Bee Mating Optimization (HBMO), respectively. In addition, increasing the number of intervals for flood routing significantly improved the accuracy of the results. The results indicated that the Sum of Absolute Deviations (SAD) using IBA for the Karahan flood was 117, which had reduced by 83%, 88%, 94%, and 12%, compared to the PSO, GA, HS, and BA, respectively. Furthermore, the achieved results for the Myanmar flood showed that SSD for IBA relative to GA, BA, and PSO was reduced by 32%, 11%, and 42%, respectively. In conclusion, the proposed Muskingum Model integrated with IBA considering the existence of lateral flow, outperformed the existing applied simple Muskingum models in previous studies. In addition, the more the number of intervals used in the model, the better the accuracy of flood routing prediction achieved. Water 2018, 10, 1130; doi:10.3390/w10091130 www.mdpi.com/journal/water Water 2018, 10, 1130 2 of 24


Introduction
Flood routing is fundamental to the design of structural, as well as nonstructural, flood control measures [1].Routing involves the calculation of changes in the magnitude, velocity, and shape of a flood wave, as a function of time at one or several points of the river [2].There are two types of flood routing methods: Hydraulic and hydrological.Hydraulic methods have complex computations and are more data-intensive, but describe the complete flood wave profile, whereas hydrological methods are much simpler, but yield the flood hydrograph at the end of a reach [3][4][5][6][7].The hydrologic methods need only the inflow hydrograph for a river reach.A common hydrologic method is the Muskingum method, which has several versions with parameters ranging from two to five.The two and three parameter versions of the Muskingum methods are more popular.In recent years, optimization methods, especially evolutionary algorithms, for estimating the Muskingum parameters have been popular [3].A brief background of such algorithms is now given.

Background
The flood routing models mainly include two different types of modeling: the hydraulic and hydrologic models.The hydraulic model is usually developed in a one or two-dimensional domain.Full three water shallow models and two diffusive models were used for an urban site, and the results had the same difference with each other because of different representation of a numerical and hydraulic method in the model algorithm process [5].Hunter et al [6] successfully set three explicit hydraulic models based on the inertia, diffusive, and shallow water models for flood simulation.The results indicated that the models with the shallow water equation were simple and could provide good accuracy, for the prediction of depth and velocity of the flood [5].Dottori and Todini [7] evaluated two-dimensional models based on the diffusive wave for urban floods, and the results indicated that the model could simulate the overall phenomenon well.Kim et al. [8] evaluated the different meshes in diffusive models, to investigate the effect of different meshes on the flood hydrograph.Prestininzi [9] applied the diffusive models based on the impulsive wave for inundation areas, and the results showed the model could simulate the flood conditions even in complex topography, based on a good match of simulated results with the observed data.Aricò et al. [10] applied diffusive wave equation based on 2-D numerical models for a slow varying flood, and the results showed that the simulated depth of the flood had a relatively good match with observed data.Classical, explicit finite differences in the hydraulic models and simple Muskingum model were used to investigate the flood routing [11].It has been reported that the applied numerical methods had numerical instability, for some case studies.As a result, the Muskingum model showed superior performance compared to the same applied hydraulic models [11].It has been reported that the 2-D models could be developed, based on the availability of enough information about topography and topology.With this information, a 2-D hydraulic model could successfully simulate the flood characteristics for different urban conditions.In fact, digital maps helped to identify all the required information about the boundary conditions, and to differentiate between numerous transitions within the urban hydraulic modes [12].
In addition, it has been reported that the 2-D models had had trouble in application to particular cases under small water depths, especially when the status comes close to wet/dry boundary conditions, so that there is a need for specific algorithms for simulation [6].Costabile et al. [13] reported that the main advantages of the one-dimensional model over the 2-D models for flood routing, are a simpler run process and low computational time where topographic data was unnecessary.
Fassoni-Andrade et al. [14] considered the development of a one-dimensional model based on the equation of hydrological models, which include the continuity equation and mass equation, such as the equation of the Muskingum models [14].It was observed that one-dimensional flow routing inertial models, based on the explicit solution were superior to the other models.These models simplified the Saint-Venant equation, and the main advantages of these models were good simulated results with a simple structure.Singh and Arvamuthan et al. [15] applied two hydraulic models that were developed based on the Kinematic and diffusion waves, in addition, the results were compared to the Muskingum model for the flood routing.The results showed that the simulation of hydraulic models was dependent on the kinematic wave number, so that when the value of this parameter was not considered based on accurate computation, the results for the hydraulic model could be worse than hydrologic models.Costabile et al. [16] reported that 2-D models could overcome the limitation of 1-D models, when the case study characterized as unsteady flow in irregular topography.The reports showed that if the flow was not one dimensional for the urban hydraulic, the one-dimensional channel network should be used instead of a one-dimensional model.In addition, the results showed the significant difference between 1-D and 2-D models to simulate the velocity and depths.
However, the results showed that the complex nonlinear form characteristic, numerical stability, high computational time, and complexity in the run process of hydraulic models, meant that the simpler and more accurate models have high importance [13].In fact, the hydraulic models need to measure the flow depth and discharge based on applying stream gaging.These models are known as complex models and difficult to use, whilst the hydrologic models need only to use the discharge data.In addition, the hydrologic models can be effective for the initial planning level, where the measuring system is undeveloped for accurate measurement [13].For example, Chatila [17] simulated flood routing based on the Muskingum model and EXTRAN hydraulic model.The hydraulic model developed was based on finite difference.Both hydrologic models and hydraulic models, were applied on simple and compound channels for flood routing.The results revealed that the Muskingum model had achieved higher accuracy compared with the hydraulic model because of its flexibility in calibration, where even the river bed geometry was not considered for this model.It has been demonstrated that the Muskingum model could simulate the peak discharge, achieving a close fit with the actual one, compared to the hydraulic model.Furthermore, it has been reported that hydraulic models are dependent too many assumptions, such as reach geometry, channel slope, and flow velocity, which causes the application of some hydraulic models to be limited to the specific case studies.
The Muskingum model is a useful and important hydrological model, due to its high accuracy and simplicity.Hydrological models could be accomplished after estimating the value of parameters, on the other hand, hydraulic models are required to simulate the complex boundary hydraulic conditions that causes an increase in the computational time [17].
Therefore, this model was used as a model with free access, fast computation, highly accuracy, and low cost.Furthermore, it can be used as a good tool, instead of complex hydraulic models, for flood simulation.Additional background of the application of the Muskingum model and its integration with an evolutionary algorithm, will be presented and discussed hereinafter [18].
Under a two parameter Muskingum method, Luo and Xie [19] applied the immune clonal selection algorithm (ICSA) for flood routing in a river in China, and found that the algorithm had faster convergence than the GA and PSO; and routed discharges had a high correlation with observed discharges.
Geem [20] obtained the two parameters of Muskingum method using a harmony search algorithm (HAS) for the Wilson flood in the USA, and obtained less root mean square between the predicted and observed discharges, than for GA and PSO, and less computational time.
Nelder-Mead simplex algorithm (NMSA) was considered for flood routing, and a case study in the USA [21].The parameters of the Muskingum model were considered as decision variables, and the results indicated that the RMSE (root mean square error) based on NMSA decreased by 20%, compared to the genetic algorithm [21].
Karahan et al. [22] applied a hybrid of GA, HAS, and nonlinear programming to a three-parameter Muskingum method for flood routing in a river, and found the hybrid algorithm more accurate.
Orouji et al. [23] used a genetic programming algorithm (GPA) for flood routing by the Muskingum method, and showed that GPA was more accurate than GA and PSO.The Muskingum model by 4 parameters was considered under different case studies [23].
Four parameters were considered as decision variables, and the results indicated that the model based on the considered parameters with the genetic algorithm decreased the RMSE and mean absolute error (MAE) by 20% and 25%, respectively, compared to the nonlinear programming methods [24].
The hybrid PSO and harmony search algorithm was considered for flood routing [25].The results indicated that the new hybrid method could increase the convergence velocity of the harmony algorithm, and decreased the error indexes, RMSE, and MAE, compared to the simple harmony and particle swarm algorithm.
Ouyang et al. [26] applied a hybrid of PSO and GA for Muskingum flood routing and showed that the hybrid algorithm was faster, and more accurately predicted the peak discharge and time to peak.
Under the three-parameter Muskingum method, Geem [27] found the harmony algorithm (HA) to have higher convergence than PSO and GA.Under the four-parameter Muskingum method, the Frog Leaping Algorithm (FLA) was found to have a lower computational error, than PSO and GA [28].Niknazar and Afzali [29] used an improved Honey Bee Algorithm (IHBA) to optimize three parameters of the Muskingum method, and found it to be superior to GA and PSO.
Using an Invasive Weed Optimization Algorithm (IWOA) for parameter optimization, Hamedi et al. [2] found the five-parameter Muskingum method to be more accurate than the four-parameter version.With PSO for parameter estimation, Moghadam et al. [30] found the fourparameter Muskingum method to be more accurate, than the three-parameter Muskingum method with GA and linear programming.
Using the gravitational search algorithm for parameter optimization, Kang et al. [31] found the four-parameter Muskingum method to be accurate for flood routing.Flood routing for a case study in China based on real code genetic algorithm was considered [3], and the results indicated that the four-parameter Muskingum model based on genetic algorithm decreased the RMSE and MAE of the two-and three-parameter Muskingum models.
Barati et al. [21] applied different kinds of GPA to the four-parameter Muskingum method for flood routing, and found the fixed genetic programming to be more accurate.For flood routing using the Muskingum-Cunge method, Wang et al. [23] found PSO to yield better results than GA.
In 2018, Lee [32] developed and applied an advanced Muskingum flood routing model by considering continuous stream flow, utilizing weighted inflow.Several statistical indicators have been used to evaluate the performance of the suggested model.The proposed model provided acceptable results, compared to those obtained from previous studies.The results showed that the vision corrected algorithm (VCA) had experienced a relatively small error index compared to the GA, PSO, and nonlinear programming models, for different flood case studies.
The literature review showed that the evolutionary algorithm has a high ability for obtaining the parameter values of the Muskingum model, but some algorithms have limitations, such as trapping in local optimums, slow convergence velocity, or insufficient accuracy for the simulation [33].For example, the GA can trap in the local optimums or the PSO may have an immature solution due to fast convergence [1,27].Some evolutionary algorithms have many random parameters, and accurate determination of these parameters is difficult.Thus, the improvement of previous algorithms or definition of new algorithms is necessary.

Problem Statement
Flood routing is nonlinear and multimodal, with noise.The most widely used model for flood routing is the Muskingum model.One of the main components of the Muskingum model procedure, is the inclusion of a particular number of parameters ranging between two and five, which have to be estimated based on the case study characteristics to be able to accurately route the flood.There are different procedures to estimate the Muskingum model parameters.Actually, the better the estimation of the parameters' value, the better the prediction of the flood routing characteristics achieved [33][34][35][36][37]. Having an optimization algorithm that is able to optimize the value of these parameters, whether two or three, is needed to enhance the ability of the Muskingum model to identify the flood routing characteristics.Thus, a robust algorithm that leads to the global optimal solution is needed.Bats exhibit a mysterious behavior that has long been attractive.They are able to orientate themselves to their surroundings and food acquisition, without depending on their eyesight.Bats consistently emit echolocation signals.Through analyzing the returning echoes in the auditory system, bats can distinguish their environment and find preys.By continuously watching and concentrating on the abilities of bats, scientists have suggested different bat-inspired algorithms (i.e., bat algorithm (BA) and bat intelligence (BI) algorithm) for the solution of optimization problems [33,34].Bozorg-Hadad et al. [35] used BA to optimize the use of a repository to reduce hydroelectricity energy shortage.This algorithm has been shown to have a faster convergence, than GA and PSO.Ahmadinafar et al. [36] used a hybrid BA to exploit a 10-repository system to increase energy production, which reduced computational time compared to other evolutionary algorithms.Studies have shown that BA is a powerful method but it has weaknesses, such as trapped in the local optimum or premature convergence [37,38].Thus, it needs to be improved.

Objective
In the light of the above, the use of Muskingum model showed its success when applied in flood routing prediction, but it has a few limitations.For example, it cannot be used for the complex boundary condition.In addition, one of the weaknesses of the model is the consideration of the lateral flow.In fact, most of the previous researches ignore the lateral flow while using the simple Muskingum model.Therefore, if lateral flow exists and has high volume, the Muskingum model simulates the flood without consideration of lateral flow [23][24][25][26][27][28][29][30][31], and hence, the achieved accuracy is relatively low.In fact, simplification of the Muskingum equation was a reason to omit the effect of lateral flow for flood routing, whilst there is lateral flow in the reach when the flood happens in nature.Although, there are a few number of references considering the lateral flow, some of them are limited to specific case studies with the low flow lateral condition and using Muskingum models with a simple structure [39][40][41][42].
Limited studies showed that the consideration of the lateral flow, with the help of hydraulic models, could simulate the actual situations that close to the real conditions [12].Osolivan et al. [17] reported the disadvantages of the hydraulic models for flood routing in the floodplain.Their reports showed that the initial and boundary conditions, and the resistance characteristics of main channels in the floodplain, are neccessary for the hydraulic modeling.
The present study develops the new Muskingum model for flood routing considering lateral flow.In addition, the Muskingum Model has been coupled with a new Improved Bat Algorithm (IBA), to optimize the estimation of the three-parameter Muskingum model.
The objective of this study, was to couple the three-parameter Muskingum method (TPMM) with an improved bat algorithm (IBA) for flood routing, with a multi-reach method and consideration of lateral flow.Citing easy trapping in the extremum for bat algorithm, an improved bat method is suggested and used to simulate flood routing.Chaos search tool is defined to enhance the uniformity and ergodicity of the population.Adapting weight is defined to balance the local and global search tools, for the bat algorithm.One of the advantages of the new bat algorithm is related to decreasing the search range, based on dynamic contraction.
The innovation of the study is the use of a new bat algorithm for flood routing lateral flow.Moreover, whilst previous studies usually considered one reach for flood routing without lateral flow, the present study compared the effect of dividing a river into different reaches on flood routing, and the prediction of peak discharge.

Flood Routing
The Muskingum method of flood routing is based on the continuity equation and a storage-discharge relation [30].The present study considers the lateral flow for the flood routing for the case studies based on a ration of inflow rate O lateral = βI t , while the other studies do not consider the effect of lateral flow for flood regime, and thus it adds one term to the equation of the Muskingum model with three parameters.If β coefficient equals to zero, the lateral flow has not been considered for the flood routing.
where O t is the output flow at time t, I t is the input flow at time t, S t is the storage at time t, ds t dt is the storage time variation at time t, K is the time coefficient of storage, and X is the weighting factor showing the effect of input and output flows on storage.
Equation ( 2) expresses a linear relation between storage, and input and output flows.However, a nonlinear has also been presented as: where m is an exponent.Using Equations ( 1) and ( 3), one can obtain [21]: Using S t , ∆S t , the storage later can be expressed as: Flood routing can be done using the following steps: 1. Consider initial values for parameters K, X, β, and m and enter them into the optimization algorithm, in the form of initial population.2. Calculate the storage based on Equation (3), assuming the equality of input and output flow.3. Calculate the change in storage relative to time, based on Equation (5). 4. Calculate the storage based on t + 1, according to Equation (6). 5. Calculate the output flow at t + 1, based on Equation (4).6. Repeat steps 2 to 5.

Optimization of Multi-Reach Muskingum Coefficients
The multi-reach Muskingum method is introduced to enhance the accuracy of the Muskingum method.The river under study was divided into several smaller reaches, and for each reach, routing was done separately.In other words, for each reach, parameters X, K, β, and m were calculated separately, and the output hydrograph was obtained based on the input flood hydrograph and the assumed values of X, K, β, and m for the first reach [3,21,28,31].This output hydrograph was considered as the input hydrograph for the second reach, and so on.For the second reach, the assumed values of X, K, β, and m were used, and the output hydrograph was calculated.This process was repeated for all the reaches, until the output hydrograph of the last was obtained.By comparing the computed hydrograph with the observed hydrograph, the error was calculated, and to reach the minimum error, Muskingum coefficients were optimized at all reaches.Figure 1 shows the division of the river into several reaches.
Water 2018, 10, x FOR PEER REVIEW 7 of 25 considered as the input hydrograph for the second reach, and so on.For the second reach, the assumed values of X, K, , and m were used, and the output hydrograph was calculated.This process was repeated for all the reaches, until the output hydrograph of the last was obtained.By comparing the computed hydrograph with the observed hydrograph, the error was calculated, and to reach the minimum error, Muskingum coefficients were optimized at all reaches.Figure 1 shows the division of the river into several reaches.

BAT Algorithm
Bat algorithm (BA) is based on bat sound reproduction and sound reflection.The difference in loudness that comes from the surrounding environment, allows the bat to identify the barrier from food.Bats produce very high sound pulses and listen to their return from the objects around.Each pulse remains only for a few milliseconds.BA is based on the following assumptions [30,[39][40][41]: 1.All bats have a high ability to receive sound, so that they can detect food after producing loud sounds.
2. Bats fly randomly at a velocity at place , capable of producing sound with frequency and λ wavelength.The sound produced by bats also has loudness 0 A .
3. The loudness of sound, of the bats ranges from 0 A to min A .
Each sound produced by the bat has a pulse rate (r), between 0 and 1.The sound frequency speed and position of the bats are updated as: ( ) where l f is the frequency of sound of bats, min f is the minimum frequency, max f is the maximum frequency, β is the random coefficient between 0 and 1, ( ) The following equation is used for local search in the bat algorithm:

BAT Algorithm
Bat algorithm (BA) is based on bat sound reproduction and sound reflection.The difference in loudness that comes from the surrounding environment, allows the bat to identify the barrier from food.Bats produce very high sound pulses and listen to their return from the objects around.Each pulse remains only for a few milliseconds.BA is based on the following assumptions [30,39-41]: 1.All bats have a high ability to receive sound, so that they can detect food after producing loud sounds.2. Bats fly randomly at a velocity at place y l , capable of producing sound with f min frequency and λ wavelength.The sound produced by bats also has loudness A 0 .3. The loudness of sound, of the bats ranges from A 0 to A min .
Each sound produced by the bat has a pulse rate (r), between 0 and 1.The sound frequency speed and position of the bats are updated as: where f l is the frequency of sound of bats, f min is the minimum frequency, f max is the maximum frequency, β is the random coefficient between 0 and 1, v l (t) is the velocity of the bat, Y * is the best position of the bat, y l (t) is the position of the bat, and T is the number of periods evaluated.
The following equation is used for local search in the bat algorithm: where ε is the random variable between −1 and 1, and A(t) is the loudness of sound.Loudness and pulse rates are updated, according to each stage of the iteration.For example, zero sound loudness means that the bat has found its prey and has temporarily stopped the search.
where α and γ are fixed as constant coefficients.For any value of α between 0 and 1, and γ greater than zero, A t l → 0 and r t l → t 0 l are true.Figure 2 shows the mathematical procedure of the bat algorithm.In addition, it should be noted that the random walk is considered as a parameter for the local search, for the bat algorithm.
and pulse rates are updated, according to each stage of the iteration.For example, zero sound loudness means that the bat has found its prey and has temporarily stopped the search.
( ) where α and γ are fixed as constant coefficients.For any value of α between 0 and 1, and γ greater than zero, are true.Figure 2 shows the mathematical procedure of the bat algorithm.In addition, it should be noted that the random walk is considered as a parameter for the local search, for the bat algorithm.

Improved Bat Algorithm (IBA)
The initial arrangement for the initial version of the BA is defined randomly, and it can be a reason for the uneven distribution that causes premature convergence.The chaos is a technique for improving different algorithms, where the basic idea is related to the exchange of members in the range of (−1,1).The logic mapping function is used for modulating the algorithms.Then individuals are inserted into the chaos sequence, so that it should satisfy the chaos variable space.Then, linear transformation is used to return the members to the corresponding position.The convert space is shown based on following mathematical equation: where, a i x : the initial position of the members.
The following equation is used to show the logic mapping function: Then, the elements and values are returned to the corresponding position by the following linear transformation: Furthermore, adapting weight is applied to the bat algorithm to have a good balance between

Improved Bat Algorithm (IBA)
The initial arrangement for the initial version of the BA is defined randomly, and it can be a reason for the uneven distribution that causes premature convergence.The chaos is a technique for improving different algorithms, where the basic idea is related to the exchange of members in the range of (−1,1).The logic mapping function is used for modulating the algorithms.Then individuals are inserted into the chaos sequence, so that it should satisfy the chaos variable space.Then, linear transformation is used to return the members to the corresponding position.The convert space is shown based on following mathematical equation: where x a i : the initial position of the members.The following equation is used to show the logic mapping function: Then, the elements and values are returned to the corresponding position by the following linear transformation: Furthermore, adapting weight is applied to the bat algorithm to have a good balance between global search ability and local search ability.
Water 2018, 10, 1130 The weight is computed based on the following equation: w max : the initial weight, w min : the final weight, and k: the current iteration number.The final level is related to the application of dynamic contraction, for adjusting the convergence speed: where y min,i : the lower bound position, and y max,i ; the upper bound position.The algorithm functions using the following steps: 1. Adjust the random parameters for the algorithm, such as loudness, pulsation rate, frequency, and other parameters.2. The individual position is computed using Equations ( 13)-( 15), and then the objective function is computed for each member, and the best solution is considered as Y * .3. The frequency and velocity are updated using Equations ( 7) and ( 8), and the position is computed using Equation ( 17). 4. The randomness value is compared with r l , and if r l is less than the randomness value, the distribution of the best position is acted based on 0.01 times the random disturbance.5.The local search is considered for this level.If the loudness is less than rand, the loudness should be updated and the pulsation rate should be improved using Equation ( 12). 6. Compute the objective function and change the range using Equation ( 16). 7. The convergence criterion is checked and if it is satisfied, the algorithm finishes or else the algorithm goes to step 2.

Genetic Algorithm (GA)
In GA, the initial version of the population is composed of different solutions [40].During an iterative process, subsequent populations are generated to improve the objective function.At each stage, some members from the current population are selected to generate individuals or children of the next generation, based on the fact that the likelihood of selecting people with better performance than others is more likely [41].The selected individuals produce the next population based on two genetic operators, composition and mutation.The following equations can be used for the composition operator [33].
where Pop new i is the i-th child, Pop old i is the i-th parent, Pop old j is the j-th parent, Pop j new is the j-th child, and α is a coefficient between 0 and 1.Moreover, mutation is based on the following equation: where Var low i,j is the lower limit of the i-th gene in the j-th chromosome, Var hi j,i is the upper limit of the i-th gene in the j-th chromosome, and β is a random coefficient between 0 and 1.In composition, the production of both new individuals is done by changing the gene.The mutation operator is used for the change in chromosomes and transforming their genes to create diversity in the population.

Particle Swarm Algorithm (PSO)
Initially, the process starts with a particle set.Each particle is considered as a random solution.In the next step, searches are performed sequentially to achieve the optimal answer.The i-th particle is associated with a position in an s-dimension space, where the value of s shows decision-making variables of the problem [42].The values of s variables, which determine the positions of particles are a possible solution for the optimization problem.Each particle i is completely determined by three vectors.Vector X i is the current position of the particle, Y i is the best position where the particle is iterated, and the vector of the particle velocity is shown by V i .Then the particle position and particle velocity vector are updated as: where is the new velocity of the particle, the personal learning coefficient c 2 is the global learning coefficient, Y iter * is the best solution among the solutions, and X iter+1 i is the new position of the particle.Moreover, w is the coefficient of inertia.

Indices of Error Measurement
1.The sum of squared deviations (SSD): SSD index is used as the objective function in the present study.The index calculates the total of squared deviations between observed and real discharges [28,[42][43][44][45][46]: where O obt is the observed discharge, O st is the simulated discharge, and n is the number of data.
4. Error of time to peak (ETP): The ETP index measures the difference between predicted and observed time differences of discharge [24,37,38].
T peak observed is the observed discharge, and T peak routed is the time related to the routed discharge.5. Mean absolute relative error (MARE): The mean of the relative error between observed and predicted discharges: N is the number of data.
Water 2018, 10, 1130 11 of 24 6.Varex Q (Variance index): This indicator shows the proximity of predicted and observed hydrographs with each other.
O observed mean is the observed average discharge.The closer the coefficient is to one, the more accurate the ability to predict the flood will be.
7. The agreement index (d) based on follow equation, shows the performance of the model well, so that the value of index can change from 0 to 1 [45,46].
O observed i : average of observed data.

Results and Discussion
This paper considered three case studies for flood routing.Two case studies were considered as bench problems, which have been used by different researchers using many methods for flood routing (Wilson and Karahan floods), and one case study was related to a river in Myanmar that had an important flood.The Wilson flood is considered as an important case study and different researchers tested different algorithms on this case study [2,3,[28][29][30], and thus a comprehensive study can be considered for this case study.The Karahan flood is considered as one of the case studies that have been investigated by different researchers as a benchmark problem [28][29][30][31]35].

Wilson Flood
Wilson Flood [44] is one of the most important benchmark applications, to investigate the performance of the Muskingum model and other hydrologic models.In fact, this flood pattern was generated under experimental conditions.Different mathematical models were examined using this flood data pattern, which received great attention from researchers in examining their models.
The data was extracted from Wilson [44].This information includes single peak inflow and outflow hydrographs with lateral flow.In addition, the applied algorithms in this research consider the lateral flow, although as shown in Table 4, previous researches have ignored the lateral flow because of low value, as shown in Table 6.
Several methods, such as the Segmented Least Square Method (SLSM), Hook and Jeeves (HJ)method, in combination with the Conjugate Gradient (HJ + CG), HJ method, in combination with Davidson Fletcher-Powell (HJ + DFP), nonlinear least squares (NONLER), Genetic Algorithm (GA), Harmony Search (HS), Particle Swarm Optimization (PSO), and Honey Beaming Optimization (HBMO), have used this flood without consideration of lateral flow.Given that evolutionary algorithms have random parameters, sensitivity analysis was used to determine the exact values of parameters.The evolutionary algorithms have random parameters, where the accurate values of these algorithms are computed based on sensitivity analysis.It means that the variation of objective function is determined versus the variation of parameter values, and when the objective function has the best value for a parameter, the value of this parameter is introduced as the optimal.The SSD was considered as an objective function for the current study.The frequency parameters were used to update the velocity and then the position of bats was computed based on velocity.When the objective function value is minimized, the value of the different parameters is considered at its optimal value.Tables 1-3 show the sensitivity analysis of IBA parameters, BA, GA, and PSO for flood routing in a single reach.The objective function was considered SSD for this study.
The best population associated with IBA was 60, with the lowest SSD.Moreover, the maximum frequency was 5, with the objective function as 5.01.The maximum loudness of sound was 0.6, with a random walk rate of 5. Furthermore, the mutation rate for GA was 0.6, and the recombination rate for GA was 0.7.In addition, the personal and global learning coefficient in PSO was 2, and the inertia coefficient was 0.6.When the inertia coefficient was 0.6, the objective function had the least value (10.82), and thus, the best value for the inertia coefficient was selected to be equal to 0.6 for the PSO.Other parameters can be seen in the Tables 1-3.Table 4 compares IBA with other evolutionary algorithms, for a single reach routing.The SSD value for IBA was 4.123, with SSD being reduced by 89%, 51%, 93%, 69%, 88%, and 97% compared to GA, PSO, PS, HS, HBMO, and SLSM, respectively.In addition, the SAD index was 7.112, reduced by 69%, 22%, 75%, 69%, 81%, and 84% relative to GA, PSO, PS, HS, HBMO, and SLSM, respectively.In addition, EP and MARE showed the superiority of IBA to other methods.VarexQ index for IBA, compared to other methods, showed more consistency with the predicted hydrograph.Furthermore, the performance of IBA was better than BA, so that SSD, SAD, and other error indexes for IBA were less than BA.For example, SSD and SAD for IBA were 4.123 and 7.112, whilst SSD and SAD for BA were 5.123 and 8.114.In addition, focusing on the peak discharge value, it could be depicted that the computed peak discharge based on IBA (85.11) was the most nearest estimated value to the real observed one (85), and in general, showed a worthy match with observed discharge during the whole period.In addition, for the estimation of the peak time which was predicated based on IBA, it was same to the observed Water 2018, 10, 1130 13 of 24 one 60 h.For further assessment, Table 5 shows the performance of all the models, and it could be observed that the proposed IBA in the present study outperformed all the other models.

Multi-Interval Flood Routing (Wilson Flood)
It can be seen from Table 5 that SSD, SAD, and MARE for IBA, BA, GA, and PSO for three reaches were less than for two reaches, and their values for two reaches were less than for the single reach.The best results by IBA were obtained for three reach flood routing.For flood routing using single, two, and three reaches, all evolutionary algorithms predicted the peak discharge accurately, such that the difference with the observed value was 0. For example, SSD for IBA, for one reach was 4.123, whilst it was 3.988 for three reaches.SSD and SAD for IBA for one, two, and three reaches were 1 and 85%, compared to other algorithms, and thus, IBA better performed than other algorithms.The results based on Varex Q showed that the generated hydrograph based on IBA with consideration of more value for Varex Q had a better performance than the other algorithms, and it had a high match with the observed hydrograph.The investigation of Ep showed the value of the index had decreased from 50% to 97% based on IBA, compared to the other methods for two and three intervals.In addition, the value of MARE had decreased from 6% to 56% based on IBA, for the two and three reaches.Furthermore, the agreement index (d) showed a better performance for IBA, achieving a value closer to 1 compared to the other algorithms.
An increase in the number of reaches for evolutionary algorithms increased the accuracy of flood routing.Figure 3 shows the performance of IBA for flood routing with one, two, and three reaches, which is more consistent with observed discharges.Moreover, the computational time showed better performance for the IBA, compared to the other algorithms.For example, the computational time for the IBA based on one reach was 5 s, whilst it was 7, 8, and 9 s for the BA, PSO, and GA, respectively.

Karahan Flood
Karahan et al. [22] routed a flood using various algorithms.This study considered the 1960 flood on the River Wye, Uk.River Wye is 69.75 km from Erwood to Belmont with consideration of lateral inflow, which was ignored in the previous studies, as shown in Table 7.In this study, the proposed model was evaluated considering the flood routing based on lateral inflow, which occurred in this

Karahan Flood
Karahan et al. [22] routed a flood using various algorithms.This study considered the 1960 flood on the River Wye, UK.River Wye is 69.75 km from Erwood to Belmont with consideration of lateral inflow, which was ignored in the previous studies, as shown in Table 7.In this study, the proposed model was evaluated considering the flood routing based on lateral inflow, which occurred in this Karahan flood.The data was extracted from Karahan [22].
This flood was also used in the present study.The population used for IBA was 50, the maximum frequency was 5, the minimum frequency was zero, and the maximum sound loudness was 95%.In addition, the number of chromosomes for GA was 50, the probability of mutation was 0.6, and the recombination rate was 0.7.Furthermore, the number of particles used in the particle swarm algorithm was 50, the inertia coefficient was 0.7, and the personal and global learning coefficients were 2. Table 6 shows a comparison of algorithms used in flood routing in a single reach.The value of SSD for IBA was 17,120.21,which had reduced by 81%, 87%, 84%, and 10% compared with PSO, GA, HS, BA, respectively.In addition, SAD for IBA was 117, which had reduced by 83%, 88%, 94%, and 12% compared to PSO, GA, HS, and BA, respectively.The other error indices also showed a more favorable performance of IBA, compared to other algorithms.The predicted peak discharge difference with observed discharge was 0.002 cm, which was less than the other algorithms (Table 6).The time difference between predicted and observed discharge peak time for IBA was one hour, whilst this time for other algorithms was 6 hours, so IBA had a better performance.In addition, VarexQ for IBA had a larger value than other methods, which indicated a better performance.The difference of peak discharge based on IBA with observed peak was 69 cms (the closest value to the observed one), whilst it was 135 cms, 134 cms, and 70 cms for HS, PSO, and BA, respectively.
Table 7 compares IBA, BA, GA, and PSO in flood routing (Karahan flood) for single, two, and three reaches.SSD for IBA for three reaches was 16,098.21,which was 6 percent lower than for a single reach with IBA.Moreover, SAD for IBA for three reaches was 102, which had decreased by 12% compared to a single reach.Thus, IBA had an improved performance in routing with three reaches relative to single and two intervals.This was also true for the other algorithms, as shown in Table 5. Figure 4 also shows the superior performance of IBA based on three reaches.SSD using IBA for three reaches was reduced by 46%, 51%, and 5.3% compared to GA, PSO, and BA, respectively.Furthermore, SAD using IBA for three, two, and one reaches was reduced by 23-89%, compared to the other algorithms.As a result, IBA had a better performance than BA, GA, and PSO because the error indexes had the lowest value using the IBA, compared to the other algorithms.Examining the computation time showed that IBA achieved the optimal value of the objective function "under any number of multi-reach interval" faster than the other algorithms.The lowest value of the MARE index was the one associated with the IBA, on the other hand, this value decreased by 12% and 91% when using two and three reaches, respectively.Furthermore, the d index showed the IBA algorithm and simulated hydrographs using IBA had a better match with the observed hydrograph.Many studies did not consider investigate whether parameters generated by calibration could show different floods on the same river reach.In this article, this issue was considered.The optimal parameters based on IBA and three intervals for the flood event in December 1960, were used to obtain the output hydrograph for a flood on the same river on the another flood event in January 1969, based on the inflow hydrograph on January 1969, and it was considered as the first scenario.Then, the output hydrograph was extracted based on application of obtained parameters of Muskingum models for January 1969.It meant that the model was calibrated with the 1960 storm "first scenario".

Chindwin River
One of the most important Myanmar Rivers is the Chindwin River, shown in Figure 6.This river, over the past twenty years, has experienced many floods.Due to severe flood conditions, the construction of hydraulic structures and living conditions in the downstream areas are difficult and there is high lateral inflow for this basin.Therefore, it is important to predict river flood conditions.The length of the river is 114 km and the surface of the basin is 7485 km 2 .A historical flood considered in this study occurred in July 2004, and the Mawliak Station recording the flood is shown in Figure 5. Table 8 shows a comparison of different methods for the first flood with flood routing using a single reach.Comparing IBA, BA, GA, and PSO based on routing with a single reach, SSD for IBA relative to GA, BA, and PSO was reduced by 32%, 11%, and 42%, respectively.For routing with two reaches, SSD for IBA was reduced by 33%, 34%, and 12% compared to PSO, GA, and BA, respectively.For three reaches, SSD for IBA was reduced by 44%, 50%, and 37% compared to PSO, GA and BA, respectively.The SAD index for the three routings indicated the superiority of IBA.Comparison of routing with three, two, and one reaches indicated that three reaches improved routing.SSD in   The value of SSD for the first scenario was 878 and the peak value was 280 cms, whilst the observed peak value was 285 and there was a small difference between the observed and the simulated peak value (Figure 5).The SSD for the second scenario was 869 and the peak value was 282 cms.Although, the second scenario acted better than the first scenario, the results for the first scenario were acceptable.The results were so similar because of the accurate sensitivity analysis considered for the Muskingum model for the objective function and different parameters, such as in pervious sections.Although the proposed model structure successfully predicted the flood routing, further research could consider different Muskingum models based on more parameters to investigate the performance of IBA.In addition, application of these models helps the hydraulic designers to have the accurate value for the peak discharge and the flood characteristics, to be able to optimally design the target structure.

Chindwin River
One of the most important Myanmar Rivers is the Chindwin River, shown in Figure 6.This river, over the past twenty years, has experienced many floods.Due to severe flood conditions, the construction of hydraulic structures and living conditions in the downstream areas are difficult and there is high lateral inflow for this basin.Therefore, it is important to predict river flood conditions.The length of the river is 114 km and the surface of the basin is 7485 km 2 .A historical flood considered in this study occurred in July 2004, and the Mawliak Station recording the flood is shown in Figure 5. Table 8 shows a comparison of different methods for the first flood with flood routing using a single reach.Comparing IBA, BA, GA, and PSO based on routing with a single reach, SSD for IBA relative to GA, BA, and PSO was reduced by 32%, 11%, and 42%, respectively.For routing with two reaches, SSD for IBA was reduced by 33%, 34%, and 12% compared to PSO, GA, and BA, respectively.For three reaches, SSD for IBA was reduced by 44%, 50%, and 37% compared to PSO, GA and BA, respectively.The SAD index for the three routings indicated the superiority of IBA.Comparison of routing with three, two, and one reaches indicated that three reaches improved routing.SSD in

Chindwin River
One of the most important Myanmar Rivers is the Chindwin River, shown in Figure 6.This river, over the past twenty years, has experienced many floods.Due to severe flood conditions, the construction of hydraulic structures and living conditions in the downstream areas are difficult and there is high lateral inflow for this basin.Therefore, it is important to predict river flood conditions.The length of the river is 114 km and the surface of the basin is 7485 km 2 .A historical flood considered in this study occurred in July 2004, and the Mawliak Station recording the flood is shown in Figure 5. Table 8 shows a comparison of different methods for the first flood with flood routing using a single reach.Comparing IBA, BA, GA, and PSO based on routing with a single reach, SSD for IBA relative to GA, BA, and PSO was reduced by 32%, 11%, and 42%, respectively.For routing with two reaches, SSD for IBA was reduced by 33%, 34%, and 12% compared to PSO, GA, and BA, respectively.For three reaches, SSD for IBA was reduced by 44%, 50%, and 37% compared to PSO, GA and BA, respectively.The SAD index for the three routings indicated the superiority of IBA.Comparison of routing with three, two, and one reaches indicated that three reaches improved routing.SSD in relation to routing with two and one reaches had reduced by 39% and 20%, respectively.ETp showed that there was no time lag in forecasting of peak flow by different algorithms.VarexQ showed that three-reach routing was better than two and one reach routings.Additionally, MARE for IBA-based routing had reduced by 64% and 62%, relative to one and two reach routing, respectively.Figure 7 shows that the three reach flood routing was superior to that based on two and single reaches.The results indicated that the IBA has the less computational time compared to the other algorithms, as shown in Table 6.In addition, considering the attained value of the d index, it was obvious that IBA had the best performance compared to the other algorithms.packages, such as HEC-RAS model.The development procedure could be carried out by utilizing HEC-RAS software as the model input for hydrographs and geometric characteristic, to estimate the travel times and attenuations peak, whilst the weighted coefficient (X) value could be achieved based on the optimization model and Muskingum model, as shown in Figure 8.As a result, the proposed model in this research showed the potential to be suitable and appropriate for studying and analyzing flood propagation and flood mapping.To enhance further, the used structure of the proposed Muskingum model coupled with IBA showed the highest abilities, for studying the flood plain with back water.Furthermore, such enhancement for the Muskingum model could be successfully applied as an easy and inexpensive method for computing the time and shape, for an overbank flood, when there is back water and inertia influences along a river channel.In fact, with the new nature-inspired optimization algorithms, the traditional Muskingum method could be integrated with the hydrodynamic software packages, such as HEC-RAS model.The development procedure could be carried out by utilizing HEC-RAS software as the model input for hydrographs and geometric characteristic, to estimate the travel times and attenuations peak, whilst the weighted coefficient (X) value could be achieved based on the optimization model and Muskingum model, as shown in Figure 8.As a result, the proposed model in this research showed the potential to be suitable and appropriate for studying and analyzing flood propagation and flood mapping.

Conclusions
The present study investigated the potential of utilizing a three-parameter Muskingum model coupled with Improved Bat Algorithm (IBA) to accurately predict flood routing.Three different case studies have been used in this study, to evaluate the performance of the proposed model.These three case studies were the Wilson flood, Karahan flood, and Myanmar River.Seven different performance indices were used to examine and compare the performance of the proposed model over other algorithms.In addition, discretization of the river stream was considered to improve model accuracy.The results showed that IBA outperformed all other algorithms and was able to reduce the SSD by percentages ranging between 20% and 84%, compared with the other algorithms.In addition, the achieved results using the IBA could predict the peak discharge accurately, with a value very close to the observed one.Under the Karahan flood, IBA considerably achieved the minimum level of error indices for a single reach, compared to other algorithms.Finally, IBA in flood routing with three intervals had a better performance than with single and two reaches.The division of the river into different reaches increased the accuracy of flood routing.The performance of IBA for the river in Myanmar, also showed that the simulated hydrograph with three reaches was more accurate.For

Conclusions
The present study investigated the potential of utilizing a three-parameter Muskingum model coupled with Improved Bat Algorithm (IBA) to accurately predict flood routing.Three different case studies have been used in this study, to evaluate the performance of the proposed model.These three case studies were the Wilson flood, Karahan flood, and Myanmar River.Seven different performance indices were used to examine and compare the performance of the proposed model over other algorithms.In addition, discretization of the river stream was considered to improve model accuracy.The results showed that IBA outperformed all other algorithms and was able to reduce the SSD by percentages ranging between 20% and 84%, compared with the other algorithms.In addition, the achieved results using the IBA could predict the peak discharge accurately, with a value very close to the observed one.Under the Karahan flood, IBA considerably achieved the minimum level of error indices for a single reach, compared to other algorithms.Finally, IBA in flood routing with three intervals had a better performance than with single and two reaches.The division of the river into different reaches increased the accuracy of flood routing.The performance of IBA for the river in Myanmar, also showed that the simulated hydrograph with three reaches was more accurate.For example, the computational time for the IBA based on three intervals was 2, 4, and 6 s less than, BA, PSO, and GA, respectively.Furthermore, the EP for the IBA was 33%, 96%, and 97% less than BA, PSO, and GA, respectively.As a result, the proposed Muskingum model coupled with the IBA, could be considered as a strong alternative method for predicting flood routing characteristics.

Figure 1 .
Figure 1.Discretization of the river stream.
is the velocity of the bat, Y * is the best position of the bat, the position of the bat, and T is the number of periods evaluated.

Figure 1 .
Figure 1.Discretization of the river stream.

Figure 2 .
Figure 2. The flow chart of the bat algorithm.

Figure 2 .
Figure 2. The flow chart of the bat algorithm.

Figure 3 .
Figure 3.The simulated Hydrographs for Wilson flood using different methods.

Figure 3 .
Figure 3.The simulated Hydrographs for Wilson flood using different methods.

Figure 5 .
Figure 5.The first and second scenario for calibration of model.

Figure 5 .
Figure 5.The first and second scenario for calibration of model.

Figure 5 .
Figure 5.The first and second scenario for calibration of model.

Figure 6 .
Figure 6.The location of the River in Myanmar.Figure 6.The location of the River in Myanmar.

Figure 6 .
Figure 6.The location of the River in Myanmar.Figure 6.The location of the River in Myanmar.Water 2018, 10, x FOR PEER REVIEW 22 of 25

Figure 8 .
Figure 8. Advanced Muskingum model for flood plain with complex boundary condition.

Figure 8 .
Figure 8. Advanced Muskingum model for flood plain with complex boundary condition.

Table 1 .
Analysis for the performance of Improved Bat Algorithm (IBA) algorithm (Wilson flood).

Table 2 .
Analysis for the performance of the Genetic Algorithm (GA) (Wilson flood).

Table 3 .
Analysis for the performance of the Particle Swarm Optimization (PSO)(Wilson Flood).

Table 4 .
Computed error indexes for Wilson flood.

Table 5 .
Flood routing for Wilson flood according to three different intervals.

Table 6 .
Inflow and outflow for Karahan flood.

Table 7 .
Results of four models (i.e., IBA, BA, PSO and GA) for Karahan flood according to three different intervals.

Table 8 .
The computed error indexes for Myanmar River.