Hydrodynamic Removal of Sediments in Open Irrigation Channel

Abstract

ABSTRACT

The enhancement of irrigation infrastructure is essential for boosting agricultural productivity, particularly in developing nations. Nevertheless, the accumulation of sediment in irrigation channels presents a major obstacle, reducing water flow efficiency and escalating maintenance expenses. This study introduces a numerical model designed for the hydrodynamic removal of bed sediment in an open irrigation channel located in Guanajuato, Mexico. Fluid dynamics are resolved using a finite difference method to simulate three‐dimensional velocity fields, while sediment transport is analysed through the particle‐in‐cell method (PICM). The model incorporates key factors influencing particle behaviour, including velocity field, turbulent dispersion, and particle‐specific properties such as size, shape, and settling velocity. Validation experiments were performed in a laboratory setting, utilizing a 1:20 scale irrigation channel with three distinct configurations of submerged structures. The results confirmed the model's capability to predict sediment transport, pinpoint erosion and deposition zones, and evaluate alterations in the hydraulic cross‐section. These outcomes underscore the critical role of effective sediment management in ensuring the optimal performance and maintenance of irrigation channels.

Introduction

Abrupt water level rises in irrigation canals often result from intense or excessive rainfall over short durations, ranking among the most hazardous weather-related natural disasters globally. Such events pose substantial threats to human lives and infrastructure.

To regulate water flow in irrigation canals, flow control gates and related structures are typically installed. These gates generally facilitate unobstructed flow acceleration as water exits regulation works, bridges, or gauging zones [1]. However, sediment deposition frequently occurs in regions where flow velocity diminishes, such as elongated canal sections or areas with directional changes [2]. This necessitates regular external maintenance measures, including dredging or sediment removal, which are time-intensive and costly.

Research on efficient sediment removal methods in irrigation canals is limited [3–5]. Sediment deposition can alter canal bed characteristics, leading to downstream erosion and scouring. Some studies have examined hydrodynamic behaviours near gates and surface structures and employed numerical models to simulate channel bottom scouring [6, 7]. Tools like FLOW 3D v11.1.0 have been utilized to analyse hydrodynamic conditions downstream of structures under various design scenarios. These simulations provide accurate insights into velocity distributions, maximum scour depths, and water levels under specified conditions [8, 9].

Improperly designed flow-regulating structures can compromise upstream and downstream hydraulic systems, potentially resulting in structural failures [10, 11]. Experimental investigations into flow dynamics through gates and control structures have been complemented by open-source software such as OpenFOAM. This software has been applied to simulate flow through gates under diverse flow and blockage conditions, with validations highlighting the significant influence of blockage scenarios on flow distribution during drawdowns. For instance, extending the lock length alters velocity profiles compared to gates alone, with higher blockage rates causing more pronounced flow redistribution [12].

Numerical modelling of sediment transport is a critical tool in understanding and predicting sediment dynamics in various water bodies, such as rivers, channels, and coastal areas. These models help address issues such as erosion, reservoir silting, and channel bed changes, which are vital for the management and development of hydraulic structures. The models vary in complexity, from two-dimensional to three-dimensional, and employ different computational techniques to simulate sediment transport processes. Below are key aspects of numerical modelling in sediment transport.

Two-dimensional models, such as the one developed by Hansen, simulate sediment transport by dividing it into suspended load and bed load components. These models use a fractional step approach to solve the advection-diffusion equation, employing high-resolution algorithms for advection and semi-implicit schemes for diffusion. Parallel computing enhances their efficiency, making them suitable for large-scale simulations [13].

Three-dimensional models provide detailed simulations of flow and sediment transport, accounting for complex interactions between water flow and sediment dynamics. These models solve governing differential equations using methods such as the finite volume method on unstructured grids, allowing for accurate representation of sediment transport in curved channels and river systems [14].

Advanced numerical models incorporate optimization algorithms, such as chicken swarm optimization, to improve prediction accuracy. These models achieve high accuracy rates, with significant improvements over traditional techniques, as demonstrated by Kumar et al., who reported an accuracy of 96.87% and a low RMSE value of 0.12 [15].

Despite advancements, challenges remain in accurately predicting processes such as scouring around hydraulic structures. Models like SedFoam, which use a two-fluid approach, offer detailed simulations but are computationally expensive. There is a need for intermediate-scale models that balance accuracy and computational cost [16].

Numerical models are indispensable for detailed sediment transport analysis; empirical and soft computing methods offer less computationally intensive alternatives for cumulative sediment load predictions. These methods, however, may oversimplify complex processes, highlighting the need for continued research and development in numerical modelling techniques [17].

Other studies have explored the hydrodynamics of vegetation in irrigation canals using three-dimensional (3D) computational fluid dynamics (CFD) models to simulate flows with varying roughness [18, 19]. While these models replicate key flow characteristics, discrepancies in velocity profiles often arise due to grid resolution and wall boundary configurations. Nonetheless, they effectively highlight roughness-induced low-velocity zones and emphasize the importance of hybrid approaches combining numerical and experimental methods for optimizing submerged structure designs.

Comprehensive validation of CFD models requires detailed experimental datasets for systematic comparisons between numerical and physical observations. Although submerged structures are not widely employed in Mexico, they are valued for their hydraulic capacity and ability to support large embankments, making them suitable for irrigation systems in low-resistance soil regions [20]. Misaligned surface or submerged structures, where flow paths intersect at skewed angles, further complicate flow dynamics.

This study investigates flow behaviour in open irrigation canals with submerged structures designed to enhance bottom velocity in specific sections. The primary aim is to evaluate the feasibility of natural hydrodynamic dredging to remove sediment from silting-prone areas and ensure efficient water delivery to irrigation modules. Hydrodynamic simulations reveal stable and time-converged velocity and turbulence fields [21], demonstrating that while these parameters remain temporally steady, they exhibit significant spatial variability.

Materials and Methods

Physical Laboratory Model

The experimental channel utilized in the study (Figure 1) is characterized by a variable slope and a constant rectangular cross-section, measuring 2.32 m in length, 0.25 m in width, and 0.18 m in height. The system includes a 22-L feed tank, a 2-horsepower pump, and a 90-L storage tank.

[IMAGE OMITTED. SEE PDF]

The channel walls are constructed from 0.012 m thick acrylic sheets, while the gate is made from the same material with a thickness of 0.006 m. The gate dimensions are 0.24 m in width and 0.15 m in height, positioned at the mid-length of the channel to mitigate oscillations induced by the inlet and outlet weirs. The gate assembly (Figure 2), including the plastic chains for servo motor attachment during opening and closing operations, weighs 0.346 kg.

[IMAGE OMITTED. SEE PDF]

A stainless steel weir with a sharp crest is installed at the channel inlet, featuring a height of 0.076 m and a base width of 0.25 m, to regulate flow from the feed tank. At the channel's outlet to the storage tank, flow can be discharged with or without a spillway. The channel bed is fabricated from 0.006 m thick stainless steel sheets, with a roughness coefficient of 0.020 assigned.

In the physical scale experiments conducted on the open channel, the system operated at a flow rate of 0.0019 m³/s, with a slope of 0.002. The experimental setup included three submerged structures of identical shape but varying widths (Figure 3).

[IMAGE OMITTED. SEE PDF]

The Figure 4 shows the Ls5 configuration submerged structure inside the laboratory channel; this configuration was the first to be proved under flow conditions and initial conditions on sediment bottom concentration.

[IMAGE OMITTED. SEE PDF]

Governing Equations

The velocity fields and free surface variations were modelled using the shallow-water equations [20]. To simplify these equations, two primary assumptions were applied. The first is the hydrostatic approximation, which assumes that the horizontal length scale is significantly greater than the vertical length scale. This condition, essential for the shallow-water hypothesis, relates pressure to the reference density, gravitational acceleration, and the reference vertical plane. The second assumption is the Boussinesq approximation, which treats density as constant in all momentum conservation terms except the gravitational term.

