Development of a distributed nonlinear Muskingum model by considering snowmelt effects for flood routing in the Red River

This research paper presents the development of a nonlinear Muskingum model which achieves precise flood routing through river reaches while considering lateral inflow conditions. Fourteen pairs of flood hydrograph found at two specific United States Geological Survey (USGS) stations located along the Red River of the North, namely Grand Forks and Drayton, are used for the calibrations and validations of the Muskingum model. To enhance the accuracy of the procedure, a reach is divided into multiple sub-reaches, and the Muskingum model calculations are performed individually for each interval using the distributed Muskingum method. Notably, the model development process incorporates the use of the Salp Swarm algorithm. The obtained results demonstrate the effectiveness of the developed nonlinear Muskingum model in accurately routing floods through the very gentle river with a bed slope of (0.0002–0.0003). The events were categorized into three groups based on their dominant drivers: Group A (Snowmelt-driven floods), Group B (Rain-on-snow-induced floods), and Group C (Mixed floods influenced by both snowmelt and rainfall). For the sub-reaches in Group A, single sub-reach (NR = 1), the Performance Evaluation Criteria (PEC) yielded the highest value for SSE, amounting to 404.9 × 106. In Group B, when NR = 2, PEC results the highest value were SSE = 730.2 × 106. The number of sub-reaches in a model has a significant influence on parameter estimates and model performance, as demonstrated by the analysis of hydrologic parameters and performance evaluation criteria. Optimal performance varied across case studies, emphasizing the importance of selecting the appropriate number of sub-reaches for peak discharge predictions.

The travel time (K) is the time required for the water to travel through the reach, which is dependent on the channel geometry, roughness, and other hydraulic characteristics.The reach weighting factor (x) is the proportion of the discharge that enters the reach from the upstream section, which is also known as the weighting coefficient.These parameters can be determined using various techniques, including trial and error, optimization algorithms, and regression analysis.The Muskingum model can be represented mathematically as follows 4 : where O is the discharge at the downstream end of the reach (m 3 /s), and I is the discharge at the upstream end of the reach (m 3 /s).x is the weighting factor for the reach (ranges between 0 and 0.5 for reservoir storage and between 0 and 0.3 for stream channels 5 ), K is the travel time for the reach(s), and S is the storage volumes of the reach (m 3 ).By combining Eq. (1) with the continuity equation an explicit equation can be obtained to calculate the outflow at the next time step: The subscripts 1 and 2 on I and O represent the values at time t 1 and t 2 respectively.C 0 , C 1, and C 2 are the coefficients.
The traditional linear Muskingum model seeks a method of parameter estimation to determine the values of K and x.However, the linear Muskingum model leads to considerable inaccuracy in the forecast of flood behavior throughout its propagation along a river because natural channel reaches often have a nonlinear storagedischarge connection.To address this limitation, models such as the Muskingum model have been modified to account for the nonlinearity of flow movement processes.Gill 6 introduced a nonlinear storage equation using the exponent of the Muskingum storage equation as the third parameter, and later models such as the Nonlinear Muskingum model (NLMM) have been developed to include lateral inflows and better simulate the nonlinear processes of flood movements in rivers.As stated by Perumal et al. 7 there exists no "truly physically based" flood routing model which does not require any calibration.Although the roughness factor is a property of the natural conditions, it is considered as a model tuning parameter in the routing process of Muskingum-Cunge models 8,9 .
Furthermore, as pointed out by Koussis, nonlinear routing models like the nonlinear Muskingum model possess an advantage in their ability to accurately replicate the rapid surge of a flood wave, a task that linear models often struggle with 10,11 .It should be noted that the Muskingum-Cunge is not a "linear model".However, it is a "time-variant linear model", which means that it is "locally linear" in time, but the overall behavior is nonlinear.Every flood routing model necessitates specific input parameters and data.In some river segments, all the required inputs are readily available, and the choice of a model can be based on personal expertise and computational capacity.When there are no constraints on these factors, the use of a dynamic wave model is the most suitable option.For the Muskingum-Cunge model, essential input parameters include initial conditions, upstream boundary conditions, Manning's roughness coefficient, length of the routing reach, river cross-sections, and the bed slope 12 , while nonlinear Muskingum model requires the initial condition, upstream boundary condition and the hydrologic parameters.One of the important motivations of the authors is to suggest alternative hydrological flood routing model for using in modeling software of hydrologic processes of watershed systems such as HEC-HMS (hydrologic modeling system) and SWAT (soil and water assessment tool).The proposed and rigorously validated distributed hydrological Muskingum model can serve as a valuable addition to hydrological software, effectively mitigating uncertainties in flood modeling.It's important to note that this study primarily emphasizes the nonlinear Muskingum routing models.

