## Abstract

This research introduces a novel nonlinear Muskingum model for river flood routing, aiming to enhance accuracy in modeling. It integrates lateral inflows using the Whale Optimization Algorithm (WOA) and employs a distributed Muskingum model, dividing river reaches into smaller intervals for precise calculations. The primary goal is to minimize the Sum of Square Errors (SSE) between the observed and modeled outflows. Our methodology is applied to six distinct flood hydrographs, revealing its versatility and efficacy. For Lawler's and Dinavar's flood data, the single-reach Muskingum model outperforms multi-reach versions, demonstrating its effectiveness in handling lateral inflows. For Lawler's data, the single-reach model (NR = 1) yields optimal parameters of *K* = 0.392, *x* = 0.027, *m* = 1.511, and *β* = 0.010, delivering superior results. Conversely, when fitting flood data from Wilson, Wye, Linsley, and Viessman and Lewis, the multi-reach Muskingum model exhibits better overall performance. Remarkably, the model excels with the Viessman and Lewis flood data, especially with two reaches (NR = 2), achieving a 21.6% SSE improvement while employing the same parameter set. This research represents a significant advancement in flood modeling, offering heightened accuracy and adaptability in river flood routing.

## HIGHLIGHTS

Flooding can have a significant impact on people and environment.

Distributed Muskingum model was introduced to improve the accuracy and efficiency of the model's calculations.

This study suggests WOA to improve the local optimum.

Considering lateral flow is the major feature developed for the new approach in our paper.

The results will be used in Red River of the North to anticipate floods in ND and MN in the US.

## INTRODUCTION

*K*, and weight coefficient of discharge,

*x*. However, over the years, researchers and hydrologists have made several modifications and enhancements to improve the accuracy and applicability of the Muskingum model. Some of these modifications include:

**Incorporate lateral flow:** During the actual flood events, lateral inflows were seen along the river reach under inquiry. The idea of considering linear lateral inflows concerning inflows in each time interval in the linear Muskingum model of O'Donnell was extended to Gill's nonlinear model (Karahan *et al.* 2015). Karahan *et al.* (2015) created a Cuckoo Search Algorithm (CSA) for parameter calibration and verification in the modified nonlinear Muskingum model. The improved model accounted for lateral inflows in an actual flood event (Karahan *et al.* 2015). Modeling flood events with a major contribution from lateral inflows have not been studied thoroughly to establish the prediction accuracy of the Non-Linear Muskingum model (NLMM). To fulfill this desire, the Non-Linear Muskingum model with the lateral flow (NLMM-L) is presented in the current study.

**Develop a nonlinear distributed Muskingum model:** To realistically simulate the nonlinear processes of flood movements in rivers, models such as the Muskingum model are modified to account for the nonlinearity of flow movement processes (Barati 2013; Karahan *et al.* 2013; Niazkar & Afzali 2015). Gill (1978) first introduced a nonlinear storage equation using the exponent of the Muskingum storage equation as the third parameter (Gill 1978). In this study, a distributed nonlinear Muskingum model that incorporates lateral inflows has been developed. This model represents a cascade of nonlinear Muskingum reaches, allowing for the flexibility of using one, two, three, or more sub-reaches to characterize flow behavior. Calibration of a single set of hydrological model parameters (*K*, *x*, and *m*) is required for all nonlinear routing calculations (Farzin *et al.* 2018). In a recent study, Atashi *et al.* (2023) introduced a novel flood routing method using a spatially variable exponent Muskingum model and the sine cosine optimization algorithm. This method demonstrates superior accuracy and efficiency compared to traditional Muskingum models and holds potential for enhancing hydraulic structure design and flood mitigation strategies (Atashi *et al.* 2023).

**Explicit numerical solution methods:** The calculation of the storage value included another modification to the flood routing process. The computation of new storage value has historically been carried out using the Euler's method (Tung 1985), however alternative explicit time integration methods, including the modified Euler's method (Vatankhah 2014), the fourth order Runge-Kutta method (Vatankhah 2014; Wang *et al.* 2014), and the Runge-Kutta-Fehlberg method (Vatankhah 2014) have also been employed. Vatankhah revealed that Euler's and improved Euler's methods are not accurate solution methods for any sized time intervals (Vatankhah 2014). The fourth order Runge-Kutta method and the Runge-Kutta-Fehlberg method, on the other hand, yield identical results and acceptable values in terms of accuracy evaluation criteria.