By incorporating these approximations and applying turbulent averaging, the governing equations for the velocity components (Equations 1–3) are derived [20]. 1 $$\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }t}+\overline{u}\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }x}+\overline{v}\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }y}+\overline{w}\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }z}=-\frac{1}{{\rho }_{0}}\frac{\mathit{\partial }\overline{P}}{\mathit{\partial }x}-\frac{\mathit{\partial }\overline{{u}^{\mathit{\prime }}{u}^{\mathit{\prime }}}}{\mathit{\partial }x}-\frac{\mathit{\partial }\overline{{u}^{\mathit{\prime }}{v}^{\mathit{\prime }}}}{\mathit{\partial }y}-\frac{\mathit{\partial }\overline{{u}^{\mathit{\prime }}{w}^{\mathit{\prime }}}}{\mathit{\partial }z}+f\overline{v}$$ 2 $$\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }t}+\overline{u}\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }x}+\overline{v}\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }y}+\overline{w}\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }z}=-\frac{1}{{\rho }_{0}}\frac{\mathit{\partial }\overline{P}}{\mathit{\partial }y}-\frac{\mathit{\partial }\overline{{v}^{\mathit{\prime }}{u}^{\mathit{\prime }}}}{\mathit{\partial }x}-\frac{\mathit{\partial }\overline{{v}^{\mathit{\prime }}{v}^{\mathit{\prime }}}}{\mathit{\partial }y}-\frac{\mathit{\partial }\overline{{v}^{\mathit{\prime }}{w}^{\mathit{\prime }}}}{\mathit{\partial }z}-f\overline{u}$$ 3 $$\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }x}+\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }y}+\frac{\mathit{\partial }\overline{w}}{\mathit{\partial }z}=0$$ Here, $$t$$ represents time, $$\overline{u}$$, $$\overline{v}$$ and $$\overline{w}$$ are the time-averaged velocity components, and $${u}^{\prime }$$, $${v}^{\prime }$$ and $${w}^{\prime }$$ are their fluctuating counterparts in the $$x$$ -, $$y$$ - and $$z$$ -directions, respectively. The Coriolis parameter $$\mathit{f}=2{\Omega }\,\mathit{sin}\,\phi $$, $${\Omega }$$ is the Earth's rotation rate and $$\phi $$ is the latitude, is also included in the model. Integrating the continuity equation [20] over the water depth and applying a kinematic boundary condition at the free surface leads to the free surface equation (Equation 4). 4 $$\frac{\mathit{\partial }\overline{\eta }}{\mathit{\partial }t}+\frac{\mathit{\partial }}{\mathit{\partial }x}\left(\int \nolimits_{-h}^{\eta }\overline{u}dz\right)+\frac{\mathit{\partial }}{\mathit{\partial }y}\left(\int \nolimits_{-h}^{\eta }\overline{v}dz\right)=0$$ where the water depth $$h$$ and the time-averaged water surface elevation $$\overline{\eta }$$ are measured from a reference level or from the undisturbed water surface.

The hydrostatic condition (Equation 5) gives the following definition of pressure [20]: 5 $$\overline{P}=g\int \nolimits_{z}^{\eta }\rho dz+{P}_{\mathrm{atm}}$$ where $${P}_{\text{atm}}$$ is the atmospheric pressure. By applying the Leibnitz integration rule in (5), the pressure terms of Equations (1) and (2) can be written as [22] in Equations (6) and (7): 6 $${-}\frac{1}{{\rho }_{0}}\frac{\mathit{\partial }\overline{P}}{\mathit{\partial }x}=-\frac{\rho g}{{\rho }_{0}}\frac{\mathit{\partial }\overline{\eta }}{\mathit{\partial }x}-\frac{g}{{\rho }_{0}}\int \nolimits_{z}^{\eta }\frac{\mathit{\partial }{\rho }^{\mathit{\prime }}}{\mathit{\partial }x}dz-\frac{1}{{\rho }_{0}}\frac{\mathit{\partial }{P}_{\text{atm}}}{\mathit{\partial }x}$$ 7 $${-}\frac{1}{{\rho }_{0}}\frac{\mathit{\partial }\overline{P}}{\mathit{\partial }y}=-\frac{\rho g}{{\rho }_{0}}\frac{\mathit{\partial }\overline{\eta }}{\mathit{\partial }y}-\frac{g}{{\rho }_{0}}\int \nolimits_{z}^{\eta }\frac{\mathit{\partial }{\rho }^{\mathit{\prime }}}{\mathit{\partial }y}dz-\frac{1}{{\rho }_{0}}\frac{\mathit{\partial }{P}_{\text{atm}}}{\mathit{\partial }y}$$

The anomalous density $${\rho }^{\prime }=\rho -{\rho }_{0}$$ plays a role in describing the horizontal pressure gradients (6) and (7). These gradients arise from free surface variations (barotropic term), horizontal density differences (baroclinic term), and atmospheric pressure effects. When density is constant and atmospheric pressure contributions are negligible, the gradients simplify to the barotropic term alone.

A state equation (Equation 8) is used to compute density variations and is defined as [23]: 8 $$\rho \left(\overline{P},\overline{T},\overline{S}\right)=\frac{{\rho }_{0}}{(1-\frac{\overline{P}}{{k}_{P}})},$$ where $${k}_{P}$$ is a coefficient function of pressure, temperature and salinity. The formulation to compute the value of $${k}_{P}$$ can be found in [23].

In summary, the hydrodynamic system is defined by four variables, $$\overline{u},\overline{v},\overline{w}$$ and $$\overline{\eta }$$, computed through Equations (1–4), while thermodynamic variables $$\rho ,\overline{T}$$ and $$\overline{S}$$, and are determined using Equations (6–8). The pressure in Equation (8) is hydrostatic and can be evaluated at any time [20].

Turbulence Modelling

The Reynolds stress correlations $$\overline{{u}^{\prime }{u}^{\prime }},\overline{{u}^{\prime }{v}^{\prime }},\overline{{u}^{\prime }{w}^{\prime }},...,\overline{{v}^{\prime }{w}^{\prime }}$$ in Equations (1) and (2) are computed using a zero-equation turbulence model. These models approximate the correlations by expressing them as the product of the kinematic eddy viscosity, denoted as $${\nu }_{E}$$, and the mean strain-rate tensor, as presented in Equations (9) and (10): 9 $${-}\frac{\mathit{\partial }\overline{{u}^{\mathit{\prime }}{u}^{\mathit{\prime }}}}{\mathit{\partial }x}-\frac{\mathit{\partial }\overline{{u}^{\mathit{\prime }}{v}^{\mathit{\prime }}}}{\mathit{\partial }y}-\frac{\mathit{\partial }\overline{{u}^{\mathit{\prime }}{w}^{\mathit{\prime }}}}{\mathit{\partial }z}=\frac{\mathit{\partial }}{\mathit{\partial }x}\left(2{\nu }_{E}\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }x}\right)+\frac{\mathit{\partial }}{\mathit{\partial }y}\left[{\nu }_{E}\left(\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }y}+\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }x}\right)\right]+\frac{\mathit{\partial }}{\mathit{\partial }z}\left({\nu }_{E}\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }z}\right)$$ 10 $${-}\frac{\mathit{\partial }\overline{{v}^{\mathit{\prime }}{u}^{\mathit{\prime }}}}{\mathit{\partial }x}-\frac{\mathit{\partial }\overline{{v}^{\mathit{\prime }}{v}^{\mathit{\prime }}}}{\mathit{\partial }y}-\frac{\mathit{\partial }\overline{{v}^{\mathit{\prime }}{w}^{\mathit{\prime }}}}{\mathit{\partial }z}=\frac{\mathit{\partial }}{\mathit{\partial }y}\left(2{\nu }_{E}\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }y}\right)+\frac{\mathit{\partial }}{\mathit{\partial }x}\left[{\nu }_{E}\left(\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }y}+\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }x}\right)\right]+\frac{\mathit{\partial }}{\mathit{\partial }z}\left({\nu }_{E}\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }z}\right)$$

The kinematic eddy viscosity $${\nu }_{E}$$ is composed of both turbulent and molecular components, such that $${\nu }_{E}={\nu }_{t}+\nu $$. The turbulent viscosity coefficient $${\nu }_{t}$$ is determined using a mixing length model [19], which can be expressed by Equation (11): 11 $${\nu }_{t}={\left({l}_{h}^{4}\left[2{\left(\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }x}\right)}^{2}+2{\left(\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }y}\right)}^{2}+{\left(\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }x}+\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }y}\right)}^{2}\right]+{l}_{v}^{4}\left[{\left(\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }z}\right)}^{2}+{\left(\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }z}\right)}^{2}\right]\right)}^{1/2}$$

