A gravity-driven runoff and erosion model for sediment transfers at the catchment scale

The aim of this work is to simulate hydraulic transfers and upstream erosion sources in steep and erodible mountain watersheds with a physicallybased hydraulic model. In such environments, immature debris flows and shallow landslides can be the largest sources of sediments transported at the outlet. To simulate these phenomena, a gravity-driven erosion model and a 1D vertical infiltration model have been developed in the TELEMAC 2D numerical code. In this new erosion model, the motion of the granular flow is described with a fully dynamic system and a Coulomb-like bottom friction treatment, more adapted to the properties of the flow. The new model is first qualitatively evaluated on a theoretical test case: a steep plot with a slope break is used to evaluate the erosion and deposition dynamics of a single immature debris flow. Then, the model is confronted to field data on a real catchment (Draix, in the Southern French Alps). First, the infiltration model is successfully calibrated in order to simulate two different rain events. Then, the new erosion model is applied at the catchment scale. The numerical results show a very realistic behavior compared to the field observation, providing erosion in the upper parts of the hillslopes and deposition in the hydraulic network. This work opens promising perspectives, for example by coupling this new model with a classical and complementary velocity-driven model for the erosion, deposition and transfers in the hydraulic network.


Introduction
In mountain watersheds, extreme events can lead to the exportation of a large amount of sediment at the outlet.Indeed, with the combination of steep slopes and highly erodible soil in badlands the erosion can be very important.The presence of a large mass of sediment in the river flows can fill the dam reservoirs which affects the hydro-electricity production and causes safety issues.At the outlet of small watersheds in the head of river basin, the prediction of the sediment transfers is important to design hydraulic structures or to prevent erosion by hillslope management.
There are several processes involved in the sediment dynamics at the catchment scale.The erosion and deposition processes being heavily reliant on the flow characteristics, the choice of the hydraulic model is crucial.The choice of the Chen and Noelle's numerical scheme [1] is justified in [2] to solve the shallow water equations in the particular case of small water depth runoff on steep slopes.The surface flow and groundwater exchanges in a watershed needs to be taken into account to represent the correct hydraulic transfers.Concerning the sedimentary processes, there is a clear distinction in the literature between velocity-driven erosion in the rivers and gravity-driven erosion in the hillslopes.The velocity-driven models are easily coupled with the hydraulic models (see [3], [4], [5]).Concerning the gravity-driven erosion, plenty of models can be found to represent the detachment, deposition and dynamic of muddy or granular flows (see [6], [7], [8], [9]).However, none of them is integrated to a global model with a hydraulic and hydrological description.
The goal is to build a model capable of representing the hydrology of a watershed, the upstream gravity-driven erosion sources that supply a classical erosion and sedimentation model in the hydraulic network.Then, this model will be able to give explicit values of the hydraulic and sedimentary parameters distributed in time and space.It offers the possibility of identifying the critical erosion zones but also anticipating watershed responses to a management strategy like slopes revegetation or dike construction.
In this paper, the modeling of all the processes above-mentioned and the coupling method are presented.Then, a theoretical test case and a field application are used to assess qualitatively the model.

Presentation of the test cases
The first test case is a straight channel with a length of 5 m and a width of 1 m.This channel is divided in two parts: an upstream part with a 50% slope and a downstream part with a 5% slope.The break in slope separating these two parts is in the middle of the channel.The figure 1 represents the geometry of the domain.On this domain, a steady rain, with an intensity of 100 mm/h, is applied during 10 s.The spatial discretization of the channel is a triangular mesh with a mean space step of 10 cm.The upstream and lateral boundaries are considered as walls.The flow being supercritical at the outlet, the boundary is treated with a free condition.
The second test case is the Laval watershed, which is a sub-catchment of the Bouinenc watershed, located on the Draix site in the Southern French Alps.Its total area is about 86.4 ha and the mean slope is 58%.Its bottom elevation is presented in the figure 1 The soil is mostly constituted of black marls and 68% of the surface is a bare soil.Thanks to data of the Draix-Bleone observatory [10], it is possible to evaluate the model by comparing the results to measured discharges and sediment fluxes from field campaigns on this catchment.At the outlet, the discharges are available for many rainfall events with a time step of 60 s.The rainfall intensity associated to the discharge is also measured every 60 s.Two different events are selected, one spring rain following three rainy days with a high soil moisture initial condition, and a summer storm with a dry initial soil state.

