Modelling suspended-sediment propagation and related heavy metal contamination in floodplains : a parameter sensitivity analysis

Fine sediments represent an important vector of pollutant diffusion in rivers. When deposited in floodplains and riverbeds, they can be responsible for soil pollution. In this context, this paper proposes a modelling exercise aimed at predicting transport and diffusion of fine sediments and dissolved pollutants. The model is based upon the Telemac hydro-informatic system (dynamical coupling Telemac-2DSysiphe). As empirical and semiempirical parameters need to be calibrated for such a modelling exercise, a sensitivity analysis is proposed. An innovative point in this study is the assessment of the usefulness of dissolved trace metal contamination information for model calibration. Moreover, for supporting the modelling exercise, an extensive database was set up during two flood events. It includes water surface elevation records, discharge measurements and geochemistry data such as time series of dissolved/particulate contaminants and suspended-sediment concentrations. The most sensitive parameters were found to be the hydraulic friction coefficients and the sediment particle settling velocity in water. It was also found that model calibration did not benefit from dissolved trace metal contamination information. Using the two monitored hydrological events as calibration and validation, it was found that the model is able to satisfyingly predict suspended sediment and dissolve pollutant transport in the river channel. In addition, a qualitative comparison between simulated sediment deposition in the floodplain and a soil contamination map shows that the preferential zones for deposition identified by the model are realistic.


Introduction
Recent years have seen a growing awareness of the central role that fine-sediment loads play in transport and diffusion of pollutants by rivers and streams (Walling, 2005).Suspended sediment can potentially carry important amounts of nutrients and contaminants, such as trace metals of which some are recognized as potentially harmful elements (PHEs).These threaten water quality in rivers and wetlands and soil quality in floodplains (Carter et al., 2006;Hissler and Probst, 2006).Contemporary data on sediment loads of rivers provide clear evidence of significant recent changes in sediment fluxes of several rivers in response to human activities (Walling, 2006).Although fine-sediment deposition in floodplains is not necessarily responsible for important topographical evolution, it can play a central role from a contamination point of view (Stewart et al., 1998;Benjankar and Yager, 2012).
Currently, many studies focusing on sediment transport modelling deal with marine and estuarine areas (e.g.Le Normant, 2000;Nguyen et al., 2009).Some studies evaluate sediment transport on basin scales and often evaluate yearly sediment fluxes using hydrologic and simplified hydraulic models (e.g.van Griensven et al., 2013).Some more theoretical studies develop and improve numerical models on the basis of physical model experiments (e.g.Belleudy, 2000Belleudy, , 2001;;Bui and Rutschmann, 2010).As a matter of fact, sediment transport modelling in small rivers on a reach/floodplain scale is a rather new research field (Simpson and Castelltort, 2006).Among the most recent studies on sediment transport modelling in river systems, one can cite Villaret et al. (2013) and González-Sanchis et al. (2014).While the study of Villaret et al. (2013) concentrates mostly on the evaluation of a specific modelling system in various test cases, González-Sanchis et al. (2014) present a study with objectives closer to the one we aim at.
In this paper, we aim at simulating sediment transport on the floodplain scale and the flood event scale in order to predict sediment spreading on alluvial soils.This simulation will help in the estimation of the potential pollution of soil due to the transport of PHEs by suspended sediments.As argued by Benjankar and Yager (2012), only a few studies have focused on fine-sediment deposition in floodplains.Moreover, Hardy et al. (2000) explain that, in this context, it is necessary to make use of a model able to consider advection and diffusion processes and to carry out unsteady simulation (physically variable state varying in time).
Numerical models are used more and more used by water resource planners, water quality managers, engineers and scientists (Singh and Woolhiser, 2002;Simpson and Castelltort, 2006).Hydrodynamic and sediment transport models, built up using in situ measurements, arguably represent a useful tool for predicting natural and man-induced environmental impacts on sediment dynamics, especially due to the complexity of physical processes involved in sediment transport (Belleudy, 2000).They are based on an approximate representation of complex natural systems.Therefore, the evaluation of such models with respect to their ability to reproduce multiple processes of a real system is still problematic.In this context, Matgen et al. (2007), Pappenberger et al. (2007), Schumann et al. (2007) and Hostache et al. (2009) demonstrated that the calibration of such models is not straightforward and needs particular attention as the data sets used in the calibration process determine the optimal parameter set.Hostache et al. (2009) came to the conclusion that data sets other than conventional stream gauge measurements are particularly useful for model calibration because they constrain model parameters better and improve the identifiability of the model parameters.According to Beven (2006), equifinality occurs when various parameter sets yield similar model performance with respect to a given observation due, for example, to compensating effects or unsuitable process understanding.
In this context, this article describes a modelling exercise using a rather unique field data set, which includes not only water surface elevation records but also geochemistry data such as temporal variations of contaminants and suspendedsediment concentrations.This extensive data set offers new opportunities for evaluating and analysing in an objective way the performance of a model, which is applied to a small river system, with respect to different physical processes.
The aim of the modelling exercise is twofold.The first objective is to set up a model capable of accurately predicting flood wave dynamics, dissolved-contaminant dispersion (tracers) and sediment propagation during flood events.The second objective is to identify the most sensitive parameters of the model using different kinds of observations and to evaluate how tracers enable a better calibration of the model.More specifically we aim at determining whether tracer concentration information can be used to calibrate a hydrodynamic model in a context similar to the study of Fenicia et al. (2010), who used tracer data to calibrate a hydrologic model.The calibration of such models is far from trivial especially because there are, in addition to physical parameters, sitespecific, empirical and semiempirical parameters that need to be calibrated (Hardy et al., 2000).In this context, we propose to carry out a sensitivity analysis of the model.
The paper is organized in three parts.First we detail the study area and the available observation data set.Next, we present the model and the method adopted for the sensitivity analysis.Finally, we show the results of the study and discuss its main outcomes.