The vertical length scale $${l}_{v}=\kappa \left(z-{z}_{b}\right)$$ for $$\left(z-{z}_{b}\right)/\delta < \lambda /\kappa $$ and $${l}_{v}=\lambda \delta $$ for $$\lambda /\kappa < \left(z-{z}_{b}\right)/\delta < 1$$. Here, $$\kappa $$ represents the von Kármán constant, typically valued at 0.41, $$\left(z-{z}_{b}\right)$$ is the distance from the boundary wall, $$\delta $$ denotes the boundary-layer thickness, and $$\lambda $$ is an empirical constant, usually set to 0.09. For shallow-water flows subjected to steady currents, the boundary-layer thickness is often approximated by the water depth, $$h$$ [24].

The horizontal length scale, $${l}_{h}$$, generally differs from the vertical length scale. A common approach assumes proportionality, defined $${l}_{h}=\beta {l}_{v}$$, where the proportionality constant $$\beta $$ is determined experimentally [24]. In scenarios involving parallel or nearly parallel flows, the eddy viscosity simplifies to its classical boundary-layer representation. If $${l}_{h}={l}_{v}$$, it transitions to its appropriate three-dimensional form, assuming vertical velocity is negligible.

Free Surface and Bottom Conditions

The shear stress at the free surface is specified by the prescribed wind stresses $${\tau }_{\omega x}$$ and $${\tau }_{\omega y}$$ Equations (12) and (13) [20]: 12 $${\nu }_{t}{\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }z}\vert }_{\text{surface}}={\tau }_{\omega x}$$ 13 $${\nu }_{t}{\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }z}\vert }_{\text{surface}}={\tau }_{\omega y}$$ where $${\tau }_{\omega x}={C}_{\omega }{\rho }_{a}{\overline{\omega }}_{x}\vert {\overline{\omega }}_{x}\vert $$ and $${\tau }_{\omega y}={C}_{\omega }{\rho }_{a}{\overline{\omega }}_{y}\vert {\overline{\omega }}_{y}\vert $$, $${\rho }_{a}$$ is the air density, and $${\overline{\omega }}_{x}$$ and $${\overline{\omega }}_{y}$$ are the mean wind velocity components. The wind-drag coefficient $${C}_{\omega }$$ is computed by the Garratt's formula $${C}_{\omega }=\left(0.75+0.067\overline{\omega }\right)\times 1{0}^{-3}$$ [25], where $$\overline{\omega }$$ is the main magnitude of wind velocity vector. In practice, the wind velocity is typically evaluated at a height of 10 m above the free surface.

The bottom shear stress is expressed using the Manning-Chezy formulation, which relates it to the velocity components through 14 $${\nu }_{t}{\frac{\mathit{\partial }\overline{u}}{\mathit{\partial }z}\vert }_{\text{bottom}}=\frac{g\sqrt{{\overline{u}}^{2}+{\overline{v}}^{2}}}{C{z}^{2}}\overline{u}$$ 15 $${\nu }_{t}{\frac{\mathit{\partial }\overline{v}}{\mathit{\partial }z}\vert }_{\text{bottom}}=\frac{g\sqrt{{\overline{u}}^{2}+{\overline{v}}^{2}}}{C{z}^{2}}\overline{v}$$ where $$Cz$$ is the Chezy friction coefficient. The velocity components are derived from values corresponding to the layer directly adjacent to the sediment-water interface.

The Governing Equation Sediment Particle Model (PICM)

The numerical model for particle transport is given under a Lagrangian approach; the particles are placed following an exponential law of concentrations or by an initial position in three-dimensional space [26]. For the movement of particles, a stochastic model is considered and discretized in three dimensions (Figure 5), considering the specific weight of each particle as well as the fall velocity of the same [27] and it is verified if these are within the domain study for a single time step (Dt) from (n) to (n + 1) is given by: 16 $\begin{align*}{x}_{i}^{n+1}={x}_{i}^{n}+{u}_{i,j,k}({\increment}t)\mathit{\pm }(2\text{rand}(\text{iseed})-0.5)\sqrt{\left(2{\upsilon t}_{i,j,k}{\increment}t\right)}\\ {y}_{i}^{n+1}={y}_{i}^{n}+{v}_{i,j,k}({\increment}t)\mathit{\pm }(2\text{rand}(\text{iseed})-0.5)\sqrt{\left(2{\upsilon t}_{i,j,k}{\increment}t\right)}\\ {z}_{i}^{n+1}={z}_{i}^{n}+{w}_{i,j,k}({\increment}t)\mathit{\pm }(2\text{rand}(\text{iseed})-0.5)\sqrt{\left(2{\upsilon t}_{i,j,k}{\increment}t-{w}_{s}{\increment}t\right)}\end{align*}$

[IMAGE OMITTED. SEE PDF]

At time t, the particle's position is represented by $${x}_{i}^{n},{y}_{i}^{n},{z}_{i}^{n}$$, while (u, v, w) denote the average velocities in the (i, j, k) directions, respectively. The parameter ($${\nu }_{t}$$) refers to the turbulent viscosity coefficient, (Dt) represents the Lagrangian time step, and $${w}_{s}$$ corresponds to the sediment settling velocity. Particle tracking is governed by Equation (16), where each particle undergoes a spatial displacement of magnitude $$\pm (2\text{rand}(\text{iseed})-0.5)\sqrt{\left(2{\upsilon t}_{i,j,k}{\increment}t\right)}$$, in any direction within the domain. The displacement's sign, positive or negative, depends on the particle's spatial orientation. At each time step, the velocity field exerts influence on each particle, resulting in motion that aligns with the dominant flow direction dictated by the velocity fields. The turbulent viscosity coefficient ($${\nu }_{t}$$) is distributed throughout the domain as a scalar field, encapsulating information about turbulent intensities.

The shear forces present in a turbulent flow, varying along its depth (z), can be expressed as follows: 17 $${\tau }_{i}=\rho \nu \frac{dU}{dz}-\rho \overline{{u}_{i}w}$$ where ($${\tau }_{i}$$) represents the bottom shear stress, ($$\rho $$) is the fluid density, ($$\upsilon $$) is the kinematic viscosity coefficient, (U) is the mean fluid velocity, and ($$\overline{{u}_{i}w}$$) denotes the near-bottom Reynolds stress term.

The critical shear stress acting on particles is defined by the following equation: 18 $${\tau }_{\text{critical}}=0.03\left({\rho }_{s}-\rho \right)g{d}_{50}$$ where ($${\rho }_{s}$$) is the solid density, ($$\rho $$) is the water density, (g) is the gravitational acceleration, and (d50) is the median particle diameter (50th percentile).

The probability function for the deposition of the particles is determined by the following equation: 19 $${P}_{\text{fall}}=\left\{\begin{array}{@{}c@{}}0\\ \left(1-\frac{{\tau }_{x,y}}{{\tau }_{\text{critical}}}\right)\end{array}\right\},\begin{array}{c}{\tau }_{x,y}\mathit{\ge }\text{critical}\\ {\tau }_{x,y}< \text{critical}\end{array}$$ and the probability function for the resuspension of the particles is established in the following way: 20 $${P}_{\text{resuspension}}=\left\{\begin{array}{@{}c@{}}0\\ \left(1-\frac{{\tau }_{\text{critical}}}{{\tau }_{x,y}}\right)\end{array}\right\},\begin{array}{c}{\tau }_{x,y}\mathit{\le }\text{critical}\\ {\tau }_{x,y} > \text{critical}\end{array}$$

Velocity Fall of Sediment

To determine the velocity fall of sediment particles, these are considered to have non-spherical shape, so the effect of the shapes has a considerable influence on their velocity, mainly on relatively large particles (> 300 μm), the expressions that determine the magnitude of velocity fall [28] are expressed below. 21 $\begin{align*}{w}_{s}=\frac{(S-1)g{d}^{2}}{0.8\nu };\ 1< d\mathit{\le }100{\upmu }\mathrm{m}\\ {w}_{s}=\frac{10\nu }{d}\left[{\left(1+0.01\frac{(S-1)g{d}^{3}}{{\nu }^{2}}\right)}^{0.5}-1\right];100< d\mathit{\le }1000{\upmu }\mathrm{m}\,\\ {w}_{s}=1.1{\left((S-1)g{d}^{2}\right)}^{0.5};\quad d > 1000{\upmu }\mathrm{m}\end{align*}$ where (d) diameter of particle, (S) specific gravity, ($$\upsilon $$) kinematic viscosity coefficient y (g) gravitational constant.

