Explicit Scheme for a Hydrological Channel Routing: Mathematical Model and Practical Application

Abstract

The computation of hydrographs in large watersheds necessitates utilizing channel routing, which calculates the movement of hydrographs along channel branches. Routing methods rely on an implicit scheme to facilitate numerical resolution, which requires more computational time than the explicit scheme. This study presents an explicit scheme channel routing model that offers a versatile approach to open channel flow analysis. The model is based on mass conservation principles and Manning equations, and it can accommodate varying bed slopes, making it highly adaptable to diverse hydraulic scenarios. In addition, the proposed model considers backwater effects, which enhances its applicability in practical scenarios. The model was tested in a practical application on a rectangular channel with a width of 7 m, and the results showed that it can accurately predict outflow hydrographs and handle different flow conditions. Comparative analyses with existing models revealed that the proposed model’s performance in generating water flow oscillations was competitive. Moreover, sensitivity analyses were performed, which showed that the model is highly responsive to parameter variations, such as Manning’s coefficient, bed slope, and channel width. The comparison of peak flows and peak times between the proposed model and existing methods further emphasized the model’s reliability and efficiency in simulating channel routing processes. This research introduces a valuable addition to the field of hydrology by proposing a practical and effective channel routing model that integrates essential hydraulic principles and parameters. The results of the proposed model (lumped routing) are comparable with the solution provided by the Muskingum–Cunge method (distributed routing). It is of utmost importance to note that the proposed model applies to channel branches with bed slopes below 6°.

1. Introduction

Open channel routing is a procedure aimed at determining the timing and magnitude of flow at specific points along a channel branch, utilizing a measured or estimated inflow hydrograph [1,2,3]. The routing is carried out using lumped (hydrological) or distributed (hydraulic) methods [4].

When an outflow hydrograph is primarily calculated using the conservation of mass equation, it is referred to as lumped routing, with the Lag and Muskingum methods among the developed methodologies.

Distributed channel routing theories are grounded in continuity and momentum equations, recognized as the Saint-Venant equations for dimensional flow. The Saint-Venant equations are the more complex solutions for one-dimensional routing methods in the current literature. These partial differential equations enable the calculation of flow and water level as functions of space and time. From the Saint-Venant equations [5], the numerical resolution has been calculated using Kinematic [6,7], Diffusive, and Dynamic waves theories, which are derived considering or neglecting some terms. Dynamic wave theory has given rise to models such as the DWOPER (Dynamic Wave Operational Model) from the U.S. National Weather Service Hydrological Research Laboratory [8], where the Saint-Venant equations were solved using a weight four-point nonlinear implicit finite difference scheme. This method has been used to model storage in rivers.

Additionally, other models like the FLDWAV [9], a combination of the NWS DAMBRK and DWOPER models, were developed to analyze natural floods and dam break analysis. Most rainfall-runoff models used to employ the Kinematic Wave [7] and the Muskingum–Cunge [10] methods for simulating hydraulic routing, which uses representative hydraulic parameters such as Manning’s coefficient, geometrical proprieties of a cross-sectional area, bed slope, and the length of an open channel branch.

Table 1 describes the main characteristics of lumped and distributed models that are considered for simulating rainfall-runoff processes [11].

Scholars have diligently endeavored to augment the efficacy of simulating hydrological processes [16]. The adoption of distributed models for simulating rainfall-runoff phenomena has grown progressively intricate, necessitating a plethora of computations. This intricacy arises from myriad factors, encompassing the numerical resolution of infiltration equations, analysis of temporal and spatial rainfall distribution, application of the numeric convolution method, and determination of unit hydrographs. These challenges must be addressed for individual sub-basins and many cells within a watershed [17,18]. Reducing computational times for open channel routing is imperative to diminish the duration of computations when executing distributed rainfall-runoff models.

Herein, the efficiency of these methodologies holds paramount significance. Lumped routing expends less computational time than distributed routing, as the latter necessitates computations for selecting spatial and space steps during simulations. In this vein, Tseng (2010) [7] propounded an efficient finite-difference implicit MacCormack scheme for simulating the Kinematic Wave method. Incorporating the longitudinal hydrostatic pressure gradient term into the Muskingum method was implemented [15]. Several adaptations of the Muskingum–Cunge method have been proposed, considering continuous flow and modified inflows [19,20]. Recently, machine learning techniques have been integrated to enhance distributed routing in rivers [21,22]. Linear and nonlinear analysis using the lumped routing method of Muskingum was analyzed [23,24]. The analysis of compound cross-sections has been implemented to account for varying water velocities, especially in floodplain zones [25,26]. Additionally, the influence of cross-sectional orientation has been considered for simulating distributed channel routing [27]. This analysis has been developed by considering the unsteady governing equations in open channels [28].

This study presents a novel model designed for open channel routing, predicated on a hydrological channel routing framework (lumped routing) that employs straightforward mathematical methodologies. The formulations adopt an explicit scheme, requiring minimal computational resources, and are founded on integrating the mass conservation equation with the Manning equation. Existing distributed channel routing methods typically employ an implicit scheme as the numerical resolution method. They also divide the channel branch into various parts to apply the partial or complete resolution of the Saint-Venant equations. However, these approaches result in longer computational times than the proposed model (an explicit scheme method).

