*Article* **Development of a New 8-Parameter Muskingum Flood Routing Model with Modified Inflows**

**Eui Hoon Lee**

School of Civil Engineering, Chungbuk National University, Cheongju 28644, Korea; hydrohydro@chungbuk.ac.kr; Tel.: +82-043-261-2407

**Abstract:** Flood routing can be subclassified into hydraulic and hydrologic flood routing; the former yields accurate values but requires a large amount of data and complex calculations. The latter, in contrast, requires only inflow and outflow data, and has a simpler calculation process than the hydraulic one. The Muskingum model is a representative hydrologic flood routing model, and various versions of Muskingum flood routing models have been studied. The new Muskingum flood routing model considers inflows at previous and next time during the calculation of the inflow and storage. The self-adaptive vision correction algorithm is used to calculate the parameters of the proposed model. The new model leads to a smaller error compared to the existing Muskingum flood routing models in various flood data. The sum of squares obtained by applying the new model to Wilson's flood data, Wang's flood data, the flood data of River Wye from December 1960, Sutculer flood data, and the flood data of River Wyre from October 1982 were 4.11, 759.79, 18,816.99, 217.73, 38.81 (m3/s)2, respectively. The magnitude of error for different types of flood data may be different, but the error may be large if the flow rate of the flood data is large.

**Keywords:** hydrologic flood routing; Muskingum flood routing model; meta-heuristic optimization; self-adaptive vision correction algorithm

#### **1. Introduction**

Water resources from rivers are sources of hydroelectric power generation, agricultural water, and industrial water; however, owing to the large volumes of water, such rivers are prone to floods that have adverse impacts on life and property [1]. To reduce or prevent such damage, engineering measures, such as the construction of flood control dams or flood walls (levees), are necessary. Therefore, the evaluation of engineering measures for flood control is critical, and these measures are generally directly related to flood routing. Flood routing can be defined as a procedure for determining the flood hydrograph at a point downstream from the base flood hydrograph at an upstream point. In other words, flood routing is the process of determining the amount by which a flood wave is reduced and how long it takes for a flood wave to pass through an arbitrary section of a river based on the amount of storage in that section.

There are two types of flood routing methodologies: hydraulic and hydrologic [2]. Hydraulic flood routing is a method for solving the partial differential continuity and momentum equations, the governing equations of an unsteady nonuniform flow, which hydraulically represent the flow of the flood wave according to the initial and boundary conditions [3]. In contrast, the hydrologic flood routing method yields an approximate solution using the storage equation based on the continuity equation of the flood wave [2]. The hydrologic flood routing method can be divided into three categories: reservoir routing, channel routing, and watershed routing. Channel routing allows the measurement of the storage effect of natural rivers on flood waves by calculating how the discharge of a flood changes as it progresses downstream and provides a standard hydrologic quantity for river planning. The Muskingum flood routing model is a representative channel-routing model [2].

**Citation:** Lee, E.H. Development of a New 8-Parameter Muskingum Flood Routing Model with Modified Inflows. *Water* **2021**, *13*, 3170. https://doi.org/10.3390/w13223170

Academic Editor: Fi-John Chang

Received: 4 September 2021 Accepted: 5 November 2021 Published: 10 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

The first Muskingum flood routing model proposed was the linear Muskingum flood routing model (LMM) with two parameters [4]. However, the LMM did not include lateral inflow, and a new Muskingum flood routing model with three variables (LMM-L) was thus proposed [5]. Additionally, a study using the nonlinear Muskingum flood routing model (NLMM) considered the nonlinear relationship between storage and outflow as means of improving upon the LMM [6]. Two types of NLMMs determined by the location of nonlinear factor have been proposed to calculate the storage [7]. The Broyden– Fletcher–Goldfarb–Shanno technique based on a mathematical gradient was applied to the NLMM [8]. NLMM incorporating lateral flow (NLMM-L) was developed to complement the existing NLMM [9]. In 2018, a Muskingum flood routing model called the advanced NLMM (ANLMM-L) was used for calculating a continuous inflow [2]. ANLMM-L is a type of nonlinear Muskingum flood routing model that considers lateral inflow and continuous flow with time. Generalized storage equations for the NLMM have also been suggested to apply more degrees of freedom in the suggested model [10].

In addition to the aforementioned studies, various other investigations have focused on recalculating the error between the outflow from flood data and the calculated outflow. Various studies on Muskingum flood routing models were conducted before the 2000s. The two parameters of the LMM including nonlinear relation between the storage and weighted flow were determined using the least-squares method [11]. The least-squares method was used to adjust the two parameters of LMM, K and X. The Muskingum parameter estimation/flood routing system was developed for linear LMMs and NLMMs and their results have been compared [12]. The results of the two different Muskingum flood routing models were compared. The genetic algorithm was used to estimate the parameters of the NLMMs [13]. In order to overcome the limitations of traditional methods used for Muskingum flood routing models, the genetic algorithm, a well-known meta-heuristic optimization algorithm, was applied.

