Abstract

Pollution in the environment has increased drastically in recent times. Environmentalists and scientists have faced challenges in finding ways to minimize pollution and manage its effects. Water pollution has tremendously increased, which causes human aquatic activities like fishing, transportation, sea mining, and others to be unpredictable. The influx of pollutants into the flow causes significant changes in flow density, bed level, flow depth, flow velocity, and consequently bed topography. Soluble pollutants affect the density of flow and hence its velocity and other flow parameters. Insoluble pollutants have a reasonable effect on bed level, flow depth, and bed topography. This study aims to mathematically model open-channel flow hydrodynamics and assess how Manning’s roughness coefficient and bed roughness influence flow velocity, water depth, and discharge under different flow conditions. The study investigates how the roughness of an open-channel bed affects the way water flows, particularly its velocity, depth, discharge, and energy loss. The differential equations governing the flow are first expressed in dimensionless form and discretized using the finite difference method with appropriate initial and boundary conditions. Numerical values are simulated using the Python program, and results obtained from this study are represented graphically. The results indicate that flow velocity decreases with an increase in pollutants in the flow; a consistent increase in average friction slope as the Manning roughness coefficient increases; and further demonstrate that flow depth increases with increasing Manning roughness coefficient.

Keywords

Manning’s Roughness, Finite Difference Method, Numerical Simulation, Pollution

Share and Cite:

1. Introduction

The depth-averaged equations are widely used to model hydraulic and environmental flows [1]. Their applications include river dynamics, flood propagation, sediment transport, and coastal circulation. Finite difference numerical methods have been used for solving depth-averaged flow equations. The study of open-channel flows requires an understanding of the flow in order to accurately capture its characteristics. Designing and studying engineering projects on rivers, understanding the hydraulic phenomena of the watercourse for water quality control, and predicting exceptional natural events associated with water flow constitute the main objectives of studying open-channel flow simulation. The effects of channel curvature and riparian vegetation on flow dynamics and sediment transport have been investigated using a depth-averaged two-dimensional model [2]. The objective of this paper is to determine the flow velocities and flow depth due to changes in the Manning roughness coefficient. Although considerable research has been conducted on the effects of either soluble or insoluble pollutants individually, studies that simultaneously incorporate both forms of pollution remain limited. This gap in the literature necessitates further investigation into the combined effects of pollutants on flow and sediment dynamics. Therefore, the main objective of the study is to model depth-averaged polluted flow and bed topography in an open channel. The specific objectives of this work are:

• To investigate the effects of variations in Manning’s roughness coefficient on flow velocity and flow depth in an open-channel flow system.

• To develop a sediment transport model that incorporates the influence of pollutant influx on flow and sediment dynamics.

• To establish the effects of increased solid pollution on bed height and bed topography.

2. Physical Model

Geometry

The flow parameters: flow height $h\left(x,y\right)$ , bed height $z\left(x,y\right)$ , and flow velocities u and v in both x and y directions, respectively, can be illustrated using Figure 1.

The fluid enters the computational domain at a velocity of $u\left(x\right)$ and $v\left(y\right)$ in the x and y directions, respectively, and exits through the outlet boundary.

3. Governing Equations

The mathematical model is based on the conservation equations of mass and momentum.

3.1. Continuity Equation

For two-dimensional depth-averaged shallow water flow, the continuity equation is given as,

Figure 1. Two-dimensional flow.

$\frac{\partial h}{\partial t}+\frac{\partial hu}{\partial x}+\frac{\partial hv}{\partial y}=q$ (1)

$h\left(x,y,t\right)$ is the water depth, $u\left(x,y,t\right)$ is the depth-averaged velocity in the x-direction, and $v\left(x,y,t\right)$ is the depth-averaged velocity in the y-direction.

Initial and Boundary Conditions

The continuity equation was initialized using a uniform water depth of 0.5 m throughout the 100 m computational domain. The initial longitudinal velocity was set to 1.5 m/s, giving an initial discharge of 0.75 m³/s for the 1 m-wide rectangular channel. At the upstream boundary, a constant inflow discharge of 0.75 m³/s was imposed. Where required, this was accompanied by an upstream flow depth of 0.5 m. At the downstream boundary, an open-flow condition was applied by prescribing the downstream water depth where observational data were available; otherwise, a zero-gradient condition for water depth was imposed. For the two-dimensional formulation, the lateral channel boundaries were treated as impermeable walls with zero transverse velocity. These conditions provide the required inflow, outflow, and initial water volume for updating the flow depth at each computational time step using the continuity equation.

3.2. Momentum Equation

A two-dimensional shallow water equation is expressed in Equation (2) in the x-direction and Equation (3) in the y-direction, respectively [3].