**Optimization using the Whale Optimization Algorithm (WOA):** In this research, the authors were motivated to adopt the WOA for the parameter optimization of the nonlinear Muskingum flood routing model. Mirjalili & Lewis (2016) developed the WOA to meta-heuristic algorithm (Diop *et al.* 2020; Malla *et al.* 2022). Guo *et al.* (2020) applied the WOA to forecast water consumption from 2004 to 2016 in Shaanxi Province of China. Their results indicated that the proposed algorithm for solving the three water resources forecasting model is better in comparison to other algorithms (Guo *et al.* 2020). Ezzeldin & Djebedjian (2020) used this algorithm for the least cost design of pipe networks. They find that the least cost network design of the WOA when compared to most of the optimization techniques available in the literature (Ezzeldin & Djebedjian 2020). Mohammadi & Mehdizadeh (2020) utilized the WOA for modeling daily reference evapotranspiration and found promising results (Mohammadi & Mehdizadeh 2020).

The parameters of the nonlinear distributed Muskingum models (*K, x, m,* and *β*) are estimated in this study using the WOA based on performance evaluation criteria (PEC). These case studies are Wilson's flood case study, Wye River flood case study, Lawler's flood case study, Linsley's flood case study, Viessman and Lewis's flood case study, and Dinavar river flood case study. As follows, the objective of this research is to develop a nonlinear Muskingum model for six well-known flood datasets: (1) Developing a nonlinear Muskingum model for river discharge estimation under lateral inflow conditions, (2) Estimation of the parameters of the nonlinear Muskingum models by using the WOA, and (3) Applying the distributed Muskingum to improve the accuracy of the procedure by splitting a reach into numerous periods, with the Muskingum model calculations conducted individually for each interval (Figure 1).

## MATERIALS AND METHODS

### Muskingum model

*K*and

*x*are stated, and cross-section details are not necessary. For each end of the river or channel reach, a minimum of two sections of the Muskingum Routing are needed. The equations used in the Muskingum Routing are the continuity equation (Das 2004):and the storage relationship:where

*I*indicates the inflow to the reach (m

^{3}/s),

*Q*indicates the outflow from the reach (m

^{3}/s),

*S*indicates the storage in the reach (m

^{3}),

*t*indicates the time (s),

*x*indicates the weighting coefficient, and

*K*indicates the storage constant (s).

*x*ranges between 0 and 0.5 for reservoir storage and between 0 and 0.3 for stream channels (Mohan 1997). It is notable that other hydrologic parameters (i.e.,

*K*and

*m*) have no general constraint (Barati 2013). Combining both equations an explicit equation can be obtained to calculate the outflow at the next time step:

*I*and

*O*represent the values at time

*t*

_{1}and

*t*

_{2}, respectively. According to Yoon & Padmanabhan (1993), when the connection between and

*S*is not linear, a nonlinear model may be more appropriate (Yoon & Padmanabhan 1993). Previous research has advocated a nonlinear Muskingum model for accounting for nonlinearity which is presented in Equation (5) (Gill 1978; Singh & Scarlatos 1987; Mohan 1997; Luo & Xie 2010; Easa

*et al.*2014):where

*m*takes the nonlinearity without lateral flow into the models. These models feature an extra parameter

*m*(=exponent power), which may be calculated using various parameter estimation approaches. In nonlinear models, however, unlike the linear model,

*K*with dimension does not describe the flood wave's travel time. Furthermore,

*x*does not have to be the same as it was in the linear model.

### Distributed nonlinear Muskingum model incorporating lateral flows

*β*parameter.

*t*is the index of time between zero and the ending time of the flood.

*j*is the spatial index between two and NR

*+*

*1*. The flowchart of the distributed nonlinear Muskingum model using the fourth-order Runge-Kutta method steps is shown in Figure 2. It should be noted that the NR can be determined by a trial-and-error approach. In other words, the model could be calculated for one, two, three, or more reaches and the best results between different NR could be selected by comparing objective function value and other PEC.

*Δt*, and an estimated slope. The slope will be a weighted average of the following slopes using the fourth-order Runge-Kutta method:

- 1.

### Whale optimization algorithm