Nonlinear Muskingum model
Previous research has advocated a nonlinear Muskingum model for accounting nonlinearity, which allows for a better representation of the nonlinear relationship between the inflow at the upstream end and outflow at the downstream end of the river channel, which is presented in Eq. (3) 5,6,13-15 : where m takes the nonlinearity without lateral inflow to the models.These models feature an extra parameter m (= exponent power), which may be calculated using various parameter estimation approaches.On the other hand, K with dimension of L 3(1−m) T m in nonlinear models unlike the linear model does not describe the travel time of the flood wave.In addition, x does not have to be the same as in the linear model.Equation (4) shows a modified storage equation that considers lateral inflow 16 : where β is the parameter accounting for the lateral inflow.The storage at time t + 1 is shown in equation below.
By substituting Eq. ( 4) into Eq.( 5), with consideration of lateral inflow in a nonlinear relationship between the inflow at the upstream end and outflow at the downstream end, the storage at time t + 1 is represented in Eq. ( 6): The NLMM with the lateral inflow (NLMM-L) has been suggested as an accurate solution method for addressing the nonlinear Muskingum model [17][18][19][20][21][22][23][24][25][26] .The Muskingum model can be solved using various numerical (1) methods, with the fourth-order Runge-Kutta method standing out as one of the most accurate approaches.Here, the fourth order Runge-Kutta method has been offered as an accurate and acceptable solution method among the different explicit solution methods for addressing the nonlinear Muskingum model since it is simpler than the Runge-Kutta-Fehlberg method 27,28 .
It's important to note that Cunge 12 established a vital connection between the flood routing parameters within the Muskingum approach and the channel properties, as well as flow characteristics.This connection was achieved by utilizing an approximation error derived from a Taylor series expansion of grid specifications and employing a diffusion analogy.Consequently, Cunge introduced a model known as the Muskingum-Cunge model, which has served as a cornerstone for further research and refinement 7,[29][30][31][32][33][34][35] .
This particular class of flood routing models necessitates a set of crucial input parameters, including the initial condition, upstream boundary condition, Manning's roughness coefficient, length of the routing reach, cross-sections along the river reach, and the bed slope.Within these inputs, the Manning's roughness coefficient stands out as a notable source of uncertainty for Muskingum-Cunge models.This coefficient is closely linked to surface roughness, vegetation, channel irregularities, channel alignment, silting, scouring, obstruction, channel size, and shape, as well as the magnitudes of stages and discharges 36 .Most of these factors exhibit variations from one flood event to another within a given river reach 9 .As a result, the roughness variations within a river reach are inherently three-dimensional, making them challenging to model.Therefore, there's a need to strategically select a single parameter as the roughness coefficient through a calibration process to align the flood routing results, particularly when only a single set of inflow and corresponding outflow hydrographs are available for the considered river reach 7,8 .
In contrast, the traditional linear Muskingum model primarily relies on the initial condition, upstream boundary condition, and various hydrologic parameters 10,11 .One of the principal motivations of the authors is to propose an alternative hydrological flood routing model suitable for integration into modeling software for hydrologic processes within watershed systems like HEC-HMS (hydrologic modeling system) and SWAT (soil and water assessment tool).Consequently, the primary focus of this study revolves around enhancing nonlinear Muskingum routing models.
To further improve its accuracy and convergence, optimization algorithms like the Salp Swarm algorithm have emerged as effective tools.Salp Swarm algorithm is a population-based optimization algorithm inspired by the swarming behavior of salps.The algorithm starts with a population of salps, each of which has a random position in the search space.In each iteration, the salps move towards the leader salp, which is the salp with the best fitness 37 .The salps that have the best fitness values are more likely to be selected for reproduction, and their offspring are added to the population.This process continues until the algorithm converges on an optimal solution 38,39 .SSA has been shown to outperform other optimization algorithms in terms of both accuracy and convergence speed 37 and has the potential to be used to solve a wide variety of problems 40,41 .
Through the integration of the Salp Swarm algorithm and NLMM-L Muskingum method, the research seeks to address the challenges associated with snowmelt-induced flooding and provide more precise flood routing solutions for the study area.The objective of this research is to develop a nonlinear Muskingum model for the Red River between two USGS stream gauging stations, Grand Forks, and Drayton, in the US, using flood hydrographs caused by snowmelt in spring.This area was selected due to the recurrent flooding observed in the central part of the Red River and its adjacent floodplain regions, stretching between Grand Forks, ND, and Emerson, ND.These flood events are notably characterized by the substantial seasonal water area that consistently forms during wet spring periods 42 .The flat terrain and downstream ice jams in the Red River and Lake Winnipeg contribute to the frequent flooding during spring seasons in wet years, including 1997, 2009, 2011, and 2013.The repetitive nature of flood events in this section of the Red River underscores the importance of comprehending and effectively managing flood risks in the region.
The primary objective of this study is to develop a Muskingum model that can accurately estimate river discharge, considering lateral inflow conditions.The research methodology involves estimating the parameters (K, x, m, and β) of the nonlinear Muskingum models using a distributed flood routing model, utilizing the Salp Swarm algorithm.This approach divides the river reach into multiple intervals, allowing individual Muskingum model calculations for each interval, thereby improving the overall accuracy of the estimation process.The findings of this research are expected to contribute to the enhancement of flood forecasting and warning systems in the Red River basin, enabling better preparedness and response to flood events.