$\frac{\partial \left(uh\right)}{\partial t}+\frac{\partial \left(u^{2}h+\frac{1}{2}g{h}^2\right)}{\partial x}+\frac{\partial \left(uvh\right)}{\partial y}=-gh\frac{\partial z}{\partial x}-gh{S}_{f_x}$ (2)

$\frac{\partial \left(vh\right)}{\partial t}+\frac{\partial \left(uvh\right)}{\partial x}+\frac{\partial \left(v^{2}h+\frac{1}{2}g{h}^2\right)}{\partial y}=-gh\frac{\partial z}{\partial y}-gh{S}_{f_y}$ (3)

where

$S_{f_x}=n_b^2{h}^{\frac{-4}{3}}u\sqrt{u^2+v^2}$ (4)

$S_{f_y}=n_b^2{h}^{\frac{-4}{3}}v\sqrt{u^2+v^2}$ (5)

$S_{f_x}$ and $S_{f_y}$ are the friction slope terms in the x and y directions, respectively. $n_b$ is the Manning roughness coefficient for a clean and straight natural stream. $h$ is the flow depth, $u$ and $v$ are the depth-averaged velocity components in the x and y directions, respectively, $z$ is the bed level, and $g$ is the gravitational acceleration.

The Manning roughness coefficient for clean and straight natural streams is approximately $n=0.03$ [4], which results in Equations (2) and (3). Due to pollution, the flow changes its turbidity, and this results in gradual changes in the Manning roughness coefficient. This leads to a change in the friction slope, which is reflected in the following flow equations. Substituting Manning roughness values of 0.04, 0.05, and 0.06 into Equations (4) and (5), respectively, we obtain the following equations of friction slope.

Manning’s roughness coefficient is employed as a surrogate scenario parameter to represent alternative channel-bed roughness conditions and associated hydraulic resistance. The different values are used to assess the sensitivity of the modeled hydrodynamic characteristics to variations in channel roughness and are not intended to represent pollutant-specific parameters or direct pollutant concentration limits. Any influence of roughness on pollutant transport is therefore interpreted as an indirect effect mediated through changes in the hydrodynamic and dispersion characteristics of the flow.

Manning roughness: 0.04.

$S_{fx}^1=\frac{16}{9}{n}_b^2{h}^{\frac{-4}{3}}u\sqrt{u^2+v^2}$ (6)

$S_{fy}^1=\frac{16}{9}{n}_b^2{h}^{\frac{-4}{3}}v\sqrt{u^2+v^2}$ (7)

Manning roughness: 0.05.

$S_{fx}^2=\frac{25}{9}{n}_b^2{h}^{\frac{-4}{3}}u\sqrt{u^2+v^2}$ (8)

$S_{fy}^2=\frac{25}{9}{n}_b^2{h}^{\frac{-4}{3}}v\sqrt{u^2+v^2}$ (9)

Manning roughness: 0.06.

$S_{fx}^3=4n_b^2{h}^{\frac{-4}{3}}u\sqrt{u^2+v^2}$ (10)

$S_{fy}^3=4n_b^2{h}^{\frac{-4}{3}}v\sqrt{u^2+v^2}$ (11)

Substituting Equations (6) and (7) into Equations (2) and (3), respectively, we obtain Equations (12) and (13), respectively.

$\frac{\partial uh}{\partial t}+\frac{\partial \left(u^{2}h+\frac{1}{2}g{h}^2\right)}{\partial x}+\frac{\partial \left(uvh\right)}{\partial y}=-gh\frac{\partial z}{\partial x}-\frac{16}{9}gh{S}_{fx}$ (12)

$\frac{\partial \left(vh\right)}{\partial t}+\frac{\partial \left(uvh\right)}{\partial x}+\frac{\partial \left(v^{2}h+\frac{1}{2}g{h}^2\right)}{\partial y}=-gh\frac{\partial z}{\partial y}-\frac{16}{9}gh{S}_{f_y}$ (13)

Substituting Equations (8) and (9) into Equations (2) and (3), respectively, we obtain Equations (14) and (15), respectively.

$\frac{\partial uh}{\partial t}+\frac{\partial \left(u^{2}h+\frac{1}{2}g{h}^2\right)}{\partial x}+\frac{\partial \left(uvh\right)}{\partial y}=-gh\frac{\partial z}{\partial x}-\frac{25}{9}gh{S}_{fx}$ (14)

$\frac{\partial \left(vh\right)}{\partial t}+\frac{\partial \left(uvh\right)}{\partial x}+\frac{\partial \left(v^{2}h+\frac{1}{2}g{h}^2\right)}{\partial y}=-gh\frac{\partial z}{\partial y}-\frac{25}{9}gh{S}_{f_y}$ (15)