The WOA is a metaheuristic optimization algorithm based on humpback whales' behavior (Mirjalili & Lewis 2016). The algorithm mimics the foraging behavior of humpback whales, where the whales search for food by diving and surfacing. The WOA uses this behavior to search for the global optimum of a given function. It is a population-based algorithm that uses a combination of stochastic and deterministic search techniques to explore the search space and find the global optimum. This study suggests this algorithm to update the individuals trapped in the local optimum to apply to the parameter estimation of the Muskingum model (Wang *et al.* 2013).

*N D*-dimensions real-valued parameter vectors in the population. In a generation, the position of the individual is the vector , and the global optimal location in a population so far is expressed as . Also, the main mathematical equation in WOA is as follows:where

*A*

*=*

*2a.r–a*,

*C*

*=*

*2.r*,

*a*linearly decreases from 2 to 0 over the course of a generation, and

*r*is the random vector between 0 and 1, is a randomly selected whale position vector. When is more than 1, the algorithm randomly selects a search agent to update the location of other whales according to the randomly selected whale position, to find more suitable prey, which can strengthen the exploration ability of the algorithm for global search. Moreover, and indicates the distance of

*i*individual,

*b*is a constant for defining the shape of the logarithmic spiral, and

*L*is a random number between −1 and 1.

The WOA parameters used in the present study are given as follows, WOA Source Codes (2016):

where and are the iteration number and a maximum number of iterations. The values of *a*, _{,} and are [2, 0], [−1, 1], and [−1, −2], respectively, when *iter* increases from zero to *iter _{max}*.

### Statistical PEC

To compare between the results of different methods, several PEC were proposed (Yoon & Padmanabhan 1993; Mohan 1997; McCuen *et al.* 2006; Fuat Toprak & Savci 2007; Toprak & Cigizoglu 2008; Luo & Xie 2010; Hosseini *et al.* 2016; Alizadeh *et al.* 2017; Kazemi & Barati 2022). Recent research by Wang *et al.* provides an overview of six accuracy evaluation criteria and nine research case datasets while also discussing current challenges and emerging trends in Muskingum model research (Wang *et al.* 2023). The following evaluation criteria are adopted for verifying the efficacy of the Excel solver and the other parameter estimation procedures.

^{6}T

^{−2}] between computed and observed outflows as follows:where and , respectively, are the observed and calculated outflow rates at each time step, and

*N*is the number of data.

To evaluate and compare the performances of the models, three PEC are used, as follows:

^{3}T

^{−1}] and deviation of peak time (DPOT) [L

^{3}T

^{−1}] is considered as the accuracy of the amount of peak outflow (Yoon & Padmanabhan 1993). DPO is given as:

SAD and SSE are measures of the bias and accuracy of the techniques, respectively (Mohan 1997). The performance criteria of the nonlinear Muskingum model for various parameter estimation methods and the hydrologic parameter values are presented in Tables 1–6. Because SSE is the objective function, it is the premier measure in the mentioned tables. Other measures derived from the data may perform differently.

Flood case study . | Period of flood . | NR in this study . |
---|---|---|

Wilson | 126 h | 1, 2, 3, 4 |

River Wye | 196 h | 1, 2, 3, 4 |

Lawler | 27 days | 1, 2, 3 |

Linsley | 10 h | 1, 2, 3 |

Viessman, Lewis | 24 h | 1, 2, 3, 4 |

Dinavar | 103 h | 1, 2, 3 |

Flood case study . | Period of flood . | NR in this study . |
---|---|---|

Wilson | 126 h | 1, 2, 3, 4 |

River Wye | 196 h | 1, 2, 3, 4 |

Lawler | 27 days | 1, 2, 3 |

Linsley | 10 h | 1, 2, 3 |

Viessman, Lewis | 24 h | 1, 2, 3, 4 |

Dinavar | 103 h | 1, 2, 3 |

Parameters . | Number of reaches . | ||||
---|---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | 5 . | |

x | 0.294 | 0.173 | 0.043 | 4.76 × 10^{−11} | 3.24 × 10^{−09} |

K (h) | 0.053 | 0.752 | 0.865 | 0.751 | 0.599 |

m | 2.381 | 1.601 | 1.478 | 1.444 | 1.443 |

