A hydro-morphodynamic model integrating extended sediment particle size distribution and flocculation processes for better simulating hydro-sedimentary fluxes

Many studies focusing on suspended sediment transport modelling in river systems only consider one class of sediment grain size. Rather recently, the SISYPHE sediment transport model has integrated sand-mud mixture transport processes using two classes of sediment. However, this new modelling framework still suffers from limitations, and increasing the number of sediment classes would arguably improve sediment transport and therefore riverbed evolution simulations. Moreover, current sediment transport models do not simulate sediment particle aggregation and disaggregation processes while these can play an important role in sediment transport. Integrating these new concepts would then contribute to significant improvements to river bed morphodynamics and sediment transport modelling. In this study, we further develop the SISYPHE model by extending the sediment particle size distribution to ten classes and integrating flocculation processes (coupling with the flocculation FLOCMOD model). The preliminary results we present in this paper are based on a large-scale flood event, which occurred in river Orne, north-eastern France. We clearly show that the proposed developments of SYSIPHE improves qualitatively and quantitatively the predictions of sediment transport and riverbed morphodynamics.


Introduction
In Europe, many riverine environments suffer from sediment contamination due to past intensive metallurgical activities.In this context, restoring and maintaining water quality remains challenging due to the presence of contaminated sediments in many areas.The stock of contaminants that has historically been buried in sediments at different depths can potentially be re-suspended during floods [1][2].As a consequence, there is a need for monitoring and predicting the potential resuspension of river sediments and relative persistent pollutants in contaminated rivers.
One of the main limitation of the actual version of SISYPHE is to consider only two classes of sediment as cohesive/non-cohesive sediment mixtures, whereas the sediment transported in rivers is composed of many different sediment particle sizes with distinct physical properties.In presence of heterogeneous sediment cover, the use of two grain sizes for hydro-morphodynamic model simulation is quite questionable.Indeed, the two main factors influencing erosion fluxes and suspended sediment loads are the sediment erodibility and the sediment particle fall velocity both closely related to the sediment diameter.As a consequence, too simplified representations of sediment particle size distribution arguably leads to over-and/or under-estimated erosion fluxes.To avoid huge errors in sediment transport simulations, models need to consider more realistic distributions of sediment particle size, i.e. including more sediment particle classes.
Moreover fine particles are not only transported as primary particle but are strongly impacted by flocculation processes due to the combined effects of turbulence and biophysico-chemical forces.Consequently, the suspended sediments are principally transported as aggregates (i.e.flocs).The flocculation process also influences the flux and the fate of particles because the formation of flocs influences the transport of sediment via variations of floc density, fall velocity and dispersion in the water column.
In this proceeding, we propose a further development of the module SISYPHE, which is now able to simulate sediment transport based on : i) a 10 classes sediment size distribution with specific density values and ii) the consideration of flocculation processes in contrasted hydrological conditions.The performances of the two new model configurations are discussed and evaluated based on the comparison with the standard version of SISYPHE.