Hydraulic Similarity

To evaluate the reliability and reproducibility of the hydrodynamic behaviour observed in the laboratory flume, a comprehensive cross-validation was conducted between experimental observations and numerical simulations. The physical model employed a 1:10 geometric scale in a variable-slope channel, allowing direct visualization and measurement of velocity fields and sediment transport patterns under controlled boundary conditions. This experimental setup aimed to replicate key flow phenomena, including turbulence intensities and sediment particle dispersion, encountered in full-scale irrigation canals.

Hydraulic similarity between numerical and physical models was ensured through the application of dimensionless criteria, such as Froude similarity, which preserves the dynamic relationship between inertial and gravitational forces governing open-channel flow. The use of consistent Reynolds and Froude numbers across both models minimized scale distortion and allowed meaningful comparisons of flow structure and sediment behaviour in regions of concern, especially downstream of submerged structures.

Numerical simulations were executed using a custom CFD model developed in MATLAB, incorporating identical geometry, boundary conditions, and sediment properties as those used in the laboratory tests. The results from both platforms were compared using statistical performance indicators, including Absolute Error (EA), Mean Absolute Error (MAE), and Mean Squared Error (MSE). These metrics provided a quantitative assessment of the deviation between the simulated and measured velocity values, confirming that the CFD model captured the essential flow dynamics with acceptable accuracy.

Notably, velocity profiles generated from the numerical simulations closely matched those measured experimentally, particularly in regions of flow contraction and acceleration. Additionally, spatial patterns of turbulent dispersion and sediment concentration aligned with physical observations, confirming that both modelling approaches were capable of characterizing the hydrodynamics near submerged elements. This high degree of hydraulic similarity supports the integration of numerical tools with laboratory data for the robust design of sediment control strategies in irrigation channels.

Validation and Analysis Methods

To ensure the accuracy and reliability of the proposed numerical and experimental modelling framework, a cross-validation process was implemented by comparing two independent datasets: (i) numerical simulations conducted with MATLAB, and (ii) experimental data obtained from scaled physical models tested in a variable-slope channel. This validation enabled the systematic evaluation of velocity fields and hydrodynamic thrust under controlled and reproducible conditions.

Cross-Validation Procedure

The validation framework was structured as follows:

• Numerical Model Comparison: CFD simulations using MATLAB were conducted for identical geometrical and boundary conditions. The numerical model incorporated geometrical configuration, water depths, and flow rate. Velocity profiles were calculated and used to estimate the corresponding hydrodynamic particle dispersion.

• Experimental Scale Model Testing: Physical simulations were carried out using 1:10 scaled hydraulic channels. Flow velocities were measured under different geometrical configurations, and the resulting data were used to reproduce realistic sediment particle movement scenarios.

Quantitative Comparison Metrics for Validation

The agreement between numerical simulated and measured data was evaluated using the following statistical performance indicators:

• Absolute Error (EA): To determine, on average, how much the predictions deviate from the actual value. It is easy to interpret because it is in the same units as the data.

22 $$EA=\vert {V}_{\text{measured}}-{V}_{\text{simulated}}\vert $$

• Mean Absolute Error (MAE): To determine, on average, how much the predictions deviate from the actual value. It is easy to interpret because it is in the same units as the data.

23 $$\text{MAE}=\frac{1}{n}\sum\limits _{i=1}^{n}\vert {V}_{\text{measured},i}-{V}_{\text{simulated},i}\vert $$

• Mean Squared Error (MSE): To penalize large errors more heavily. It's useful when you want large errors to be even less frequent.

24 $$\text{MSE}=\frac{1}{n}\sum\limits _{i=1}^{n}{\left({V}_{\text{measured},i}-{V}_{\text{simulated},i}\right)}^{2}$$

• Root Mean Square Error (RMSE): quantifies the average magnitude of error between observed (sensor-acquired) and predicted (modelled) velocities.

25 $$\text{RMSE}=\sqrt{\text{MSE}}$$

• Percentage Error (ε%): reflects the normalized difference between the average simulated and measured velocities.

26 $$\varepsilon \mathit{\%}=\left(\frac{\vert {V}_{\text{measured}}-{V}_{\text{simulated}}\vert }{{V}_{\text{measured}}}\right)\ast 100$$

• Willmott's Concordance Index (d-index): To evaluate the overall concordance of the model.

$$\text{d\mbox{-}index}=1\to \text{perfect agreement}$$ $$\text{d\mbox{-}index}=0\to \text{There is no agreement}$$

In order to evaluate the accuracy of the developed numerical model, the flow rates obtained through simulation were compared with those measured experimentally using a flow metre installed in a test channel. For this validation, three statistical metrics were employed: root mean square error (RMSE), Willmott's Concordance Index (d-index), and mean absolute percentage error (ε %).

Multivariate ANOVA Analysis

The objective of the ANOVA analysis is to assess whether there are significant differences in sediment concentration based on three categorical factors:

• Channel hydraulic configuration (Current, L3, L4, L5)

• Measurement points (P1 to P4)

• Sediment diameter type (D15, D50, D90)

The analysis was performed in MATLAB for a three-way ANOVA with full interaction on sediment concentration data, considering hydraulic configuration, measurement point, and diameter type as factors. Post hoc multiple comparisons between configurations are then performed, and visualizations are generated to support statistical interpretation. This approach allows for detecting whether channel design, sampling position, or sediment size significantly influences the concentration of transported particles.

Irrigation Channel Ing. Antonio Coria, Guanajuato Mexico

The Ing. Antonio Coria irrigation channel is an important waterway located in the city of Guanajuato, Mexico. This canal was built during the colonial era to help supply water to the city and operate the hydraulic systems for mining. Currently, the canal is used for agriculture, especially as an irrigation canal that provides water to irrigation modules in different municipalities of the State of Guanajuato.

The Antonio Coria irrigation channel is of interest for this work, especially Section 72 + 253 (Figure 6). Of which there is geometric information and annual average flows for the study of hydrodynamics in both physical and numerical models.

[IMAGE OMITTED. SEE PDF]

The Antonio Coria channel generally has trapezoidal sections throughout its length, with the exception of rectangular sections at the location of regulation gates; for this reason, for numerical modelling purposes, it was determined to use an equivalent rectangular section to optimize the laboratory physical model and the numerical model.

Three configurations of submerged structures described in the laboratory physical model section were made. These were designed in a trapezoidal shape with gradual concentric reduction, with the purpose of conducting the flow with the least possible turbulence at the inlet and searching for a velocity field that Based on the evaluation of the shear stresses produced at the bottom between the water and the sediments, the natural movement will be generated to produce hydrodynamic dredging.

Hydraulic and Hydrological Analysis

Steady-Flow Hydrodynamic Analysis

This analysis was performed using the design flow rate provided by the hydrological study prepared by Herrera (2021) [29] for the micro-basin associated with the Coria Canal as input, using HEC-HMS software. The updated flow rate for 2024 was 40.53 m3/s. The purpose of this analysis was to obtain the depths in the different cross-sections of the Coria Canal and the channel velocities in the centre of the channel using HEC-RAS software for the 25-year return period. Furthermore, current Mexican regulations for these analyses were followed, and a permanent and stable flow was defined as the initial condition, obtaining the velocities shown in Figure 7.

[IMAGE OMITTED. SEE PDF]

The maximum velocity occurring in the centre of the channel is 1.63 m/s.

Hydrodynamic Analysis in Unsteady Flow

The behaviour of a 1-year simulation of the channel was analysed, focusing on areas where the channel direction changes, flooding, and flood durations. This analysis was carried out using the adopted hypotheses and data provided by the characteristics of the analysed channel [29]. This allowed obtaining the maximum depths and velocities that can be reached by a flood occurring during the 25-year return period.