β | − 0.022 | − 0.012 | − 0.008 | − 0.008 | − 0.008 |

SSE (cm^{2}) | 35.466 | 8.601 | 7.013 | 20.680 | 52.301 |

SAD (cm) | 23.921 | 11.330 | 9.142 | 18.184 | 29.183 |

DPO (cm) | 1.418 | 0.576 | 0.467 | 0.532 | 1.39 |

DPOT | 0 | 0 | 0 | 0 | 1 |

Parameters . | Number of reaches . | ||||
---|---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | 5 . | |

x | 0.294 | 0.173 | 0.043 | 4.76 × 10^{−11} | 3.24 × 10^{−09} |

K (h) | 0.053 | 0.752 | 0.865 | 0.751 | 0.599 |

m | 2.381 | 1.601 | 1.478 | 1.444 | 1.443 |

β | − 0.022 | − 0.012 | − 0.008 | − 0.008 | − 0.008 |

SSE (cm^{2}) | 35.466 | 8.601 | 7.013 | 20.680 | 52.301 |

SAD (cm) | 23.921 | 11.330 | 9.142 | 18.184 | 29.183 |

DPO (cm) | 1.418 | 0.576 | 0.467 | 0.532 | 1.39 |

DPOT | 0 | 0 | 0 | 0 | 1 |

Bold value of SSE indicates optimum Number of Sub-reaches (NR).

Parameters . | Number of reaches . | ||||
---|---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | 5 . | |

x | 0.30 | 0.33 | 0.50 | 0.30 | 0.11 |

K (h) | 0.01 | 0.46 | 1.55 | 1.29 | 1.13 |

m | 2.17 | 1.45 | 1.20 | 1.19 | 1.18 |

β | 0.06 | 0.02 | 0.00 | 0.01 | 0.01 |

SSE × 10^{3} | 84.00 | 31.91 | 27.86 | 43.03 | 54.89 |

SAD (cm) | 893.17 | 724.20 | 767.46 | 828.04 | 903.63 |

DPO (cm) | 141.95 | 88.30 | 55.81 | 130.10 | 162.16 |

DPOT | 0.30 | 0.33 | 0.50 | 0 | 0 |

Parameters . | Number of reaches . | ||||
---|---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | 5 . | |

x | 0.30 | 0.33 | 0.50 | 0.30 | 0.11 |

K (h) | 0.01 | 0.46 | 1.55 | 1.29 | 1.13 |

m | 2.17 | 1.45 | 1.20 | 1.19 | 1.18 |

β | 0.06 | 0.02 | 0.00 | 0.01 | 0.01 |

SSE × 10^{3} | 84.00 | 31.91 | 27.86 | 43.03 | 54.89 |

SAD (cm) | 893.17 | 724.20 | 767.46 | 828.04 | 903.63 |

DPO (cm) | 141.95 | 88.30 | 55.81 | 130.10 | 162.16 |

DPOT | 0.30 | 0.33 | 0.50 | 0 | 0 |

Bold value of SSE indicates optimum Number of Sub-reaches (NR).

PEC parameters . | Number of reaches . | ||
---|---|---|---|

1 . | 2 . | 3 . | |

x | 0.027 | 1.25 × 10^{−09} | 2.45 × 10^{−07} |

K (day) | 0.392 | 0.207 | 0.224 |

m | 1.511 | 1.485 | 1.327 |

β | 0.010 | 0.001 | − 0.001 |

SSE (cm^{2}) | 0.365 | 1.952 | 3.345 |

SAD (cm) | 2.333 | 5.885 | 7.977 |

DPO (cm) | 0.061 | 0.304 | 0.400 |

DPOT | 0 | 0 | 0 |

PEC parameters . | Number of reaches . | ||
---|---|---|---|

1 . | 2 . | 3 . | |

x | 0.027 | 1.25 × 10^{−09} | 2.45 × 10^{−07} |

K (day) | 0.392 | 0.207 | 0.224 |

m | 1.511 | 1.485 | 1.327 |

β | 0.010 | 0.001 | − 0.001 |

SSE (cm^{2}) | 0.365 | 1.952 | 3.345 |

SAD (cm) | 2.333 | 5.885 | 7.977 |

DPO (cm) | 0.061 | 0.304 | 0.400 |

DPOT | 0 | 0 | 0 |