Substituting Equations (10) and (11) into Equations (2) and (3), respectively, we obtain Equations (16) and (17), respectively.

$\frac{\partial uh}{\partial t}+\frac{\partial \left(u^{2}h+\frac{1}{2}g{h}^2\right)}{\partial x}+\frac{\partial \left(uvh\right)}{\partial y}=-gh\frac{\partial z}{\partial x}-4gh{S}_{fx}$ (16)

$\frac{\partial vh}{\partial t}+\frac{\partial \left(uvh\right)}{\partial x}+\frac{\partial \left(v^{2}h+\frac{1}{2}g{h}^2\right)}{\partial y}=-gh\frac{\partial z}{\partial y}-4gh{S}_{fy}$ (17)

The Manning roughness coefficient increases with increasing turbidity and sediment concentration, rising from 0.03 under clean-water conditions to approximately 0.06 in polluted flows [5]. The increase in roughness reflects the enhanced hydraulic resistance associated with pollutant and sediment influx. Consequently, the flow variables, including the velocity components (u) and (v), flow depth (h), and bed elevation (z), exhibit a gradual and systematic variation over time at a fixed spatial location, indicating the dynamic response of the flow field and bed morphology to increasing levels of pollutants.

4. Initial and Boundary Conditions

4.1. Initial Conditions

$u\left(x,y,0\right)=10,v\left(x,y,0\right)=10,h\left(x,y,0\right)=0.5$ (18)

We define the non-dimensional variables. The velocity component is non-dimensionalised using the characteristic velocity $U_o$ ; water depth is non-dimensionalised using the characteristic depth $H_o$ , such that

$u^{\text{*}}=\frac{u}{U_o},v^{\text{*}}=\frac{v}{U_o},h^{\text{*}}=\frac{h}{H_o}$

Consequently, the prescribed non-dimensional initial conditions

$u^{\text{*}}\left(x,y,0\right)=10,v^{\text{*}}\left(x,y,0\right)=10,h^{\text{*}}\left(x,y,0\right)=0.5$

Corresponding to the dimensional initial conditions

$u\left(x,y,0\right)=10U_o,v\left(x,y,0\right)=10U_o,h\left(x,y,0\right)=0.5H_o$

Thus, the initial u and v velocity components are each ten times the selected reference velocity, while the initial depth is one-half of the selected reference depth; therefore, the corresponding dimensional values are,

Since $U_o=1.5$ and $H_o=1.2$ , then $u=15,v=15,h=0.6$ .

The initial velocity components $u=10$ and $v=10$ indicate that the flow possesses momentum in both the x-direction and y-direction at the beginning of the simulation. This is appropriate for modeling rivers, channels, or pollutant flows where flow is continuously moving. Assigning non-zero velocities avoids the trivial stationary solution that would result from zero initial velocities and allows the model to capture the evolution of flow dynamics immediately after the simulation starts. The initial water depth $h=0.5$ ensures that a finite volume of flow exists throughout the computational domain. A positive depth is required for the shallow water equations to remain physically meaningful and numerically stable, since zero depth may lead to computational instability.

4.2. Boundary Conditions

$u\left(0,y,t\right)=10,u\left(100,y,t\right)=10$ (19)

The computational domain was defined over a channel length of 100 m, discretised using a uniform spatial grid $\text{Δ}x=0.5025\text{m}$ . An explicit time-stepping scheme with a time step of $\text{Δ}t=0.01\text{s}$ was employed for a total simulation time of 2.0 s. Numerical stability was maintained by satisfying the Courant-Friedrichs-Lewy (CFL) condition with the CFL number restricted to 0.5 or less. A prescribed velocity boundary condition was applied at the inlet, while a zero-gradient outflow condition was imposed at the downstream boundary. With g = 9.81 and wave speed $c=\sqrt{gh}=4.43\text{m}/\text{s}$ , CFL = 0.108, which is below 0.5 and therefore satisfies a typical explicit CFL limit.

The condition $u\left(0,y,t\right)=10$ represents a prescribed inflow velocity entering the computational domain. This simulates a continuously supplied flow, such as water entering a channel from an upstream source. The condition $u\left(100,y,t\right)=10$ maintains a consistent downstream flow velocity at the outlet boundary, preventing unrealistic velocity gradients and numerical reflections that may distort the solution.

Using equal inlet and outlet velocities allows the study to focus on changes in flow depth, sediment transport, bed evolution, or pollutant effects without introducing additional disturbances caused by varying boundary inflow and outflow conditions. These boundary values are particularly suitable when investigating the influence of factors such as bed roughness, sediment concentration, or pollutant influx on flow behavior under approximately uniform flow conditions.

