One Dimensional Numerical Simulation of Aggradation-Degradation in a Channel Using Finite Difference Method Case Study Chashma Right Bank Canal ( CRBC )

A one-dimensional MacCormack explicit finite difference model is developed for simulating hydraulics and bed changes in irrigation channels. The Saint-Venant equations describing unsteady flow in open channels and the continuity equation for the conservation of sediment mass are numerically solved. These equations are highly nonlinear and therefore do not have analytical solutions. For this purpose the MacCormack scheme is used. The scheme is second order accurate; it is a coupled solution as it is a two-step predictor corrector method. Model gives results in terms of bed level changes, flow depth and discharge provided physical boundaries of the system are valid for simulation time. Model execution and accuracy is very sensitive to time step and stability. The simulated results show a good agreement with previous studies in the downstream section, and can predict an average value of measured profile in the upstream section. The application of this model to Chashma Right Bank Canal in Pakistan and results of the model are compared with the published results gave very convincing result. Citation: Riaz MZB, Shakir AS, Masood M (2017) One Dimensional Numerical Simulation of Aggradation-Degradation in a Channel Using Finite Difference Method Case Study Chashma Right Bank Canal (CRBC). Fluid Mech Open Acc 4: 178. doi: 10.4172/2476-2296.1000178


Introduction
The design of hydraulic structures and water resources system, in the analysis of channel mechanics problems and in the development of channel control works is of vital importance in gradually varied unsteady flow in open channels described by solution of equations. Transport of sediment by flow and its scour or depositions are involved in many of these problems of rivers and canals. Even when the flows are steady, sediment transport phenomena are often time variant. Appropriate quantification of aggradation and/or degradation and variations in channel form of such channel still has been a subject of considerable research. From the viewpoint of improved planning and designing of various water development projects, a systematic study on the aggradation-degradation phenomenon of these channels of Pakistan is, therefore, of utmost importance.
To study the effect of the long term and short term bed level changes in alluvial channels with different flow conditions, a number of experimental had been conducted. Begin et al. [1] experimentally studied degradation of alluvial networks in response to lowering of the bed level. Soni et al. [2] conducted an experiment that covered a extensive range of flow and sediment loading conditions. They observed that after a long time, when the hydraulic conditions became compatible with the increased sediment load, the aggradation in downstream of the section of increased sediment supply was stopped. Yen et al. [3] found out that both the aggradation wave speed and the mean sediment transport velocity increase with the initial equilibrium bed slope and with decreasing loading ratio after performing a series of overloading experiments with uniform coarse sediment. To study the reversibility of an alternation aggradation-degradation process, a series of experiments were conducted by Yen et al. [4]. The results show that as standard deviation of sediment gradation increases the recovery rate decreases. MacCormack scheme was applied by Alam [5] to the study of aggradation-degradation in alluvial channels. Kassem and Chaudhry [6] developed a two dimensional numerical model to predict the time variation of bed deformation in alluvial channel bends. In this model, Beam and Warming alternating direction implicit scheme are used to solve the depth-averaged unsteady water flow equations along with the sediment continuity equation. Deng and Li [7] used the one dimensional unsteady flow and sediment transport model to study the river channel equilibrium profile. For sediment transport, erosion and deposition in straight channels and rivers some of the widely used one dimensional models are MIKE11 (DHI [8] and HEC-6 (USACE). Current study uses the explicit finite difference method to determine the rate and extent of bed level-changes due to transient conditions of flow and sediment (routing of sediment laden water). The finite difference explicit method of MacCormack is used to numerically solve the continuity and momentum equations along with the sediment mass balance equation.
The present work is aimed at investigating more closely where relatively rapid changes in both fluid and sediment discharge are imposed at the upstream boundary, with coupled solution through MacCornaack.

Computational Procedures
Flow computations, and sediment routing are two parts in each step of the simulation method. Firstly flow computation is executed then sediment routing is executed to compute the quantity of channel bed variations. (Figure 1) shows the computational flow.

Numerical simulations
The effective application of a model basically depends on the

Governing equations
The one-dimensional partial differential equations describing unsteady flow in a channel with movable bed (Figure 2) are given by:

1)
Continuity equation for water