Study area
The Red River Basin, spanning both the United States and Canada, encompasses an area of 117,185 km 2 , with most of its expanse situated in North Dakota, South Dakota, and Minnesota 43 .Figure 1a shows the location of the basin.Characterized by a semi-arid climate, the region experiences cold winters and hot, dry summers, with the primary streamflow occurring during spring and early summer due to snowmelt or heavy rainfall 44 .The Red River itself is a prominent watercourse within the basin, flowing northward along the border between Minnesota and North Dakota, as well as the Canadian provinces of North Dakota and Manitoba.Notably, the river is prone to frequent flooding, owing to its gentle flow and flat topography, and its broad and shallow floodplain exacerbates the vulnerability to heavy rain or spring snowmelt, resulting in historically devastating floods [45][46][47][48][49] .
A recent study by Atashi et al. 50investigated various forecasting methods for water levels in flood warning systems.The study found that the Long Short-Term Memory (LSTM) method demonstrated superior accuracy and precision compared to classical statistical and machine learning approaches, making it a reliable choice for flood prediction, particularly for downstream stations lacking discharge information 50 .

Grouping of dataset
For conducting flood routing analysis, it is essential to have observed flow hydrographs at specific upstream and downstream cross-section pairs.For this study, we selected the existing USGS streamflow gauging stations at Drayton (Station No. 05092000) and Grand Forks (Station No. 05082500).These stations were chosen because they offer vital streamflow data necessary for hydrograph analysis in the region extending from Grand Forks to the US-Canada border.As mentioned, this area has experienced recurrent flooding, particularly in the central part of the Red River and the nearby floodplain regions, which extend from Grand Forks, ND, to Emerson, ND.The locations of the gauging stations are shown in Fig. 1b.
A total of fourteen flood events occurring between 1990 and 2022 were utilized to calibrate and validate the proposed model, with twelve events designated for calibration and two events for validation.These validation events specifically corresponded to the flood events in 2020 and 2022.
To form the groups, specific criteria were considered based on the distinct routing characteristics observed in different flood types.For instance, the 1997 flood primarily resulted from snowmelt with minimal rain-on-snow impact, while the 2022 flood predominantly comprised rain-on-snow conditions.The selection of these criteria was crucial to ensure accurate calibration results across all events.The summarized data for the fourteen flood occurrences between 1990 and 2022 is presented in Table 1, with the events categorized into Group A, Group B, and Group C. The criteria used for forming the groups include: 1. Snowmelt dominant events: flood events where the primary driver was snowmelt with minimal rain-onsnow impact were categorized into Group A. These events typically exhibit specific routing characteristics associated with snowmelt-dominated hydrological processes.2. Rain-on-snow dominant events: flood events predominantly characterized by rain-on-snow conditions were categorized into Group B. These events show distinct routing behavior resulting from the combined effects of rain and snowmelt on the hydrological system.3. Mixed events: flood events that had a combination of snowmelt and rain-on-snow conditions were categorized into Group C. These events exhibit routing characteristics influenced by both snowmelt and rainfall contributions.www.nature.com/scientificreports/By categorizing the flood events into these distinct groups, it was possible to calibrate the model separately for each flood type, considering the specific routing behaviors associated with each category.Group A and Group B exhibit distinct variations in terms of precipitation amounts.Specifically, Group A demonstrates a lower range of monthly precipitation, ranging from 57. 2 to 141.0 mm during the selected years.In contrast, Group B showcases a higher range of monthly precipitation, ranging from 187.3 to 276.4 mm, observed specifically in the years 1999, 2004, 2013, and 2022.

