COURLIS : a 1 D suspension and bedload code

COURLIS is a 1D sedimentology module coupled with MASCARET, 1D hydraulic code of the TELEMAC-MASCARET open source system. The code has been developed for more than 10 years, mainly for suspension sediment transport. Recently, the need of a 1D bedload code has been identified to model the long term evolution of rivers and reservoirs (several decades). New numerical schemes were implemented, some improvements were done in the geometry evolution algorithms. In terms of performance and robustness, the best scheme implemented is a finite volume upwind/downwind scheme. Several solutions are implemented to reduce calculation time. This new version of COURLIS for bedload transport was validated successfully on test-cases (Soni and Newton experiments). A real case has been simulated during an 11 year period. The calculation time is very similar to those obtained with codes tested in the benchmark and the results are in a good agreement with measurements and other code results. COURLIS (suspension and bedload transport) will be released in the next version of the TELEMAC-MASCARET open source system and so, it will be freely available for sedimentology community. Further developments are planned in 2018.


Introduction
Modeling sediment is crucial to deal with environmental and safety issues such as saving biodiversity and flood risks.Even now, with the increase of computational capacity and development of 2D and 3D CFD codes, one dimensional sediment transport codes are efficient engineering tools to predict reservoir and river bed evolution.They are more efficient in terms of computational times and generally, 1D models are simpler to build and require less field data.
COURLIS has been developed for more than 10 years.The first objective of this code was to understand the processes of fine sediments in reservoirs, i.e. deposition, erosion and advection-diffusion of fine sands and mud.The code has been used to prepare the emptying of several French reservoirs [1,2].The main particularity of COURLIS as 1D code is that output results include not only bed evolutions on the river longitudinal profile but also on initial cross-sections, preserving real geometry of the river.
Recently, due to the emergence of bedload continuity issues, a 1D bedload code, which would be robust and with fast computing time, is needed.It should be able to model long term evolution of rivers and reservoirs (several decades).This tool must help engineers to better understand and predict sediment deficit or deposit, in order to improve sediment management in hydroelectricity reservoirs and in rivers impacted by dams.
Based on a benchmark of several 1D codes (commercial and research, including HEC-RAS), it was decided to develop a bedload option in COURLIS.The Exner equation has been implemented in the code and, an efficient and well-adapted numerical scheme has been used to solve it.
In the followings, COURLIS suspension (which stands for suspension option) will be described and some applications will be mentioned; then, COURLIS bedload (which stands for bedload option) new developments and validation test cases are presented; finally, some perspectives and future developments are given.

COURLIS Suspension
COURLIS suspension is coupled with MASCARET [3], the hydraulic 1D software of the opensource TELEMAC-MASCARET system (www.opentelemac.org).Both codes are internally coupled.That is to say, during one time step, first MASCARET solves the 1D shallow water equations with finite volume scheme [3], then, with variables calculated by the hydraulic model, COURLIS calculates sediment transport processes and bottom evolutions.(1) Where is the cross sectional flow area, the suspended concentration, the flow rate, the dispersion coefficient, and respectively the Erosion and Deposition linear rates.
Deposition and erosion are calculated according to the type of sediment:  For cohesive sediment, empirical formulae are used : Krone (1962) formula (deposition) [4]: Partheniades (1965) formula (erosion) [5]: Where is the erosion rate coefficient,   the settling velocity,   and   deposition and erosion critical shear stresses.The local shear stress is written: Where is the gravity acceleration, the local water depth (variable in the cross section), the mean velocity,   the Strickler coefficient and  ℎ the hydraulic radius.