For the unsteady flow, the data were updated to the year 2024, and the associated hydrograph was obtained (Figure 8). The simulation intervals are every 24 h, thus obtaining the maximum and minimum levels reached in each section of the Coria Canal.

[IMAGE OMITTED. SEE PDF]

Figure 9 shows the variation in velocity at the centre of the channel for the simulation year, where a small fluctuation is observed, taking as an average value 1.62 m/s.

[IMAGE OMITTED. SEE PDF]

Below is a summary of the maximum conditions for the A. Coria channel (Table 1).

TABLE 1 Summary of HEC-RAS hydraulic simulation.

Steady flow	Unsteady flow (Max. condition)
W.S. Elev (m)	1723.93	W.S. Elev (m)	1722.76
Q total (m3/s)	40.53	Q total (m3/s)	39.02
Top width (m)	15.78	Top width (m)	15.70
Vel total (m/s)	1.63	Vel total (m/s)	1.58
Max Chl Dpth (m)	1.93	Max Chl Dpth (m)	1.88

Numerical Model and Grid

A numerical model for sediment particle simulation has been developed based on the finite difference method to replicate three-dimensional velocity fields. The transport dynamics of sediment particles are determined by the velocity field and turbulent dispersion mechanisms, incorporating stochastic Brownian-type random motion through the Particle-In-Cell Method (PICM). The processes of particle dispersion and resuspension are modelled stochastically using probability functions, which account for parameters such as settling velocities, bed-shear stress, bed-shear velocity, particle size, and shape factors. Model validation was conducted in prior studies. Subsequently, the model was applied to a scaled laboratory channel for irrigation with a variable slope at a 1:20 ratio.

For each scenario, a computational grid was constructed using finite elements (Figure 10). The hydrodynamic model was employed to generate velocity fields over a simulation duration of 1500 s with a time step (∆t) of 0.1 s. The flow rate magnitude served as the driving force across all configurations.

[IMAGE OMITTED. SEE PDF]

Initial Conditions

The initial hydrodynamic condition is a constant flow rate of 0.0019 m3/s for the three scenarios, where the velocity field is obtained to set the sediment particles planted in the channel in motion.

The initial sediment particle conditions for each configuration, detailed in Table 2, involved a total of 700,000 particles, distributed according to a granulometric curve comprising 87% sand and 13% coarse and bulky material. Sediment particles were introduced along the bottom of the laboratory channel using an approximation method based on the PICM approach to establish an initial sediment concentration at time zero. The concentration within each cell was calculated by dividing the total sediment mass within the cell by its corresponding volume.

TABLE 2 Sediment particle distribution.

Grid	Sand particle diameter (mm)	Bulky particle diameter (mm)
D15	0.19	1.25
D50	0.28	1.65
D90	0.42	1.80

Results and Application

The analyses were based on steady-state flow conditions. This scenario corresponds to the maximum value of the Coria Canal hydrograph for 2024 under unsteady flow. It is considered that this flow rate presents the most favourable conditions for sediment transport.

The velocity profiles derived from the physical modelling in the experimental channel were measured at the centreline of the cross-section. Figures 11–13 illustrate the specific points where velocity profiles were obtained for the three analysed scenarios: configurations Ls3, Ls4, and Ls5. These configurations represent varying hydrodynamic conditions modelled to evaluate sediment transport behaviour.

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Quantitative Assessment of Model Accuracy: Numerical Versus Experimental Data

A rigorous statistical analysis was conducted to assess the accuracy and predictive capability of the numerical model compared to the experimental measurements. The following statistical indicators were evaluated: the Mean Absolute Error (MAE), Mean Squared Error (MSE), Root Mean Squared Error (RMSE), Willmott's Index of Agreement (d-index), and the mean relative percentage error (ε %).

The computed values are as follows:

• MAE = 0.0055

• MSE ≈ 0.0000

• RMSE = 0.0069

• Willmott's (d-index) = 0.3507

• Relative error (ε %) = 3.38%

The low MAE and RMSE values indicate that the pointwise numerical predictions closely approximate the experimental data, exhibiting low absolute and squared errors. These low error magnitudes suggest a high degree of numerical precision in isolated data points (Figure 14).

[IMAGE OMITTED. SEE PDF]

Furthermore, Willmott's index of agreement (d = 0.3507) reinforces this observation. Although d-ranges from 0 (no agreement) to one (perfect agreement), a value around 0.35 indicates poor overall concordance between observed and predicted values.

Despite a reasonable relative error (ε % = 3.38%), which might indicate local accuracy in terms of magnitude, d-index points to a low capability in the model's ability to reproduce some underlying physical or dynamic trends observed in the experimental results.

The numerical model demonstrated acceptable pointwise accuracy in terms of small absolute and percentage errors. Nevertheless, the model falls short in quantitatively reproducing the experimental trend and variability. This discrepancy suggests the presence of systemic deviations between the numerical approach and the physical process being modelled; however, we can say that it is generally accepted for the purposes of modelling sediment particles.

Particle Sediment Transport Numerical Modelling

Following the computation of the velocity fields, sediment particles were seeded along the channel bottom. The characteristic sediment diameters used in the simulations were derived from granulometric analyses, as detailed in Table 2. These particle sizes served as inputs for the sediment transport model, ensuring a realistic representation of sediment behaviour in the study area.

Figure 15 presents the temporal evolution of sediment transport for the Ls3 configuration over a simulation period of 1500 s. The results demonstrate the spatial distribution of sediment particles at distinct time intervals (500, 1000, and 1500 s). Notably, sediment accumulation was observed in areas adjacent to hydraulic contractions, highlighting the influence of bottom geometry on particle deposition patterns.

[IMAGE OMITTED. SEE PDF]

For the Ls4 configuration, the results (Figure 16) reveal a different distribution of sediment particles. Enhanced sediment accumulation near the channel inlet was observed, attributed to variations in the hydrodynamic field and the flow velocity gradients near the bottom boundary.

[IMAGE OMITTED. SEE PDF]

In the Ls5 configuration, the particle transport exhibited a reduced concentration near the channel bottom (Figure 17). The greater distance between submerged obstacles mitigated particle deposition, resulting in a more uniform sediment distribution throughout the channel.

[IMAGE OMITTED. SEE PDF]

To evaluate the influence of particle size on transport behaviour, the distribution of particles with diameters D15, D50, and D90 was analysed for all configurations. Figure 18 illustrates the time evolution of particle concentration for these sizes across the different scenarios. The results underscore the variability in sediment transport dynamics based on particle size, with finer particles (D15) exhibiting prolonged suspension, while coarser particles (D90) displayed more rapid settling.

[IMAGE OMITTED. SEE PDF]

To further analyse sediment dynamics, the particle distribution was evaluated at four cross-sectional points (A, B, C, and D) located before and after two bottom contractions (as shown in Figure 19). These sections were chosen to assess the influence of flow constrictions on sediment concentration in the water column.

[IMAGE OMITTED. SEE PDF]

For the Ls3 configuration, the concentration profiles at each section (Figure 20) reveal a notable increase in particle suspension near sections A and B, where the flow initially interacts with the bottom contractions. The velocities generated at these points contribute to the upward transport of particles into the water column.

[IMAGE OMITTED. SEE PDF]

In the Ls4 configuration (Figure 21), significant changes were observed in the concentration of finer particles (D15) at sections B and C. This behaviour is attributed to the stochastic nature of the particle motion and the local hydrodynamic conditions, which amplify turbulence and enhance suspension in these regions.

[IMAGE OMITTED. SEE PDF]

For the Ls5 configuration, Figure 22 shows the particle concentration profiles for each section. Notably, section D exhibited distinct behaviour for finer particles (D15), likely influenced by the velocity field downstream of the second contraction. This observation highlights the role of flow deceleration in promoting sediment deposition further from the initial contractions.

[IMAGE OMITTED. SEE PDF]

Overall, the results demonstrate that the channel configurations significantly impact sediment transport dynamics. Sections C and D, located further downstream, exhibited lower sediment concentrations compared to sections A and B, which are directly influenced by the initial flow interactions with the bottom contractions. These findings emphasize the critical role of channel geometry and flow dynamics in determining sediment transport and deposition patterns.

ANOVA Analysis