Since the 2000s, studies applying various meta-heuristic optimization algorithms to the Muskingum flood routing models have continued. An immune clonal selection algorithm was suggested to improve the convergence speed and it was applied to estimate parameters of the NLMM [14]. The Nelder–Mead simplex algorithm was introduced to improve the usability, and it was used to estimate parameters of the NLMM [15]. Furthermore, the simulated annealing and shuffled frog leaping algorithms were used to estimate the parameters of the Muskingum flood routing model in two benchmark/real case studies, and they were compared with the results of Tung's method [16]. The honeybee mating optimization algorithm with past convergence speed has also been applied for the parameter estimation of the NLMM [17]. The elitist-mutated particle swarm optimization and improved gravitational search algorithm were applied to estimate the parameters of LMMs and NLMMs [18]. Particle swarm optimization was applied to the parameter estimation of the NLMM with four parameters to fit the multiple-peak hydrographs [19], and various NLMMs with different storage calculations such as parameterized initial storage have been proposed using a weed optimization algorithm [20]. Various meta-heuristic optimization algorithms, such as the genetic algorithm, evolution, particle swarm, and a harmony search have been used for parameter estimations of the nonlinear Muskingum model and the variable parameter McCarthy–Muskingum model [21]. The adaptive genetic algorithm was used to estimate the various exponent parameters of the NLMM and it was applied to Wilson's flood data [22]. In addition, genetic expression programming with faster convergence speed than existing genetic programming was developed for parameter estimation in the Muskingum flood routing model [23]. The water cycle algorithm was applied to estimate the parameters of the NLMM and compared with the genetic algorithm, particle swarm optimization, harmony search, and imperialist competitive algorithm [24]. Although various meta-heuristic optimization algorithms have been tested, studies focusing on comparing the results of each algorithm have indicated limited improvements for the Muskingum flood routing model.

Studies have also been conducted on the application of hybrid meta-heuristic optimization algorithms, combining a charged system search and particle swarm optimization, for parameter estimation of the Muskingum flood routing models [25]. For example, a hybrid meta-heuristic optimization algorithm combining particle swarm optimization and the Nelder–Mead simplex method was used to estimate the parameters of the Muskingum flood routing model [26]. Parameter estimation was conducted using a hybrid metaheuristic optimization algorithm combining the shuffled frog leaping algorithm and the Nelder–Mead simplex method [27]. The improved real-coded adaptive genetic algorithm and the Nelder–Mead simplex algorithm were combined for the parameter estimation of two improved NLMMs [28]. The hybrid meta-heuristic optimization algorithm applied in particle swarm optimization and the bat algorithm were used to reduce the computational time of the Muskingum flood routing model [29]. Although good results were obtained in some studies, it is difficult to accurately compare them with the results of other existing Muskingum flood routing models because they were calculated using additional variables.

General improvements of the Muskingum flood routing model have also been considered. For example, a modified Muskingum flood routing approach, in conjunction with the HEC-RAS model, was implemented to determine floodplain flows [30]. A new NLMM with four parameters has been suggested [31]. A parameter estimation method of the Muskingum flood routing model in ungagged channel reaches has also been suggested [32].

Studies have been conducted to apply the hybrid method to Muskingum flood routing models. A hybrid harmony search combined with local search algorithm such as Broyden–Fletcher–Goldfarb–Shanno technique was developed and was applied to estimate parameters in NLMM [33]. The new hybrid optimization technique was suggested by combining the modified honeybee mating optimization and generalized reduced gradient algorithm for the application of the new Muskingum model with six parameters [34]. The particle swarm optimization hybridized with Nelder–Mead simplex method was proposed to improve precision and convergence speed in Muskingum model [26]. The hybrid algorithm combining the shuffled frog leaping algorithm and Nelder–Mead simplex was applied to NLMM with four parameters and NLMM with five parameters, and it was compared with the genetic algorithm-generalized reduced gradient [27]. The improved NLMM was suggested for flood prediction using the hybrid algorithm of particle swarm optimization and bat algorithm and it was compared with particle swarm optimization and bat algorithm [29]. The parameters in the two types of NLMM were estimated to improve precision using the hybrid algorithm combining the improved real-coded adaptive genetic algorithm and the Nelder–Mead simplex [28]. Most of the previously proposed hybrid methods combine optimization algorithms, but the hybrid method of this study is a method combining the inflow at the previous time and the inflow at the next time in the Muskingum flood routing. The honey bee mating optimization algorithm was combined with the generalized reduced gradient algorithm to estimate parameters of improved Muskingum flood routing model and applied to the single and multi-peak flood hydrographs [35].

In this study, a new Muskingum flood routing model was suggested. The new Muskingum flood routing model, which considers continuous inflow at previous and next time in the storage and inflow calculations, can enable accurate flood routing in various flood data. The self-adaptive vision correction algorithm (SAVCA), a recently developed meta-heuristic optimization algorithm, was applied to calibrate various parameters in the new Muskingum flood routing model. SAVCA can overcome the disadvantages of the previously developed vision correction algorithm (VCA). When applied to mathematical benchmark functions and water distribution problems, SAVCA has previously displayed good performance [36]. Various meta-heuristic optimization algorithms as well as SAVCA can be applied to the new Muskingum flood routing model to show good results. In the previous study, the type of meta-heuristic optimization algorithm did not significantly affect the results of Muskingum flood routing models [37].

#### **2. Materials and Methodologies**

#### *2.1. Overview*

The two primary methods used in this study were the SAVCA and new Muskingum flood routing model. The error between the flood outflow data and the calculated outflow in the new Muskingum flood routing model was used as an objective function in the SAVCA, which applied an iterative calculation to minimize the error. The iterative calculation progresses as follows:


The solution refers to the values of the parameters for Muskingum flood routing models. The calculation in the SAVCA is as follows: If all initial solutions are calculated according to the new Muskingum flood routing model process, the errors of the initial solutions are calculated and sorted in ascending order. SAVCA consists of two types of parameters: self-adaptive and fixed. Division rate 1 (DR1), division rate 2 (DR2), and the compression factor (CF) are self-adaptive parameters. The modulation transfer function rate (MR) and astigmatic rate (AR) are fixed parameters.