For non-cohesive sediment, the flow is assumed to be saturated by non-cohesive sediment at an equilibrium concentration.The Engelund-Hansen (1967) formula [6] is used to calculate a transport capacity : where Where   is the grain friction coefficient and could be written   = 26  90 ,  is the sediment mean diameter,   is the sediment density,  is the water density and  is the reduced density ( = (  − )  ⁄ ).This transport capacity is used to determine the equilibrium concentration: Finally, with the equilibrium concentration, deposition and erosion rates are determined: Then the bottom evolution depends on erosion and deposition rates: Where   is the bottom elevation,   the concentration of deposition layer and   the concentration of the layer where the sediment is eroded.
Besides, the bank deformation is taken into account using a simple model.The bank slope is compared to a stability slope (submerged or emerged): if the critical slope is exceeded, sediment deposit is supposed to collapse immediately.
COURLIS suspension is described with more details in [7].
COURLIS suspension has been used in EDF on several cases.It has been used to prepare flushing or emptying reservoir operations, it allows to predict output sediment concentrations and to limit the environmental impacts [1, 2, 8].COURLIS has also been used to predict flood level to study possible dike overflows considering bottom erosion.Settling basins with a 1D configuration can be modeled also with COURLIS suspension.Propagation of cohesive sediment in rivers has also been studied with COURLIS suspension, for example in the case of Arc and Isère rivers [9].

COURLIS bedload is coupled with MASCARET, using Exner equation:
(, ) +     (, ) = 0 (9) with Az (m 2 ) sediment volume in a cross-section and Qs (m 3 /s) volumetric sediment flux integrated along the transversal direction.The flow mean velocity is calculated by MASCARET and used by COURLIS to calculate the mean shear stress on the river bed.The bottom evolution is calculated with the Exner equation ( 9) and a relation between Az and Z.This new bottom elevation is sent to MASCARET, and so on.To close the Exner -shallow water equation system a bedload formula is required.Four different bedload formulas are coded in COURLIS bedload (Meyer-Peter and Müller (1948) [10], Engelund and Hansen (1967) [6], Recking (2011) [11], Lefort (2015) [12]).
To solve the Exner equation, a 1D finite volume discretization is adopted: numerical flux evaluated at x i+1/2 interface, Δt a timestep and Δx a space step (X i n = X(nΔt, iΔx) ).Three different schemes were tested to solve numerical fluxes: -a Roe scheme [13] (named ROE in the following), -a staggered scheme [14,15] (named STAG in the following), -an uncentered scheme [16,17], upwind or downwind according to the flow regimecritical or subcritical -before and after the cells interface (named UNCENTERED in the following) The third scheme was identified during the development as the most stable, efficient and robust.In order to save computational time, in addition to the implementation of bedload formula and numerical schemes to solve the Exner equation, a new method has been developed to optimize the update of the bottom elevation.Initially, MASCARET subroutines calculating 1D geometric quantities from real geometry (integration of 1D cross-sections) were used to update the bottom from COURLIS calculations.However, these routines slowed dramatically the simulations.The new method includes two steps.
First, a new threshold condition has been defined.The bottom is now modified only when sediment erosion or sediment deposition is higher than a percentage of the water depth defined in the code.
Secondly, a method to calculate deposition was implemented.Initially, deposition was done according to water depth.This definition is coming from COURLIS suspension.It allows to take into account a constant deposition on each node of the cross-section under the water level (named Delta Constant Bottom Evolution method, DCBE in the following).A constant elevation deposition was implemented (named Level Constant Bottom Evolution method, LCBE).This method is a priori more suitable for bedload sediment transport and moreover, it is saving computational time.It allows to calculate the bottom evolution easily and the 1D quantities required by MASCARET like wetted areas are directly calculated without the 1D cross-section integration subroutines of MASCARET.

Validation test cases
A set of two experiments was selected to test the bedload developments in COURLIS.Specific attention is given to computational time because the code will be used in engineering studies to simulate long periods of time.

Soni Experiment
The aim of this experiment, carried out by J.P. Soni in 1980 [18,19], is to study aggradation of sediment.The sediment supply is increased above the equilibrium sediment transport capacity of the flume.The evolution of the deposit is monitored during the experiment at 30, 60 and 90 min.    Cs=4.88 kg.m -3   Results are presented in Fig. 1 for the uncentered scheme and the Meyer-Peter and Müller formula.The comparison shows a fairly good agreement with measurements.Table 2 summarize the calculation times.The uncentered scheme drastically reduces the calculation time.The subcritical kernel of MASCARET and the bottom evolution threshold reduce also this time.The new bottom evolution method, LCBE, does not increase performance in this case.