A three-way analysis of variance (ANOVA) was performed to evaluate the effect of hydraulic configuration, measurement point, and sediment diameter type on the relative concentration of particles transported in the channel. The results obtained (Figure 23) showed that the three main factors Configuration, Point, and Diameter had statistically significant effects on concentration (p < 0.05). In particular, the Configuration factor showed highly robust significance (F = 42.58, p < 0.001), indicating that there are substantial differences in the behaviour of the different hydraulic configurations evaluated (Current, L3, L4, and L5) with respect to their capacity to mobilize sediment.

[IMAGE OMITTED. SEE PDF]

Likewise, significant interactions were detected among the factors, specifically Configuration × Point, Configuration × Diameter, and Point × Diameter, all with p-values < 0.001. These interactions suggest that the effect of hydraulic configuration on sediment concentration is not uniform across all measurement points or for all particle sizes, highlighting the importance of considering the spatial and particle size context in sediment transport analysis.

Figure 24 shows the distribution of normalized sediment concentration by configuration. Configuration L4 has the highest mean concentration (≈ 0.75) and a wide interquartile range, suggesting greater particle mobilization efficiency compared to the other configurations. In contrast, the Current configuration has a lower mean concentration (≈ 0.60) and greater dispersion, which could indicate less favourable flow conditions for sediment transport. Configuration L3, although presenting a similar average concentration to that of L5, shows a high number of outliers towards low concentrations, which may be indicative of unstable behaviour or punctual accumulations.

[IMAGE OMITTED. SEE PDF]

These results confirm that the hydraulic design of the channel significantly influences sediment transport, and that selecting the appropriate configuration can optimize the hydraulic performance of the system, minimize siltation and improve the channel's operational efficiency.

Spatial Distribution of Sediment Concentration

Figure 25 shows the statistical distribution of suspended sediment concentrations along four measurement points (P1 to P4) within the experimental channel. This visualization using box plots allows for identifying data variation and the presence of potential outliers at each location.

[IMAGE OMITTED. SEE PDF]

Median concentrations are observed to remain relatively constant at all four points, with values ranging from 0.58 to 0.67 units. However, data dispersion presents notable differences. In particular, point P3 exhibits a greater number of outliers below the interquartile range, suggesting greater fluctuations in flow conditions or local sedimentation behaviour in that section of the channel. On the other hand, points P1 and P2 exhibit more symmetrical distributions, while P4 shows a slight concentration toward higher values, possibly related to a decrease in flow velocity and greater entrainment capacity for fine particles.

The overall analysis of these boxplots suggests that, while there are no extreme differences in the medians, the variability and distribution of the data allow us to infer patterns in the spatial dynamics of sediment transport. These results reinforce the importance of channel geometric design and the strategic placement of hydraulic structures in controlling the mobility and deposition of suspended particles.

Influence of Particle Size on Sediment Concentration

Figure 26 presents the statistical distribution of sediment concentration as a function of characteristic particle size, represented by the particle size percentiles D15, D50, and D90. Box plots show how the concentration of fine, medium, and coarse particles varies along the experimental channel.

[IMAGE OMITTED. SEE PDF]

It is observed that the finest particles (D15) have the highest mean concentration and a wide dispersion, indicating greater stability in suspension within the flow. This is due to their low weight and lower terminal fall velocity, which facilitates their transport in the active zone of the channel. Furthermore, the D15 concentration shows a high interquartile range, reflecting high variability depending on local hydraulic conditions.

In contrast, the intermediate particles (D50) have a visibly lower median concentration accompanied by a more uniform distribution without any notable outliers. This behaviour is characteristic of sediments in partial transit, where transport depends on specific energy conditions that vary along the channel.

On the other hand, coarse particles (D90) exhibit a median concentration similar to that of D15, but are accompanied by a greater number of outliers below the lowest percentile. This suggests that, while some larger particles manage to remain suspended in certain areas of the channel, in most cases they tend to settle rapidly or concentrate in specific sections with greater turbulence. The high frequency of outliers in D90 indicates uneven transport conditions for larger sediments, possibly influenced by changes in flow morphology and internal structures.

These results confirm the strong dependence between sediment transport behaviour and particle size distribution, a key aspect for the design of hydraulic structures focused on erosion and sedimentation control in open channels.

Effect of Geometric Configuration on Mean Sediment Concentration

Figure 27 shows the mean sediment concentration for each of the geometric configurations evaluated in the experimental channel: the current geometry and the modified variants L3, L4, and L5. The error bars represent one standard deviation, allowing the variability associated with each configuration to be visualized.

[IMAGE OMITTED. SEE PDF]

Configuration L4 had the highest average concentration (≈ 0.71), followed by L5 and L3, while the current geometry yielded the lowest values (≈ 0.57). This result suggests that modifications to the channel geometry directly influence sediment transport and suspension, possibly by inducing more turbulent flow regimes that favour the resuspension of particles from the bottom.

The differences observed between configurations are also reflected in the width of the error bars, with greater dispersion in configurations L3 and L4, which may be associated with local interactions between internal hydraulic structures and the flow regime along the measurement points. In contrast, the current geometry shows less variability, suggesting a more uniform but less efficient flow pattern for maintaining sediment suspension.

These findings support the hypothesis that appropriate geometric modification can optimize sediment transport in artificial channels, a relevant element for the hydraulic design of agricultural and transport channels in areas with the presence of fine sediments.

Interaction between hydraulic configuration and measurement point on sediment concentration

Figure 28 illustrates the variation in average sediment concentration as a function of channel configuration (Current, L3, L4, L5) and longitudinal measurement point (P1 to P4). This graph visualizes the interaction between both variables and their influence on sediment transport behaviour within the hydraulic system.

[IMAGE OMITTED. SEE PDF]

It is observed that configuration L4 presented the highest average concentration values at most points, especially at P3, where it reached the highest recorded value (∼0.77). This behaviour suggests that the geometry of configuration L4 facilitates flow conditions that favour particle suspension, possibly due to a combination of contractions and expansions that generate zones of effective turbulence. L5 also shows high and consistent efficiencies, with concentrations higher than the overall average at all measured points, primarily at P3 and P4.

In contrast, the current configuration presents lower concentrations at intermediate points (P2 and P3), evidencing less efficient sediment transport in those sections of the channel. In particular, at P2, this configuration obtained the lowest reported value (∼0.47), which could be related to a decrease in local flow velocity or to deposition processes due to gravitational settling.

Comparatively, configurations L3 and L4 stand out in P2 and P3, significantly increasing suspended sediment concentrations, reflecting a significant improvement over the original channel design. In P4, although the differences are less marked, both L4 and L5 maintain higher concentrations compared to the base geometry, suggesting greater sustained entrainment capacity in the final reaches of the channel.

These results reinforce the hypothesis that the geometric modifications applied in configurations L3, L4, and L5 induce more favourable hydraulic conditions for the continued transport of suspended sediment, especially at critical points where geometry influences flow energy.

Relationship between channel configuration and particle size on sediment concentration

Figure 29 shows the average suspended sediment concentration for each channel geometry configuration (Current, L3, L4, and L5), broken down by characteristic particle diameter: D15 (fine), D50 (medium), and D90 (coarse). This representation allows analysis of the combined effect of channel geometry and particle size on sediment transport efficiency.

[IMAGE OMITTED. SEE PDF]

In the current geometry, the highest concentration was recorded for fine particles (D15), with a progressive decrease toward larger sizes, with D50 presenting the lowest value (∼0.47). This pattern is consistent with the low hydraulic energy of the base configuration, which favours the entrainment of light sediments and the deposition of coarser fractions.

In contrast, the L4 configuration exhibits the highest average concentrations for all three sediment sizes, reaching a maximum value (> 0.75) for fine particles (D15) and maintaining equally high levels for D50 and D90 (∼0.67 and ∼0.68, respectively). This suggests a more turbulent and uniform flow regime capable of supporting a higher proportion of suspended particles regardless of their particle size distribution.

Configurations L3 and L5 show intermediate behaviour. In configuration L3, the concentration is more evenly distributed across the three diameters, although intermediate particles (D50) have a slight advantage. In L5, coarse particles (D90) exceed D50 in concentration, indicating areas with possibly higher local velocities that manage to maintain the entrainment of heavier fractions.