Continuity equation for sediment
Terms used in these equations are defined as, Where, n=Manning roughness coefficient.
The sediment discharge is predicted by an empirical function of the flow velocity: Where, a, b are empirical constants to be adjusted. Values of a, b depend upon sediment properties and channel geometry.

Numerical scheme
The eqn. (1) through eqn. (3) characterizes as nonlinear hyperbolic partial differential equations and closed form solutions of these equations are not conceivable with the exception of rare simplified cases. Therefore, they are solved by numerical schemes. In this model, a finite difference method is used to solve MacCormack scheme [9]. Two possibility of this scheme are conceivable for one dimensional flow. In one alternate, forward finite differences are used in the predictor part and backward finite differences are used in the corrector part. In the second alternative, backward finite differences are used in the predictor part and forward finite differences are used in the corrector part. The values of the variables so determined during the predictor part are used during the corrector part. The finite difference approximations for the derivatives of a variable 'f' are as follows: Predictor step: Where, i=the node in spatial grid, k=the node in time grid, Δx=the distance step and Δt=the time step. An asterisk refers to the predicted values of the variable.
Corrector step:  Where, the two asterisks denote the value of the variables calculated in the corrector step.

Boundary conditions
The above equations determine the values of h, q, and z at the new time level k+1 at each interior grid points (i=2, 3… N-1). The values at the boundary grid points (i=and N) cannot be calculated by using these equations. Therefore, they are calculated by using boundary conditions. For one-dimensional flows, characteristic method gives acceptable results and is employed here in. For the proper functioning of the model two boundary conditions at upstream boundary and one boundary condition at downstream end are required. In addition to these, every computational grid needs initial conditions. Uniform unit discharge (q 0 ), uniform flow depth (h 0 ) and bed levels as calculated from initial bed slope (S 0 ) were provided as initial conditions, i.e., The flow depth, h at node 1 was determined by extrapolation from the already calculated values at the interior nodes using the MacCormack scheme. The downstream boundary condition was the constant depth, which was specified by h(x n ,t)=h 0 for t ≥ 0. This was based on the assumption that the channel was long and the bed transients would not reach the downstream end within the period for which conditions were computed and that the variation in flow depth would be negligible. The values at the interior nodes are used for extrapolation to determine to discharge and the bed elevation at the downstream end.

Model Performance Assessment
The Coding of the model has been checked by comparing it results with the already computed by the implicit model [10] namely, Chen, Chang and Richard and LUM [10]. The results are shown in Figure 3. It indicates that the algorithm performs well and the overall agreement is satisfactory.
Channel section is wide rectangular; the length of the reach, L=14 km; bed slope, S 0 =0.0003; and Manning's n=0.01.
The value of porosity p was assumed as 0.528 and the sediment load was calculated by using eqn. (5) as given below: The inflow hydrograph was defined as a steady flow of 6.2m 3 per second per unit width. The following equation was described by Chen [11] for second upstream boundary condition.
( ) Q s1 is the sediment input discharge at the upstream node, Q s2 is sediment discharge at node 2 calculated by eqn. (10), n represent the time step.
The initial flow along the reach was assumed to be 6.2 m 3 per second per unit width. The initial upstream and downstream depths used are given at different nodes ranges from 2.15 to 4.0 meters respectively. An assumed constant water stage was defined. The equation used is; h n and z n are flow depth and total bed elevation change at downstream end of reach since start of simulation, respectively. The sum of h n and z n remains constant. Initially z n at downstream is zero.
The deviation is anticipated to arise from different steps in linearization, solution methods choose (Implicit and explicit) and handling of insignificant terms in the governing equations. Term concerning for rate of bed level change with time is neglected in the governing equations used in this study while same is explicitly coupled in the governing equations used in the models developed by Ghumman [10]. It might be due to this insignificant term that may become significant at point of maximum deposition. Linearization of friction force, pressure force and convective acceleration terms may also be one reason. Wetted perimeter used for unit width of channel is a valuable variable that affects results.
Lustrous discrepancy trends of convective acceleration, pressure force and friction force terms of momentum equation during normal execution of model are shown in Figure 4. While, tremendous and senseless variation of these terms due to the instability when model stops execution is given in Figure 5.
Mostly stability is affected due to hydraulic parameters, physical boundary system and sensitivity to input data e.g. slightly change any of the input parameter and see how it goes. Inappropriate initial conditions also introduce false transients into the simulation which may lead to improper results.