The computational domain was defined as a 100 m long and 1 m wide rectangular channel with an initial uniform water depth of 0.5 m. The channel bed was assigned a uniform longitudinal slope of 0.05, with the upstream bed elevation taken as the reference datum of 0 m. Accordingly, the initial bed elevation was defined by $z_b\left(x,0\right)=-0.05$ , resulting in a downstream bed elevation of −5.0 m at x = 100 m. The initial longitudinal flow velocity was set to 1.5 m/s, corresponding to an initial discharge of 0.75 m³/s, while the initial transverse velocity was set to zero.

At the upstream boundary, a constant water depth of 0.5 m and longitudinal velocity of 1.5 m/s were prescribed. The downstream boundary was treated as an open boundary, with the longitudinal and transverse velocity gradients set to zero and the downstream water depth specified where required. The two lateral boundaries were treated as impermeable channel walls, imposing zero transverse velocity. The bed elevation at all boundaries followed the prescribed bed geometry. For a fixed-bed simulation, the bed elevation was maintained throughout the simulation by imposing $\partial z_b/\partial t=0$ . The model was subsequently run for Manning roughness coefficients of 0.03, 0.04, 0.05, and 0.06 to examine their influence on the longitudinal velocity distribution shown in Figure 1.

5. Dimensionless Parameters

Non-dimensional numbers, also known as dimensionless numbers, are fundamental in fluid mechanics and engineering, as they provide insights into the relative importance of different physical quantities and parameters. They are mathematical ratios that express the relationship between various forces, velocities, and dimensions in a system, helping to simplify and generalize the analysis of complex phenomena. According to Icha (2012), non-dimensional numbers are particularly useful for comparing and categorizing flows, such as characterizing the relative influence of inertial forces to viscous forces in a fluid. One of the most well-known dimensionless numbers is the Reynolds number (Re), which is used to predict whether a flow is laminar or turbulent. The Reynolds number is expressed as the ratio of inertial forces to viscous forces and is essential in understanding the flow behavior in different scenarios.

5.1. Reynold Number, Re

Reynolds number (Re) is a dimensionless quantity used in fluid dynamics to characterize the flow of fluids, such as liquids or gases, around objects. It was named after the British engineer Osborne Reynolds, who popularized its use in the late 19th century. The Reynolds number is an essential parameter in understanding the transition between laminar and turbulent flow. The Reynolds number (Re) is defined as the ratio of inertial forces to viscous forces within a fluid flow. It helps predict the flow regime (laminar or turbulent) and provides insight into flow behavior. Reynolds number is given by Equation (20).

$\mathrm{Re}=\frac{\text{Inertia force}}{\text{Viscous force}}=\frac{\rho uL}{\mu }$ (20)

5.2. Froude Number, Fr

The Froude number (Fr) is a dimensionless value representing the ratio of inertia forces to gravitational forces in fluid dynamics, given by Equation (21)

$\text{Fr}=\frac{\text{Inertia force}}{\text{gravitational force}}=\frac{V}{\sqrt{gL}}$ (21)

The Froude number is essential for analyzing open channels, waves, and ship design, and it defines free surfaces: Fr < 1 subcritical, Fr = 1 critical, and Fr > 1 supercritical.

6. Numerical Method

6.1. Continuity Equation in Finite Difference Form

At grid point $\left(i,j\right)$ and time level n, using forward time and central spatial differences

$\frac{h_{i,j}^{n+1}-h_{i,j}^n}{\text{Δ}t}+\frac{{\left(hu\right)}_{i+1,j}^n-{\left(hu\right)}_{i-1,j}^n}{2\text{Δ}x}+\frac{{\left(hv\right)}_{i+1,j}^n-{\left(hv\right)}_{i-1,j}^n}{2\text{Δ}y}=q_{i,j}^n$ (22)

Hence, the depth at the next time step is

$h_{i,j}^{n+1}=h_{i,j}^n-\frac{\Delta t}{2\Delta x}\left[{\left(hu\right)}_{i+1,j}^n-{\left(hu\right)}_{i-1,j}^n\right]-\frac{\Delta t}{2\Delta y}\left[{\left(hv\right)}_{i,j+1}^n-{\left(hv\right)}_{i,j-1}^n\right]+\Delta t{q}_{i,j}^n$ (23)

For q = 0,

$h_{i,j}^{n+1}=h_{i,j}^n-\frac{\Delta t}{2\Delta x}\left[{\left(hu\right)}_{i+1,j}^n-{\left(hu\right)}_{i-1,j}^n\right]-\frac{\Delta t}{2\Delta y}\left[{\left(hv\right)}_{i,j+1}^n-{\left(hv\right)}_{i,j-1}^n\right]$ (24)

u and v are depth-averaged velocities, so the discharge per unit width is $q_x=hu$ , $q_y=hv$ .

6.2. Momentum Equation in Finite Difference Approximation