Modelling framework
The hydrodynamics are simulated using the TELEMAC-3D model solving the Navier Stockes equations in a hydrostatic mode.The modelling of the turbulence is carried out using constant horizontal and vertical viscosities (1.10 -6 m 2 /s).The bottom friction coefficient is calculated using the Nikuradse's effective roughness-length depending on the median grain size in the river bed (in absence of bedforms).In the riverbanks, the apparent roughness is fixed to 4cm to take into account of the presence of dense vegetation.The bed friction estimated using the bottom sediment repartition has already been successfully applied in estuary system with a highly heterogeneous sediment cover [4].
The morphodynamic and sediment transport model SISYPHE [5][6] is dynamically coupled with TELEMAC3D.SISYPHE computes the morphological bottom evolution using the Exner equation [7][8], calculates the bedload flux using the Meyer-Peter-Müller equation [9] with a compensation for slope effect (magnitude and direction).It calculates the advection and diffusion of the suspended sediment concentration (SSC), and computes the sediment erosion and deposition fluxes based in particular on the sediment particle fall velocity.The sediment fall velocity is computed based on sediment diameter and density according to Van Rijn [10] formulation.The model of Van Ledden [11] is applied to calculate the erosion and the deposition rates of a cohesive/non-cohesive sediment mixture, based on the proportion of cohesive sediment in the bottom sediment.The sediment mixture is eroded as purely cohesive sediment if the mud fraction is beyond 50% and as purely non-cohesive sediment if the mud fraction is below 30%.Between 50% and 30% of mud fraction, a linear interpolation is applied.When the sediment is cohesive the formulation of Krone and Partheniades [12] is used and the initiation of motion is based on the critical shear strength of the mud, which is measured in situ using a scissometer.The erosion rate for cohesive sediment and the critical shear stress for deposition are set to 2.4 10 -5 kg.(m.s) -2 and 0.005 Pa respectively.When the sediment is not cohesive, the formulation of Smith and Mc Lean [13]  3 1 3 ⁄   the dimensionless particle diameter,  the kinematic viscosity of the fluid, g the gravitational constant and  =    ⁄ the relative density.As initial condition, the bed is structured in ten horizontal layers with thicknesses defined as follows: with ES(N) the thickness of the layer number N and D50(N) the D50 of the layer N.Each layer thickness is characterised by the volume percentile of each sediment class.In case of non-cohesive mixture, the second layer is eroded only if the coarsest sediment (i.e.sediment class with the highest particle diameter) has disappeared from the top layer.As the riverbed is made of a cohesive/non-cohesive sediment mixture the mud concentration and the critical shear stress is defined for each layer.Using a scissometer, the critical shear strength of mud erosion was estimated at 0.48 Pa for the top layer and 0.84 Pa at 15 cm depth (linear interpolation is performed to attribute to each bottom's layer the appropriate critical shear stress).In case of cohesive mixture, there is no limitation of eroded thickness.
The flocculation model FLOCMOD [16], coupled with SISYPHE, simulates the aggregation/disaggregation of mineral particles via a population model [17] based on a finite number of floc size classes with a logarithmic distribution.15 floc diameter classess can be considered from 1.5m to 1600m.This limits the number of flocs classes and hence the computational time.FLOCMOD is constraint with a 5 m primary aggregate, which can be or aggregated to coarser flocs (diameter > 5m) or disaggregated to finer flocs (diameter < 5m).The relationship between the floc diameter and the floc density is assumed to be based on a fractal dimension [18].The shear rate G (s -1 ) that represents the flow turbulence known as the main mechanisms controlling the floc formations is given by the following equation:  is the dissipation rate with K is the Von Karman constant (0.41), z the distance above the bed and  * = √    ⁄ is the friction velocity.

Study area and available data and model set up
The Orne River, located in the north-eastern part of France, drains about 1270 km 2 at the confluence with the Moselle River.Within the Orne basin, a 4 km-long control section has been chosen for the modelling exercise of suspended sediment transport.In this area, the riverbed has an average width of 30 m.This river reach is composed of two large meanders.Its downstream boundary is equipped with a dam that contributes to the formation of mud banks due to suspended load deposition as a result of reduced streamflow velocities.Pebbles, coarse gravels, sand and a small silt portion mainly compose the streambed.The riverbanks are mainly composed of a sand/mud mixture covered by dense vegetation.The riverbanks are locally made of concrete or silted up rockfill.
Since January 2014 monitoring efforts have been concentrated on continuously recording streamflow and water turbidity, as a proxy of SSC, and more punctually on measuring bottom sediment particle size distribution, river bathymetry and sediment deposition at several locations on selected banks and in the riverbed.The continuous data used in this study refer to the important flood event, which occurred in the Orne River basin in June 2016.During this event a peak discharge of 206 m 3 .s - has been recorded.An ISCO © automatic water sampler was used to evaluate the turbidimeter's ability to perform an indirect measurement of the SSC.Because of the difficulty of keeping the equipment in such hydrological conditions, the ISCO © water sampler did not sample this particular flood event.Consequently, the experimental data are suffering from large uncertainties and an overestimation of turbidity values are expected.The figure 3 shows the monitored turbidity and measured SSC performed synchronously at the downstream boundary of the studied river section since 2014.The corresponding turbidity values range from 0 to 300 NTU.This linear regression between the two datasets exhibits a Pearson's correlation coefficient of 0.96.Unfortunately, those turbidity values previously recorded remained well below of those recorded in June 2016.Fig. 3. Linear regression between the river water turbidity and the SSC measured at the downstream boundary of the modelled Orne River section.
During this flood event, 7 water samples were collected and the corresponding suspended sediment grain size distribution was estimated using a laser scattering SYMPATEC that combines two lenses R2 (0.45 µm to 87.5 µm), and R5 (4.5 µm to 850 µm) after 20s of ultra-sonication.The grain size distribution in the bottom sediment (Fig. 1) has been estimated on collected sediment samples in representative locations of the river section.The sediment density has been measured for each defined grain size class (Fig. 2).The TELEMAC-MASCARET is based on a finite element unstructured mesh allowing for representing complex geometry.For the study area, the unstructured mesh is composed of 16492 nodes and the distance between two adjacent nodes ranges between 7m and 25m.It has been generated using the mesher POLYMESH © (developed by A. Roland, T.U.Darmstatd).The six bridge piles (2 bridges) located in the domain are also represented in the mesh model.
Three model configurations have been designed in order to assess the influence of both the sediment size distribution and the flocculation on the model predictions.The parameterisation of the SISYPHE and TELEMAC3D models is the same for the different modelling configurations.The first configuration (2CL) corresponds to the standard SYSIPHE configuration, which only considers two classes of sediment particle sizes having distinct sediment densities.In the second configuration (10CL) ten classes of sediment particle size are considered and the input SSC time series are defined for each class.The third configuration (10CLF) combines the latter (10CL) and the added FLOCMOD module.
The initial riverbed sediment composition (Fig. 1) received two distinct repartitions: one for the riverbed and one for the riverbanks.The two grain size classes of the configuration 2CL are determined as follow: the coarser grain is the median diameter of the non-cohesive sediment for the riverbed and the smaller grain is the median diameter of the cohesive sediment of the riverbank.Instead of imposing an initial SSC over the domain at the beginning of the flood event simulation, it has been decided to carry out a hot start over six days.The hydrodynamic model is forced by the recorded streamflow time series at the upstream boundary and the downstream recorded water level time series (a few meters upstream of the dam).As SSC measurements were unavailable at the upstream boundary, SISYPHE is forced upstream by those derived from the turbidity measurements at the downstream boundary.

Sensitivity to the number of sediment class (10CL vs 2CL)
Increasing the number of sediment classes in the riverbed from two to ten leads to very different SSC simulations, bed evolution, and particle size distribution in resulted bottom riverbed sediments (Fig. 4 to 7).
As can be seen in Fig. 4 the 2CL configuration significantly overestimates the SSC during the first small sediment peaks.In particular, this configuration exhibits at the beginning of the simulation additional small peaks that are not visible in the other configuration simulations and during the field measurements.In addition, the 2CL simulation results in very limited riverbed and riverbanks evolution (Fig. 5c).This is likely due to the fact that: i) the small particles are systematically eroded and never settle down due to their really small mass and ii) the larger particles are not eroded due to their higher mass.Therefore, the finest particle class is "washed out" at the beginning of the simulation and do not contribute to the sediment deposited on the riverbed at the end of the flood event.The main peaks of SSC of the 10CL configuration is well in phase with the observations (Fig. 4).However, the intensity of each peak is underestimated.This is likely due to the fact that the upstream boundary condition is imposed based on downstream measurements and we might expect a higher upstream concentration during such a flood event.In addition, we found that the 10CL simulation exhibits a significant spatial rearrangement of sediment particle size distributions in the riverbed (Fig. 5b).From one sediment class to another, preferential erosion and deposition zones change and their spatial patterns present a better accordance to the real morphodynamic of the river.Especially at the end of the simulation, a gradient in the sediment size is visible on the banks: the settled sediment is finer at the top of the bank and vice versa.In addition, the comparison between the simulated particle size distribution in the bottom sediments at the end of the simulation and the one observed in a sample taken from a deposition area after the flood (composed mainly of carbonate sediments of 500-200μm and sand of 200-100μm, data not shown) shows a good agreement.
Same agreement can be see for the suspended sediment size distribution.Fig. 6 compares measured (as described in 4.2) and simulated suspended sediment size distributions.As the spectral bands of the two lenses overlap between 4.5μm and 87.5μm the two combined measurements provide information on measurement uncertainties.In Fig. 6, the maximum deviation between the cumulative density functions (CDF) of the two measurements (R2 and R5) is 20%, this corresponds to a largest difference between the two measurements about 25μm.The 10 CL configuration correctly simulated the suspended sediment size distribution, as it underestimated only slightly (less than 4%) the CDF for particle diameters between 5 and 30μm.The 2CL configuration is unable to provide such a CDF because it only can transport one class of sediment, making the CDF starting directly at 100% with the class 5μm (75% of error on this granulometric class).This result lends more weight to the importance of using a more realistic particle size distribution in morphodynamic models.
Comparing the simulated bathymetry evolution maps is not straightforward for such a small river (Fig. 5c).The use of bathymetry evolution during a flood event instead of absolute bathymetry offers the possibility of differentiating simulated erosion and deposition in the river bed and also of comparing their respective quantity of eroded and deposited material estimated during the simulation over the entire area.For this purpose, the evolutions coming from the 10CL simulation is used as a reference in a scatter plot (Fig. 7).The areas corresponding to deposition, erosion and that stay stable have been plotted using different colours to ensure an easy interpretation of the results.Therefore, the evolution between the configuration 10CL and 2CL represented in the scatter plot clearly shows a very limited correlations (respectively 0.21 and 0.001 for the erosion and deposit) between the two configurations.

