Temporal variability of partially-contaminated sediments in a strongly regulated reservoir of the upper Rhine River

The upper Rhine River is a highly harnessed and regulated river. EDF (a French electricity company) is in charge of eight dams on the upper Rhine River for producing hydro-electricity. In order to increase the safety and the competitiveness of the installations, but also to reduce their environmental impact, the sediment dynamics in these reservoirs has become a key factor to control and predict. In this study, we focused on the Marckolsheim reservoir, which is located 50 kilometers upstream the city of Strasbourg. Since its construction in 1961, this reservoir has been filled continuously with cohesive sediments, partially contaminated. Two field campaigns were performed in 2015 and 2016 under two different discharge conditions, with the objectives of quantifying the complex velocity fields on this site. The numerical codes TELEMAC-2D and SISYPHE were used to simulate in 2D the hydrodynamic and the suspended sediment transport of the reservoir. A ten kilometers long model was built and calibrated with the measured data of the 2015 and 2016 field campaigns, but also with measurements of sediment parameters that have been done separately. The originality of this model consists in an explicit 3D representation of the dam gates. An algorithm was implemented in TELEMAC in order to adapt the gates position at each time step, in conformity with the real regulation rules followed by the dam operator. By using upstream measured data of discharge and suspended sediment concentration, a four months period was simulated. The comparison of the simulated results with bathymetric surveys shows good agreements if specific properties of sediments related to settling processes are taken into account. Finally, the dynamics of the contaminated sediments was simulated. A 3D spatial distribution of the contaminated sediments in the reservoir was defined at the initial state by using in situ measurements. The fully coupled hydraulic-sediment-pollutant simulation performed over a single flood event gives first interesting highlights on the resuspension conditions of the contaminated sediments.


Introduction
The upper Rhine River is a highly harnessed and regulated river. Its main channel is navigable and its water is used for agriculture, drinking water supply and electricity production. EDF (a French electricity company) is in charge of eight dams on the upper Rhine River for producing hydro-electricity. In order to increase the safety and the competitiveness of the installations, but also to reduce their environmental impact, the sediment dynamics in these reservoirs has become a key factor to control and predict.
In this study, we focused on the Marckolsheim reservoir, which is located 50 kilometers upstream the city of Strasbourg. Since its construction in 1961, this reservoir has been filled continuously with cohesive sediments, partially contaminated. Figure 1 shows a general overview of the installation. To keep the water level suitable for navigation and to optimize electricity production, the dam is regulated with a high-frequency repositioning of its gates. Under usual discharge conditions, the dam gates are closed and the main part of the Rhine discharge is transferred to the hydro power plant along the navigation channel. When the Rhine discharge is higher than the maximum admissible discharge of the power plant, the dam gates are open in order to transfer the additional discharge downstream the dam. This regulation rules, combined with the bifurcation configuration of the channel, lead to a complex and unsteady behavior of the hydrodynamics in the reservoir. Additionally, the high temporal variations of the suspended sediment supply makes the deposition in the reservoir even more difficult to predict.
In order to quantify the hydrodynamics of the reservoir under different discharge conditions, two field campaigns were performed in 2015 (usual discharge condition) and 2016 (during a flood event). During these campaigns, aDcp measurements were carried out in order to evaluate the flow velocities along seven cross sections. These field campaigns are presented in details in section 2.
The numerical codes TELEMAC-2D and SISYPHE (www.opentelemac.org) were used to simulate in 2D the hydrodynamics and the suspended sediment transport on this site. A ten kilometers long model was built and calibrated using the measured data from the 2015 and 2016 field campaigns, but also with measurements of sediment parameters that have been done separately (density, critical erosion shear stress). The numerical results are presented in section 3.
Finally, the dynamics of the contaminated sediments was simulated. A 3D spatial distribution of the contaminated sediments in the reservoir was implemented in the SISYPHE model, inspired by in situ measurements. The 2016 flood event was simulated in order to evaluate the resuspension conditions of the contaminated sediments. This work is detailed in section 4.