The finite difference discretization is done on a Cartesian grid with indices $\left(i,j\right)$ and time levels n for the current grid and n + 1 for the next grid. ∆x, ∆t are the grid and time spacing.

${\left(uh\right)}_{i,j}^n$ denotes $u^{\text{*}}{h}^{\text{*}}$ at cell $\left(i,j\right)$ and time n, and face flux values at x and y are denoted as $i+\frac{1}{2}$ and $j+\frac{1}{2}$ , respectively. In discretizing, we utilize the time derivative for the forward difference and the spatially centered derivative.

Time derivative for forward difference for both x and y grids:

$\frac{{\left(uh\right)}_{i,j}^{n+1}-{\left(uh\right)}_{i,j}^n}{\text{Δ}t}$ for x grid

and for the y-grid

$\frac{{\left(vh\right)}_{i,j}^{n+1}-{\left(vh\right)}_{i,j}^n}{\text{Δ}t}$

We also apply spatially centered derivatives.

${\frac{\partial \left(uvh\right)}{\partial x}|}_{i,j}\approx \frac{{\left(uvh\right)}_{i+1,j}-{\left(uvh\right)}_{i-1,j}}{2\text{Δ}x}$

${\frac{\partial \left(vh+\frac{1}{2}{h}^2\right)}{\partial y}|}_{i,j}\approx \frac{{\left(vh+\frac{1}{2}{h}^2\right)}_{i,j+1}-{\left(vh+\frac{1}{2}{h}^2\right)}_{i,j-1}}{2\text{Δ}y}$

Discretizing Equations (4) and (5), respectively, by the finite difference method on a Cartesian grid, we obtain Equations (25) and (26), respectively.

$S_{fx,i,j}^{\text{*}}={\left(n_{b,i,j}^{\text{*}}\right)}^2{h}_{i,j}^{\text{*}-\frac{4}{3}}{u}_{i,j}^{\text{*}}\sqrt{{\left(u_{i,j}^{\text{*}}\right)}^2+{\left(v_{i,j}^{\text{*}}\right)}^2}$ (25)

$S_{fy,i,j}^{\text{*}}={\left(n_{b,i,j}^{\text{*}}\right)}^2{h}_{i,j}^{\text{*}-\frac{4}{3}}{v}_{i,j}^{\text{*}}\sqrt{{\left(u_{i,j}^{\text{*}}\right)}^2+{\left(v_{i,j}^{\text{*}}\right)}^2}$ (26)

The above equations are second-order accurate in space only.

We discretize Equation (2) in the x-direction to obtain Equation (27).