The proposed model facilitates the assessment of hydrological routing within a channel, explicitly incorporating trapezoidal, rectangular, and triangular cross-sectional geometries. The proposed model undergoes practical application to validate its efficacy and is juxtaposed with conventional lumped and distributed models (such as Lag, Muskingum, Kinematic Wave, and Muskingum–Cunge methods), demonstrating its suitability as a new routing approach. The results confirm that the Muskingum–Cunge method, a more comprehensive distributed method, yielded similar results to the proposed model.

2. Materials and Methods

The methodology of this study comprises two main components: (i) the development of the proposed model (Section 2.1), which is the primary focus of this investigation as it introduces a novel approach related to lumped routing; (ii) the evaluation of existing channel routing methods (Section 2.2), which was conducted to facilitate a comparative analysis with the proposed model.

2.1. Proposed Model

This mathematical model falls into the category of lumped routing methods. It is distinguished by considering key hydraulic characteristics of channels, such as bed slope, cross-sectional area, roughness coefficient, and total length of a channel branch. Notably, these attributes are overlooked by existing lumped channel methods like the Lag and Muskingum models. Furthermore, the proposed model eliminates the need to select different spatial steps in channels, a requirement in methods such as the Kinematic Wave and Muskingum–Cunge approaches. As a result, in watersheds with multiple sub-basins and routing analyses, computational times using the proposed model are shorter than those of existing distributed methods. The proposed model is applicable to mild slopes (less than 6°). The primary sources of uncertainty in the proposed model stem from the friction slope and the Manning’s coefficient. The assumptions and governing equations of the proposed model are presented below.

2.1.1. Assumptions

The proposed model has the assumptions as follows:

• The conservation of mass equations is used for the lumped routing.

• The Manning’s equation is utilized to evaluate the friction losses of a channel branch, and it is applicable for bed slopes lower than 6° (hydrostatic pressure distribution).

• Bed slope and cross-sectional area are considered constant over time.

• Backwater effects are considered for channel routing.

2.1.2. Governing Equations of the Proposed Model

Considering the flow hydrograph depicted in Figure 1, the average flow rate for times $t$ and $t+∆t$ is determined by the following:

(1)$\overline{Q}=\frac{Q_t+Q_{t+∆t}}{2}$

where $\overline{Q}$ = average flow rate for times $t$ and $t+∆t$, $Q_t$ = flow rate at $t$, and $Q_{t+∆t}$ = flow rate at $t+∆t$.

The volume ($V$) that is passing the time interval $∆t$ can be computed as follows:

(2)$V=\left(\frac{Q_t+Q_{t+∆t}}{2}\right)∆t$

Now, considering a channel with two cross-sectional areas perpendicular to the flow (as shown in Figure 2), located at chainages $i$ and $i+1$, the following relationship can be established by considering the conservation of mass equations over the time interval $∆t$. Equation (3) is derived from a lumped routing method, which considers both the inflow and outflow hydrographs and the storage volume.

(3)$\left(\frac{Q_{I,t}+Q_{I,t+∆t}}{2}\right)∆t-\left(\frac{Q_{O,t}+Q_{O,t+∆t}}{2}\right)∆t=∆y\overline{T}L$

where $Q_I$= inflow hydrograph at chainage $i$, $Q_{O,t}$ = outflow hydrograph at chainage $i+1$, $\overline{T}$ = average top width of the cross-sectional areas at chainages $i$ and $i+1$, $L$ = length between sections $i$ and $i+1$, and $∆y$ = variation in depth of flow between sections $i$ and $i+1$.

Expanding upon Equation (3), is possible to determine the following:

(4)$Q_{I,t}+Q_{I,t+∆t}-Q_{O,t}=\left(\frac{2L}{∆t}\right)\overline{T}∆y+Q_{O,t+∆t}$

Considering the following:

(5)$a=\left(\frac{2L}{∆t}\right) \overline{T}$

and

(6)$C=Q_{I,t}+Q_{I,t+∆t}-Q_{O,t}$

Thus, reorganizing terms in Equation (4):

(7)$Q_{O,t+∆t}=C-a∆y$

Considering a backwater downstream occurrence for a time $t$, then the conservation of mass is expressed as follows:

(8)$\left(\frac{Q_{I,t}-Q_{I,t+∆t}}{2}\right)∆t-\left(\frac{Q_{O,t}-Q_{O,t+∆t}}{2}\right)∆t=∆y\overline{T}L$

For this condition is possible to demonstrate the following:

(9)$Q_{O,t+∆t}=-(C-a∆y)$

Given a trapezoidal channel as shown in Figure 3, then the flow area is given by the following:

(10)$A=By+\frac{1}{2}{y}^2\left(z_1+z_2\right)$

where $A$ = flow area, $B$ = width of a cross-sectional, $y$ = depth of flow, $z_1=$ left side slope, and $z_2$= right side slope.

The analysis of a trapezoidal section involves the examination of values for $z_1$ and $z_2$. When $z_1$ and $z_2$ both equal 0, a rectangular cross-sectional area is depicted. Conversely, when $B$ equals 0, a triangular section is formed. This suggests that the proposed model has been developed considering these cross-sectional areas.

Now, if $K=(z_1+z_2)/2$, then:

(11)$A=By+K\ast y^2$