Study area
Due to important urban and industrial developments since 1900, the southern part of the Grand Duchy of Luxembourg suffers from substantial PHE contamination from various origins (Hissler et al., 2008).The upper Alzette River, which drains this historical steel basin in Luxembourg presents all the characteristics of a good test site for evaluating small river system and alluvial plain contamination.The study area of about 2.2 km 2 is located at the outlet of a 290 km 2 river basin (see Fig. 1).In this part of the basin, the Alzette riverbed has a mean slope of around 0.1 %, a mean depth of around 4 m and a mean top width of around 12 m.It has two main tributaries (see Fig. 1), namely the Bibeschbach on its left side and the Crauthemerbach on its right side.
As a test case, we propose to focus on two flood events that occurred in January and December 2011.Whereas the January flood event (return period of 8 years) was responsible for a rather large floodplain inundation, flows mostly remained in-channel with only sparse overbank flow during the December one (return period of 1 year).

Topographic data
The set up of the model requires accurate information about the terrain and the riverbed topography.The terrain elevation data have been derived from a lidar digital elevation model (DEM) representing the ground surface elevation with a pixel size of 2 m and a theoretical accuracy on the elevation of ±15 cm.The riverbed elevation data for River Alzette are available as ground-surveyed cross sections (26 in the study area, see Fig. 1) with a theoretical elevation accuracy of a few centimetres.Between these cross sections, the Alzette riverbed elevation has been linearly interpolated in order to draw a continuous riverbed across the study area.The riverbed topography of the two tributaries have been interpolated between two ground-surveyed cross sections (one at the upstream end and one at the downstream end of each tributary).This interpolation has been performed  along digitized riverbed lines.The banks, the roads and the other hydraulic singularities present in the study area have also been digitized.The mesh representing the model domain was drawn from the river, banks and hydraulic singularity lines and the elevation of each node was derived from the DEM in the floodplain and interpolated between observed cross sections in the riverbeds.The resulting unstructured triangular mesh contains 25 086 nodes.It is refined in the riverbed and close to the riverbanks and the hydraulic structures where the average distance between nodes is around 1 to 3 m.In the flat parts of the floodplain, the distance between nodes can reach 15 m.The model simulation time step has been set to 0.1 s in order to respect the Courant-Friedrichs-Lewy condition.