Sensitivity to the flocculation (10CL vs 10CLF)
No substantial difference between the 10CL and 10CLF configuration is observed.The two estimated total SSC present similar temporal evolutions (Fig. 4).In addition, in Fig. 7, the scatter plot comparing changes in bed elevation between the 10CL and 10CLF models shows a very high correlation.However, on the one hand, Fig. 6 shows a slight improvement in the SSC particle size distribution from 10CL to 10 CLF configurations.On the other hand, the dimensions of the flocs formed using FLOCMOD are in accordance to the measures, without any overestimation of aggregated flocs.The 10CLF configuration is able to represent a realistic size distribution of the aggregated flocs (>5m) during the event.But it also show that flocs can be disaggregated during peaks, due to elevated shear stress (about 20s -1 ), which tended to break the coarse flocs instead of aggregating the particles.Higher flocculation efficiency can be expected in smaller floods, which would present lower shear stress.This hypothesis will be tested in future flood events.Moreover, due to a lack of continuous field observations on floc size distribution, the actual 10CLF model configuration only considers a primary aggregate size of 5m at the boundary conditions defined in the model and hence do not take into account a real floc size distribution.The next step would be to constrain such a floc size distribution in the developed 10CLF model in order to improve the dynamics of floc aggregation/disaggregation along the studied river section during flood events.