Then, the top width can be expressed as follows:

(12)$T=B+2Ky$

The wetted perimeter ($P$) for a trapezoidal section, featuring distinct left and right side slopes, is expressed as follows:

(13)$P=B+y\left(\sqrt{z_1^2+1}+\sqrt{z_2^2+1}\right)$

Considering that $J=\sqrt{z_1^2+1}+\sqrt{z_2^2+1}$, then:

(14)$P=B+Jy$

According to Manning’s equation, friction losses in a channel branch can be computed as follows:

(15)$Q=\frac{1}{n}\frac{A^{5/3}}{P^{2/3}}{S}^{1/2}$

where $Q$ = water flow, $n$ = Manning’s coefficient, and $S$ = bed slope of a channel.

Substituting Equations (11) and (14) into Equation (15), then:

(16)$Q=\frac{1}{n}\frac{{(By+K{y}^2)}^{5/3}}{{(B+Jy)}^{2/3}}{S}^{1/2}$

Using $b=S^{1/2}/n$, then Equation (16) is expressed for a time $t$ as follows:

(17)$Q_t=b\frac{{(By+K{y}^2)}^{5/3}}{{(B+Jy)}^{2/3}}$

For a time $t+∆t$, the depth of flow is giving by $y+∆y$, thus:

(18)$Q_{t+∆t}=b\frac{{\left[B\right(y+∆y)+K{\left(y+∆y\right)}^2]}^{5/3}}{{[B+J\left(y+∆y\right)]}^{2/3}}$

Considering a variable $M={\left[B\right(y+∆y)+K{\left(y+∆y\right)}^2]}^{5/3}$ and applying the binomial theorem for ${(y+∆y)}^2$, then:

(19)$M={\left[By+B∆y+K{y}^2+2Ky∆y+{K∆y}^2\right]}^{5/3}$

Neglecting the term ${∆y}^2$:

(20)$M={\left[By+K{y}^2+(B+2Ky)∆y\right]}^{5/3}$

If $D=By+K{y}^2$ and $E=B+2Ky$, then Equation (20) is expressed as follows:

(21)$M={\left[D+E∆y\right]}^{5/3}$

In Equation (18), the denominator is giving by a variable $N={[B+J\left(y+∆y\right)]}^{2/3}$. If $F=B+Jy$, then is possible to demonstrate the following:

(22)$N=F^{2/3}+\frac{2}{3}{F}^{-1/3}J∆y$

Inserting Equations (21) and (22) into Equation (18):

(23)$Q_{t+∆t}=b\left[\frac{D^{5/3}+\frac{5}{3}{D}^{2/3}E∆y}{F^{2/3}+\frac{2}{3}{F}^{-1/3}J∆y}\right]$

Replacing Equation (23) into Equation (7), thus:

(24)$C-a∆y=b\left[\frac{D^{5/3}+\frac{5}{3}{D}^{2/3}E∆y}{F^{2/3}+\frac{2}{3}{F}^{-1/3}J∆y}\right]$

Expanding Equation (24):

(25)$b{D}^{5/3}+\frac{5}{3}{D}^{2/3}bE∆y=C{F}^{2/3}+\frac{2}{3}{F}^{-1/3}CJ∆y-a∆y{F}^{2/3}-\frac{2}{3}{F}^{-1/3}Ja{∆y}^2$

If the term ${∆y}^2\approx 0$, then:

(26)$∆y\left(\frac{5}{3}{D}^{2/3}bE-\frac{2}{3}{F}^{-1/3}CJ+a{F}^{2/3}\right)=C{F}^{2/3}-b{D}^{5/3}$

Finally, for channel routing without accounting for backwater effects:

(27)$∆y=\left[\frac{C{F}^{2/3}-b{D}^{5/3}}{\frac{5}{3}{D}^{2/3}bE-\frac{2}{3}{F}^{-1/3}CJ+a{F}^{2/3}}\right]$

Similarly, for the condition of backwater downstream occurrence, considering Equations (9) and (23), then:

(28)$∆y=\left[\frac{C{F}^{2/3}+b{D}^{5/3}}{a{F}^{2/3}-\frac{2}{3}{F}^{-1/3}CJ-\frac{5}{3}{D}^{2/3}bE}\right]$

Considering that the formulation of the proposed model (see Section 2.1.2) is based on the Manning equation, it can be applied to calculate composite roughness, wherein the wetted area is divided into $N_p$ parts. Table 2 presents the various formulas for computing the equivalent roughness, depending on whether water velocities remain constant or vary within an analyzed cross-sectional area [29].

Table 2

Calculation of composite roughness in open channels.

Composite Roughness Formula	Application	Equation
$n={\left(\frac{\sum _1^{N_p}{P}_{N_p}{n}_{N_p}^{1.5}}{P}\right)}^{2/3}$	This formula applies to channels of suitable dimensions with varying Manning’s coefficients, provided no flooding exists. Consequently, each segment of the cross-sectional area experiences an equivalent average velocity.	(29)
$n=\frac{P{R}^{2/3}}{{\sum }_1^{N_p}\left(\frac{P_{N_p}{R}_{N_p}^{5/3}}{n_{N_p}}\right)}$	This equation is employed for channels containing floodplain zones, where different Manning coefficients and disparate water velocity distributions are observed.	(30)