Distributed nonlinear Muskingum model incorporating lateral inflows
Nonlinear Muskingum models consist of a series of nonlinear Muskingum reaches, which are further subdivided into equal nonlinear Muskingum sub-reaches.The distributed nonlinear Muskingum model, as depicted in Fig. 2, provides an illustrative representation of this arrangement.Only one set of hydrological model parameters (K, x, and m) needs to be calibrated and used in the nonlinear routing calculations.The flood hydrograph is routed from the main inflow hydrograph at the upstream section to the downstream section of the first sub-reach.The outcome is treated as the inflow for the second sub-reach and is routed subsequently to the downstream section of the second sub-reach 52 .To get the flood hydrograph at the downstream section of the final sub-reach, this process is repeated sequentially.The number of sub-reaches (NR) can be determined by trying different options and selecting the one that gives the best results.An objective function value and other performance evaluation criteria can be used to compare the different NR options.The continuity and storage equations used in the distributed nonlinear Muskingum model that includes lateral inflows are presented as follows: where the lateral inflows varied linearly along the river reach and could be represented as a ratio of the inflow rate by considering the β parameter.β allows for the consideration of lateral inflow or outflow from the main channel during flood events.It represents the ratio of the inflow or outflow to the main channel flow within the reach.One of the assumptions used in the modeling process is that β is constant in time that means hydrograph shape of the lateral inflow wave is proportional to the upstream hydrograph inflow.t is the measure of time between zero and the flood's finish time.The spatial index between 2 and NR + 1 is called j.The following stages are used in the routing strategy for the distributed nonlinear Muskingum model utilizing the fourth order Runge-Kutta method: 1. Choose one, two, three, or more sub-reaches as NR and assume random values for the hydrological model parameters K, x, and m, as well as β. 2. Use Eq. ( 8) to estimate the starting storage.The starting flow rate at each sub-reach's downstream part is the same as the initial flow rate at the sub-reach's upstream section.Calculate the next storage.3. The next storage is computed by the present value plus the product of the size of the interval, Δt, and an estimated slope.The slope will be a weighted average of the following slopes using the Fourth order Runge-Kutta method: (7)

Salp swarm algorithm (SSA)
The Salp swarm algorithm (SSA) is a population-based swarm intelligence algorithm developed in 2017 by Mirjalili et al. 37 .The food source, which is the objective of the swarm, is represented by F. The leader of the swarm updates its position using a specific equation below: where x 1 j is the position of leader in jth dimension, Ub j are the upper and lower boundary at jth dimension, F j is the food source position.The coefficient c 1 plays an important role in SSA balancing exploration and exploita- tion.During the process of optimization, exploration refers to searching the search space thoroughly to find better solutions, while exploitation refers to utilizing the information present in the local region to improve the current solution.The parameter c 1 is gradually decreased over iterations and can be calculated using the follow- ing formula.
where l is the current iteration and L is the maximum number of iterations.The parameters c 2 and c 3 are random numbers generated within the interval [0,1].c 3 is responsible for indicating whether the next position of current leader salp should be toward + ∞ or -− ∞.The other members of the salp swarm update their positions based on Newton's law of motion, which is expressed using the following equation: where i ≥ 2 , x i j is the position of the ith follower in the jth dimension, t is the time, v 0 is the initial speed, and a = v final v 0 where v = (x − x 0 )/t.Since the time is considered as iterations and v 0 = 0 , Eq. ( 15) can be reformulated as the equation below: where i ≥ 2 , x i j is the position of the ith follower in the jth dimension.The main steps of the SSA can be summarized as follows (see Fig. 3): • Parameter initialization: the algorithm starts by initializing the parameters such as the population size N, number of the iterations t and the maximum number of iterations max itr .• Initial population: we generate the initial population x i , i = {1, ..., n} randomly in the range [u,l], where u and l are the upper and lower boundaries, respectively.
Vol.:(0123456789) • Individual evaluations: every individual (solution) within the population is assessed by determining its value using the objective function, and the best overall solution is designated as F. • Exploration and exploitation: to balance between the exploration and the exploitation of the algorithm, the value of the parameter c 1 is updated as shown in Eq. ( 16).• Update the position of the solutions: the position of the leader solution and the other follower solutions are updated as shown in Eqs. ( 15) and ( 18), respectively.• Boundary violations: boundary violations occur when a solution goes beyond the allowable range of the search space while updating, and it is then adjusted to fall within the problem's range.• Termination criteria: the number of iterations t is increased gradually until it reaches to the maximum number of iterations max itr .Then the algorithm terminates the search process and produces the overall best solution found.