Conclusion
This study evaluates the sensitivity of the TELEMAC-MASCARET modelling suite to the number sediment size classes considered in the sediment transport module SYSIPHE and to the flocculation process by combination between SYSIPHE and FLOCMOD.
First we compared two model configurations with respectively 2 and 10 sediment size classes (without flocculation) based on simulations of a rather high magnitude flood event in a small river show that using only two sediment class leads to spurious simulated erosion fluxes.Furthermore, the deposition/erosion areas are highly impacted by the chosen configuration (10CL and 2CL) as no correlation where found between these two models.
The integration of a flocculation module in the sediment transport module of the modelling suite was further analysed.During the simulated flood event, the influence of the flocculation process on the total SSC and on the riverbed evolution is rather limited, but it allows for better addressing the particle size distribution of the suspended sediment.It is important to highlight that simulating the flocculation process is not straightforward as this ideally implicates to integrate the impact of the bio-chemistry, the nature of the particles that are aggregating each other (ex: aggregates instead of primary particles, humic substance …) and the level of degradation of the aggregates.

Fig. 1 .
Fig. 1.Grain size distribution for the bottom sediments collected in the (a) riverbed and (b) the riverbanks of the Orne River.The sediment size distribution is divided in 10 classes according to the development operated on SYSIPHE.

Fig. 2 .
Fig. 2. The Orne River bottom sediment relative density measured for each of the 10 grain size classes.

Fig. 4 .
Fig. 4. Simulated and observed (estimated based on turbidity) SSC time series at the downstream boundary for the three model configurations.

Fig. 5 :
Fig. 5: Initial bathymetry of the entire river section before the simulations (a) and mapping of the three model simulations for the final simulated median grain size (b) a and the final bathymetry evolution after the flood event (c) in selected part of the studied river section.

Fig. 7 .
Fig. 7. Cross comparison of simulated bathymetry evolutions between the reference model configuration (10CL) and the two other configurations (2CL and 10CLF).The colours are corresponding to the 10CL's bathymetry evolution: the yellow is the deposition; the grey is the stable bathymetry and the blue is the erosion.