Bold value of SSE indicates optimum Number of Sub-reaches (NR).

PEC parameters . | Number of reaches . | |||
---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | |

x | 0.225 | 0.149 | 2.24 × 10^{−08} | 1.80 × 10^{−08} |

K (h) | 0.527 | 0.344 | 0.262 | 0.227 |

m | 1.115 | 1.024 | 1.000 | 1 |

β | 0.018 | 0.010 | 0.013 | 0.016 |

SSE (cm^{2}) | 1.717 | 1.509 | 3.616 | 14.660 |

SAD (cm) | 3.666 | 3.347 | 6.375 | 14.657 |

DPO (cm) | 0.183 | 0.257 | 0.140 | 0.883 |

DPOT | 0 | 0 | 0 | 0 |

PEC parameters . | Number of reaches . | |||
---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | |

x | 0.225 | 0.149 | 2.24 × 10^{−08} | 1.80 × 10^{−08} |

K (h) | 0.527 | 0.344 | 0.262 | 0.227 |

m | 1.115 | 1.024 | 1.000 | 1 |

β | 0.018 | 0.010 | 0.013 | 0.016 |

SSE (cm^{2}) | 1.717 | 1.509 | 3.616 | 14.660 |

SAD (cm) | 3.666 | 3.347 | 6.375 | 14.657 |

DPO (cm) | 0.183 | 0.257 | 0.140 | 0.883 |

DPOT | 0 | 0 | 0 | 0 |

Bold value of SSE indicates optimum Number of Sub-reaches (NR).

PEC parameters . | Number of reaches . | |||
---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | |

x | 0.230 | 0.177 | 0.019 | 6.20 × 10^{−09} |

K (h) | 0.021 | 0.069 | 0.098 | 0.170 |

m | 1.611 | 1.345 | 1.248 | 1.146 |

β | 0.016 | 0.005 | 0.004 | 0.003 |

SSE × 10^{−3} | 71.62 | 52.04 | 61.10 | 102.13 |

SAD (cm) | 925.600 | 828.514 | 909.494 | 1,252.53 |

DPO (cm) | 71.547 | 25.470 | 49.911 | 80.872 |

DPOT | 0 | 0 | 0 | 0 |

PEC parameters . | Number of reaches . | |||
---|---|---|---|---|

1 . | 2 . | 3 . | 4 . | |

x | 0.230 | 0.177 | 0.019 | 6.20 × 10^{−09} |

K (h) | 0.021 | 0.069 | 0.098 | 0.170 |

m | 1.611 | 1.345 | 1.248 | 1.146 |

β | 0.016 | 0.005 | 0.004 | 0.003 |

SSE × 10^{−3} | 71.62 | 52.04 | 61.10 | 102.13 |

SAD (cm) | 925.600 | 828.514 | 909.494 | 1,252.53 |

DPO (cm) | 71.547 | 25.470 | 49.911 | 80.872 |

DPOT | 0 | 0 | 0 | 0 |

Bold value of SSE indicates optimum Number of Sub-reaches (NR).

### Application case studies

The performances of the proposed model are evaluated in six case studies making use of the inflow–outflow hydrograph. The first case study provided by Wilson (1990), a smooth single-peak hydrograph (Wilson 1990). The second case study is the River Wye, UK, with no tributaries linked from Erwood to Belmont; its total length is 69.75 km. The third case study is flood data for 27 days presented by Lawler (1964). The fourth case study presented a 10-day flood by Linsley *et al.* (1975). Records for the flood of March 15–31, 1936 were used. The time interval was 12 h. Local inflow was included by increasing the gaged tributary flow at E. Liverpool in proportion to the difference between inflow and outflow to the routing reach. A multiple-peak hydrograph studied by (Viessman *et al.* 2003) is demonstrated as the fifth case study. Also, the last case study is the data related to the Dinavar River in Iran, Kermanshah, which is presented in this study. Moradi *et al.* (2022) conducted a study focused on the Dinavar River in Iran, Kermanshah, where they proposed a new flood routing technique using a nonlinear Muskingum model and the artificial gorilla troop optimizer algorithm. This innovative approach demonstrated superior accuracy and efficiency compared to existing Muskingum models (Moradi *et al.* 2022). Information regarding the flood data examples, including their duration and the NR considered for the current study, is provided in Table 1. All flood data used in this study have been applied in several Muskingum flood routing models in existing studies.