DR1 determines whether a new solution should be generated in the range of each variable (global search) or if one solution should be selected from the existing solution group (local search). If the generation of a new solution is determined in DR1, the positive and negative direction searches are determined by DR2. The decision variables of the new solution, generated by global search or selected by local search, are adjusted in detail by MR, CF, and AR. After generating a new solution, the calculation process of the new Muskingum flood routing model is applied. The error (SSQ) is calculated for the inflow, storage, and outflow during each time period. If the error of the new solution is lower than that of the worst solution among the existing solution groups, the new solution is included in the existing solution group. DR1 and DR2 are adjusted according to the calculation process for the new solution. All processes are repeated until a certain number of iterations of SAVCA.

The SSQ was used to calculate the first error value for the Muskingum flood routing models in this study. The SSQ between the observed and calculated outflows was used as the objective function in the optimization process. In the new Muskingum flood routing model, eight parameters were used as decision variables, and the objective function is shown in Equation (1).

$$\text{Minimize SSQ} = \sum \left( \text{O}\_{\text{o}} - \text{O}\_{\text{s}} \right)^{2} \tag{1}$$

where Oo is the observed outflow (m3/s), and Os is the calculated outflow (m3/s). The NSE was used to calculate the second error value for the Muskingum flood routing models in this study. The equation of NSE is shown in Equation (2).

$$\text{NSE} = 1 - \frac{\sum\_{i=1}^{n} \left(\text{O}\_{\text{o}} - \text{O}\_{\text{s}}\right)^{2}}{\sum\_{i=1}^{n} \left(\text{O}\_{\text{o}} - \overline{\text{O}}\_{\text{o}}\right)^{2}} \tag{2}$$

where Oo is the average of observed outflow (m3/s), and n is the number of data. The RMSE was used to calculate the third error value for the Muskingum flood routing models in this study. The equation of RMSE is shown in Equation (3).

$$\text{RMSE} = \sqrt{\frac{\sum\_{i=1}^{n} \left(\text{O}\_{\text{o}} - \text{O}\_{\text{s}}\right)^{2}}{n}} \tag{3}$$

#### *2.2. New Muskingum Flood Routing Model*

The initial LMM was calculated by assuming the amount of storage in the channel as the sum of prism storage and wedge storage. Prism storage is proportional to the outflow, and wedge storage is proportional to the difference between the inflow and outflow. In the LMM, storage is calculated using Equation (4).

$$S\_t = K[XI\_t + (1 - X)O\_t] \tag{4}$$

where *St*, *It*, and *Ot* are the storage, inflow, and outflow at time t, respectively, and *X* is the weighted factor. In the NLMM, the nonlinear factor is added in the exponential form of Equation (4). Equation (5) represents the storage in an NLMM.

$$S\_t = K[XI\_t + (1 - X)O\_t]^m \tag{5}$$

where *m* is a nonlinear factor, namely different from 1. The storage calculation in the new Muskingum flood routing model is applied by considering an existing generalized storage. The storage is calculated by considering not only the inflow at the current time (*t*) but also the inflow at the next time point (*t* + 1). The reason for considering the inflow at *t* + 1 instead of *t* − 1 in the storage calculation is as follows. In the nonlinear Muskingum flood routing models, it is assumed that the storage at time *t* depends on the upstream storage *Sin*, and the downstream storage *Sout*. Inflow (*I*), outflow (*O*), *Sin*, *Sout* were organized as follows according to water depth [38]. *I* and *Sin* are shown in Equations (6) and (7).

$$I = a\_1 y^{c\_1} \tag{6}$$

$$S\_{\rm in} = b\_1 y^{d\_1} \tag{7}$$

where *y* is the water depth and *a*1, *b*1, *c*1, *d*<sup>1</sup> are coefficients. If *c*<sup>1</sup> and *d*<sup>1</sup> are equal, then *Sin* is shown in Equation (8). 

$$S\_{\rm ini} = b\_1 \left(\frac{I}{a\_1}\right) \tag{8}$$

*O* and *Sout* are shown in Equations (9) and (10).

$$O = a\_2 y^{c\_2} \tag{9}$$

$$S\_{out} = b \mathfrak{z} y^{d\_2} \tag{10}$$

where *a*2, *b*2, *c*<sup>2</sup> and *d*<sup>2</sup> are coefficients. If *c*<sup>2</sup> and *d*<sup>2</sup> are equal, then *Sout* is shown in Equation (11).

$$S\_{\rm out} = b\_2 \left(\frac{O}{a\_2}\right) \tag{11}$$

In NLMM, inflow, outflow and storage are assumed to be water depth related. The storage can be summarized in Equation (12) [38].

$$S = XS\_{in} + (1 - X)S\_{out} \tag{12}$$

