One Dimensional Numerical Simulation of Bed Changes in Irrigation Channels using Finite Volume Method

Irrigation channels that transport water in order to irrigate croplands are of significant importance in agricultural studies. The failure of irrigation channel in proper transport of water may cause substantial costs. Numerical prediction of irrigation channel erosion is of practical significance in agricultural practice. Mean while, fertilizers and nutrients essential to crops are attached to the soil surface and transported along the furrows by the irrigation flow. To optimize nutrient delivery to crops, one requires accurate predictions of flow and sediment transport in irrigation furrows [1].


Introduction
Irrigation channels that transport water in order to irrigate croplands are of significant importance in agricultural studies. The failure of irrigation channel in proper transport of water may cause substantial costs. Numerical prediction of irrigation channel erosion is of practical significance in agricultural practice. Mean while, fertilizers and nutrients essential to crops are attached to the soil surface and transported along the furrows by the irrigation flow. To optimize nutrient delivery to crops, one requires accurate predictions of flow and sediment transport in irrigation furrows [1].
The hydraulics of irrigation channel flows has been commonly calculated by kinematic wave and zero inertia models in early studies. A zero inertia model of the stream flow was presented in the context of negligible accelerations by Strelkoff and Katopodes [2]. Schwankl and Wallender [3] adopted the zero inertia furrow model with a combination of spatially-varying infiltration and steady infiltration. Clemmens [4] also used the zero inertia model for simulating irrigation channels. Elliott et al. [5] modified the zero-inertia model for continuous furrow irrigation and showed that the zero inertia model can effectively simulate the hydraulics of furrow irrigation. Oweis and Walker [6] developed a surge flow furrow irrigation model based on the zero inertia concept. Chen [7,8] used kinematic wave theory to solve the problem of irrigation flow over a wide porous bed and concluded that the kinematic wave method may only be valid for supercritical flow. Woolhiser [9] showed that the method would also be applicable to subcritical flows but may lead to poor results if water pounds at the downstream boundary and a moving backwater extended over an appreciable length of the field. Smith [10] adopted a kinematic-wave method of characteristics and the Lax-Wendroff scheme to simulate irrigation flow. Walker and Humpherys [11] developed and verified a kinematic-wave model of furrow irrigation under both continuous and surged flow management. Lu et al. [12] developed a numerical kinematic wave model to simulate the dynamic erosion process in a ridge-furrow system. Based on the kinematic wave method, Reddy and Singh [13] used moving control volume approach to model the flow in irrigation systems and used fixed control volume to model the nearly stationary phase and the runoff rate, and showed a good agreement with the experimental data.
These studies showed that the zero inertia model is applicable for surface irrigation flow with negligible errors, the kinematic wave model is more applicable for steeper slopes in irrigation furrows with a larger Reynolds number. Nonetheless, the accuracy of these models may be limited for unsteady irrigation flow due to the neglection of flow acceleration.
One-dimensional models solving St. Venant equations are considered to be accurate and rigorous for hydraulics of irrigation flow. Sakkas et al. [14] used a mathematical model based on the complete hydrodynamic equations describing open-channel flow developed for simulation of surface irrigation systems. Abbasi et al. [15] presented a system combining an overland waterflow and solute transport model and applied it to furrows/borders. The integrated model numerically solves the 1D dispersion-advection and simplified form of the Saint-Venant equations as the governing equations in a decoupled fashion. Burguete et al. [16] adopted a coupled system including solute transport and implicated the model on furrow irrigation and obtained reasonable results.
Several sediment transport models have been examined for furrow irrigation flows. Lu et al. [12] examined two concepts of the erosion process of the loose, deposited sediment and concluded that the best concept for describing erosion of this loose layer of sediment was one where the deposited sediment was eroded uniformly across the width of the channel. Strelkoff and Bjorneberg (2001) implemented Yang's and Yalin's sediment transport method on furrow erosion and concluded that furrow erosion by upstream inflow may be largely different from hillside erosion by rain-fed overland flow and predictive approaches and formulas for the latter problem may not be satisfactory for the former.