## RESULTS AND DISCUSSIONS

The six cases involving the specific datasets are evaluated in this study and will be discussed in the following sections. The proposed flood routing model (NLMM-L) was solved for more than 100 runs using WOAs, to compare the performance of different numbers of lateral reaches.

**1.****A case study of the dataset given by Wilson**

*et al.*2013; Vatankhah 2014). For Wilson flood data, observed inflow and outflow hydrographs are displayed in Figure 3 for 126 h (Wilson 1990). The number of reaches for each run is indicated by NR. Table 2 illustrates the routing parameters associated with the Muskingum model and PEC for the nonlinear Muskingum model parameters discussed in the material and methods section. Using the current study algorithm, the optimized parameters for the Wilson flood data were determined to be

*K*= 0.865,

*x*= 0.043,

*m*= 1.478,

*β*= −0.008 for NR = 3. The largest outflow discharge, which occurred at the 60th hour of flood data, has been changed by −1.67, 1.01, 0.13, 1.18, and 0.96% for NR = 1 to 5, respectively, with NR = 3 displaying the lowest difference in peak discharge with the Wilson flood data. The findings show that SSE, SAD, and DPO for NR = 3 have the lowest values (SSE = 7.013, SAD = 90.142, and DPO = 0.467) compared to other numbers of reach. Vatankhah (2014) obtained SSE = 62.59 using the same storage equation and the fourth order Runge-Kutta method (Vatankhah 2014). Therefore, the proposed distributed hydrological Muskingum model incorporating lateral flow yielded 792% better results in terms of SSE by using only two more calibration parameters.

**2.****Case study of the dataset of the Wye River**