Statistical performance evaluation criteria
Statistical performance evaluation criteria are metrics used to assess the accuracy and reliability of mathematical models, such as hydrological or hydraulic models.Some common statistical performance evaluation criteria used in the papers referenced 5,14,[53][54][55][56][57][58][59] are applied to assess the performance of the SSA-based routing results.These criteria are described below.where Q i and Q ι respectively are the observed and calculated outflow rates at the ith time, and N is the number of data.
The sum of absolute differences (SAD): the SAD measures the sum of the absolute differences between predicted and observed values.It is a measure of the overall deviation of the model from the observed data and is useful for evaluating the model's performance under conditions where large errors may have a significant impact.Difference of peak observed (DPO): the DPO measures the difference between the predicted discharge values from the observed peak discharge values.It is a measure of the model's ability to accurately predict extreme events, such as floods or droughts, that may have a significant impact on the environment or society 53 .

The deviation of peak time of routed and actual outflows (DPOT).
T pobs and T pest denote the observed and estimated times to peak discharge, respectively.All the criteria pre- sented are measurements of the accuracy of a routing model, with the optimum value at 0.

Table 2 presents estimates of hydrologic parameters and performance evaluation criteria (PEC) values for Group
A, considering different numbers of sub-reaches.The table encompasses sub-reach configurations ranging from 1 to 4, each characterized by Muskingum parameters (k, x, m, and β).The performance evaluation criteria include SSE, SAD, DPO, and DPOT.
The table clearly demonstrates that the number of sub-reaches employed has a substantial influence on both the estimates of hydrologic parameters and the corresponding PEC values.Analysis of the performance criteria indicates that the optimal performance is achieved when the Red River at Drayton station is considered as a single sub-reach.This is primarily due to the single sub-reach model exhibiting the lowest SSE values.Using the current study algorithm for Group A flood, the optimized parameters were determined to be K = 0.54, x = 0.24, m = 1.478, and β = 0.19 for NR = 1.
The hydrologic parameters estimates and performance evaluation criteria (PEC) values for Group B were investigated using various numbers of sub-reaches, ranging from 1 to 4. Each sub-reach was characterized by Muskingum parameters, namely K, x, m, and β.These findings are presented in Table 3, along with performance  www.nature.com/scientificreports/metrics including SSE, SAD, DPO, and DPOT.The analysis of performance criteria revealed that the optimal performance was achieved when employing two sub-reaches (NR = 2).This was primarily attributed to the fact that NR = 2 yielded the lowest SSE values.It is noteworthy that the Muskingum parameters for NR = 2 were determined as follows: K = 0.06, x = 0.06, m = 1.46, and β = 0.16.The calculation of performance evaluation criteria (PEC) values necessitates the availability of at least one year of observed flood data for the corresponding validation period.Unfortunately, for Group C, no year of observed data was available, rendering the calculation of PEC values impossible.However, the hydrologic parameters estimate provided in Table 4 can still be employed to evaluate the performance of the model.These estimates allow for comparisons between the model's predictions and those of other models.The Muskingum parameters determined using the SSA technique are presented in Table 4.
For evaluation of the developed model in real field condition validation step has been considered.Tables 5  and 6 present examples of calibration simulation for Group A and B, respectively.Table 5 presents the measured inflow data from the Grand Forks USGS station, measured outflow data from the Drayton USGS station, and the corresponding routed outflow values for Group A, specifically when using a single sub-reach (NR = 1).The calibration years considered for Group A include 1997, 2001, 2005, 2006, 2009, 2010, and 2018.Similarly, Table 6 displays the measured inflow data from the Grand Forks USGS station, measured outflow data from the Drayton USGS station, and the associated routed outflow values for Group B. For Group B, the calibration years selected are 1999, 2004, and 2013, and these results correspond to the scenario where two sub-reaches are utilized (NR = 2).
These tables provide essential data for the calibration and validation of the respective models and facilitate the comparison between the simulated and observed flow values during the specified calibration years for each group.
Table 7 displays the validated inflow data from the Grand Forks USGS station, measured outflow data from the Drayton USGS station, and the corresponding routed outflow values for Group A in the year 2020 in Drayton, specifically when utilizing a single sub-reach (NR = 1).
Likewise, Table 8 presents the validated inflow data from the Grand Forks USGS station, measured outflow data from the Drayton USGS station, and the associated routed outflow values for Group B in the year 2022 in Drayton.For Group B, the model configuration involved two sub-reaches (NR = 2).
Furthermore, during the evaluation of the model's performance, it was observed that the model accurately predicted the maximum outflow discharge for both Group A in 2020 (NR = 1) and Group B in 2022 (NR = 2).Specifically, for Group A, the model's prediction of the maximum outflow discharge on April 14th exhibited only a 3.7% difference compared to the measured value.Similarly, for the 2022 flood classified as Group B with NR = 2, the difference was 2.84%.These minimal differences indicate a close agreement between the model's predictions and the actual measured values, suggesting a high level of accuracy in the model's ability to forecast future floods.
Figures 4 and 5 complement the information found in Tables 7 and 8 by presenting a visual depiction of the validation results for the flood data of Group A and B. The figures provide a graphical representation that enhances the understanding and analysis of the validation outcomes for the respective datasets.
Figures 4 and 5 serve as visual representations of the validation results for the flood data of Group A and B, respectively.These figures complement the information presented in Tables 7 and 8 by providing a graphical depiction of the validation outcomes.By utilizing visual representations, it becomes easier to comprehend and analyze the results of the validation process for each dataset.These figures display important information such as observed and simulated flood hydrographs, peak flow values, and timing of peak flows.They also show how well the simulated hydrographs match the observed data, indicating the accuracy of the routing model.
Upon careful examination of Tables 2, 3 and 4 it is evident that during the validation phase, the maximum deviation in peak time between the routed and actual outflows was only 2 units.This observation aligns with our previous experiences and findings, where we have encountered similar discrepancies when employing a dynamic wave model that considers all the necessary inputs for flood modeling 1,34 .It is worth noting that when soft computing models are brought into the equation, such errors tend to be more pronounced.
It is essential to note that the absence of attenuation in the flood peak on the downstream side could be due to various factors.These include the absence of significant floodplain storage, a sufficiently large channel capacity, the lack of hydraulic controls such as dams or levees, and specific hydrological conditions such as continuous heavy rainfall or rapid snowmelt.