$\begin{array}{l}\frac{{\left(u^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^{n+1}-{\left(u^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^n}{\text{Δ}{t}^{\text{*}}}+\frac{1}{2\text{Δ}{x}^{\text{*}}}\left[{\left(u^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i+1,j}^n-{\left(u^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i-1,j}^n\right]\\ +\frac{1}{2\text{Δ}{y}^{\text{*}}}\left[{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i,j+1}^n-{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i,j-1}^n\right]=-\frac{z_{i+1,j}^{\text{*}}-z_{i-1,j}^{\text{*}}}{2\text{Δ}{x}^{\text{*}}}-S_{fx,i,j}^{\text{*}n}\end{array}$ (27)

Similarly, we discretize Equation (3) in the y-direction to obtain Equation (28).

$\begin{array}{l}\frac{{\left(v^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^{n+1}-{\left(v^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^n}{\text{Δ}{t}^{\text{*}}}-\frac{1}{2\text{Δ}{y}^{\text{*}}}\left[{\left(v^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i,j+1}^n-{\left(v^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i,j-1}^n\right]\\ +\frac{1}{2\text{Δ}{x}^{\text{*}}}\left[{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i+1,j}^n-{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i-1,j}^n\right]=-\frac{z_{i,j+1}^{\text{*}}-z_{i,j-1}^{\text{*}}}{2\text{Δ}{y}^{\text{*}}}-S_{fy,i,j}^{\text{*}n}\end{array}$ (28)

We discretize Equations (6) and (7) for Manning roughness 0.04 to obtain standard Equations (29) and (30).

$S_{fx,i,j}^{\text{*}}=\frac{16}{9}\left[{\left(n_{b,i,j}^{\text{*}}\right)}^2{h}_{i,j}^{\text{*}-\frac{4}{3}}{u}_{i,j}^{\text{*}}\sqrt{{\left(u_{i,j}^{\text{*}}\right)}^2+{\left(v_{i,j}^{\text{*}}\right)}^2}\right]$ (29)

$S_{fy,i,j}^{\text{*}}=\frac{16}{9}\left[{\left(n_{b,i,j}^{\text{*}}\right)}^2{h}_{i,j}^{\text{*}-\frac{4}{3}}{v}_{i,j}^{\text{*}}\sqrt{{\left(u_{i,j}^{\text{*}}\right)}^2+{\left(v_{i,j}^{\text{*}}\right)}^2}\right]$ (30)

Substituting Equation (29) into Equation (12), we obtain Equation (31).

$\begin{array}{l}\frac{{\left(u^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^{n+1}-{\left(u^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^n}{\text{Δ}{t}^{\text{*}}}+\frac{1}{2\text{Δ}{x}^{\text{*}}}\left[{\left(u^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i+1,j}^n-{\left(u^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i-1,j}^n\right]\\ +\frac{1}{2\text{Δ}{y}^{\text{*}}}\left[{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i,j+1}^n-{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i,j-1}^n\right]=-\frac{z_{i+1,j}^{\text{*}}-z_{i-1,j}^{\text{*}}}{2\text{Δ}{x}^{\text{*}}}-\frac{16}{9}{S}_{fx,i,j}^{\text{*}n}\end{array}$ (31)

Similarly, substituting Equation (30) into Equation (13), we obtain Equation (32).

$\begin{array}{l}\frac{{\left(v^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^{n+1}-{\left(v^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^n}{\text{Δ}{t}^{\text{*}}}-\frac{1}{2\text{Δ}{y}^{\text{*}}}\left[{\left(v^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i,j+1}^n-{\left(v^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i,j-1}^n\right]\\ +\frac{1}{2\text{Δ}{x}^{\text{*}}}\left[{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i+1,j}^n-{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i-1,j}^n\right]=-\frac{z_{i,j+1}^{\text{*}}-z_{i,j-1}^{\text{*}}}{2\text{Δ}{y}^{\text{*}}}-\frac{16}{9}{S}_{fy,i,j}^{\text{*}n}\end{array}$ (32)

We discretize Equations (8) and (9) for Manning roughness 0.05 to obtain Equations (33) and (34).

$S_{fx,i,j}^{\text{*}}=\frac{25}{9}\left[{\left(n_{b,i,j}^{\text{*}}\right)}^2{h}_{i,j}^{\text{*}-\frac{4}{3}}{u}_{i,j}^{\text{*}}\sqrt{{\left(u_{i,j}^{\text{*}}\right)}^2+{\left(v_{i,j}^{\text{*}}\right)}^2}\right]$ (33)

$S_{fy,i,j}^{\text{*}}=\frac{25}{9}\left[{\left(n_{b,i,j}^{\text{*}}\right)}^2{h}_{i,j}^{\text{*}-\frac{4}{3}}{v}_{i,j}^{\text{*}}\sqrt{{\left(u_{i,j}^{\text{*}}\right)}^2+{\left(v_{i,j}^{\text{*}}\right)}^2}\right]$ (34)

Substituting Equation (33) into Equation (14), we obtain Equation (35).

$\begin{array}{l}\frac{{\left(u^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^{n+1}-{\left(u^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^n}{\text{Δ}{t}^{\text{*}}}+\frac{1}{2\text{Δ}{x}^{\text{*}}}\left[{\left(u^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i+1,j}^n-{\left(u^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i-1,j}^n\right]\\ +\frac{1}{2\text{Δ}{y}^{\text{*}}}\left[{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i,j+1}^n-{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i,j-1}^n\right]=-\frac{z_{i+1,j}^{\text{*}}-z_{i-1,j}^{\text{*}}}{2\text{Δ}{x}^{\text{*}}}-\frac{25}{9}{S}_{fx,i,j}^{\text{*}n}\end{array}$ (35)

Similarly, substituting Equation (37) into Equation (15), we obtain Equation (36).

$\begin{array}{l}\frac{{\left(v^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^{n+1}-{\left(v^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^n}{\text{Δ}{t}^{\text{*}}}-\frac{1}{2\text{Δ}{y}^{\text{*}}}\left[{\left(v^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i,j+1}^n-{\left(v^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i,j-1}^n\right]\\ +\frac{1}{2\text{Δ}{x}^{\text{*}}}\left[{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i+1,j}^n-{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i-1,j}^n\right]=-\frac{z_{i,j+1}^{\text{*}}-z_{i,j-1}^{\text{*}}}{2\text{Δ}{y}^{\text{*}}}-\frac{25}{9}{S}_{fy,i,j}^{\text{*}n}\end{array}$ (36)

We discretize Equations (10) and (11) for Manning roughness coefficient 0.06 to obtain Equations (37) and (38).

$S_{fx,i,j}^{\text{*}}=4\left[{\left(n_{b,i,j}^{\text{*}}\right)}^2{h}_{i,j}^{\text{*}-\frac{4}{3}}{u}_{i,j}^{\text{*}}\sqrt{{\left(u_{i,j}^{\text{*}}\right)}^2+{\left(v_{i,j}^{\text{*}}\right)}^2}\right]$ (37)

$S_{fy,i,j}^{\text{*}}=4\left[{\left(n_{b,i,j}^{\text{*}}\right)}^2{h}_{i,j}^{\text{*}-\frac{4}{3}}{v}_{i,j}^{\text{*}}\sqrt{{\left(u_{i,j}^{\text{*}}\right)}^2+{\left(v_{i,j}^{\text{*}}\right)}^2}\right]$ (38)

Substituting Equation (37) into Equation (17), we obtain Equation (39).

$\begin{array}{l}\frac{{\left(u^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^{n+1}-{\left(u^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^n}{\text{Δ}{t}^{\text{*}}}+\frac{1}{2\text{Δ}{x}^{\text{*}}}\left[{\left(u^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i+1,j}^n-{\left(u^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i-1,j}^n\right]\\ +\frac{1}{2\text{Δ}{y}^{\text{*}}}\left[{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i,j+1}^n-{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i,j-1}^n\right]=-\frac{z_{i+1,j}^{\text{*}}-z_{i-1,j}^{\text{*}}}{2\text{Δ}{x}^{\text{*}}}-4S_{fx,i,j}^{\text{*}n}\end{array}$ (39)

Similarly, substituting Equation (38) into Equation (17), we obtain Equation (40).

$\begin{array}{l}\frac{{\left(v^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^{n+1}-{\left(v^{\text{*}}{h}^{\text{*}}\right)}_{i,j}^n}{\text{Δ}{t}^{\text{*}}}-\frac{1}{2\text{Δ}{y}^{\text{*}}}\left[{\left(v^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i,j+1}^n-{\left(v^{\text{*}}{h}^{\text{*}}+\frac{1}{2}{h}^{\text{*}2}\right)}_{i,j-1}^n\right]\\ +\frac{1}{2\text{Δ}{x}^{\text{*}}}\left[{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i+1,j}^n-{\left(u^{\text{*}}{v}^{\text{*}}{h}^{\text{*}}\right)}_{i-1,j}^n\right]=-\frac{z_{i,j+1}^{\text{*}}-z_{i,j-1}^{\text{*}}}{2\text{Δ}{y}^{\text{*}}}-4S_{fy,i,j}^{\text{*}n}\end{array}$ (40)

The governing equations were solved numerically using the finite difference method. The numerical simulations were performed using Python.

7. Results and Discussion

7.1. Results and Discussion for Continuity Equation

Flow Depth Update during Each Time Step

Simulating Equation (24) at each time step, the following results are obtained.

Figure 2 shows how the flow depth h is updated from one time step to the next. At time step n, the model starts with a known flow depth h. The continuity equation is then used to calculate the change in water stored in the channel. The curve illustrates how flow depth increases with an increase in time step.

Figure 2. Flow depth vs time.

7.2. Results and Discussion for Momentum Equation

7.2.1. Comparison of Average Friction Slope for Different Manning Roughness Coefficients

Equation (27) was numerically simulated for a range of Manning roughness coefficients, and the resulting average friction slopes were evaluated on a common computational grid. A comparative analysis of the average friction slope corresponding to each roughness coefficient yielded the results presented below.

Figure 3. Average friction for different Manning roughness.

Figure 3 presents a bar chart illustrating the relationship between the Manning roughness coefficient and the average friction slope. The results show a consistent increase in average friction slope as the Manning roughness coefficient increases from 0.03 to 0.06. In particular, the highest roughness coefficient (n = 0.06) corresponds to the greatest average friction slope, indicating increased resistance to flow. This trend demonstrates the strong influence of channel roughness on frictional losses within the flow system. Since elevated pollutant concentrations contribute to higher Manning roughness values, the results suggest that pollutant-induced increases in channel roughness lead to greater average friction slopes and consequently enhanced energy dissipation within the flow.

7.2.2. Variation of Velocity with Time at Three Fixed Observation Points for Different Manning Roughness Coefficients

By numerically simulating Equation (27) at three fixed spatial locations, (x = 200), (x = 600), and (x = 800), for different Manning roughness coefficients, the temporal evolution of velocity was investigated. The simulations were performed on various computational grids to facilitate a comparative analysis of the influence of grid resolution and channel roughness on velocity variation with time. The corresponding results are presented below.

Figure 4. Velocity vs time for three different points.

Figure 4 illustrates the temporal variation of flow velocity at three fixed spatial locations (x = 200), (x = 600), and (x = 800), for different Manning roughness coefficients. At all three observation points, a similar trend is observed, whereby velocity decreases with time for Manning roughness values of 0.03, 0.04, 0.05, and 0.06. Furthermore, an increase in the Manning roughness coefficient results in a corresponding reduction in flow velocity, indicating the significant influence of channel roughness on flow resistance.

Among the roughness values considered, a Manning roughness coefficient of 0.06 consistently produces the lowest velocity throughout the simulation period at all three locations. The similarity of the velocity profiles at (x = 200), (x = 600), and (x = 800) suggests that the effect of Manning roughness on flow velocity is spatially consistent across the computational domain. Consequently, the velocity-distance relationship is expected to exhibit a similar pattern, characterized by lower velocities for higher Manning roughness coefficients.

7.2.3. Effects of Distance on Velocity for Different Manning Roughness Coefficients

Simulating Equation (27) for velocity with distance for various Manning roughness coefficients, the following results are obtained.

Figure 5. Velocity vs distance.

Figure 5 illustrates the variation of flow velocity with distance for different Manning roughness coefficients. The results show a progressive decrease in velocity as the flow advances downstream. A comparison of the different Manning roughness coefficients reveals that an increase in Manning’s roughness coefficient leads to a reduction in flow velocity. Specifically, a roughness coefficient of n = 0.06 produces lower velocities than a coefficient of n = 0.03. This trend indicates an inverse relationship between flow velocity and channel roughness, whereby lower roughness values correspond to higher velocities and higher roughness values correspond to lower velocities.

7.2.4. Effects of Manning Roughness Coefficient on Water Depth

Simulating Equation (27) for water depth with distance for various Manning roughness coefficients, the following results are obtained.

Figure 6. Water depth vs distance.

Figure 6 illustrates the variation of water depth with distance for different Manning roughness coefficients. The results indicate that the flow depth increases progressively along the channel length. In the present study, it is assumed that the bed slope exceeds the friction slope, resulting in subcritical flow conditions ($F_r<1$ ). Under such conditions, the flow depth is not constant but adjusts continuously in response to the balance between gravitational forces, flow resistance, and energy dissipation along the channel. The results further demonstrate that flow depth increases with increasing Manning roughness coefficient. Specifically, a roughness coefficient of n = 0.06 produces greater flow depths than a coefficient of n = 0.03.

8. Conclusions

To satisfy the continuity of flow and maintain the discharge, the water depth adjusts by increasing, thereby compensating for the loss in flow velocity. Consequently, rougher channels tend to sustain higher flow depths than smoother channels. The interaction between flow depth and channel roughness has important implications for hydraulic engineering design and river management. In the design of artificial channels, both parameters must be carefully considered to ensure that the desired flow capacity and velocity are achieved. Higher channel roughness may require greater flow depths to convey the same discharge effectively. In natural rivers and streams, roughness varies spatially due to factors such as bed material composition, sediment accumulation, vegetation growth, and channel irregularities.

The graph of flow depth in the continuity equation shows that the flow depth is updated progressively at each time step, with the depth changing toward a stable value as the simulation advances. This indicates that the continuity equation is being applied consistently, allowing the model to account for the variation in water volume over time. Overall, the results demonstrate a stable and physically consistent flow-depth update process during the simulation.

The relationship between flow velocity and the Manning roughness coefficient is of considerable hydraulic significance, as the roughness coefficient directly influences flow resistance and, consequently, the overall flow characteristics within an open channel. The numerical simulations demonstrate that an increase in the Manning roughness coefficient leads to a reduction in flow velocity and promotes greater sediment deposition along the channel bed. This behavior arises because higher roughness values increase the resistance exerted on the flowing water, resulting in greater energy dissipation through friction.

The observed results are consistent with Manning’s equation, which predicts an inverse relationship between flow velocity and the roughness coefficient. As the Manning roughness coefficient increases, the frictional resistance along the channel boundary intensifies, causing a decrease in flow velocity. The reduction in flow velocity lowers the sediment transport capacity of the flow, thereby enhancing the tendency for sediment particles to settle and accumulate on the bed.

Flow depth and Manning’s roughness coefficient are closely interrelated hydraulic parameters that govern flow behavior in open channels. Flow depth influences the hydraulic radius, which, together with the Manning roughness coefficient, determines the flow velocity as described by Manning’s equation. In general, an increase in flow depth enhances the hydraulic radius and tends to increase flow velocity. However, the extent of this increase is moderated by the channel roughness, which imposes resistance to flow and dissipates energy through frictional effects. The numerical results indicate that channels with higher Manning roughness coefficients exhibit greater flow depths than channels with lower roughness coefficients under similar hydraulic conditions. This behavior can be attributed to the increased hydraulic resistance associated with larger roughness values. As roughness increases, frictional losses become more significant, leading to a reduction in flow velocity.

Acknowledgements

The authors would like to thank Chuka University for their support of this research.

Nomenclature

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this paper.

References