The storage with *K* = *b*1/*a*<sup>1</sup> = *b*2/*a*<sup>2</sup> can be expressed as Equation (13).

$$S = KXI + K(1 - X)O\tag{13}$$

Nonlinear parameter *m* was applied to Equation (11) and it can be expressed as Equation (14).

$$S = K[XI + (1 - X)O]^m \tag{14}$$

Additionally, it is assumed that there is an interdependence between the storage at time *t* and storage at time *t* + 1 in generalized storage [10].

$$S\_t = X\_1 S\_{in, \, t} + X\_2 S\_{in, \, t+1} + (1 - X\_1 - X\_2) S\_{out, \, t} \tag{15}$$

Equation (16) can be rearranged as in the process from Equations (12)–(14) and it represents storage in the new Muskingum flood routing model.

$$S\_t = K[X\_1 I\_t + X\_2 I\_{t+1} + (1 - X\_1 - X\_2)O\_t]^m \tag{16}$$

where *X*<sup>1</sup> is the weighted factor at time *t*, and *X2* is the weighted factor at time *t* + 1. In addition, *It* is the inflow at time *t*, and *It*+1 is the inflow at time *t* + 1. Based on Equation (4), the outflow calculation is summarized in Equation (17).

$$O\_t = \frac{1}{\left(1 - X\_1 - X\_2\right)} \left(\frac{S\_t}{K}\right)^{\frac{1}{m}} - \frac{X\_1}{\left(1 - X\_1 - X\_2\right)} I\_t - \frac{X\_2}{\left(1 - X\_1 - X\_2\right)} I\_{t+1} \tag{17}$$

A weighted inflow, including a continuous inflow has been proposed previously [2]. In this study, the inflow at time *t* + 1 was included to consider the additional continuous inflow, and the weighted inflow could be calculated as shown in Equation (18).

$$\mathcal{W}\_{l} = \left[ (1 - \theta\_1 - \theta\_2 - \theta\_3)I\_l + \theta\_1 I\_{t-1} + \theta\_2 I\_{t-2} + \theta\_3 I\_{t+1} \right] \tag{18}$$

where *Wt* is the weighted inflow at time t, and *θ<sup>1</sup>* is the weighted factor of the previous inflow at time *t* − 1. In addition, *θ<sup>2</sup>* is the weighted factor of the previous inflow at time *t* − 2, and *θ<sup>3</sup>* is the weighted factor of the next inflow at time *t* + 1. In previous studies, the inflow at time *t* − 1 (*It*−1) and the inflow at time *t* − *2* (*It*−2) were considered [2,9]. In this study, all inflows before and after the current time were considered by including the inflow at time *t* +1(*It*+1). If the weighted inflow of Equation (18) is substituted into Equation (17), the outflow is calculated using Equation (19).

$$O\_l = \frac{1}{\left(1 - X\_1 - X\_2\right)} \left(\frac{S\_l}{K}\right)^{\frac{1}{m}} - \frac{X\_1}{\left(1 - X\_1 - X\_2\right)} W\_l - \frac{X\_2}{\left(1 - X\_1 - X\_2\right)} W\_{l+1} \tag{19}$$

where *Wt*+1 is the weighted inflow at time *t* + 1. The storage at time *t* + 1 can be calculated from the outflow in Equation (19) and the observed inflow. The general storage equation is calculated as Equation (20).

$$\frac{dS}{dt} = I\_t - O\_t \tag{20}$$

Equation (21) shows the modified storage equation that considers a change in lateral flow.

$$\frac{dS}{dt} = \frac{\Delta S}{\Delta t} = (1 + \beta)I\_t - O\_t \tag{21}$$

where *β* is the parameter accounting for the lateral flow. The storage at time *t* + 1 is shown in Equation (22).

$$S\_{t+1} = S\_t + \Delta S \tag{22}$$

Equation (23) represents the storage at time *t* + 1, and it can be obtained by substituting Equation (21) into Equation (22).

$$\mathbf{S}\_{t+1} = \mathbf{S}\_{t} + [(1+\beta)I\_{t} - O\_{t}]\Delta t \tag{23}$$

where *St*+1 is the storage at time *t* + 1. The outflow in the new Muskingum flood routing model is calculated using eight variables, i.e., *K*, *X*1, *X*2, *m*, *β*, *θ*1, *θ*2, and *θ*3. The initial Muskingum flood routing model is based on mass conservation. Therefore, the new Muskingum flood routing model was calculated based on the mass conservation.

#### *2.3. Self-Adaptive Vision Correction Algorithm*

SAVCA has a total of six parameters: DR1, DR2, MR, CF, AR, and AF. Among these, DR1, DR2, and CF are self-adaptive, and MR, AR, and AF are fixed. Table 1 shows the parameter types of SAVCA [36].

**Table 1.** Parameter types of SAVCA.


In SAVCA, the initial decision variables and decision variables generated by the global search are randomly generated within the range between the upper and lower boundaries. Decision variables in the global search are between the current optimal value and the upper boundary or between the lower boundary and the current optimal value based on the probability of DR2. The decision variable generated between the current optimal decision variable and the upper boundary by a global search is shown in Equation (24).

$$\mathbf{x}\_{n} = \mathbf{x}\_{b} + random(0, \ 1) \times (b\_{u} - \mathbf{x}\_{b}) \tag{24}$$

where *xn* is the new decision variable, and *xb* is the current optimal value. In addition, *random*(0, 1) is a random value generated from 0 to 1, and *bu* is the upper boundary. The decision variable generated between the lower boundary and the current optimal decision variable through a global search is shown in Equation (25).

$$\mathbf{x}\_{\text{il}} = b\_{\text{l}} + random(\mathbf{0}, \mathbf{1}) \times (\mathbf{x}\_{\text{b}} - b\_{\text{l}}) \tag{25}$$

where *bl* is the lower boundary. In SAVCA, each decision variable is adjusted by the parameters used in the local search, such as MR, CF, AR, and AF. Equation (26) shows the calculation of the new decision variable.

$$\mathbf{x}\_n = \mathbf{x}\_n \times \left\{ 1 + MTF \times random(-1, 1) \times \left( 1 - \frac{current \ iteration}{total \ iteration} \right)^{\rm CF} \right\} \tag{26}$$

where *MTF* is the calculated modulation transfer function value. The term *random*(−1, 1) is a random value between −1 and 1. In addition, *CF* is a parameter for lens compression. The calculation of MTF is based on the distance (*dx*) between the current best decision variable and the selected decision variable. *dx* can be calculated using Equation (27).

$$d\mathbf{x} = \frac{rank(\mathbf{x}\_s) - rank(\mathbf{x}\_b)}{rank(\mathbf{x}\_l) - rank(\mathbf{x}\_b)} \tag{27}$$

where *dx* is the relative distance between each decision variable, *rank(xs)* is the fitness rank of the selected decision variable (*xs*), *rank(xb)* is the fitness rank of the best decision variable (*xb*), and *rank(xl)* is the fitness rank of the worst decision variable (*xl*). The fitness rank is the order in which the values of the objective function are sorted. The worst decision variable has the lowest fitness rank. Figure 1 shows the relative distance of *dx*.

**Figure 1.** Relative distance of *dx*.

The calculation of the MTF by applying *dx* is shown in Equation (28).

$$MTF\_s = \left(\frac{dx\_s}{\left(\sum\_{i=1}^k dx\_i^2\right)^{0.5}}\right)^{0.5} \tag{28}$$

where *MTFs* is the *MTF* value of the selected decision variable, and *k* is the total number of decision variables. In addition, *dxs* is the relative distance of the selected decision variable, and *dxi* is the relative distance of the *i*-th decision variable. The *CF* in SAVCA can be calculated as shown in Equation (29).

$$\mathcal{CF} = 10 \times \left\{ \frac{standard \, deviation(\mathbf{x}\_i)}{average(\mathbf{x}\_i)} \right\} \tag{29}$$

where *xi* is the *i*-th decision variable. The probability of applying the astigmatism correction process is determined using the AR. The new decision variable adjusted by the application of the astigmatism correction process on a local search is shown in Equation (30).

$$\mathbf{x}\_n = \mathbf{x}\_n \times \left\{ 1 + random(-1, 1) \times \sin^2(AF) \right\} \tag{30}$$

where *AF* is the angle of the astigmatic axis. The application process of SAVCA is summarized as follows: (1) generation of an initial solution group, (2) calculation of the fitness of the initial solution groups, (3) generation of a new solution, (4) application of MR and AR, (5) decision to replace after comparing the new solution with the worst solution in the existing solution group, and (6) repeating (2)–(5) until the total number of iterations is reached. The worst solution is the solution with the lowest fitness rank among the existing solution group. Figure 2 shows the application process.

**Figure 2.** Application process.

An initial solution group was created to apply SAVCA to the new Muskingum flood routing model. The weighted inflow, storage, outflow, and SSQ were then calculated. According to the probability of DR1, a new solution was generated by a global search, or one of the existing solution groups was selected through a local search. When creating a new solution with a global search, a new decision variable was created in the positive and negative directions based on the current best decision variable. Each new decision variable was corrected using the MR and AR. In addition, the weighted inflow, storage, outflow, and SSQ were calculated using the new solution with new decision variables. Whether a new solution should be added to the existing solution group was determined by comparing the SSQ of the new solution with the SSQ of the worst solution among the existing group of solutions.

#### *2.4. Flood Data*

Five types of flood data were applied to the Muskingum flood routing models: Wilson's flood data, Wang's flood data, flood data for River Wye December in 1960, Sutculer flood data, and flood data for River Wyre October in 1982 [5,39–41]. All flood data used in this study have been applied in several Muskingum flood routing models in existing studies. The most important aspect of the Muskingum flood routing model is the range of each parameter. The range of each parameter in the new Muskingum flood routing model applied to the five flood datasets is presented in Table 2.


**Table 2.** Range of parameters in the new Muskingum flood routing model.

Various Muskingum flood routing models, i.e., LMM, LMM-L, NLMM, NLMM-L, ANLMM-L, and the new Muskingum flood routing model, were compared herein. The parameters used in each Muskingum flood routing model are listed in Table 3.

**Table 3.** Parameters used in each Muskingum flood routing model (-: applied, X: not applied).


Differences were observed in the results of the Muskingum flood routing models proposed in other studies. However, the simulation used in this study was conducted according to the parameters listed in Table 3. The data values from existing studies were considered only up to two decimal points when using them as an input for the models in this study. The parameters of Muskingum flood routing models without results from previous studies were obtained by applying SAVCA. However, the values of each parameter were all calculated differently.

#### **3. Application and Results**

The first flood dataset used for the application of all Muskingum flood routing models was Wilson's flood data. The parameters of the LMM for Wilson's flood data were determined to be 29.164640 for K, and 0.118200 for X1. The parameters of the new Muskingum flood routing model for Wilson's flood data were determined to be 0.943442 for K, 0.340333 for X1, −0.00102 for X2, 1.744439 for m, −0.02166 for β, 0.758873 for θ1, 0.230779 for θ2, and 0.047773 for θ3. The results, including those obtained using the new Muskingum flood routing model, are compared in Table 4.


**Table 4.** Results when using Wilson's flood data.

Among the results in Table 4, those of the LMM-L, NLMM, NLMM-L, and ANLMM-L were calculated in previous studies [2,5,9,42]. The results of LMM and new Muskingum flood routing model were calculated using Wilson's flood data by applying SAVCA. Notably, the results of the LMM are better than those of the LMM-L. This is because the results of the LMM-L are the results of a previous study wherein the optimization method was not used, while the results of the LMM were obtained using SAVCA. It should be noted that the results of the new Muskingum flood routing model were better than those of the LMM, LMM-L, NLMM, NLMM-L, and ANLMM-L because the new Muskingum flood routing model showed the smallest error in the initial part from 0 to 24 h and showed the smallest error in the overall results. Because the errors of the NLMM-L and ANLMM-L were small, the new Muskingum flood routing model did not lead to a substantial improvement.

Among the existing Muskingum flood routing models, ANLMM-L showed the best results (smallest error). The difference in SSQ between the ANLMM-L and new Muskingum flood routing model was 0.43 (m3/s)2. The differences between the results of the two models were most notable from 6 to 18 h and from 108 to 126 h. The new Muskingum flood routing model showed the closest outflow to the observed outflow from 6 to 18 h. Among other Muskingum flood routing models, the outflows obtained from LMM-L at 6 h, LMM at 12 h, and NLMM-L at 18 h were close to the observed outflow. The difference between the observed outflow and the outflow obtained from the new Muskingum flood routing model was smaller than that attained using other Muskingum flood routing models.

The second flood dataset was Wang's flood data. The parameters of the LMM-L for Wang's flood data were determined to be 1.075331 for K, −0.762101 for X1, and −0.003024 for β. The parameters of the new Muskingum flood routing model for Wang's flood data were determined to be 0.079266 for K, −1.49742 for X1, 0.011592 for X2, 1.360300 for m, −0.000450 for β, 0.421275 for θ1, 0.044483 for θ2, and 0.261537 for θ3. The results using Wang's flood data, including those for the new Muskingum flood routing model, are compared in Table 5.

The results of LMM, NLMM, NLMM-L, and ANLMM-L were calculated in previous studies [2,8,9,39]. The results of LMM-L and new Muskingum flood routing model were calculated by applying SAVCA. Notably, the results improved dramatically when the lateral inflow was considered, indicated by the difference between the results of the LMM and LMM-L and between those of the NLMM and NLMM-L. A difference between the results of the ANLMM-L and new Muskingum flood routing model was also observed, although it was not due to a lateral inflow but to differences in the calculation equations of the weighted inflow and storage. The results of the new Muskingum flood routing model were overwhelmingly better than those of other Muskingum flood routing models because the new Muskingum flood routing model showed the smallest error in the latter part from 19 to 29 h and showed the smallest error in the overall results.

Among the existing Muskingum flood routing models, the ANLMM-L showed the best results (smallest error). The difference in SSQ between the ANLMM-L and new Muskingum flood routing model was 149.56 (m3/s)2. The difference between the two results was clear from 228 (19) to 288 (24) h. The differences in the outflow obtained using the Muskingum flood routing models and the observed outflow was small. The outflow calculated by the new Muskingum flood routing model was closest to the observed outflow. Among other Muskingum flood routing models, the outflow obtained using NLMM at 21 h was close to the observed outflow. The difference between the observed outflow and the outflow obtained using the new Muskingum flood routing model was smaller than that when using other Muskingum flood routing models.

The third flood dataset was the flood data of River Wye December in 1960. The parameters for the LMM using these data were determined to be 23.877307 for K, and 0.153174 for X1. The parameters for the new Muskingum flood routing model using these data were determined to be 2.318963 for K, 0.499999 for X1, 0.000390 for X2, 1.359406 for m, 0.057839 for <sup>β</sup>, 0.805567 for <sup>θ</sup>1, 0.233550 for <sup>θ</sup>2, and 6.88×10−<sup>11</sup> for <sup>θ</sup>3. The results of all the models for these data are compared in Table 6.

LMM-L, NLMM, NLMM-L, and ANLMM-L results were calculated in previous studies [2,5,9,42]. The results of LMM and new Muskingum flood routing model were calculated by applying SAVCA. Table 6 displays a notable difference between the results of the LMMs and NLMMs; the difference was large because the error in the flood data of River Wye from December 1960 was large. The new Muskingum flood routing model results were better than those of other Muskingum flood routing models because the new Muskingum flood routing model showed the smallest error from 102 to 198 h including the peak value and showed the smallest error in the overall result.

Among other Muskingum flood routing models, ANLMM-L showed the best results (smallest error). The difference in SSQ between the ANLMM-L and new Muskingum flood routing model was 1677.99 (m3/s)2. The greatest difference between the models occurred at 138–186 h. Among other Muskingum flood routing models, the outflow from ANLMM-L was close to the observed outflow at 180 and 186 h. The difference between the observed outflow and the outflow for the new Muskingum flood routing model was smaller than that for other Muskingum flood routing models.


**Table 5.** Results when using Wang's flood data.


**Table 6.** Results when using flood data of River Wye December in 1960.

The fourth flood dataset used was Sutculer flood data which is a flood data with a double-peak. The parameters of the LMM for Sutculer flood data were determined to be 1.0 for K, and −0.006097 for X1. The parameters of the LMM-L for Sutculer flood data were determined to be 1.0 for K, −0.025914 for X1, and −0.041042 for β. The parameters of the NLMM for Sutculer flood data were determined to be 1.0 for K, −0.053787 for X1, and 1.002498 for m. The parameters of the new Muskingum flood routing model for Sutculer flood data were determined to be 0.931599 for K, −0.092988 for X1, 0.009066 for X2, 1.000013 for m, −0.036144 for β, 0.817639 for θ1, 0.214801 for θ2, and 0.745272 for θ3. The results of all models are shown in Table 7.

Among the results in Table 7, those of NLMM-L and ANLMM-L were calculated in previous studies [2,9]. The results of LMM, LMM-L, NLMM, and new Muskingum flood routing model were calculated by applying the SAVCA. Notably, the models that consider the lateral inflow show good results. In addition, the results of LMM-L, NLMM-L, and ANLMM-L were better than those of LMM and NLMM because LMM-L, NLMM-L, and ANLMM-L showed relatively small errors in the overall results. Moreover, the difference between the results of the new Muskingum flood routing model and other Muskingum flood routing models was substantial.

The ANLMM-L showed the best results (smallest error) among the considered models. The difference in SSQ between the ANLMM-L and new Muskingum flood routing model was 63.22 (m3/s)2. The time required to show the difference between the two results ranged from 1 to 3 h. At 1 h, the outflow of most Muskingum flood routing models was close to the observed outflow. At 2 and 3 h, except for the new Muskingum flood routing model, the outflow of most Muskingum flood routing models showed a difference from the observed outflow.

The last flood dataset analyzed was the flood data of River Wyre October in 1982. The parameters of the LMM for the flood data of River Wyre from October 1982 were determined to be 3.950351 for K and 0.295668 for X1. The parameters of the NLMM for the flood data of River Wyre from October 1982 were determined to be 8.248204 for K, 0.284338 for X1, and 0.812821 for m. The parameters of the new Muskingum flood routing model for the flood data of River Wyre from October 1982 were determined to be 0.931599 for K, −0.092988 for X1, 0.009066 for X2, 1.000013 for m, −0.036144 for β, 0.817639 for θ1, 0.214801 for θ2, and 0.745272 for θ3; the results, and their comparison with those of the other models, are shown in Table 8.

Of the results given in Table 8, the LMM-L, NLMM-L, and ANLMM-L results were calculated in previous studies [2,5,9], while those of the LMM and NLMM were calculated by applying SAVCA; the results of the latter two models were the same. The errors were the greatest for the results of both LMM and NLMM, and some of the calculated outflow values obtained were negative. The new Muskingum flood routing model results in Table 8 were also calculated by applying SAVCA.

Notably, the results obtained for the models considering the lateral inflow were good; those of the LMM-L, NLMM-L, and ANLMM-L were significantly better than those of the LMM and NLMM. The difference in the results of all Muskingum flood routing models occurs from 26 to 31 h. The outflows of NLMM-L, ANLMM-L, and new Muskingum flood routing model were the closest to the observed outflow from 26 to 31 h. Although the results of the NLMM-L and ANLMM-L were similar to the observed outflow at 26 h, these differed over time. However, the new Muskingum flood routing model results did not differ significantly from the observed outflow from 26 to 31 h. ANLMM-L showed the best results (smallest error) among the existing Muskingum flood routing models. The difference in SSQ between the ANLMM-L and new Muskingum flood routing model was 1.35 (m3/s)2. Although varying results were obtained when using different flood data for the various models, the overall results of new Muskingum flood routing model were better than those of other Muskingum flood routing models because the error of new Muskingum flood routing model was relatively small in the latter part from 16 to 31 h.


**Table 7.** Results when using Sutculer flood data.


**Table 8.** Results when using flood data of River Wyre October in 1982.

The new flood data in Daechung were applied to calibrate and validate the new Muskingum flood routing model. The flood data in April, 2014 were used for calibration and the flood data in April, 2018 were used for validation. Among various Muskingum

flood routing models, LMM-L considering lateral inflow to LMM, NLMM considering nonlinearity to LMM and new Muskingum flood routing model were applied to Daechung flood data and compared.

The parameters of LMM were 3.989981 for K and −0.034950 for X. The parameters of LMM-L were 3.865970 for K, −0.043293 for X, and −0.020080 for β. The parameters of NLMM were 2.225924 for K, −1.5 for X, and 1.0 for m. The parameters of NLMM-L were 2.225076 for K, −1.406061 for X, 1.000000 for m, −0.007072 for β, 1.000000 for θ. The parameters of ANLMM-L were 2.123526 for K, −1.500000 for X, 1.000000 for m, −0.010223 for β, 1.000000 for θ1, and 0.017871 for θ2. The parameters of new Muskingum flood routing model were 2.220910 for K, −1.498329 for X1, 0.094832 for X2, 1.000008 for m, −0.012616 for β, 0.999660 for θ1, 0.000093 for θ2, and 0.000288 for θ3. Each parameter was applied equally for 2014 and 2018 data. A total of 100 simulations were conducted for each Muskingum flood routing model, yielding the best results. Table 9 shows the results of Daechung flood data in 2014.

Based on the variable values determined from the 2014 flood data, it was applied to the 2018 flood data. Table 10 shows the results of Daechung flood data in 2018.

Figure 3 shows the results of calibration and validation in Daechung flood data.

**Figure 3.** Results of calibration and validation in Daechung flood data (**a**) 2014 flood data; (**b**) 2018 flood data..

In the 2014 flood data, the SSQs of LMM, LMM-L, NLMM, NLMM-L, ANLMM-L and new Muskingum flood routing model were 88.23, 73.81, 43.79, 42.32, 40.16, and 39.55, respectively. In the 2014 flood data, the error of new Muskingum flood routing model was relatively small and the calibration of new Muskingum flood routing model was relatively accurate. The SSQs of new Muskingum flood routing model in the 2018 flood data was relatively small. In the 2018 flood data, the SSQs of LMM, LMM-L, NLMM, NLMM-L, ANLMM-L and new Muskingum flood routing model were 221.92, 180.41, 171.45, 161.99, 159.84, and 157.64, respectively. In the results of Figure 2, more accurate flood routing was performed by applying the new Muskingum flood routing model compared to LMM, LMM-L, NLMM, NLMM-L and ANLMM-L.


**Table 9.** Results of Daechung flood data in 2014.


**Table 10.** Results of Daechung flood data in 2018.

#### **4. Discussion**

Because the calculation process differs for each Muskingum flood routing model, the time required to find the parameters when applying a meta-heuristic optimization algorithm is different for each method. The time required to apply SAVCA was summarized to determine the parameters of each Muskingum flood routing model for Wilson's flood data. The parameters of SAVCA were set at a constant, and the simulation was conducted 10 times. In addition, the number of iterations was set to 100,000. Table 11 shows the time taken for SAVCA when using Wilson's flood data.

**Table 11.** Time required by Muskingum flood routing models when using Wilson's flood data.


Depending on the parameters of the SAVCA and flood data, the results over time using the Muskingum flood routing models differ from the results in Table 9. As the number of parameters of the Muskingum flood routing models increased, the time required also increased. The time required by LMMs, LMM and LMM-L, was the shortest. NLMM and NLMM-L required a greater amount of time compared to the LMMs. The new Muskingum flood routing model required more time than the ANLMM-L, which in turn required more time than the NLMM-L. The time required by the new Muskingum flood routing model was approximately 1.5-times that required by LMM and LMM-L. In conclusion, the new Muskingum flood routing model produced more accurate results but took more time owing to the greater number of parameters and calculations.

An analysis was conducted on how each parameter of new Muskingum flood routing model affects the results. Daechung flood data in April, 2014 was applied to analyze the sensitivity of each parameter in the new Muskingum flood routing model. Figure 4 showed the results of the sensitivity analysis for the parameters in the new Muskingum flood routing model.

**Figure 4.** Results of sensitivity analysis for the parameters in the new Muskingum flood routing model.

SSQ decreases and then increases as parameter K increases, and SSQ increases as parameter X1 increases. SSQ increases as parameter X2 increases and SSQ increases as parameter m increases. SSQ decreases and then increases as parameter β increases, and SSQ decreases as the parameter θ<sup>1</sup> increases. SSQ increases as parameter θ<sup>2</sup> increases, and SSQ increases as parameter θ<sup>3</sup> increases. As each parameter changed, the change of the results was not constant. It is difficult to find a uniform pattern in all the results. However, what can be confirmed from the results of sensitivity analysis is that a meta-heuristic optimization algorithm such as SAVCA is required to produce results with a low SSQ.

#### **5. Conclusions**

The Muskingum flood routing model is a representative hydrologic flood routing model that is widely used owing to its easy applicability. The proposed Muskingum flood routing model in this study is a simple model that can be applied by researchers that use the existing Muskingum models for accurate flood routing.

In this study, the new Muskingum flood routing model was applied to various flood data, and the results obtained were compared with those of previously developed Muskingum flood routing models. As an index for comparison, the error was calculated between the observed and simulated outflows using the SSQ, NSE and RMSE. In addition, SAVCA, a meta-heuristic optimization algorithm, was applied to adjust the parameters of the new Muskingum flood routing model.

In the sensitivity analysis, the changes of the eight parameters in the new Muskingum flood routing model are different. There are parameters whose results are improved as the value (θ1) increases, some parameters (m, θ2, θ3, x1, x2) whose results are improved as the value decreases, and some parameters (K, β) whose results are changed (improved and then deteriorated) as the value increases. The eight parameters of the new Muskingum flood routing model are decision variables of SAVCA and are calculated through the optimization process.

Muskingum flood routing models considering the lateral inflow are capable of relatively sophisticated simulations, which corroborates that the influence of lateral inflow on the outflow can be considered. Among the existing models, the ANLMM-L showed the highest accuracy, although the difference between its results and those of the other Muskingum flood routing models was insignificant. In the new Muskingum flood routing model, the improved calculation method of the inflow at previous time and next time reflected the trend of the observed outflow.

Since the Muskingum flood routing model proposed in this study has eight parameters, the calculation process is more complicated than the existing Muskingum flood routing models. Accurate flood prediction is possible due to the complicated calculation process, but the calculation time is long.

Many studies, including this study, have been performed by applying various metaheuristic optimization algorithms to the Muskingum flood routing models. Since various optimization algorithms cannot show advantages in all problems, the results can be improved by appropriately selecting the operators of each meta-heuristic optimization algorithm. In addition, deep learning techniques have been widely used to apply various flood prediction methods including Muskingum flood routing models. By replacing the optimizer in deep learning techniques with meta-heuristic optimization algorithms, it would be possible to produce improved results in flood prediction.

**Funding:** This research was funded by the National Research Foundation (NRF) of Korea (NRF-2019R1I1A3A01059929).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** This work was supported by a grant from The National Research Foundation (NRF) of Korea (NRF-2019R1I1A3A01059929).

**Conflicts of Interest:** The author declares no conflict of interest.

#### **Abbreviations**


#### **References**