Field campaigns
In order to quantify the hydrodynamics of the reservoir under different Rhine discharge conditions, two field campaigns were performed in 2015 and 2016. During these campaigns, aDcp measurements were done along seven cross sections from upstream to downstream the bifurcation. During the 2015 field campaign, the Rhine discharge was equal to 1200 m 3 /s and the dam gates were closed (a compensation discharge of 15 m 3 /s is transferred to the downstream part). This campaign will be designated UD (Usual Discharge) for thereafter. For the 2016 field campaign, measurements were done during a single flood event. The maximum discharge value for this flood event was 3000 m 3 /s, but measurements were done at the end of the decreasing part of the hydrogram, with an almost constant value of 2000 m 3 /s. For this discharge value, approximately 300 m 3 /s is transferred downstream the dam by opening the gates. This campaign will be designated HD (High Discharge) for the rest of the document.
The velocity measurements were performed with an aDcp profiler Rio Grande, developed by the company Teledyne. For the two field campaigns, the acquisition frequency was fixed to 600 kHz. The size of the vertical acquisition cell was respectively fixed to 10 cm for the UD campaign and to 50 cm for the HD campaign. For each of the seven cross-sections, five measurement crossings were performed. The final results were computed with the software VMT (developed by USGS) by calculating the average of the five measurements, projected on the target profile. Figure 2 shows the depth integrated velocities measured along profiles T3 to T7 (profiles T1 and T2 were far away upstream), for the UD (green light) and the HD (red) campaigns. Figure 2 shows that for most of the profiles, the measured velocity magnitudes are higher for the HD campaign than for the UD campaign, which is in line with the difference of discharge level. The profiles T1 to T3, T6 and T7 show similar horizontal velocity profiles, with velocity magnitudes almost homogeneous along the cross section. The profile T5 highlights a recirculation shape of the flow in the reservoir, which was accentuated during the HD campaign. For the T4 profile, placed at the bifurcation, the shape of the velocity profile is different for each campaign. The spatial distribution of velocities measured during the UD campaign shows that the main part of the flow was oriented toward the navigable channel, whereas a significant part of the velocities was oriented toward the dam during the HD campaign. Figure 3 shows the vertical and transversal distribution of the velocity magnitude along the profile T4 for the UD and HD campaigns. We can see on this figure that for the UD campaign (up), the velocities measured in the navigable channel (left part of the cross section) were up to 1 m/s, whereas very small velocities (smaller than 0.2 ms/) were measured in the reservoir (right part of the cross section). Velocities slightly higher than 0.2 m/s were measured along the left side of the reservoir, highlighting the initialization of a large and slow recirculation. For the HD campaign (down), velocities were much higher in the navigable channel (up to 1.6 m/s) than during the UD campaign, and the velocities measured in the left side of the reservoir were significant, up to 0.6 m/s. This highlights a fully developed recirculation in the reservoir, accelerated by the opening of the dam gates in flood conditions.