Newton Experiment
A degradation is often observed downstream of dams due to the interruption of sediment flow from upstream.The aim of the experiment carried out by T. Newton in 1951 [20], is to study the erosion process in a flume in order to better understand this phenomenon.The evolution of the bottom of the flume following the interruption of the sediment supply is observed during 24 hours.Results are given on Fig. 2 for the uncentered scheme and the Meyer-Peter and Müller formula.The comparison show a fairly good agreement with the measurements.4 gives calculation times.The conclusion are the same as for the Soni case except for the bottom evolution threshold that does not increase performance here.

Real reservoir test case
This test case is based on an EDF reservoir.The study area is 3.9 km long (1.8 km of reservoir and 2.1 km of river upstream the reservoir).The objective is to simulate the bottom evolution during 12 years.This period includes reservoir flushing events.Consequently, in different zones of the model, the bed alternates between periods of deposition and erosion of sediments generated by high flows.The river bed presents strong discontinuities (significant increases of the slope, variations of width).These discontinuities could create hydraulic instabilities.They are a real challenge and the main difficulty of this test case for the code.The first step of the study addresses the stability of the simulation and the calculation time.The second step, calibration and comparison with the measurement results will be carried out in the future.
Two meshes are used: a fine mesh, with one cross section every 20 meters, and a coarse mesh, with one cross section every 100 m.The geometric and physical parameters of the model are presented in Table 5. Results are presented in Fig. 5 (the lowest points of the cross-sections are plot).There are few differences between the two meshes (less than 10 cm).
The difference between the two bottom evolution methods is predictable.The LCBE method tends to suppress the trenches in the cross-section.Calculation times are given in the table 6.The bottom evolution method does not decrease the calculation time.Important parameters are the mesh and the bottom evolution threshold.Computational time improvements may be possible in the future with the use of the subcritical scheme of MASCARET.The computational time (and results) obtained with the coarse mesh, the UNCENTERED scheme, the supercritical hydraulic kernel and the LCBE method is acceptable and allow long term operational application.

Conclusion & Perspectives
COURLIS is a one dimension code allowing the calculation of suspension sediment transport and bedload sediment transport for rivers and reservoirs.It has been developed in order to be stable and fast.
Laboratory test cases (Soni and Newton experiments) show that COURLIS bedload gives results in a fairly good agreement with measurements.Moreover, several numerical schemes have been tested and the uncentered scheme showed good properties (stability, robustness, efficiency, etc.).Some developments have been done to reduce calculation times.A real test case shows that developments carried out COURLIS allow the software to simulate long term evolution of gravel rivers with reasonable calculation time.Future work will focus on -Implementation of other bedload formulas, -Integration of an extended grain-size model, -Management of the transition from a rocky bottom to an alluvial bottom (adaptation of the friction law).COURLIS will be fully open source in TELEMAC-MASCARET platform (www.opentelemac.org).Better stability and better efficiency are expected with the next updates.
Finally the bottoms calculated by Courlis are used in MASCARET to calculate the next hydraulic state.COURLIS suspension is based on an advection-diffusion equation for the sediment concentration transported by the flow for two classes of fine sediments, mud (cohesive sediment) and sand (non-cohesive sediment):

Fig. 2 .
Fig. 2. Bottom evolution during the Newton experiment: measurement and calculations (MPM and UNCENTERED) Table4gives calculation times.The conclusion are the same as for the Soni case except for the bottom evolution threshold that does not increase performance here.

Fig. 3 .
Fig. 3. Initial and final bottom elevation of the real test case calculation with fine and coarse mesh and

Table 2 .
Soni experiment computational time

Table 3 .
Newton experimental parameters

Table 4 .
Newton experiment computational time

Table 5 .
Real test case parameters

Table 6 .
Real test case computational time