2.2. Current Channel Routing Methods

The proposed model underwent comparison with traditional methods as delineated in Table 3. Lumped and distributed channel routing methods were performed in this research.

Table 3

Current channel routing methods.

Channel Routing Method	Formulation	Equation No.	Observations
Lag model	$Q_{O,t}=\left\{\begin{array}{lll}{Q}_{I,t}& \mathrm{f}\mathrm{o}\mathrm{r}& t<t_{lag}\\ Q_{I,t+t_{lag}}& \mathrm{f}\mathrm{o}\mathrm{r}& t\ge t_{lag}\end{array}\right.$where $t_{lag}$ = lag time.	(31)	The simplest method for hydrological routing.
Kinematic Wave	$Q_{i+1}^{j+1}=\frac{\frac{∆t}{∆x}{Q}_i^{j+1}+\alpha \beta Q_{i+1}^j{\left(\frac{Q_{i+1}^j+Q_i^{j+1}}{2}\right)}^{\beta -1}}{\frac{∆t}{∆x}+\alpha \beta {\left(\frac{Q_{i+1}^j+Q_i^{j+1}}{2}\right)}^{\beta -1}}$ where $\alpha ={\left(\frac{n{P}^{2/3}}{S^{1/2}}\right)}^{3/5}, \mathrm{and} \beta =0.6$.	(32)	$\mathrm{The} \mathrm{numerical} \mathrm{resolution} \mathrm{is} \mathrm{conducted} \mathrm{using} \mathrm{a} \mathrm{x}\text{–}\mathrm{y} \mathrm{plane}, \mathrm{where} \mathrm{an} \mathrm{unknown} \mathrm{flow} Q_{i+1}^{j+1}$$\mathrm{is} \mathrm{computed} \mathrm{based} \mathrm{on} \mathrm{known} \mathrm{flows} (Q_i^j$$, Q_{i+1}^j$$, \mathrm{and} Q_i^{j+1}$$) \mathrm{considering} \mathrm{appropriate} \mathrm{time} (∆t$$) \mathrm{and} \mathrm{space} (∆x$) steps.
Muskingum model	$S_t=K{Q}_{O,t}+K_{t}X(Q_{I,t}-Q_{O,t})$where $S_t$ = storage at time $t$, $K_t$ = travel time in a channel (or storage constant), and $X$ = dimensionless weight.	(33)	This practical method requires the calibration of $K_t$$\mathrm{and} X$ parameters for performing accurate simulations.
Muskingum–Cunge model	$\frac{\partial Q}{\partial t}+c\frac{\partial Q}{\partial x}=\mu \frac{{\partial }^{2}Q}{\partial x^2}+c{q}_L$where $c$ = wave celerity, $\mu$ = hydraulic diffusivity, an $q_L$ = lateral inflow (if applicable)	(34)	The hydraulic diffusitivity depends on hydraulic parameters.

It is of utmost importance to emphasize that the proposed model utilizes the variables of distributed routing methods—Manning’s coefficient, bed slope, cross-sectional area, and length of the channel branch—for channel routing computation despite being classified as a lumped model.

3. Results for a Practical Application

In order to examine the application of the proposed model, an open channel routing was carried out, focusing on a rectangular cross-sectional area. The following data were considered: $n=$ 0.015 (complete rectangular section), $S=$ 0.12%, $B=$ 7.0 m, $z_1=z_2=0$, and $L=$ 1200 m. A time step ($∆t$) of 0.25 min was used for simulations. The inflow hydrograph is presented in Figure 4, showing a peak flow of 15.56 m³/s occurring at 40.75 min. The total duration of the inflow hydrograph is 83.75 min.

The channel routing is executed employing the proposed model outlined in Section 2.2. Figure 5 illustrates the outcomes of the channel routing. The peak value of the outflow hydrograph is 14.10 m³/s, observed at 44.0 min, shifting by 3.5 min compared to the inflow hydrograph. Under these conditions, the outflow hydrograph lasts around 120 min, presenting a maximum flow depth of 1.00 m. The depth of the flow pattern aligns with the outflow hydrograph as these variables are interconnected.

A sensitivity analysis was conducted, varying the critical parameters of the proposed model. Table 4 outlines the parameters considered for the analysis, including Manning’s coefficient ($n$), bed slope ($S$), and width of the open channel ($B$). Figure 6 illustrates the outcomes. A comparison is made of the results of the outflow hydrograph. The initial data (Figure 5) is labeled as the baseline solution in Figure 6.

The impact of Manning’s coefficient on the outflow hydrograph is evident in Figure 6a. A lower Manning’s coefficient corresponds to higher peak flow values. For example, with $n$ = 0.013, a peak value of 14.29 m³/s is reached, while with $n$ = 0.017, a value of 13.91 m³/s is obtained. Generally, within a reasonable range, variations in Manning’s coefficient do not significantly influence the pattern of the outflow hydrograph.

The bed slope ($S$) can markedly alter the pattern of the outflow hydrograph, depending on the value considered, as demonstrated in Figure 6b. A greater bed slope leads to higher peak flow values and earlier peak times. For a slope of 4.00%, a maximum flow value of 15.24 m³/s is attained at 41.5 min, closely resembling the inflow hydrograph (Figure 4), since steep slopes cannot retain sufficient water volume within the 1200 m length considered. Conversely, mild bed slopes tend to store significant water volume during open channel routing, as observed with a slope of 0.06%, where the resulting peak value is 13.52 m³/s at 44.5 min, indicating a reduction of around 13.0% compared to the peak value of the inflow hydrograph (15.56 m³/s).

The width of the open channel varied from 4 to 20 m. Higher width values correspond to lower peak flow values. For instance, with a width ($B$) of 4 m, a peak flow value of 14.15 m³/s is observed, whereas with $B$ = 20 m, it decreases to 13.23 m³/s.

4. Discussion

In this section, an examination of the proposed model has been undertaken. Initially, the time step was varied by adding 1.00 and 5.00 min increments. Figure 7 illustrates the findings, demonstrating that employing a time step of 0.25 (refer to Section 3) and 1.00 min yields comparable peak flows and analogous water pulses. Conversely, a slight beat disparity emerges when a time step of 5 min is utilized. Notably, a maximum water flow of 13.99 m³/s is calculated for this time step, in contrast to a peak value of 14.10 m³/s (for $∆t$ = 0.25 min). Furthermore, the peak time is advanced, occurring at 40.0 min (for $∆t$ = 5.00 min) as opposed to the timing of 43.5 min (for $∆t$ = 0.25 min).

The results of the proposed model were subject to a comparative evaluation against the lumped and distributed methods outlined in Table 1 for practical application. Table 5 delineates the factors considered in calculations, encompassing Lag, Kinematic Wave, Muskingum, and Muskingum–Cunge models, along with the parameters utilized, simulated via the HEC-HMS software (Version 4.12).

Figure 8 compares outflow hydrographs between the proposed model and the lumped and distributed models. The Muskingum method yields values akin to the proposed models, with peak water flows of 14.5 and 14.1 m³/s computed, respectively, exhibiting a discrepancy in peak time of 3 min. The Muskingum–Cunge method yields the closest values related to the peak flow, with 14.5 m³/s occurring at 44 min. Across all simulations, the Kinematic Wave and Lag methods yield higher peak water flows than those predicted by the proposed model. It is essential to mention that the proposed model (a lumped routing method) provides values similar to those of the Muskingum–Cunge (a distributed routing method).

A comparative analysis was conducted between the proposed model and other lumped and distributed models, considering the attributes outlined in Table 1. The proposed model operates as a lumped channel routing technique akin to Lag and Muskingum models. It is a pragmatic approach, leveraging typical hydraulic parameters such as Manning’s coefficient, bed slope, geometric cross-sectional area characteristics, and channel branch length to facilitate the routing of an inflow hydrograph while accommodating backwater effects. Unlike the Muskingum model, the proposed model does not rely on intricate parameters for determination. In contrast, it necessitates hydraulic parameters like those required by distributed models such as the Kinematic Wave and Muskingum–Cunge models. This renders it applicable for both mild and steep bed slopes without the need for adherence to a Courant condition, as it is not a distributed model. The Muskingum–Cunge method has yielded more precise outflow hydrograph computations than other hydrological and hydraulic methods [30].

The Root Mean Square Error (RMSE) was calculated for the outflow hydrographs, considering both lumped and distributed models (refer to Figure 8). The Muskingum–Cunge method was employed as the reference values. Subsequently, the RMSE was determined as follows:

(35)$RMSE=\sqrt{\frac{1}{N}{\left(C-O\right)}^2}$

where $C$ represents the computed values obtained from the proposed model, Lag, or Muskingum methods, $O$ corresponds to the observed values (in this instance, using the Muskingum–Cunge method), and $N$ denotes the number of observations modelled.

Table 6 displays the RMSE results. The proposed model exhibited the best accuracy compared to the Lag, Muskingum, and Kinematic Wave methods, underscoring its significance as a novel mathematical model for channel routing since the lowest value of RMSE was obtained (0.53 m³/s).

The bed slope is the most sensitive parameter for channel routing, as depicted in Figure 6b. Considering this, the results of the proposed model were contrasted with those of the Muskingum–Cunge method across bed slopes ranging from 0.06 to 4.0%. Figure 9 illustrates the outcomes, indicating that the proposed model (a lumped routing model) yielded comparable values to the Muskingum–Cunge method (a distributed routing model). This demonstrates the applicability of the proposed model as an explicit scheme, with the significant advantage of requiring less computational time than the Muskingum–Cunge method.

5. Conclusions

Channel routing methods have developed significantly to determine the transformation of an inflow hydrograph along a channel branch to ascertain the resultant conditions (outflow hydrograph). These methods typically fall into two categories: lumped and distributed models. Some models overlook critical terms, while others fail to account for backwater effects. Additionally, specific models rely on complex parameters that necessitate calibration for accurate determination.

This study devised a novel channel routing method rooted in the principles of mass conservation and Manning equations. This method is versatile and can accommodate both steep and mild bed slopes. It employs standard hydraulic parameters—including Manning’s coefficient, bed slope, geometric cross-sectional area characteristics, and channel branch length—throughout its calculations and incorporates considerations for backwater effects. The proposed model is applicable for channels with bed slopes lower than 6°.

To showcase the effectiveness of the proposed model, a practical application was undertaken, demonstrating its capabilities and yielding tangible results for a rectangular channel measuring 1200 m in total length and 7 m in width. This application considered a peak flow ranging from 15.6 to 14.1 m³/s, reflecting inflow and outflow hydrographs, respectively. A sensitivity analysis was conducted, varying Manning’s coefficient, bed slope, and the width of the open channel. The robustness of the proposed model was evident as variations in these parameters consistently influenced the generation of different outflow hydrographs.

The proposed model was compared with the Lag, Kinematic Wave, Muskingum, and Muskingum–Cunge models. The results of the proposed model reveal that it produces water flow oscillations within a similar range to existing models.

Future research should be focused on the assessment of the time-saving benefits offered by the proposed model compared to distributed routing models such as Kinematic Wave and Muskingum–Cunge. Furthermore, the implementation of an explicit scheme, akin to the proposed model, should be conducted while considering two-dimensional routing methods.

Author Contributions

Conceptualization, A.A.-P.; methodology, A.A.-P. and O.E.C.-H.; validation, A.A.-P.; formal analysis, A.A.-P., O.E.C.-H. and J.R.C.-H.; writing—original draft preparation, A.A.-P. and O.E.C.-H.; supervision, J.R.C.-H. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data may be acquired by contacting the corresponding author (A.A.-P.).

Conflicts of Interest

The authors declare no conflicts of interest.

Nomenclature

The following abbreviations were used in this manuscript:

$A$:	flow area (m2)
$a$:	value calculated as $\left(\frac{2L}{∆t}\right) \overline{T}$ (m2/s)
$B$:	width of a cross-sectional (m)
$b$:	relationship of $S^{1/2}/n$ (m1/3/s)
$C$:	value computed as $Q_{I,t}+Q_{I,t+∆t}-Q_{O,t}$ (m3/s)
$C_{index}$:	index celerity (m/s)
$c$:	wave celerity (m/s)
$D$:	value calculated as $By+K{y}^2$ (m2)
$E$:	value computed as $B+2Ky$ (m)
$F$:	value computed as $B+Jy$ (m)
$J$:	relationship between side slopes (-)
$K$:	average of side slopes (-)
$K_t$:	travel time in a channel (s)
$L$:	length between sections $i$ and $i+1$ (m)
$M$:	value computed as ${\left[B\right(y+∆y)+K{\left(y+∆y\right)}^2]}^{5/3}$ (m10/3)
$N$:	value computed as ${[B+J\left(y+∆y\right)]}^{2/3}$ (m2/3)
$n$:	Manning’s coefficient (s/m1/3)
$P$:	wetted perimeter (m)
$Q$:	water flow (m3/s)
$\overline{Q}$:	average flow rate (m3/s)
$Q_t$:	flow rate at $t$ (m3/s)
$Q_{t+∆t}$:	flow rate at $t+∆t$ (m3/s)
$Q_{index}$:	index flow (m3/s)
$q_L$:	lateral inflow (m3/s)
$S$:	bed slope (m/m)
$S_t$:	storage at time $t$ (m3)
$\overline{T}$:	average top width of the cross-sectional areas at chainages $i$ and $i+1$ (m)
$t_{lag}$:	lag time (s)
$V$:	water volume passing in an interval time (m3)
$v$:	bankfull water velocity (m/s)
$y$:	depth of flow (m)
$z_1$:	left side slope (-)
$z_2$:	right side slope (-)
$∆t$:	time step (s)
$∆x$:	space step (m)
$∆y$:	variation in depth of flow between sections $i$ and $i+1$ (m)
$\alpha$:	relationship of ${\left(\frac{n{P}^{2/3}}{S^{1/2}}\right)}^{3/5}$ (m2/5)
$\beta$:	coefficient with a value of 0.6 (-)
$\mu$:	hydraulic diffusivity
Subscript
$I$:	refers to an inflow hydrograph
$i$ and $i+1$	refers to chainages along a channel branch.
0:	refers to an outflow hydrograph
$t$:	time (s)

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Flow hydrograph representation.

Figure 2. Free surface scheme.

Figure 3. Diagram of a trapezoidal cross section.

Figure 4. Inflow hydrograph.

Figure 5. Results of the open channel routing.

Figure 6. Sensitivity analysis for the channel routing: (a) Manning’s coefficient; (b) bed slope; (c) width of the channel.

Figure 7. Variation of time step in the proposed model.

Figure 8. Comparison between the proposed model and the other channel routing methods.

Figure 9. Comparison between the proposed model versus Muskingum–Cunge method considering different bed slopes.

Table 1

Comparison of lumped and distributed channel routing methods.

Channel Routing Method	Type of Routing		Observations	References
	Lumped	Distributed
Lag model	X		Fails to incorporate the attenuation of an inflow hydrograph. Computes the outflow hydrograph solely by adding the lag time of the inflow hydrograph’s occurrence time.	[11]
Kinematic wave		X	Suitable for steep bed slopes. Require an appropriate selection of space and time steps for numerical stability. Neglects convective and local acceleration and pressure forces when solving the Navier–Stokes equations. Fails to account for backwater effects.	[1,7]
Muskingum model	X		Regards the total storage in a channel branch as the summation of a prism and wedge. Entails complex parameter calculations for routing, including travel time through the channel and dimensionless weight parameters.	[12,13,14]
Muskingum–Cunge model		X	Require an appropriate selection of space and time steps for numerical stability. Utilizes variable coefficients necessitating recalculations at each time step. Requires the adoption of an index method for calculations.	[10,15]

Table 4

Considered parameters for the sensitivity analysis.

Parameter	Units	Range
		From	To
$\mathrm{Manning}\prime \mathrm{s} \mathrm{coefficient} \left(n\right)$	-	0.013	0.017
$\mathrm{Bed} \mathrm{slope} \left(S\right)$	%	0.06	4.00
Width of the open channel $\left(B\right)$	m	4	20

Table 5

Characteristics of the different channel routing methods.

Channel Routing Method	Considerations	Used Parameters
Lag model	The lag time was computed based on the consideration of the storage constant in the Muskingum model.	Lag time ($t_{lag}$ = 13.2 min)
Kinematic Wave	The celerity index method was utilized, considering an index celerity ($C_{index}$) of 1.52 m/s. The number of branches was 4.	$n=$$0.015, S=$$0.12\%, B=$$7.0 \mathrm{m}, L=$$1200 \mathrm{m}, \mathrm{and} C_{index}$ = 1.52 m/s
Muskingum model	The calculation of the storage constant ($K_t$) is founded upon a bankfull water velocity $\left(v\right)$ of 0.91 m/s, computing wave celerity as $c$ = $1.67v$ = 1.52 m/s.	$K_t$ = 0.22 h and $X$ = 0.2
Muskingum–Cunge model	For the simulations, the values of $∆t$ and $∆x$ were automatically selected by the software. An index flow ($Q_{index}$) of 15.6 m3/s is established as the maximum flow.	$n=$$0.015, S=$$0.12\%, B=$$7.0 \mathrm{m}, L=$$1200 \mathrm{m}, \mathrm{and} Q_{index}$ = 15.6 m3/s

Table 6

Computation of Mean Square Error.

Method	RMSE (m3/s)	Observation
Proposed model	0.53	-
Muskingum	1.34	-
Kinematic Wave	0.64	-
Lag	1.45	-
Muskingum–Cunge	-	This method was considered as reference values

References

1. Chow, V.T.; Maidment, D.R.; Mays, L.W. Applied Hydrology; McGraw-Hill: New York, NY, USA, 1988.

2. Fenton, J.D. Flood Routing Methods. J. Hydrol.; 2019; 570, pp. 251-264. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2019.01.006]