Hydraulic and suspended sediment numerical modeling
The numerical codes TELEMAC-2D and SISYPHE were used to simulate in 2D the hydrodynamics (Shallow Water Equations) and the suspended sediment transport (Advection- Dispersion Equation). A ten kilometers long model was built and calibrated with the measured data of the UD and HD field campaigns. The originality of this model consists in an explicit 3D representation of the dam gates. An algorithm was implemented in TELEMAC-2D in order to adapt the vertical gates position at each time step, in conformity with the real regulation rules followed by the dam operator. This numerical procedure allowed to simulate the right water level in the reservoir, which is maintained constant by the operator. Figure 4 shows a 3D representation from downstream of the model geometry and mesh, including the numerical dam gates. In order to calibrate the numerical model, several zones of constant Strickler coefficients have been defined in the reservoir (2 zones) and the navigable channel (2 zones). Figure 5 shows plan-views comparing measured and simulated velocity fields for the two field campaigns. The simulated velocity fields show very good agreements in shape and magnitude with the measured data, and highlight the recirculation patterns in the reservoir for the two different discharge conditions. Figure 6 shows a comparison between the magnitudes of the measured and simulated velocities along the T4 profile for the HD campaign. The simulated velocities fit well with the measured ones, validating both calibration and numerical dam gates procedure. For the sediment transport modeling, the suspension module of SISYPHE was used. In this module the suspension sediment transport is simulated with an advection-dispersion equation and two classical source terms for erosion (Partheniades law) and deposition (Krone law) of cohesive sediments. The parameters of these laws were defined with measurements performed in laboratory tests.
Mass concentration measurements and erosion tests were performed by Westrich in 2010 [2] on several cores sampled in the sediment deposits of the reservoir. According to these measurements, two layers of deposited sediments were defined in the numerical model of the reservoir: a first layer with a homogeneous thickness of 20 cm was defined, with a critical erosion shear stress of 1 Pa and a mass concentration of 1000 kg/m 3 ; a second layer starting from the bottom of the first layer and ending at the initial altitude of the reservoir was defined with a critical erosion shear stress of 3 Pa and a mass concentration of 1500 kg/m 3 . SCAF measurements (Wendling et al., 2015 [1]) were done with resuspended sediments in order to calibrate the settling velocities value of the Krone law. Highest values of the measured spectra of the settling velocities (1 mm/s) were used. This value corresponds to an aggregated state of the primary particles which is realistic regarding the resuspension processes that may occur in the upstream reservoirs during flood events. The critical shear stress for deposition was fixed to 0.2 Pa.
By using upstream measured data of discharge and suspended sediment concentration as boundary conditions, a four months period was simulated (starting in mid-May 2016 and ending in mid-September). During the first 60 days of this period, several successive flood events occurred on the Rhin River, with discharge values measured between 2000 and 3500 m 3 /s. So during this first period of the hydrograph, the dam gates were open continuously, providing velocities conditions in the reservoir favorable to erosion. During the next 60 days of the hydrograph, discharge values were most of the time smaller than 2000 m 3 /s, leading to deposition conditions in the reservoir. As the time step was equal to 1 second, the clusters of EDF R&D were used to perform this long-term simulation in 7 days of CPU time. The comparison of the numerical results with bathymetric surveys are presented in Figure 7. In this figure, the real bed evolution (left) is computed from two bathymetric surveys performed before and after the simulated period. The simulated bed evolution is computed at the last time of the simulation.

Contaminated sediments modeling
Starting from this hydraulic and suspended sediment transport model, the dynamics of the contaminated sediments was simulated. In 2010, a core sampling campaign has been performed, and measurements of concentration have been done on 32 cores for several chemical agents. In this study, we focused on one example of chemical agent in order to built the modeling methodology. A 3D spatial distribution of the pollutant concentration in the reservoir was calculated with a kriging algorithm. The dimensionless results are presented in Figure 8. The SISYPHE code was used and modified for allowing the simulation of two cohesive sediment classes: one class for contaminated sediments (class 1), one class for the noncontaminated ones (class 2). With this approach, both sediment classes have the same physi-cal properties and are defined by a concentration threshold of pollutant concentration. Seven threshold values have been applied to the 3D distribution of pollutant concentration, from 5% to 50% of the maximum observed concentration. The obtained 3D spatial distributions of the two sediment classes were integrated vertically to be represented in the two sediment layers of the SISYPHE model. Finally, the hydrogram of the HD campaign flood was simulated.  Figure 9 shows the simulated concentrations of the class 2 at the outlet of the dam during the flood event, for the seven considered thresholds. Over the flood period, we can see that the effect of threshold is significant for two time periods, at the beginning and the end of the event. During this two periods, these differences can represent 25% of the concentration values, leading to a similar or higher ratio in terms of sediment fluxes.

Conclusion and perspectives
In this study, we focused on the Marckolsheim reservoir. Two field campaigns were performed in 2015 and 2016 under two different discharge conditions, higlighting the corresponding velocity fileds in the reservoir. The numerical codes TELEMAC-2D and SISYPHE were used to simulate in 2D the hydrodynamics and the suspended sediment transport of the reservoir. A ten kilometers long model was built and calibrated with the measured data of the 2015 and 2016 field campaigns, but also with measurements of sediment parameters that have been done separately. The calibrated model fitted very well with the measured velocity fields of the two filed campaigns. By using upstream measured data of discharge and suspended sediment concentration, a four months period was simulated. The comparison of the simulated results with bathymetric surveys shows good agreements if specific properties of sediments related to settling processes are taken into account. Finally, the dynamics of the contaminated sediments was simulated. A 3D spatial distribution of the contaminated sediments in the reservoir was defined at the initial state by using in situ measurements. The fully coupled hydraulic-sediment-pollutant simulation performed over a single flood event gives first interesting highlights on the resuspension conditions of the contaminated sediments, depending on a threshold pollution value.