Overland flow simulation
To simulate rain induced runoff, the following shallow water equations are solved: where h is the water height in m, t the time in s, u the flow velocity in m/s, R the rain intensity in m/s, I the infiltration rate in m/s, g the gravity constant in m/s 2 , z the bottom elevation in m and S f the friction slope.
The shallow water equations [11] are used in most applications of runoff modeling ([3], [12], [13]).As it is composed of two conservation laws, a finite volume method ( [14]) is used to solve the system numerically.The challenges of the numerical resolution of these equations is to preserve the positivity of the water depths, be able to treat dry zones, be massconservative and be well-balanced in the sense of [15], meaning preserving the hydrostatic balance of the lake at rest.Few schemes verify all these properties.In order to have a robust resolution of the shallow water equations on an entire watershed, a better compromise needs to be found in the literature between precision of the numerical scheme and the good properties mentioned above.The Audusse et al.'s scheme [16] is widely used for its ability to preserve the stability of the lake at rest and verify all the required properties.However, [17] shows that this scheme can be in default for some combination of slope, water depth and mesh size.Indeed, when the water depth is smaller than the bottom difference between two neighbor cells, the slope gradient is not taken into account.The Chen and Noelle's scheme [1] has been developed to overcome this limitation by modifying the hydrostatic reconstruction in this particular case.

Infiltration model
Inspired by the the Green and Ampt's model [18], a vertical 1D model is presented in [19] and [13].This model, used to represent the infiltration, is computed at each cell of the domain.The infiltration velocity is discribed like with K the soil conductivity under less than 1 cm of hydraulic head in m/s, h f the capillarity head in m and z f the wetting front in m.Then, the wetting front is updated following the equation: with I C the cumulated height infiltrated since the beginning of the event in m, θ s the saturated soil moisture or the porosity and θ i the initial soil moisture.The soil is vertically divided in two layers, a first with an associated width Z c in m and conductivity K c and a second with a width considered as infinite and a conductivity K s .The conductivity K varies with the wetting front evolution according to this equation:

Debris flow modeling
The detachment criterion is calculated in accordance to [7].At each cell interface, the slope of the ground θ between two nodes of the mesh is calculated.Then, depending on the water depths and the soil properties of these nodes, a stability angle θ 1 is evaluated and compared to θ.The soil is considered stable if θ 1 < θ and unstable if θ 1 ≥ θ.The formula to calculate the critical angle is: where , with c the apparent cohesive strength of the sediment layer in Pa, k a numerical coefficient near unity, C * the maximal volumetric concentration allowed in the flow, ρ s the sediment density in kg/m 3 , ρ the water density in kg/m 3 and φ is the internal friction angle of the sediment.The figure 2 shows the behaviour of the detachment criterion depending on the slope and the water depth for C * = 0.65, c = 35 Pa, ρ s = 2650 kg/m 3 , ρ = 1000 kg/m 3 and k=0.85.If the rainfall intensity is large enough to generate a runoff with a water depth in the unstable zone, a layer with a thickness e is eroded.This thickness is proportional to the water depth following with C e = ρ tan θ (ρ s −ρ)(tan φ−tan θ) the equilibrium concentration.Concerning the deposition velocity, the formula is similar to the one used for the cohesive sediment deposition: where u s is the debris flow velocity in in m/s, u c is the critical deposition velocity in m/s, m is a coefficient less than 1, C is the sediment concentration in the flow, V s is the settling velocity in m/s and (.) + = max(0, .).The debris flow motion is described by its velocity u s in m/s and depth h s in m.The evolution of these variables is governed by the shallow water equation ( 1).The source term is modified with the erosion E = e ∆t and the deposition velocities D. The friction term becomes: the Chézy friction part is treated semi-implicitly and the coulomb friction is added explicitly to the scheme (see [20]).

Coupling method
The coupling between the infiltration and the runoff is made through the source term I in the mass conservation equation.Concerning the debris flow dynamics, it interacts with the runoff equations by modifying the bottom elevation like: The entire model is presented in the figure 3. The available eroded layer is limited by the wetting front calculated with the infiltration model.Indeed, only the saturated part of the soil is considered available for the debris flow simulation.