Overall, the results indicate that channel geometry not only affects the magnitude of sediment transport but also the selectivity toward different particle sizes. Configurations such as L4 demonstrate a greater capacity to move sediments of different particle sizes simultaneously, which represents a relevant criterion for hydraulic design when seeking to avoid selective sedimentation.

Distribution of Sediment Concentration by Measurement Point and Particle Size

Figure 30 represents the average sediment concentration along four measurement points in the channel (P1 to P4), differentiated by the particle size percentiles D15, D50, and D90. This graph allows for the evaluation of the spatial behaviour of sediment transport based on the size of the solid particles present.

[IMAGE OMITTED. SEE PDF]

In general terms, fine particles (D15) show the highest concentrations at all points, with a maximum value at P4 (∼0.69). This behaviour reflects the high transport capacity of light suspended particles, which tend to remain distributed throughout the turbulent flow without easily settling. This trend is consistent with what was observed in previous studies on the behaviour of small-diameter particles in open channels.

An interesting difference is observed at P3, where coarse particles (D90) reach an even higher value than D15 (∼0.68), suggesting the presence of a zone with hydraulic conditions that favour the entrainment of heavy particles. This phenomenon could be associated with local structures that increase turbulence or form vortices, temporarily favouring the suspension of larger particles.

In contrast, intermediate particles (D50) exhibit more variable behaviour. At points such as P1 and P3, their average concentration drops below D90, while at P2 and P4 they slightly exceed coarse particles, indicating that the transport of these particles strongly depends on local flow conditions, such as transverse velocity and specific kinetic energy.

Taken together, these results indicate that local flow characteristics along the channel have a significant impact on the distribution of sedimentary material, especially in relation to particle size. The system's ability to maintain different particle fractions in suspension varies spatially, which is a critical aspect for the design and evaluation of channels in agricultural or runoff management systems.

Discussion

The results highlight the significant influence of increased flow velocity near the channel bottom on sediment transport dynamics. Sediment particles respond directly to changes in flow velocity, exhibiting distinct movement patterns along the channel under different configurations of submerged structures.

The design of submerged structures, particularly their shape and dimensions, plays a critical role in modulating flow velocity and influencing sediment transport. Trapezoidal-shaped structures, as implemented in this study, were found to enhance sediment mobilization with greater efficiency while minimizing the generation of high-magnitude turbulence (Figures 11–13). The hydrodynamic effects induced by these structures create localized flow conditions that favour sediment transport along the channel, which is especially advantageous during periods of increased discharge, such as during rainy seasons.

Likewise, in addition to the hydraulic and functional analysis of the irrigation canal, it is essential to consider the economic dimension associated with its operation and maintenance. In particular, periodic dredging of the canal represents one of the main recurring costs. The average cost per cubic metre is $15.60 USD, and sediment, trash, and other objects are removed in an average annual volume of 1024 m3, representing a cost of $15,974.40 USD per year. If this maintenance is not performed, sediment accumulation can significantly reduce hydraulic capacity, generate obstructions, and affect system performance. These cleaning operations often require specialized machinery, skilled labour, and intervention times that disrupt the canal's normal operation, increasing operating costs year after year.

In this regard, it was identified that the implementation of submerged hydraulic works, such as sediment movement structures or lined pipelines, can represent a higher initial investment compared with traditional methods. However, a long-term cost-benefit analysis shows that these works are more cost-effective. The cost per cubic metre of concrete with steel is $220.00 USD. If option L4 is considered, the volume of each submerged structure is 2.10 m3. Using four elements, this gives a volume of 8.4 m3 of concrete, with a total construction cost of $1848.00 USD, significantly cheaper than current maintenance in the study section. By designing and building these structures to operate in submerged conditions, they not only contribute to reducing sedimentation but also significantly reduce the need for frequent dredging, extending the intervals between major maintenance. This results in a considerable reduction in annual operating costs and greater efficiency in the use of water resources.

The three-way analysis of variance (ANOVA) identified that the three main factors—channel configuration, measurement point, and particle size—had a statistically significant influence on suspended sediment concentration (p < 0.05 in all cases). Among these, the geometric configuration of the channel proved to be the most influential factor (F = 42.58), demonstrating that modifications to the channel design are crucial for improving sediment transport conditions.

The various visualizations and statistical analysis clearly indicate that the L4 configuration is the most effective at maintaining suspended particles along the channel. This conclusion is reinforced in Figure 3, where L4 reaches the highest mean concentration, indicating that its design promotes the greatest entrainment capacity in that section.

The significant interaction between configuration and point (F = 24.62) reinforces this observation, showing that the L4 geometry maintains optimal performance in different channel sections, unlike the base configuration (Actual), where there is a marked drop in efficiency at point P2. Similarly, the configuration: diameter interaction (F = 17.83) revealed that L4 exhibits high performance in transporting all sediment sizes (D15, D50, and D90), consistently outperforming L3 and L5, which show favouritism for specific fractions.

In comparison, the Actual configuration was statistically inferior in its ability to retain suspended sediments, especially intermediate particles (D50). Configurations L3 and L5 offered improvements over the base design, although with variable responses depending on the point and particle diameter.

The combined analysis of the figures and the ANOVA allows us to conclude that L4 is the most efficient and complete configuration, as it maximizes the average sediment concentration and maintains consistent performance across spatial and particle size variations. This finding is crucial for the hydraulic design of canals requiring high transport capacity, especially in agricultural settings where sedimentation control directly impacts the system's productivity and sustainability.

In summary, maintaining a constant hydraulic cross-section in irrigation canals is critical for ensuring the efficient and effective conduction of water from primary irrigation districts to secondary modules. This not only promotes optimal irrigation practices and maximizes crop yields but also contributes to the operational sustainability of the agricultural system. In this context, the strategic integration of submerged hydraulic structures represents a practical and technically efficient solution to the challenges associated with sedimentation. While these structures may entail a larger initial investment, their implementation allows for maintaining the canal's hydraulic capacity, significantly reducing periodic dredging costs, and decreasing maintenance frequency. Therefore, it is concluded that canal planning that incorporates both hydraulic and economic criteria should prioritize the use of durable structural solutions that ensure more efficient and cost-effective long-term operation of the irrigation system, thus optimizing investment in public and private hydraulic infrastructure.

Conclusions

Sediment transport in irrigation channels is a critical factor affecting hydraulic efficiency, operational sustainability, and long-term cost reduction. The findings of this study confirm that controlled increases in bottom flow velocity, induced by submerged structures, significantly enhance sediment removal without compromising the channel's hydraulic cross-section.

Numerical modelling, validated through scaled physical channel experiments, demonstrated that submerged structure configurations can be optimized to generate velocity patterns that promote sediment entrainment, particularly during high-flow events. Among the tested scenarios, the L4 configuration exhibited superior statistical and hydraulic performance as confirmed by both graphical analyses and a three-way ANOVA. This configuration consistently achieved the highest average sediment concentrations across different particle sizes and longitudinal points along the channel, maintaining a more uniform suspension and minimizing localized deposition.

These results support the feasibility of implementing a natural hydrodynamic dredging approach in irrigation channels through the strategic redesign of the channel beds with submerged structures. In particular, the configuration L4 stands out not only for maximizing sediment suspension but also for performing consistently across spatial and granulometric variables, outperforming other tested alternatives, including the original (unmodified) geometry.

The statistical analysis revealed that geometry (configuration) was the most influential factor (p < 0.001), followed by sediment diameter and measurement location. Significant interactions between configuration and diameter, and between configuration and channel position, further emphasize that optimal flow and sediment behavior arise from considering these variables in combination rather than isolation.

This innovative approach helps preserve hydraulic capacity, reduces the frequency of mechanical dredging, and improves water delivery to secondary channels—key aspects for enhancing agricultural productivity and strengthening the climate resilience of irrigation systems.

Future research should focus on the field-scale implementation of submerged structures in secondary or distributary irrigation channels, complemented by in situ hydraulic and sediment monitoring and continuous model validation. In parallel, the development of design and operational guidelines based on performance indicators such as sediment retention, cost-benefit analysis, and maintenance intervals is essential to promote a new generation of sustainable, efficient, and climate-adaptive irrigation infrastructure. In addition, in situ measurements of sediment transport are being conducted and several environmental variables are being monitored jointly. In the long term, with in situ measurements, the use of artificial neural networks to predict sediment transport under different hydrodynamic conditions is being considered.