3. Gąsiorowski, D.; Szymkiewicz, R. Inverse Flood Routing Using Simplified Flow Equations. Water Resour. Manag.; 2022; 36, pp. 4115-4135. [DOI: https://dx.doi.org/10.1007/s11269-022-03244-8]

4. Todini, E. Hydraulic and Hydrologic Flood Routing Schemes. Recent Advances in the Modeling of Hydrologic Systems; Bowles, D.S.; O’Connell, P.E. Springer: Dordrecht, The Netherlands, 1991; pp. 389-405. [DOI: https://dx.doi.org/10.1007/978-94-011-3480-4_18]

5. Theodor, S. Numerical Solution of Saint-Venant Equations. J. Hydraul. Div.; 1970; 96, pp. 223-252. [DOI: https://dx.doi.org/10.1061/JYCEAJ.0002262]

6. Smith, A.A. A Generalized Approach to Kinematic Flood Routing. J. Hydrol.; 1980; 45, pp. 71-89. [DOI: https://dx.doi.org/10.1016/0022-1694(80)90006-2]

7. Tseng, M.-H. Kinematic Wave Computation Using an Efficient Implicit Method. J. Hydroinform.; 2010; 12, pp. 329-338. [DOI: https://dx.doi.org/10.2166/hydro.2010.068]

8. Unver, O.; Mays, L.W.; Lansey, K. Real-Time Flood Management Model for Highland Lake System. J. Water Resour. Plan. Manag.; 1987; 113, pp. 620-638. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-9496(1987)113:5(620)]

9. Moreda, F.; Gutierrez, A.; Reed, S.; Aschwanden, C. Transitioning NWS Operational Hydraulics Models from FLDWAV to HEC-RAS. World Environmental and Water Resources Congress; ASCE: Reston, VA, USA, 2009; pp. 1-11. [DOI: https://dx.doi.org/10.1061/41036(342)292]

10. Barry, D.A.; Bajracharya, K. On the Muskingum-Cunge Flood Routing Method. Environ. Int.; 1995; 21, pp. 485-490. [DOI: https://dx.doi.org/10.1016/0160-4120(95)00046-N]

11. USACE. Hydrologic Modeling System HEC-HMS Technical Reference Manual; US Army Corps of Engineers: Davis, CA, USA, 2000.

12. Salvati, A.; Moghaddam Nia, A.; Salajegheh, A.; Shirzadi, A.; Shahabi, H.; Ahmadisharaf, E.; Han, D.; Clague, J.J. A Systematic Review of Muskingum Flood Routing Techniques. Hydrol. Sci. J.; 2024; 69, pp. 810-831. [DOI: https://dx.doi.org/10.1080/02626667.2024.2324132]

13. Tahiri, A.; Che, D.; Ladeveze, D.; Chiron, P.; Archimède, B. Network Flow and Flood Routing Model for Water Resources Optimization. Sci. Rep.; 2022; 12, 3937. [DOI: https://dx.doi.org/10.1038/s41598-022-06075-0] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35273171]

14. Cunge, J.A. On The Subject of A Flood Propagation Computation Method (Musklngum Method). J. Hydraul. Res.; 1969; 7, pp. 205-230. [DOI: https://dx.doi.org/10.1080/00221686909500264]

15. Tang, X.-N.; Knight, D.W.; Samuels, P.G. Volume Conservation in Variable Parameter Muskingum-Cunge Method. J. Hydraul. Eng.; 1999; 125, pp. 610-620. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-9429(1999)125:6(610)]

16. Lozano Sandoval, G.; Aníbal Monsalve, E.; García Reinoso, P.L.; Rodríguez Mejía, C.A.; Gómez Ospina, J.P.; Triviño Loaiza, H.J. Environmental Flow Estimation Using Hydrological and Hydraulic Methods for the Quindío River Basin: WEAP as a Support Tool. Inge Cuc; 2015; 11, pp. 34-48. [DOI: https://dx.doi.org/10.17981/ingecuc.11.2.2015.04]

17. Dotson, H.W. Watershed Modeling with HEC-HMS (Hydrologic Engineering Centers-Hydrologic Modeling System) Using Spatially Distributed Rainfall. Coping with Flash Floods; Springer: Berlin/Heidelberg, Germany, 2001; pp. 219-230.

18. Francés, F.; Vélez, J.I.; Vélez, J.J. Split-parameter structure for the automatic calibration of distributed hydrological models. J. Hydrol.; 2007; 332, pp. 226-240. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2006.06.032]

19. Lee, E.H. Development of a New 8-Parameter Muskingum Flood Routing Model with Modified Inflows. Water; 2021; 13, 3170. [DOI: https://dx.doi.org/10.3390/w13223170]

20. Lee, E.H.; Lee, H.M.; Kim, J.H. Development and Application of Advanced Muskingum Flood Routing Model Considering Continuous Flow. Water; 2018; 10, 760. [DOI: https://dx.doi.org/10.3390/w10060760]

21. Bindas, T.; Tsai, W.-P.; Liu, J.; Rahmani, F.; Feng, D.; Bian, Y.; Lawson, K.; Shen, C. Improving River Routing Using a Differentiable Muskingum-Cunge Model and Physics-Informed Machine Learning. Water Resour. Res.; 2024; 60, e2023WR035337. [DOI: https://dx.doi.org/10.1029/2023WR035337]

22. Li, L.; Jun, K.S. Review of Machine Learning Methods for River Flood Routing. Water; 2024; 16, 364. [DOI: https://dx.doi.org/10.3390/w16020364]

23. Singh Vijay, P.; Scarlatos Panagiotis, D. Analysis of Nonlinear Muskingum Flood Routing. J. Hydraul. Eng.; 1987; 113, pp. 61-79. [DOI: https://dx.doi.org/10.1061/(ASCE)0733-9429(1987)113:1(61)]

24. Meisam, B.; Reza, B.; Emrah, D.; Gokmen, T. Reverse Flood Routing in Rivers Using Linear and Nonlinear Muskingum Models. J. Hydrol. Eng.; 2021; 26, 04021018. [DOI: https://dx.doi.org/10.1061/(ASCE)HE.1943-5584.0002088]

25. Balamurugan, M.; Murty Bhallamudi, S. Flood Routing in an Ephemeral Channel with Compound Cross-Section. Sādhanā; 2016; 41, pp. 771-785. [DOI: https://dx.doi.org/10.1007/s12046-016-0511-x]

26. Costabile, P.; Macchione, F. Analysis of One-Dimensional Modelling for Flood Routing in Compound Channels. Water Resour. Manag.; 2012; 26, pp. 1065-1087. [DOI: https://dx.doi.org/10.1007/s11269-011-9947-2]

27. Jesna, I.; Bhallamudi, S.M.; Sudheer, K.P. Impact of Cross-Sectional Orientation in One-Dimensional Hydrodynamic Modeling on Flood Inundation Mapping. J. Flood Risk Manag.; 2023; 16, e12893. [DOI: https://dx.doi.org/10.1111/jfr3.12893]

28. Rahimi, H.; Yuan, S.; Tang, X.; Lu, C.; Singh, P.; Dehrashid, F.A. Study on Conveyance Coefficient Influenced by Momentum Exchange Under Steady and Unsteady Flows in Compound Open Channels. Water Resour. Manag.; 2022; 36, pp. 2179-2199. [DOI: https://dx.doi.org/10.1007/s11269-022-03130-3]

29. Chow, V.T. Open-Channel Hydraulics; McGraw-Hill: Tokyo, Japan, 1959.

30. Prawira, D.; Soeryantono, H.; Anggraheni, E.; Sutjiningsih, D. Efficiency Analysis of Muskingum-Cunge Method and Kinematic Wave Method on the Stream Routing (Study Case: Upper Ciliwung Watershed, Indonesia). IOP Conf. Ser. Mater. Sci. Eng.; 2019; 669, 012036. [DOI: https://dx.doi.org/10.1088/1757-899X/669/1/012036]