Discussion
In our recent paper, we conducted an existing case study to assess the performance of the proposed model using the inflow-outflow hydrograph data from Wilson 60 , which represents a smooth single-peak hydrograph.The distributed Muskingum method, implemented with the WOA algorithm, was employed for the routing process.Table 9 provides an overview of the Muskingum model routing parameters and performance evaluation criteria (PEC) obtained for the nonlinear Muskingum model parameters, as described in the material and methods section of the paper.The optimal parameter values for the Wilson flood data were determined using our research technique, yielding K = 0.865, x = 0.043, m = 1.478, and β = − 0.008 for NR = 3. www.nature.com/scientificreports/Furthermore, we analyzed the impact of varying the number of sub-reaches (NR) on the model's performance.The maximum outflow discharge, occurring at the 60th hour of the flood data, was examined for NR values ranging from 1 to 5. The results showed that the difference in peak discharge varied as follows: − 1.67%, 1.01%, 0.13%, 1.18%, and 0.96% for NR = 1 to NR = 5, respectively.Notably, NR = 3 exhibited the lowest difference in peak discharge compared to the Wilson flood data.
It is clear from comparing the findings of this case study with prior studies that the proposed model demonstrates consistent accuracy in predicting peak discharge across different datasets.Both the Group A and Group B case studies, along with the Wilson flood data, indicate close agreement between the model's predictions and the observed values.This reaffirms the model's reliability and its potential for effectively predicting future floods.
In comparison to the previous studies, our research offers distinctive findings.While Ayvaz and Gurarslan 61 introduced a novel partitioning approach for flood routing models, our study focuses on analyzing the impact of varying the number of sub-reaches on the performance of the nonlinear Muskingum model.In contrast to the primary emphasis of Hirpurkar and Ghare 62 on parameter estimation, our research evaluates the accuracy of peak discharge predictions using the proposed model.Furthermore, while Barbetta et al. 63 addresses river discharge estimation and rating curve development, our study concentrates on peak discharge prediction and the model's reliability across diverse datasets.Importantly, our study distinguishes itself by utilizing a specific case study with real flood data, providing a unique and valuable contribution to the field.By analyzing actual data, our findings offer practical insights, showcasing the effectiveness of the proposed model in real-world scenarios.This emphasis on practicality enhances the applicability and relevance of our research.