Abstract
A one-dimensional Saint Venant model is developed using the Einstein sediment transport formula for simulating hydraulics and bed changes in irrigation channels. The governing equations are discretized using Finite Volume Method. The Central-Upwind method is used in this study to calculate the numerical flux. The model gives a stable and good prediction of sub, super and transcritical flows and captures harp changes such as hydraulic jump without leading to numerical oscillation or diffusion. An experimental test case is adopted to examine the performance of this mode. The simulated results show a good agreement with experimental data in the upstream section, and can predict an average value of measured profile in the downstream section of a steep channel. The knick point movement is successfully predicted by the model. The model appears to be a potential and reasonably accurate tool to predict bed changes in irrigation channels due to the water flow. This study focuses on the bedload sediment transport which is defined as the type of transport where sediment grains roll or slide along the bed. A lot of sediment transport models have been developed to calculate the erosion of channel bed such as the Grass equation [17], Meyer-Peter & Muller formula [18], Van Rijn equation [19][20][21], Nielsen formula [22], Kalinske equation [23,24], Einstein formula [25,26], Smart equation [27] etc. which are generally obtained by empirical methods. Some of these models are deterministic while others are derived based on probabilistic methods. After trials and comparison, the Einstein equation is adopted to compute the channel bed erosion in this study.
This research is aimed at developing a one-dimensional numerical model based on finite volume method for simulating flow and sediment transport. The Central-Upwind scheme is used here to calculate the numerical flux. The model is developed by combining a set of one-dimensional Saint Venant equations and an accurate sediment transport formula. The experimental data from Brush and Wolman (Brush and Wolman, 1960) are used in this research as a test case to examine the performance of the numerical model.

Hydraulic Flow and Channel Bed Erosion
For soil loss cases where the flow through the breach is usually unsteady, the hydraulics can be formulated based on continuity and momentum equations. In this study, a set of 1-D St. Venant equations based on non-cohesive sediment and mobile bed are applied focusing on the incompressible flow in channel. A coordinate system is defined where the x-axis represents the horizontal direction along the channel and the z-axis represents the vertical upward direction. In this study, the width of the irrigation channel is assumed to be constant and equal everywhere. Ignoring the flow infiltration, the governing equations are adopted based on mass and momentum conservation for each phase.

Governing equations
The 1-D St. Venant equations describing the flow through a breach can be written as: The term f S is friction slope term which can be expressed using the Manning formula

Bed change model
The erosion of the soil surface can be calculated by the sediment continuity equation based on the sediment mass conservation. The continuity equation (the Exner equation) is given by q is the sediment transport rate per unit width which can be expressed in a general form as: where d is the particle size.
To control the stability and accuracy of the numerical modeling result, an explicit modified Lax scheme [28] is introduced as follows: where α is a weighting coefficient which can be adjusted from Cb is the bed celerity in a frictionless system which is also given by [29]: where u = flow velocity in the x-direction; h = water depth; b q = sediment transport rate per unit width and Fr is Froude number in a rectangular channel which can be expressed as / Fr u gh = .

Numerical Scheme
In this study, we assume the breach width to be constant (unit width). The 1-D St. Venant equations in a conservative form can be written as: The St.Venant formula can also be written in a non-conservative form: where A is the Jacobian matrix which can be written as: and c gh = s the wave velocity.
The eigen values i a corresponding to the matrix A are 1 2 , a u c a u c = + = − respectively.
The St. Venant equations are integrated over every control volume based on the finite volume method.
An explicit scheme for discretizing the time step integration term is used in this study. The time step is restricted by the Courant-Friendrichs-Lewy (CFL) condition. In a practical simulation, a constant time-step is usually used and the CFL condition is used to approximately determine the time-step before the simulation starts. The CFL condition is given by: where t ∆ = time step size; x ∆ = spatial step size; n C is the flow Courant number and its value must be less than 1 for numerical stability.
For the one dimensional case, a standard splitting algorithm is employed for discrete equation (3.8) [30] where superscript n = time step index; subscript i = spatial node index; The key problem in the above discretization is how to calculate the numerical fluxes 1/2 n i F + and 1/2 n i F + at the interface C Γ between two neighboring cells. According to Godunov [30], the normal fluxes can be obtained by an approximate Riemann solver in the direction normal to the interface. A wide variety of reliable finite-volume methods for the St. Venant system have been developed in the past few decades. In this study, the Central upwind scheme is applied to calculate the fluxes.   The sediment continuity equation can also be discretized by using the finite difference method and we obtain:

Numerical Results
In this section, two numerical experiments are performed to evaluate the performance of the model. The first case deals with a hydraulic test and the second test considers sediment transport and erosion in a channel.

Hydraulic test case
Here, we perform a numerical test to evaluate the performance of hydrodynamic solver. This test is considered to predict the water surface profile changing in an irrigation channel when the upstream flow volume rate decreases. The initial upstream flow rate is 0.6 m 3 /s and water depth is 0.6 m. The bed elevation is 1 m and the width is 1 m. The flow rate decreases to 0.3 m 3 /sover at a certain time. A time step of 0.001 s is considered here and a comparative test with a time step of 0.0001 s is also carried out with a higher resolution. Figure 2 shows the water profile changing at different time in both tests.
The computed results show that the predicted water surface profiles are smooth without any numerical oscillation or diffusion. The predicted surface changing with higher accuracy is a little flatter than the one with larger time step. The two curves are very close at different times. Thus, the changing of time step has a tiny effect on the modeling result.

Experimental Test Case
An experimental case of channel bed erosion is considered in order to test the performance of this numerical scheme.  Many different sediment transport equations were employed to predict the sediment discharge rate and the development of the shape of channel bed. In this numerical study, the sediment discharge rate was estimated using the bed load formula of Einstein [25] which is also used to predict the channel bed development by [31]  A moderate time step of 0.1 s and a spatial step of 0.1 m were used. Both channel bed and water surface profiles are smooth without any numerical oscillations. Figure 4 shows the comparison between the measured and computed channel bed profiles at a time of t = 160 min. It can be observed that the computed bed profile has a very good agreement in the upstream section with the measured data. It accurately predicted the knick point migration in the upstream section and agreed with the measured values very well. In the downstream section of the reach, the computed profile was smoother than the original measured data and it predicted an average value of experimental data.

Conclusions
The modeled results of irrigation channel erosion have a reasonable agreement with the measured data considering the large uncertainty of the sediment transport equation and the assumption of fixed width. From the report of Brush & Wolman (1960), the width of the channel increased downstream of the knick point by as much as ± 20% at the widest section. The increase of the downstream width may cause more sediment deposition and maybe one of the reasons accounting for the differences between the measured and the computed profiles in the downstream section.
The developed numerical model is found to be a convenient and adequate tool in predicting the bed erosion of irrigation channels, the main features of this model are: (1) It uses a 1-dimensional grid with finite volume method.
(2) The numerical flux is calculated using the Central Upwind method which is able to simulate complex flows including sub, super and transcritical currents with a good accuracy. It shows a good numerical stability without oscillations or diffusion under a moderate spatial step of 0.1m, which is crucial for an accurate simulation.
(3) The numerical model has a good ability to capture shock waves and is able to predict the fronts of hydraulic jump accurately.
(4) The model with Einstein bed load formula can dynamically update the elevation of the channel bed during the deposition process and is in good overall agreement with experimentally measured data.
The main limitations of current model are: (1) It ignores sediment suspension.
(2) Sediment particle size is assumed to be equal to the mean value.
(3) The 1-D model assumes a constant width of the irrigation channel.
(4) The sediment transport formula is empirical and is not suitable for all situations.
From the results of the test case, the computed values were in good overall agreement with the measured ones. Although some disagreements exist between the simulation and experimental results, these may be explained by the limitations listed above. These preliminary results show a good potential to simulate the bed erosion of irrigation channels.