The second case study suggested the application of flood routing methods in a flood on the Wye River in England in 1960. The Wye River flood in England in 1960 was a severe flood that affected the town of Roydon, Essex, and the surrounding areas along the river Wye. It was caused by heavy rainfall and the river overflowing its banks, leading to widespread damage and disruption. Wye River flood event has a sharp peak that was caused by heavy rainfall (Karahan *et al.* 2015). The shape of hydrograph is important for basins with high slope in mountain areas. The Wye River has no tributaries and very little lateral flow during its 69.75 km journey from Erwood to Belmont. This flood is an excellent test case for routing algorithms (O'donnell 1985). The input value ranges from 54 to 1,145 cm, with the flow reaching its peak in 84 h. The time step in this case study is 6 h, the number of time steps is 33, and the total time is 198 h. The highest output was 969 cm in 102 h. This case study was utilized by several researchers to estimate the Muskingum parameters (Ayvaz & Gurarslan 2017; Akbari *et al.* 2020; Lee 2021; Lu *et al.* 2021). The reason for choosing these data is to enable the comparison between this model and the other methods to verify the improvement in flood routing owing to the new method for large outflows. Table 3 displays a case study of computation and optimization for a one-to-five-reach problem using the Wye River flood data.

The study model was applied to the flood data for River Wye in December 1960 which agreed well with the method proposed in this study. The optimal parameters of the model for the flood data of river Wye were determined to be 1.55 for *K*, 0.5 for *x*, 1.20 for *m*, and 0 for *β* using WOA associated with three lateral reaches. Karahan *et al.* (2013) find SSE = 37,944.14 using the same storage equation and traditional Euler's method. Consequently, the proposed flood routing model finds 36.2% better results in terms of SSE (Karahan *et al.* 2013).

**3.****A case study of the dataset given by Lawler**

Lawler flood data has a smooth peak which is caused by a rainfall with a constant intensity during a period. This flood is important for rivers with a high capacity of river bankfull. Lawler is a flood with 27 days of inflow and outflow discharge information (Lawler 1964). The maximum outflow discharge is 18.86 for the 13th day of the flood. In this investigation, NR = 1, 2, and 3 were used as the suggested method to observe the outcomes with the lateral flow and utilize the WOA. The findings show that NR = 1 has the lowest SSE value among the different numbers of reach. Single reach also has the lowest value for SAD and DPO, implying that it performs well across the board. Notably, the value Muskingum parameters for NR = 1 are *K* = 0.392, *x* = 0.027, *m* = 1.511, and *β* = 0.010. The maximum outflow discharge, which occurred at the 10th hour of flood data, was increased by 0.32, 1.61, and 2.11% for NR = 1 to 3, with NR = 1 showing the smallest variation in peak discharge with the Lawler flood data. For *NL* = 1 to 3, Table 4 shows the optimal outflows and intermediate results for flood routing.

**4.****A case study of the dataset given by Linsley**

Linsley flood events have a partially linear relationship between weighted-flow and storage volume with a low period of flood. Simulation of such flood is important for urban flood management. Linsley flood is a 10-h flood with 0.5-h time intervals. The maximum outflow is 19.22 cm which occurred in the fourth hour of the flood (Linsley *et al.* 1975). Table 5 shows the Muskingum parameters and PEC using Linsley flood data for NR = 1, 2, 3, and 4. For NR = 1, 2, 3, and 4, the changes in peak outflow are 0.95, 0.38, −1.66, and −5.49%, respectively. The results demonstrate that among the various numbers of reach, NR = 2 has the lowest SSE value. Two reaches perform well across the board since it also has the lowest value for SAD and DPO. SSE, SAD, and DPO for NR = 2 are 1.509, 3.347, and 0.257, respectively.

**5.****A case study of the dataset given by Viessman and Lewis**

This section applied the proposed algorithm for different numbers of reaches on the flood data presented by Viessman and Lewis with a double-peak and non-smooth flow hydrograph (Viessman *et al.* 2003). The Viessman and Lewis flood dataset has a double-peak and non-smooth flow hydrograph as a response to a unique rainfall pulse. Double-peak hydrographs are important for non-urban catchments, where the first and the delayed discharge peaks are explained by different runoff mechanisms (Martínez-Carreras *et al.* 2016). Viessman and Lewis flood data is also a widely used resource for the study of hydrology and the analysis of floodplains. This case study uses *t* = 1 day and 24-h time increments. On the 11th and 18th hours, the outflow hydrograph shows two peaks. The inflow varies between 114.33 and 1,768.29 cm. The greatest discharge rate is 1,509.3 cm. Table 6 demonstrates the optimal outflows and intermediate outcomes in flood routing for NR = 1 to 4. This data has very little lateral flow contribution. Following the execution of the suggested technique for this flow data, the convergence results shown in Table 6 for NR = 1 to 4 are achieved. It can be shown that the NR of 2 has the greatest match with the Viessman and Lewis flood data for both flow peaks that occurred on days 9 and 18 of flood data. The performance criteria of the nonlinear Muskingum model for various parameter estimation methods and the hydrologic parameter values are presented in Table 6. The SSE value also indicates that the NR = 2 has a highly favorable performance, and when compared to other studies, the SSE value for this case is quite appropriate. As shown, a significant improvement is observed in terms of SSE, SAD, and DPO as the PEC parameter values when the value of NR is 2. The results demonstrate that, in comparison to other numbers of reach, the SSE, SAD, and DPO for NR = 2 have the lowest values. Additionally, while employing the WOA technique, the Muskingum parameters are as follows: *x* = 0.177, *K* = 0.069, *m* = 1.345, and *β* = 0.005. The improved nonlinear Muskingum model with variable exponent parameter of Easa (2015) resulted in SSE = 63,300. Thus, the developed model of the present study returned 21.6% better results in terms of SSE with the same number of model parameters (Easa 2015).

**6.****A case study of the dataset given by the Dinavar River**

The sixth case study is a hydrograph of the inflow of the Dinavar River in Iran's Kermanshah region, which is reported in this paper. The Dinavar River flood data has significant lateral flow. This flood event is important for a river system with multiple inflows from several branches. This sample includes 120 time steps and *t* = 1 h. The lowest and highest inflows are 14.9 and 121 cm at 51 and 54 h of the flood, respectively. The hydrograph data are connected to the Mianrahan and Heidarabad hydrometric stations, where Kermanshah Regional Water Company took measurements. The river flow is recharged by groundwater flows via this channel. Table 7 shows the optimal outflows for flood routing based on minimizing the *PEC* parameters. When the model is single-reach, the best results are obtained. The Muskingum parameters for single reach using the study algorithm are *x* = 0.3, *K* = 5.051, *m* = 1.0, and *β* = 1.083. The maximum outflow discharge, which occurred at the 38th hour of the flood data, was reduced by 1.28, 1.39, and 1.62% for NR = 1 to 3, respectively, showing that a single reach had the smallest variation in peak discharge with the Dinavar flood data.

PEC parameters . | Number of reaches . | ||
---|---|---|---|

1 . | 2 . | 3 . | |

x | 0.30 | 0.20 | 0.09 |

K (h) | 5.05 | 2.41 | 1.59 |

m | 1.00 | 1.00 | 1.00 |

β | 1.08 | 0.44 | 0.28 |

SSE × 10^{−3} | 5.36 | 5.73 | 5.84 |

SAD (cm) | 561.28 | 568.01 | 571.80 |

DPO (cm) | 2.19 | 0.82 | 0.97 |

DPOT | 1.00 | 1.00 | 1.00 |

PEC parameters . | Number of reaches . | ||
---|---|---|---|

1 . | 2 . | 3 . | |

x | 0.30 | 0.20 | 0.09 |

K (h) | 5.05 | 2.41 | 1.59 |

m | 1.00 | 1.00 | 1.00 |

β | 1.08 | 0.44 | 0.28 |

SSE × 10^{−3} | 5.36 | 5.73 | 5.84 |

SAD (cm) | 561.28 | 568.01 | 571.80 |

DPO (cm) | 2.19 | 0.82 | 0.97 |

DPOT | 1.00 | 1.00 | 1.00 |

Bold value of SSE indicates optimum Number of Sub-reaches (NR).

Figure 3 compares the convergence plots for the model for the optimized model for each flood data. Black and triangle lines indicate inflow and outflow hydrographs, respectively, based on observed data in all six data floods. The green hydrograph is the optimized routed hydrograph resulting from the model. The figures are (a) Wilson flood; (b) Wye River flood; (c) Lawler flood; (d) Linsley flood; (e) Viessman and Lewis flood; and (f) Dinavar River flood. The calibrated and observed output hydrographs in Figure 3 exhibit high coherence, demonstrating that the given model in this study worked effectively to describe the flood data based on these six samples.

## CONCLUSION

Flood routing models, as is widely known, are used to anticipate the changing magnitude, speed, and form of a flood wave as a function of time at one or more places. Hydrological river routing models are the critical components of the software that is used to mimic the whole hydrologic processes of dendritic watershed systems. The conventional linear Muskingum model, for example, has been utilized in HEC-HMS (Hydrologic Modeling System) and SWAT (Soil and Water Assessment Tool). The suggested distributed hydrological Muskingum model may be used in such software to reduce uncertainty in flood simulation.

Moreover, the Muskingum method is a simplification of the complex hydrological processes that occur in a river system, and it does not consider lateral inflow by default. However, if the lateral inflow is significant and has a noticeable impact on the overall behavior of the system, including it in the model can lead to an improved accuracy. In this paper, a novel method for nonlinear Muskingum flood routing models is provided. The suggested method involves segmenting each river reach into sub-reaches, with the Muskingum process being conducted separately for each sub-reach to increase the precision of conventional or single-reach Muskingum model results. The suggested method is then included in an optimization model that makes use of the WOA. By minimizing the differences between the observed and simulated outflows of the flood event, the WOA-based optimization model aims to establish the dimensions and related model parameters within each reach.

The results also show that when the lateral inflow is implemented in the Muskingum model for Lawler and Dinavar's flood data, the single-reach Muskingum model has significantly more accuracy than the multi-reach Muskingum model, whereas the multi-reach Muskingum model has better results for the flood data cited by Wilson, Wye, Linsley, and Viessman and Lewis. For example, in the case of the Lawler flood data, the results indicate that NR = 1 outperforms other reach numbers, with parameters *K* = 0.392, *x* = 0.027, *m* = 1.511, and *β* = 0.010. But in the Wilson flood data case, the optimized Muskingum model parameters were determined with NR = 3: *K* = 0.865, *x* = 0.043, *m* = 1.478, and *β* = −0.008. Notably, PEC for NR = 3 showed favorable values: SSE = 7.013 cm^{2}, SAD = 90.142 cm, and DPO = 0.467 cm. These results demonstrate that the proposed distributed hydrological Muskingum model, incorporating lateral flow, yielded significantly improved results compared to existing models.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.