Author Contributions

Israel E. Herrera-Díaz: conceptualization, investigation, formal analysis, writing – original draft. Rosario Martínez-Yáñez: validation and supervision. Edith A. Gamiño-Ramírez: project administration, writing – review and editing. César Gutiérrez-Vaca: methodology, software.

Conflicts of Interest

The authors declare no conflicts of interest.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

References

A. Arifjanov, A. Fatxulloev, and T. Kaletová, “Movement of the Suspended Flow in the Open Irrigation Channels,” Acta Horticulturae et Regiotecturae 2 (2019): 80–83, https://doi.org/10.2478/ahr‐2019‐0015.

K. Katayoun, S. S. Mir Jafar, M. Mirali, and S. Saeed, “Sediment Transport Modeling in Open Channels Using Neuro‐Fuzzy and Gene Expression Programming Techniques,” Water Science and Technology 79–12 (2019): 2318–2327, https://doi.org/10.2166/wst.2019.229.

K. S. Black, S. Athey, P. Wilson, and D. Evans, “The Use of Particle Tracking in Sediment Transport Studies: A Review,” Coastal and Shelf Sediment Transport 274, no. 1 (2007): 73–91: Geological Society of London, Special Publications, https://doi.org/10.1144/gsl.sp.2007.274.01.09.

R. Kosaj, R. S. Alboresha, and S. O. Sulaiman, “Sediment Incipient Motion in the Rigid Boundary Open Channels of the City of Ramadi, Iraq,” Design Engineering, no. 6 (2021): 1474–1484.

X. Liu, A. Mohammadian, and I. Sedano, “One Dimensional Numerical Simulation of Bed Changes in Irrigation Channels Using Finite Volume Method,” Irrigation & Drainage Systems Engineering 1, no. 2 (2012): 1–103, https://doi.org/10.4172/2168‐9768.1000103.

A. Syarifudin, “Erosion and Sedimentation Analysis of Rural Channels (RC) and Main Drainage Channels (MDC) on Tidal Lowland in Indonesia,” International Journal of Innovative Science and Research Technology 7, no. 1 (January 2022): 419–426.

Y. Li, P. Yang, T. Xu, et al., “CFD and Digital Particle Tracking to Assess Flow Characteristics in the Labyrinth Flow Path of a Drip Irrigation Emitter,” Irrigation Science 26, no. 5 (2008): 427–438, https://doi.org/10.1007/s00271‐008‐0108‐1.

H. A. Omran, M. Y. Fattah, and S. Q. Sulaiman, “A Model for Removing Sediments From Open Channels,” International Journal of the Physical Sciences 9, no. 4 (2014): 61–70, https://doi.org/10.5897/IJPS2013.4074.

A. A. Kamanbedast, “Investigation of Sediment Depth‐Volume and Seepage (With Using Water Stops Materials) in Canals of Irrigation Network,” World Applied Sciences Journal 11 (2010): 1353–1360.

H. Ouarriche, M. El Bouhali, F. Bouisfi, M. Chaoui, and M. Hannaoui, “Numerical Simulation of Flow in Trapezoidal Labyrinth‐Channels of Drip Irrigation,” Faculty of Mechanical Engineering, Belgrade, Transactions 49 (2012): 147–153, https://doi.org/10.5937/fme2101147O.

S. Xing, Z. Wang, J. Zhang, N. Liu, and B. Zhou, “Simulation and Verification of Hydraulic Performance and Energy Dissipation Mechanism of Perforated Drip Irrigation Emitters,” Water 13, no. 2 (2021): 171, https://doi.org/10.3390/w13020171.

L. S. de Sousa, R. M. Wambua, J. M. Raude, and B. M. Mutua, “Assessment of Water Flow and Sedimentation Processes in Irrigation Schemes for Decision‐Support Tool Development: A Case Review for the Chókwè Irrigation Scheme, Mozambique,” AgriEngineering 1 (2019): 100–118, https://doi.org/10.3390/agriengineering1010008.

N. B. Hansen, Two‐Dimensional Sediment Transport Model Using Parallel Computers (LSU University, 2022), https://doi.org/10.31390/gradschool_theses.1113.

L. S. Giang and T. T. M. Hong, “3D Numerical Modeling of Flow and Sediment Transport in Rivers and Open Channels,” Science & Technology Development Journal ‐ Science of The Earth & Environment 3, no. 1 (2019): 23–36, https://doi.org/10.32508/STDJSEE.V3I1.508.

S. Kumar, P. R. Hanmaiahgari, and H. P. Singh, “Developing a Numerical Model for Sediment Transport in Channel Network by Considering Multigrade Bed Load Sediment,” ISH Journal of Hydraulic Engineering 30 (2023): 1–16, https://doi.org/10.1080/09715010.2023.2250303.

M. Renaud, J. Chauchat, C. Bonamy, and O. Bertrand, Ongoing Development of a Three‐Dimensional Sedimenttransport Model for Local Scour Study (EGU General Assembly 2023, 2023), https://doi.org/10.5194/egusphere‐egu23‐11987

G. Tayfur, “Empirical, Numerical, and Soft Modelling Approaches for Non‐cohesive Sediment Transport,” Environmental Processes 8, no. 1 (2021): 37–58, https://doi.org/10.1007/S40710‐020‐00480‐1.

S. Choi and J. H. Hwang, “Lagrangian Coherent Structure Analysis on the Vegetated Compound Channel With Numerical Simulation,” Water 14, no. 3 (2022): 406, https://doi.org/10.3390/w14030406.

H. Romdhane, A. Soualmia, L. Cassan, and G. Belaud, “Effect of Vegetation on Flows and Sediment Transport,” E3S Web of Conferences 40 (2018): 02017, https://doi.org/10.1051/e3sconf/20184002017.

I. E. Herrera, C. Rodríguez, C. Couder, and J. R. Gasca, Modelación numérica hidrodinámico‐ hidrológica en zonas de inundación con presencia de infraestructura, Vol. 6 (Tecnología y Ciencias del Agua, 2015): Num. 1 Ene‐Feb, 139–152.

C. Rodríguez, C. Couder, E. Flore, I. E. Herrera, and R. Cisneros, “Modelling Shallow Water Wakes Using a Hybrid Turbulence Model,” Journal of Applied Mathematics 2014 (2014): 1–10, https://doi.org/10.1155/2014/714031.

P. Broomans, “Numerical Accuracy in Solutions of the Shallow‐Water Equations,” Master Thesis, TU Delf & WL, Delf Hydraulics (2003).

UNESCO, Background Papers and Supporting Data on the International Equation of State of Seawater, Vol. 38 (UNESCO technical papers in marine science, 1981), 192.

P. Stansby, “A Mixing‐Length Model for Shallow Turbulent Wakes,” Journal of Fluid Mechanics 495 (2003): 369–384, https://doi.org/10.1017/s0022112003006384.

J. R. Garratt, “Review of Drag Coefficients Over Oceans and Continents,” Monthly Weather Review 105, no. 7 (1977): 915–929, https://doi.org/10.1175/1520‐0493(1977)105<0915:rodcoo>2.0.co;2.

M. Robinson, J. Monaghan, and J. Mansour, “SPH Simulation of 2d Wall‐Bounded Turbulence,” in SPHERIC, Smoothed Particle Hydrodynamics European Research Interest Community (Second International Workshop, 2007), 107–110.

I. Hernandez, “Modelos euleriano‐lagrangeanos en flujos a superficie libre: Aplicación al transporte de partículas suspendidas y al crecimiento de microorganismos”, Tesis Maestría, DEPFI (UNAM, 2003).

L. C. Van‐Rijn, D. J. Walstra, and M. Van‐Ormondt, “Unified View of Sediment Transport by Currents and Waves IV: Application of Morphodynamic Model,” Journal of Hydraulic Engineering 133–7, no. 7 (2007): 776–793, https://doi.org/10.1061/(asce)0733‐9429(2007)133:7(776).

I. E. Herrera, Proyecto ejecutivo de rehabilitacion de camino rural la Compañia ‐ Cardenas, y puente vehicular, en el Municipio de Salamanca, Gto. Méx, Technical Report, ASESPROYCo. (2021).