Conclusion
This study aims to develop and evaluate a nonlinear Muskingum model for river modeling, with a specific focus on enhancing accuracy through the incorporation of lateral inflow.The Grand Forks and Drayton USGS stations serve as case studies for parameter estimation using the distributed Muskingum method.The results demonstrate that the developed nonlinear Muskingum model effectively routes floods through these stations.Our selection of this area is rooted in a prior spatial analysis we conducted, which highlighted a particular region's vulnerability.This area, situated between Grafton city, Grand Forks, and Emerson, demonstrated a high susceptibility to severe floods when we examined the permanent water area (PWA) and seasonal water area (SWA).Building upon these previous findings, our current research hones in on this specific locale, aiming to delve deeper into flood dynamics and comprehensively grasp the underlying factors that make it exceptionally vulnerable 42 .Through an analysis of hydrologic parameters and model performance, the study underscores the significant impact of the number of sub-reaches on parameter estimation and modeling precision.Optimal model performance varies between case studies, underscoring the importance of selecting the appropriate number of sub-reaches for precise peak discharge predictions.
For Group A, where snowmelt is the primary driver with minimal rain-on-snow impact, the findings indicate that a single sub-reach provides the best performance in accurately predicting peak discharge.This suggests that a simplified representation of the river system adequately captures the routing characteristics of this flood type.In contrast, for Group B, representing flood events primarily characterized by rain-on-snow conditions, the study reveals that employing two sub-reaches optimizes model performance.This implies that considering additional sub-reaches is essential for accurately capturing the routing dynamics and behavior of this specific flood type.
The findings hold practical significance for flood prediction and management, offering valuable insights to decision-makers and stakeholders involved in flood mitigation.The proposed model holds potential for broader application beyond the specific case studies, contributing to enhanced river system modeling and flood management practices.For future studies, non-linear Muskingum model and Muskingum-Cunge approach should be compared in real field case studies to better understand their abilities compare to each other.Moreover, β parameter of the proposed model can be consider variable to examine effects of its variation on the model prediction performance.

Figure 2 .
Figure 2. Models for distributed nonlinear Muskingum model: (a) single reach with no sub-reaches, (b) two sub-reaches within a reach, (c) three sub-reaches within a reach, and (d) multi-interval sub-reach within a reach.

Figure 4 .
Figure 4. Validation for Data of Group A, 2020 flood data.

Figure 5 .
Figure 5. Validation for Data of Group B, 2022 flood data.

Table 1 .
Characteristics of the spring snowmelt flood events at two Red River stream gauging stations.

Table 2 .
Hydrologic parameters estimates and PEC values for different numbers of sub-reaches applied for Group A. Significant values are in [bold].

Table 3 .
Hydrologic parameters estimates and PEC values for different numbers of sub-reaches applied for Group B. Significant values are in [bold].

Table 4 .
Hydrologic parameters estimates and PEC values for different numbers of sub-reaches applied for Group C.

Table 5 .
Calibration and calculations for single-reach Muskingum flood routing applied to Data of Group A.

Table 6 .
Calibration and calculations for single-reach Muskingum flood routing applied to Data of Group B.

Table 8 .
Validation and calculations for multi-reach Muskingum flood routing applied to Data of Group B.

Table 9 .
Comparison of the outflow hydrographs calculated for the Wilson flood data.