Channel test case
To fulfill the condition of instability, we set c = 1 Pa, k = 0.85, ρ s = 2650 kg/m 3 , ρ = 1000 kg/m 3 , C * = 0.65, tan φ = 0.8, u c = 1 m/s, m = 0.3 and V s = 0.01 m/s.The saturated zone is initially set at the constant value 0.1 m and the infiltration is not considered.The figure 4 presents the bottom profile along the channel at the end of the simulation.In the upstream part, the erosion starts at 20 cm of the upstream boundary, when the depth of the water is sufficient to reach the unstable state.Then, a fully eroded zone is observed until the deposition starts, when the velocity decreases with the gentlest slope.These results are consistent with what is expected from the model.

Watershed hydrology
The hydrological model is confronted to field data by comparing the measured and simulated outlet discharges on two events.The events are selected because they are the most erosive of 2012, in terms of sediment volume exported at the outlet.The first event is recorded the 29th of May succeeding six rainy events from 21st to the 27th of May.The maximal value of the rain intensity is 84 mm/h and the peak discharge is 3 m 3 /s.Concerning the other event, it is a summer storm of the 25th August with an initial dry state because the last recorded event is the 25th July.The maximal intensity of the rain is 156 mm/h and the peak discharge is 6.6 m 3 /s.
To simulate these events, we set the properties of the soil constant, modifying only the initial soil moisture (θ i in the infiltration model).Inspired by field data [21], the parameters are chosen to have a porous surface layer and a more structured base layer: Z c = 80 mm, θ 1 = 0.35, K C = 30 mm/h, K s = 1 mm/h, θ 2 = 0.25, h f = 50 mm.The initial soil moisture is set to θ i = 0.22 for the spring event and θ i = 0.03 for the summer storm.The figure 5 presents the measured and simulated outlet discharges for the spring event and the summer event.The simulated results are in good agreement with the observations.Despite some important simplifications of the model, the outlet hydrograph are well represented.Indeed, the soil properties are considered as uniform in space, as well as the friction coefficient.In addition, the exfiltration is not represented and compensated by a low conductivity of the base layer.The representation of the vertical structure of the soil is sufficient to have a satisfactory reproduction of the hydrographs.

Watershed gravity-driven erosion
The debris flow model is now applied to the Laval watershed The figure 6 shows the erosion and the deposition in the Laval catchment.The erosion is mainly localized at the upstream of the watershed, in the small gullies.Then the sediments are deposited in the main channels, distributed in the entire watershed.

Conclusion
To model sediment transfers on mountain catchment with a physically-based approach, the upstream erosion sources needs to be represented.Numerous infiltration and debris flow model exist but few can represent both phenomenon together.
A derivation of the Green-Ampt hydrological model [8] is first presented.It is a 1D vertical model which gives a precise representation of the properties of the soil.This model has been tested on the Laval watershed, on two different rainy events.The simulated and measured outlet discharges are in good agreement, keeping the same parameters of the soil properties and modifying only the initial soil moisture between the two events.
Then, a gravity-driven erosion model has been tested.It consists in simulated the debris flow formation, evolution and deposition in the model.On a simple test case, the behavior of the debris flow is expected and the form of the deposition is close to what can be observed in the experiments.Then, on the Laval watershed, it is interesting to see that the debris flows are generated in the gullies upstream and deposited in the hydraulic network.The amount deposited can be considered as an input for a velocity-driven model for erosion and sedimentation in the network.
This model still need to be validated.The concentrations at the outlet are available on several events on the Laval watershed.Data on other catchment can also be used, like tracking of deposition and erosion on the main river event by event.

Figure 1 .
Figure 1.Bottom elevation of the channel test case (left) and the Laval catchment (right)

Figure 2 .
Figure 2. Detachment criterion for the debris flow depending on the slope and the water depth

Figure 3 .
Figure 3. Schematic representation of the model

Figure 4 .
Figure 4. Bottom elevation and total evolution at the end of the simulation

Figure 5 .
Figure 5. Discharge at the outlet of the Laval catchment for the two rainy events . The goal is to observe quantitatively if the upstream sources are well represented and if the sediment are deposited in the hydraulic network.The model is tested on the summer event, because of its high rain intensity.The parameters used for the detachment and deposition are: c = 35 Pa, ρ s = 2650 kg/m 3 , ρ = 1000 kg/m 3 , k = 0.85, m = 0.3, tan φ = 0.8, C * = 0.65 and u c = 1 m/s.

Figure 6 .
Figure 6.Bottom evolution at the end on the simulation in the Laval catchment