Hydrometric data
The study area is equipped with four stream gauges (see Fig. 1) recording the water surface elevation every 15 min.During the January 2011 flood event, the stream gauge in the upstream part of the Crauthemerbach tributary was not operational.Considering that this event is really interesting due to substantial floodplain inundation, the discharge hydrograph at this section was estimated from observations of other flood events.In addition to the water surface elevation records, discharge measurements during flood events have been carried out at the following stream gauges: upstream Alzette (BCAl1), upstream Crauthemer-  bach (BCCr), upstream Bibeschbach (BCBi) and downstream Alzette (BCAl2) (Fig. 1).These measurements allowed us to estimate rating curves for each gauge (i.e. the discharge-water-surface-elevation relationship at BCAl1, BCCr, BCBi and BCAl2.Figure 2a and b present the discharge hydrographs recorded during the January and December 2011 flood events respectively.

Geochemistry data
Sites BCAl1, BCCr and BCAl2 were instrumented for multitracer monitoring of the two flood events of January and December 2011.The multitracing approach we proposed for the calibration of the model includes various physico-chemical parameters.Three distinct trace metals, namely gadolinium (Gd), lead (Pb) and zinc (Zn), were chosen to characterize the temporal evolution of the dissolved phase in the water column and were used for the calibration of the hydrodynamic model.These three trace metals were chosen as they are considered as PHE that potentially have a strong impact on water and soil quality.Moreover, they behave differently in river systems.Pb and Zn cannot be strictly considered as conservative tracers because their concentrations in the dissolved phase are known to be impacted by exchanges with the particulate phases due to redox conditions, temperature and biological activity.On the contrary, in such a contaminated river system, Gd remains in the dissolved fraction of the water.It can be considered as a conservative tracer for hydrological purposes (Möller et al., 2000).
In addition, the suspended-sediment concentration that characterizes the evolution of the particulate phase in the water column during the flood events was used to calibrate the sediment transport model.After the water sampling, using ISCO © autosamplers at each of the three sites, the filtration of the water (filters with 0.45 mm poresize) allowed separation of the dissolved (including colloidal phases) and particulate phases of each sample.The dissolved trace metal concentrations were determined by inductively coupled plasma mass spectrometry (ICP-MS).The errors due to sample preparation and analysis were negligible in comparison to the evolution of tracer concentrations between the different samples collected during the flood events.Figure 2c and d present the dissolved Gd concentrations recorded during the January and December 2011 flood events respectively and Fig. 2e and f present the suspended-sediment concentrations recorded during the same events.Additionally, the spatial distribution of the trace metal contamination at the surface of the alluvial soils was estimated in order to evaluate the fine-sediment deposition maps obtained with the Telemac-2D simulations.Zn is recognized as a tracer of the anthropogenic contamination that comes from the river to the soils of the Alzette River floodplain area (Horckmans et al., 2005).The distribution of this contaminant at the soil surface of the study area indicates the dispersion of contaminated sediment during floodings.It represents integrated information of all the contamination events within the alluvial area due to the river sediment deposition.Consequently, the most heavily contaminated areas of the floodplain may correspond to the areas most impacted by the sediment deposition.Following a regular grid of 100 m spacing that covers a large part of the study area, 100 soil surface samples (0-5 cm depth) were collected in January 2006.The Zn concentration was determined by ICP-MS after digestion of the soil samples using a HCl / HNO 3 mixture.

Material and methods
This section presents the modelling approach that has been adopted.
Following the recommendations proposed by Hardy et al. (2000) and Benjankar and Yager (2012), we made use of a modelling system able to carry out unsteady simulations, to take account of advection and diffusion processes, feedback processes between topography and hydrodynamic variables, and to integrate input suspended-sediment concentration as boundary conditions.According to these requirements, the model that has been set up is based on the open source Telemac hydro-informatic system (release 6.2) (Hervouet and Bates, 2000).The latter allows for dynamically coupling a 2D-shallow water hydrodynamic model -Telemac-2D -and a sediment transport model -Sisyphe.In this study, we made use of the parallel version of the code, using 16 parallel nodes.

The hydrodynamic model
Telemac-2D is a two-dimensional shallow-water hydrodynamic model.It solves the de Saint-Venant equations -also called the shallow-water equations (conservation of mass and momentum, see Eqs. 1-3) -and predicts, among other hydraulic variables, water depth and fluid velocity at every node of a triangular mesh representing the model domain.The following paragraphs present the main mathematical features of Telemac-2D that are relevant for our study.For more details, readers are referred to the Telemac-2D user manual (Lang, 2010).
In Eqs. ( 1)-( 3) h is the water depth, η the water surface elevation, t the time, w the fluid velocity (vector of component u and v along spatial dimensions x and y), ν t the momentum diffusion coefficient, and F x and F y represent momentum source/sink terms (e.g.friction force on the river bottom, and wind force).Turbulence could play an important role in suspended-sediment transport processes.To model the turbulence and take account of its vertical component that is neglected due to averaging over the vertical dimension, we make use of the k-epsilon formulation (Launder and Sharma, 1974).
Telemac-2D can also simulate current entrainment as well as the diffusion of tracers (e.g.dissolved PHE concentrations).For this purpose, the model solves the conservation equation for the tracer T (see Eq. 4).
In Eq. ( 4), T represents the tracer quantity (e.g.concentration in g l −1 ), S T a source/sink term and D T the tracer diffusivity coefficient.In this equation, the source can be contaminant inputs or outputs.In this context, it has to be noted that in the event of nonconservative tracers, it is possible to program a decay function.
In this study, we used the dissolved gadolinium (Gd), dissolved lead (Pb) and dissolved zinc (Zn) concentrations as tracers.
Discharge hydrographs and tracer concentration temporal evolutions are imposed at BCAl1, BCBi and BCCr (Fig. 1) as upstream boundary conditions.As downstream boundary conditions (see point BCAl2, Fig. 1), a water surface elevation hydrograph is imposed and a free exit of tracer concentrations is used.
As an initial condition we used the steady state (water surface elevation, velocity and tracer spatial distribution) reached after a 6 h-simulation using as boundary conditions the same as those at the start of the dynamic run.

The sediment transport model
Sediment transport is simulated using Sisyphe, part of the Telemac hydro-informatic system.Sisyphe does not allow directly for the simulation of contaminant transport in the particulate phase, as the chemical processes related to the adsorption/desorption of contaminants on/from suspendedsediment particles are not included.In this first study, we assume that the modelling of suspended-sediment transport might help estimate contaminant concentrations of the particulate phase in the water column.
Telemac-2D and Sisyphe were dynamically coupled with an identical simulation time step (0.1 s).After each hydrodynamic simulation time step, Telemac-2D sends its hydrodynamic state variables to Sisyphe, which carries out the sediment transport simulation and outputs the sedimenttransport-related state variables back to Telemac-2D.The dynamical coupling offered by the Telemac hydro-informatic system allows the effect of bed evolution on water propagation to be taken into account.
Sisyphe simulates erosion, deposition, bed load and suspended-sediment transport in the water column (Villaret, 2010).It is based on established semiempirical equations for sediment transport and decomposes the underlying processes into bed load and suspended load.Considering that the estimation of contaminant deposition in the floodplain is of primary interest, the target processes are the suspendedsediment transport and deposition.Sisyphe carries out the bed evolution computation using the Exner equation (Exner, 1920(Exner, , 1925)).
The transport of suspended sediment is computed using Eq. ( 5): where C is the depth-averaged suspended-sediment concentration, E and D the erosion and deposition rates at the bed load/suspended load layer interface, s the coefficient of dispersion and U conv and V conv the convection velocity components.
For the computation of erosion/deposition rates (see Eq. 5), numerous empirical equations are implemented in Sisyphe.In particular, equations used for erosion and deposition calculation depend on the cohesive properties of riverbed material.The riverbed material can then be considered as non-cohesive or cohesive.A recent development of Sisyphe offers the capability to simulate cohesive/non-cohesive sediment mixtures.The erosion rate computation follows the method proposed by Waeles (2005) (see Eq. 6).
If the bed is composed of less than 30 % of mud (cohesive sediment with grain size lower than 63 µm as defined by Villaret, 2010) the erosion rate computation is conditioned by the suspended-sediment equilibrium concentration (see Eq. 6).By default in Sisyphe, the equilibrium concentration is estimated using the formula proposed by Zyserman and Fredsoe (1994) based on the Shields parameter (Eq.6).
If the bed is composed of more than 50 % of mud, then the erosion/deposition rate computation is based on the Krone and Partheniades formulation (Partheniades, 1965), conditioned by the nominal settling velocity of mud particles in (still) water and on two critical shear velocities, namely for erosion and deposition (Eqs.6 and 7).If the mud ratio is between 30 and 50 %, a linear interpolation is performed between the two above-mentioned formulations (Eq.6).
In Eqs. ( 6) and ( 7), indices 1 and 2 refer respectively to sand and mud, M is the Partheniades constant, W s the settling velocity of particles in water, P 2 the percentage of mud in the riverbed, θ c the Shields parameter, C the suspendedsediment concentration, C eq the suspended-sediment equilibrium concentration, C Z ref the suspended-sediment reference concentration (estimated from C based on the Rouse profile), u * the local shear velocity and u * d and u * e the critical shear velocities respectively for deposition and erosion.The Shields parameter corresponds to the adimensional shear stress value beyond which erosion starts occurring.Using this formulation, erosion is computed as a global flux of sand and mud together, whereas deposition is computed separately for mud and sand particles.This reflects real behaviour in the sense that the ability of the riverbed material to resist erosion is a global property of riverbed material, whereas the deposition of individual particles depends on individual particle characteristics.Moreover, looking at Eqs. ( 6) and ( 7), it has to be noted that on the one hand parameters θ c and W s 1 and on the other hand parameters u * d and W s 2 can potentially have compensating effects on erosion and deposition.
The concept of suspended-sediment equilibrium concentration is used in Eq. ( 6) in order to characterize a tendency to erode or deposit sediment.Indeed, a simulated suspendedsediment concentration larger than the equilibrium concentration implies a tendency to deposit, whereas a simulated suspended-sediment concentration smaller than the equilibrium concentration implies a tendency to erode the riverbed.It has to be noted that the equilibrium concentration in Eq. ( 6) is a global value computed for all sediment classes.The equilibrium concentration is distributed in each sediment class proportionally to the sand grading in the riverbed.Therefore, the sand grading curve can have a non-negligible effect on erosion and deposition processes.
In our study, suspended-sediment measurements only considered particles with diameters lower than 63 µm.Particles of such diameters are defined as cohesive in Sisyphe according to Villaret (2010).This was supported by observation as the riverbed material appeared cohesive during field campaigns.Consequently, we took advantage of the recent development of Sisyphe that allows for modelling graded cohesive/non-cohesive sediment mixtures.Using this approach two classes of sediment can be defined, the first one corresponding to non-cohesive (sand) and the second to cohesive (mud) sediment.

The sensitivity analysis and the evaluation of the model results
The aim of the sensitivity analysis is twofold.On the one hand, we want to analyse the identifiability of model parameters.On the other hand, we want to evaluate which observation data set is the most powerful for the calibration.
In particular, we want to see if the concentration of selected chemical tracers can be used as an alternative to more traditional hydrometric data for the calibration of a hydrodynamic model.First, tests of an all-at-once evaluation of model parameters governing the hydrodynamic and the sediment transport processes show that it is quite difficult to interpret the results due to parameters having compensating effects on each other.It has to be noted that this problem has already been observed and discussed in many modelling exercises and is often referred to as the equifinality of model parameters (Beven, 2000).
As a consequence, we decided to carry out a two-step sensitivity analysis considering that the data sets used for the calibration of the hydrodynamic model (discharge, water surface elevation and tracer concentrations) and for the calibration of the sediment transport model (sediment concentrations) are independent.First, we focused on hydrodynamic parameters and only afterwards did we consider the sediment transport parameters.In this context, a first analysis is performed using only Telemac-2D with tracer transport simulation.Next, simulations of sediment transport are carried out and analysed considering that the hydrodynamic para-meters have already been calibrated during the first analysis.The sensitivity analysis is based on a random sampling of model parameter sets from an a priori defined range of physically plausible parameter values.To do so, each parameter value is randomly sampled from a uniform distribution having a range corresponding to the minimum and the maximum plausible value of the considered parameter.Once the parameter sets have been generated, model simulations are carried out for each generated parameter set.Finally, each model result is compared to a set of observations and subsequently, model skill scores are computed for each parameter set.
For each generated parameter set, the results are compared against observations using a slightly modified version of the Nash-Sutcliffe Efficiency (NSE) (Nash and Sutcliffe, 1970).
The NSE is a model skill score expressing the percentage of variance of the observations explained by the model results.An NSE value of 1 means a perfect fit between model results and observations.In the modified version of the NSE we used, called NSE* hereafter, we took account of observation uncertainty, meaning that the fit between observations and model results is considered perfect if model results are included in intervals representing the observation uncertainty (see Eq. 8).
In Eq. ( 8), S is the simulated variable, O the observed variable, O the time average of the observed variable and u O the uncertainty of the observation.The uncertainty is set to ±2 cm for water surface elevation, and to ±10 % of the observed values for discharge and tracer concentrations.These uncertainty bounds have been estimated from past experiments where the same measurement was repeated many times subsequently or at slightly different locations.NSE* has the main advantage of exhibiting comparable values for different kinds of observations.Moreover, model results falling within the uncertainty bounds of the observation are considered equal as further distinction is not possible with the observations.
The selected hydrodynamic parameters are the Strickler friction coefficient and tracer diffusivity coefficient.The tracer diffusivity coefficient is assumed to be spatially uniform.In a previous study, Werner et al. (2005) and more recently Hostache et al. (2010) argued that it is quite difficult, even with advanced assimilation techniques (Lai and Monnier, 2009) and spatially distributed water level data, to calibrate distributed friction coefficients in a river channel.As a consequence we decided to consider two different values of the friction coefficients: one for the river streams and one for the floodplain.150 plausible parameter sets have been randomly generated within uniform distributions ranging from 25 to 45 m 1/3 s −1 for the channel Strickler coefficient, from 15 to 35 m 1/3 s −1 for the floodplain Strickler coefficient and from 10 −7 to 10 −4 m 2 s −1 for tracer diffusivity coefficients.One could argue that the number of generated parameter sets is rather limited but considering the duration of one simulation (approximately 1 day) it does not seem easily feasible to increase this number.Moreover, we believe that, albeit limited, the number of parameter sets might be sufficient for capturing parameter sensitivity.
The aim of the second step of the sensitivity analysis is to analyse the identifiability of model parameters using observed suspended-sediment concentration.The targeted parameters for the sensitivity analysis of the sediment transport model are the critical Shields parameter (θ c ) and the settling velocities for sand (W s 1 ) and mud (W s 2 ) particles (see Eq. 6).In addition to these three parameters, as it is difficult to estimate an average sand grading curve for the whole river reach and as we saw previously that the sand grading curve could play a role in the sediment erosion/deposition processes, we decided to consider the percentage of mud in the Alzette (pma) and Crauthemerbach (pmc) riverbeds as two additional model parameters.To carry out the sensitivity analysis, we adopted the approach previously presented for the hydrodynamic model.We first randomly generated 180 parameter sets within the following ranges: pma in [.01 .10], pmc in [.10 .30],W s 1 in [10 −6 10 −3 ] m s −1 , W s 2 in [10 −7 10 −4 ] m s −1 and θ c in [.01 .45].Again, one could argue that the number of generated parameter sets is rather limited but considering the duration of one simulation (more than 1 day) it does not seem easily feasible to increase this number.Moreover, we believe that, albeit limited, the number of parameter sets might be sufficient for capturing parameter sensitivity.Next, for each generated parameter set, we carried out dynamically coupled Telemac-2D-Sisyphe simulations.For each simulation, we compared simulated and observed suspended-sediment concentrations at BCAl2.For this evaluation, we made use of NSE*, defined in Eq. ( 8).The uncertainty for the suspended-sediment concentration has been assumed to be equal to ±10 % of the observed values.This uncertainty value was estimated from measurements performed at the same time for various river discharge conditions in locations close to each other (various depth in the water column and various points on a river cross section).It is worth keeping in mind here that the suspended-sediment concentration measurements considered mud particles.

Results and discussion
In the following sections, the results of the two-step sensitivity analysis proposed in Sect.3.3 are presented and discussed using as a test case the January 2011 flood event.The January 2011 flood event has been chosen for the sensitivity analysis because the flood level and the observed suspendedsediment concentrations were higher than during the December 2011 flood event.In addition, this section proposes to qualitatively evaluate model results by comparing the simulated sediment deposition map with point samples of floodplain soil contamination.Moreover, we propose to further evaluate the calibrated model on the basis of a simulation of the December 2011 flood event.

Hydrodynamic model sensitivity analysis
The NSE* has been computed for all available observations and all model simulations.Figure 3 shows the results of the evaluation.In these dotty plots, each dot corresponds to one model simulation (one parameter set).The values on the x axis are the parameter values (channel or floodplain Strickler coefficients).The values on the y axis correspond to the NSE* values calculated by comparing simulated and observed data (water surface elevation, discharge and Gd, Pb and Zn concentrations).Each colour reflects a different observation location (stream gauge location).A clear trend in the dotty plot shape indicates a high sensitivity of the model result with respect to the considered parameter and a good identifiability of the parameter.On the contrary, a rather flat dotty plot indicates that the model is less sensitive to the considered parameter and that this parameter most likely suffers from the equifinality issue.The range of the dotty plot provides additional information, showing the overall sensitivity of the model results to the considered parameter with respect to the observed variable.A large range shows that variations of the considered parameters are responsible for significant differences in model results whereas a small range indicates that the model results are similar for all parameter values tested.
The dotty plots obtained with respect to the tracer diffusivity coefficient exhibited no significant sensitivity and are not presented in this paper.This is rather surprising since we expected this parameter to have a significant effect on simulated tracer concentrations.This shows that the tracer concentrations are mainly governed by current entrainment instead of diffusion phenomena.The relatively sharp peak exhibited in Fig. 3a shows that the water surface elevations simulated by the model are mainly sensitive to Kc, the channel friction coefficient, whereas Fig. 3b indicates a limited sensitivity to Kf, the floodplain friction coefficient, when comparing simulated and observed water surface elevation.In Fig. 3c and d, when comparing observed and simulated discharge at the BCAl2 stream gauge location it appears that the model is rather sensitive to Kc and to a lesser extent to Kf.When comparing Fig. 3b and d, the sensitivity of the model to Kf is more visible for the evaluation of the discharge at BCAl2 (downstream Alzette) than of water surface elevation at BCAl1 (upstream Alzette).This is most likely due to the fact that BCAl2 is the downstream boundary condition.As a consequence, the  downstream hydrograph is a more integrative data set with respect to both Kc and Kf.Moreover, in Fig. 3a and c the best model performances are obtained for similar values of Kc, indicating good agreement of the model results with the two kinds of data.Figure 3e-j indicate, at first view, significant sensitivities of the model to Kc when comparing simulated and observed Gd, Pb and Zn tracer concentrations at the BCAl2 stream gauge.However, one can notice that these three plots are in contradiction since the best model performances are obtained for markedly different Kc values.From this point many questions with respect to the model results arise.To better address this issue, Fig. 4 shows simulated water surface elevations and concentrations of Gd, Pb and Zn.In this figure, each black line corresponds to the results provided by Telemac-2D using one parameter set and the red stars correspond to the observations.By comparing Fig. 4a with Fig. 4b-d, one can see that the behaviour of the model ensemble is markedly different for water surface elevation than for the chemical tracer concentrations.Indeed, whereas the spread of simulated water surface elevation is rather large and encompass the observed one, the range of the simulated tracer concentrations are limited and the observations frequently lie outside the ensemble of simulated concentrations.This shows   3e-j are then dominated by a few time steps for which some parameter sets provide better results, but no parameter set clearly outperforms any other for the whole simulation period.This hampers a straightforward interpretation of the trends observed in Fig. 3e-j.In Fig. 4b-d, it is reasonable to say that the temporal dynamics of tracer concentration are globally well captured by the ensemble of model results.This shows that the model, if calibrated using only water surface elevation and/or discharge, is able to predict tracer concentration with an acceptable accuracy.However, the calibration of the hydraulic model using only observations of chemical tracer concentrations is not a suitable option.In addition, the similar results obtained for Gd, Pb and Zn (see Fig. 3) do not provide any conclusion regarding the importance of considering strictly conservative tracers.
As a matter of fact, one of the main outcomes of the sensitivity analysis of the hydrodynamic model is that these three trace metal concentrations cannot substitute traditional calibration data sets (water surface elevation and discharge).This implies moreover that the model calibrated using water surface elevation and discharge hydrographs is able to predict tracer concentrations satisfactorily, whereas a hydrodynamic model calibrated using only observations of tracer concentrations would not necessarily be able to predict water surface elevation and discharge accurately.
To define the best parameter set that will be used for the sediment transport modelling, we decided to use the parameter set providing the best performance over the hydrometric calibration data set (water surface elevation and discharge).
To do so, we first rescaled the NSE* obtained for each observed data set between 0 and 1. Next, we computed an overall performance score as the average of the previously computed rescaled values of NSE*, and the parameter set providing the best overall performance score was considered as optimal.This method for computing the overall performance score has been chosen since it gives the same weight to each observation.The parameter set identified as optimal is the following: Kc = 32 m 1/3 s −1 , Kf = 23 m 1/3 s −1 and D T = 10 −6 m 2 s −1 .

Sediment transport model sensitivity analysis
The dotty plots resulting from the model evaluation are presented in Fig. 5.The critical Shields parameter (Fig. 5e) and the mud and sand ratio in the Alzette (Fig. 5a) and Crauthemerbach (Fig. 5b) riverbeds are of rather limited importance since similar model performance values have been obtained using markedly different values of these parameters (rather flat dotty plots).As a consequence, it is worth noting in this figure that the most sensitive parameter of the model is the settling velocity of mud particles (Fig. 5d).This result is in agreement with the remark by Villaret (2010)     of riverbed material, the settling velocity of mud rather than the particle diameter might be the governing parameter for sediment transport.Figure 5c exhibits some sensitivity of the model to the settling velocity of sand particles as the model evaluation clearly penalizes very small values for the sand particle settling velocity.This result makes sense as sand particles have larger diameter than mud particles and must therefore have higher settling velocity in water.
To select the optimal parameter set we consider the parameter set yielding the best performance (NSE*).The parameter set identified as optimal is the following: pma = 5 %, pmc = 19 %, W s 1 = 6.81 × 10 −5 m s −1 , W s 2 = 1.99 × 10 −5 m s −1 , and θ c = 0.02.Table 1 summarizes the performances reached by the calibrated model.Figure 6 shows the simulated (calibrated model) and observed suspended-sediment concentrations at BCAl2.In this figure, the simulated suspended-sediment concentration captures the temporal evolution of the observed suspendedsediment concentration rather well although the first peak is  slightly overestimated and the second one is inversely underestimated.This is likely due to erosion simulated as too important during the beginning of the flood event and, by contrast, deposition simulated as too important during the following time steps.However, the overall fit between observed and simulated sediment concentration is rather good, which indicates that the calibrated model is capable of predicting the downstream suspended-sediment concentration satisfactorily.

Qualitative evaluation of calibrated models on the basis of soil contamination
In order to further evaluate model results we qualitatively compared the deposition map simulated by the model with soil contamination maps derived on the basis of Fig. 7.The latter presents the map of fine-sediment deposition simulated by the calibrated model during the January 2011 flood event.
The black circles represent the observed Zn concentration in the surficial layer (0-5 cm depth) of the soils.The larger this circle the greater the soil contamination.As this contamination is mainly due to sediment deposition in the floodplain  during the successive flood events, the soil contamination sampling can be considered as a proxy for contaminatedsediment deposition.As a matter of fact, there should be some similarities between the sediment deposition map simulated by the model and the soil contamination.An important point to mention when comparing these two sources of information is that they do not represent the exact same thing.Whereas the observed contamination is the consequence of all overbanking flood events since the beginning of the contamination, the simulated deposition only represents the result of the January 2011 flood event.The qualitative comparison of the two sources of information is possible if the extent of the flooding is fairly stable and if there is a significant link between sediment deposition and soil contamination.
In Fig. 7, most of the areas with high deposition (in red and yellow) are overlapped by circles of rather large diameters, especially in the downstream areas (north-western part).By contrast, circles of rather limited diameter (locations less impacted by river sediment deposition) are located in areas where the model predicts limited deposition (blue) to no deposition (white).In addition, in the upstream part of the domain (western part), one can remark that some large circles are located where the model does not simulate any deposition.This can be easily explained by the fact that the soil contamination data set integrates information from many flood events, with some of higher magnitude than the January 2011 flood event.These are, therefore, responsible for larger flood extents and larger zones of deposition.The overall comparison also highlights the similarity of the dispersion between the geochemical information and the model simulation.Indeed, the spatial pattern of simulated sediment deposits is in rather good agreement with the spatial distribution of Zn contamination.
Consequently, a rather good consistency between the two sets of information is obtained despite some local mismatch.This indicates that the capability of the model to identify the main sediment deposition areas is satisfying.Moreover, this indicates at the same time that the simulated deposition map can be used to infer interesting information about the contamination dispersion in the floodplain area.The main drawback in this respect is that only a qualitative characterization of the pattern of contamination can be made as it is difficult to estimate the level of contamination the deposited sediment suffers from.

Evaluation of calibrated models on the basis of the December 2011 flood event
In order to further evaluate the calibrated model, we propose in this section to carry out a simulation for the December 2011 flood event.This second flood event is interesting because it is markedly different from the January 2011 flood event that has previously been used for the sensitivity analysis and the model calibration.First, the December 2011 event was a flood event of lower magnitude, and one key feature is that limited overbank flow occurred during this event.Moreover, the contribution of the river Crauthemerbach in terms of sediment fluxes was very limited compared to the January flood event.This last point is interesting as it allows for evaluating the Alzette response more specifically.The results of the model have been evaluated against observed data using the NSE* skill score.The obtained performances are listed in Table 1 in addition to the performances obtained during the January flood event.
The performances obtained are good especially with respect to the downstream discharge and the upstream water surface elevation with NSE* values higher than 0.9.The temporal evolution of tracer concentrations are predicted satisfactorily by the model, with values higher than the ones obtained during the calibration.In particular, the model performance with respect to concentrations of Gd and Zn is good.The performance in predicting suspended-sediment concentrations is slightly lower than during the calibration.Figure 8 presents simulated and observed suspended-sediment concentration at BCAl2.In this figure, one can see that the model tends to slightly underestimate the sediment concentration, but the results appear satisfying.This underestimation could be due to the simplified representation of the flow velocity field in two dimensions in Telemac-2D.As a matter of fact the shear velocity could be underestimated especially for lower discharge.As a consequence, the overall performances of the model for a flood event different from the one used for the calibration are satisfying.This indicates that the model can be used for predicting sediment transport/deposition and the associated dispersion of PHE contamination for many flood events.

Conclusions
This study aimed at predicting sediment and contaminant dispersion in a small river system.To do so, a model was set up and a sensitivity analysis was carried out.An innovative point in this study is the rather unique measurement database including water surface elevation, discharge, dissolved trace metal concentration and suspended-sediment concentration at two point locations.
The most sensitive parameters of the hydrodynamic model were the friction coefficients.It was found that tracer concentration was not an informative data set for calibrating the hydrodynamic model as the variability of the simulated tracer concentration with changing parameter values was low.Moreover, considering that tracer concentrations were satisfyingly predicted by the model, we inferred that the model calibrated using usual hydrometric data (water surface elevation and discharge) was able to predict trace metal dispersion satisfactorily, whereas the model calibrated using only tracer concentration was not suitable.The performances reached by the calibrated hydrodynamic model were satisfying.
The most sensitive parameters of the sediment transport model were found to be the settling velocity of sediment particles, divided into two diameter classes: mud (cohesive) and sand (non-cohesive).The mud particle settling velocity was identified as the most sensitive parameter, but this result is mainly due to the fact that the model was evaluated using observed suspended-sediment concentrations that were composed of mud particles only.The comparison of the simulated sediment deposition map with floodplain spatially distributed soil contamination information exhibited rather good agreement.This result, albeit qualitative, is quite encouraging as the zones where the model tends to deposit sediment are located on the most contaminated areas.
The evaluation of the coupled model using a second flood event also yielded satisfying performances.This is also encouraging as it provides evidence that the model is robust and can be used to predict various flood events accurately.

Fig. 2 :
Fig. 1: Presentation of the study area and extension of the model domain

Figure 1 .
Figure 1.Presentation of the study area and extension of the model domain.

Fig. 3 :
Fig. 3: Dotty plots of the hydrodynamic model sensitivity analysis: evaluation of simulated (a) water surface elevation (WSE) depending on channel Strickler coefficient (Kc), (b) water surface elevation depending on floodplain Strickler coefficient (Kf ), (c) discharge (Q) depending on Kc, (d) discharge depending on Kf , (e) dissolved gadolinium depending on Kc, (f) dissolved gadolinium depending on Kf , (g) dissolved lead depending on Kc, (h) dissolved lead depending on Kf , (i) dissolved zinc depending on Kc, (j ) dissolved zinc depending on Kf .

Figure 3 .
Figure 3. Dotty plots of the hydrodynamic model sensitivity analysis: evaluation of simulated (a) water surface elevation (WSE) depending on channel Strickler coefficient (Kc), (b) water surface elevation depending on floodplain Strickler coefficient (Kf ), (c) discharge (Q) depending on Kc, (d) discharge depending on Kf , (e) dissolved gadolinium depending on Kc, (f) dissolved gadolinium depending on Kf , (g) dissolved lead depending on Kc, (h) dissolved lead depending on Kf , (i) dissolved zinc depending on Kc, (j) dissolved zinc depending on Kf .

Fig. 6 :
Fig.6: Observed and simulated temporal evolution of suspended sediment concentrations at BCAl2 during the January 2011 flood event using the best parameter set.

Figure 5 .
Figure 5. Dotty plots of sediment transport model sensitivity analysis: evaluation of the simulated suspended-sediment concentration according to (a) mud percentage in the Alzette riverbed (pma), (b) mud percentage in the Crauthemerbach riverbed (pmc), (c) settling velocity of sand particles (W s 1 ), (d) settling velocity of mud particles (W s 2 ) and (e) critical Shields parameter (θ c ).

Figure 6 .
Figure 6.Observed and simulated temporal evolution of suspendedsediment concentrations at BCAl2 during the January 2011 flood event using the best parameter set.

Fig. 7 :
Fig. 7: Comparison between the simulated fine sediment deposition map and the spatial distribution of the Zn contamination in the alluvial soils of the study area.

Figure 7 .
Figure 7.Comparison between the simulated fine-sediment deposition map and the spatial distribution of the Zn contamination in the alluvial soils of the study area.

Figure 8 .
Figure 8. Observed and simulated temporal evolution of suspendedsediment concentrations at BCAl2 during the December 2011 flood event using the calibrated model.
arguing that, in Sisyphe, due to the cohesive component of the second class