Application to CRB canal
Chashma Right Bank Irrigation Project (CRBIP) is a large irrigation system spread over two provinces (Punjab and North West Frontier Province) of Pakistan. It is the 1 st irrigation system in Pakistan, initially conceived as a crop based supply system, and then followed by the design and construction of its physical infrastructure.

Parameter examination
The total simulation time was divided into low sediment flow period and high sediment flow period. High sediment flow period was from May to Sept whereas rest of the period being low sediment flow. For the year 2001, measured hydraulic and sediment data can be used to calculate the inflow sediment load. The sediment rating was prepared for the low sediment flow and high sediment flow separately. Inflow sediment rating used the outflow sediment rating for the next reach and so on. A single representative cross section was chosen, to avoid an enormous load of calculations. The geometric properties of all the individual sections over the entire reach were averaged for this representative cross-section. The reach length, L=15 mile (24 km); bed slope, S 0 =0.000125; Manning's n=0.016 The value of porosity, p=0.3; sediment discharge was calculated using eqn. (5) in which values of empirical constants were taken as; α=0.00000006; β=3.248; frictional slope was estimated by Manning's formula.
The initial depths were assumed as normal depth is calculated by Manning's formula using the assumed cross-section as given below.

( )
Pw=2*y+T; Pw is wetted perimeter and T Substituting these terms in the Manning's eqn.; The For dx =1000 ft.; dt=c*dt1; Any increase and decrease of dx and dt have tremendous effect on stability and execution of model. If c=0.75, it stops only after 48 time steps and for c=0.65 its simulation is completed upto 90 days. Stability problem was overcome by using c=0.55 that simulated for more than 198 days and 10841 time steps when ram quota exceeded with the message "Error during write, Disk  quota failure" as output file is too large and big. Hence direct display of results was necessary to avoid excess of memory error and for model debugging. However, unnecessary output data could be avoided through managed output file.
The difference between calculated and observed bed level changes is ±0.82 ft. It is due to the limitations of the model.
The bed profile computed with the model after simulation was compared with the measured bed profile. The computed and measured bed profiles are well matching with each other as shown in Figure 8.
The sediment deposition pattern depends upon the sediment inflow, hydraulic discharge, flow velocity, roughness, bed slope and bed material gradation. The sediment deposition pattern was estimated by plotting the change in bed elevation at different time intervals. The Changes in bed elevation at different time interval are plotted in Figure 9.

Verification
Using the parameters determined in the previous analysis, the model is applied to simulate the channel bed evolutions of CRB canal from 2002 to 2003. The total simulation period is 330 days. The comparisons of the observed and computed bed profile and bed change are shown in (Figures 10 and 11) respectively. The overall accuracy is good. The computed and measured bed profiles are well matching with each other as shown in (Figure 10).
The deposited sediment resulted in raising the bed of the canal at different locations. The graph was plotted for bed changes and it is shown ( Figure 11).
Running MacCk. Couple Model with this data is humble effort of this study. Generally set trends of results are fairly responsive but need further minute model calibration and refinement which is expected in future research to make the model a generalized one and viable for sediment as well as for hydraulic routing.

Conclusions
In this study the change in bed level of channel is studied with application of a one-dimensional mathematical model based on MacCormack scheme. The following conclusions can be made after summarizing the present study:

1-
Model developed in this study has fairly good qualitative and quantitative accuracy when applied within the limitations of the model. However, quantitative estimates may be demoralized sometimes.

2-
Results and model execution was more stable and flourish for relatively smaller computational time intervals.

3-
By using different equations for similar conditions will yield transport rates that might differ by several order of magnitude. Therefore, it is difficult to reach at a universal transport equation.

4-
The MacCormack scheme can be applied to develop a twodimensional model to study the two-dimensional problems and also a three-dimensional model to study three-dimensional problem.