Large-scale morphodynamics structures in the Arc en Maurienne River (France)

This work focuses on the morphodynamic study of alternate gravel bars in a engineered montainous alpine river, the Arc en Maurienne (France). The experimental site is a selected long reach, which evolved from a natural braided pattern into a single channel formed by two straight reaches, connected by a curve acting as a forcing point on the sediment motion. The channel’s width variation also affected its morphology, leading to the formation of an alternate bar system. Numerical results shows that the Telemac-Mascaret modelling system is able to reproduce the morphodynamics behaviour observed in a gravel bar-scale reach (approx. 1 km) and the alternate bar patterns observed in a 8 km river reach of the Arc en Maurienne river.

Abstract-This work focuses on the morphodynamic study of alternate gravel bars in a engineered montainous alpine river, the Arc en Maurienne (France). The experimental site is a selected long reach, which evolved from a natural braided pattern into a single channel formed by two straight reaches, connected by a curve acting as a forcing point on the sediment motion. The channel's width variation also affected its morphology, leading to the formation of an alternate bar system. Numerical results shows that the Telemac-Mascaret modelling system is able to reproduce the morphodynamics behaviour observed in a gravel bar-scale reach (approx. 1 km) and the alternate bar patterns observed in a 8 km river reach of the Arc en Maurienne river.

I. INTRODUCTION
Rivers generally evolve from an initial single thread into a configuration which is averagely in dynamic equilibrium, also called quasi-equilibrium state [12]. The mutual interaction between hydraulics and sediment transport processes controls the shape of these dynamic systems. The evolution of the channel pattern is therefore dominated by the sediment motion, through the erosion and deposition of the alluvial material [2].
Rivers are subjected to either natural perturbations (climatic and tectonic) or anthropic pressures (water retention, channel derivation, construction of groynes, etc.), conducting to variations of the solid transport rate. Such changes trigger the modification of the river bed, towards a new equilibrium state. In this way, the system's morphological response is adapted to the given perturbation type and magnitude. Since the industrial revolution, river engineering strongly developed as a response to special needs like land reclamation, flood control or water retention and distribution. Several effects of artificial structures on the dynamics of the river can be discerned, but can also lead to side (undesired) effects. Disturbances can be originated from hydrodynamics, by controlling the upstream discharge, from morphodynamics, by dragging or limiting the sediment supply, or finally because of geometric changes of the river course (channelization, deviation, etc.). Depending on the perturbation's importance, the river morphology adjustments in space can be limited to little areas or expanded to large ones, and evolving in time from weeks to several decades.
The Arc en Maurienne (France) is an alpine mountainous river located next to Sainte-Marie-de-Cuines, approximately at 35km of the confluence with the Isère river. This river has been subjected to a large number of training works during the last centuries, where the consequences of such human activities led to considerable changes of the river morphology and dynamics. As a response of land reclamation and flood control, a wide operation of channelization started from the 19 th century, by means of embankments. A succession of hydro-power dams was built during the 20 th century, regulating the flow discharge and reducing the sediment supply. Due to artificial embankments, this river has evolved from a natural braided pattern into a single channel formed by two straight reaches, connected by a curve acting like a forcing point on the sediment motion and presenting spatial width variations. In response to these perturbations, the channel evolved from a braided natural pattern into a single-thread river confined in a straight channel. Consequently, a sequence of alternating deeps and shoals developed, with impact on the navigation, flood control and ecosystem.
In the last few years, prediction of bar formation and development has been considerably improved through the application of analytical theories and empirical and physicallybased models. Nevertheless, these models usually fail on predicting the behaviour of these macroforms for complex cases, wherein most of hypothesis become invalid. Therefore, the study of such bed forms has been promoted by the use of numerical models. Although these models have shown to successfully reproduce the formation and behavior of bars, i.e. straight channels [1], [4], [14], curved channels [3], [13], or width varing channels [5], their domain of applicability has been restricted to simple geometrical, hydrological and flow and sediment characteristics.
This work focuses on the numerical simulation of an alternate gravel bar system in a selected reach of the Arc en Maurienne river. Previous studies have been presented by Jodeau [11], Jaballah [8] and Jaballah et al. [9]. The aim of the present paper is twofold. First, the development and validation of a morphodynamic model at a gravel bar-scale reach of the order of 1 km, using data recorded during a dam flushing event; and second, the application of the model for the prediction of alternate bar patterns for a 8 km river reach. Numerical simulations were performed with the Telemac-Mascaret modelling system.

A. 2D flow model
The hydrodynamics solver Telemac-2D is internally coupled to the sediment transport and bed evolution module Sisyphe in order to investigate the behaviour of bars in the study site. The hydrodynamics module is based on the solution of shallow-water equations (SWE) obtained from several strong assumptions (hydrostatic pressure distribution, vertical acceleration negligible, etc.), wherein turbulence effects are took into account using a constant viscosity model: (1) where g is the acceleration of gravity = 9.81 ms −2 , h is the water depth [m], z s = z b + b the free surface [m], with z b the elevation of the bottom topography above datum, u (resp. v) the fluid velocity along the Cartesian x−axis (resp. y−axis) [m/s], ν t the depth-averaged eddy viscosity, and F x (resp. F y ) the source terms along the Cartesian x−axis (resp. y−axis).

B. 2D morphodynamics model
By considering only bedload transport, the bed evolution is computed from the Exner's sediment mass balance equation: where ǫ 0 = (1 − λ) with λ the bed porosity (= 0.40), (q bx , q by ) = q b (cos α, sin α) are the bedload components in the x− and y− directions, respectively, with α the deviation angle between the sediment transport and the flow direction and q b the sediment transport rate computed with the Meyer-Peter and Müller (MPM) formula: with θ the non-dimensional Shields parameter, θ c the critical Shields parameter (= 0.047), and µ a correction factor that depends on the ration between the total roughness and the skin roughness of the bed. Calibration of MPM formula for α and γ coefficients lead to other formulations popular in the literature [16]. The sediment transport model is parametrized in order to reproduce both relevant physical processes corresponding to (i) bed slope effects and (ii) secondary currents effects.

1) Bed Slope effects:
Bed slopes causes the increase (decrease) of bed-load transport rate in the downslope (upslope) direction. In this work, bed slope effects are accounted using the Soulsby's expression for magnitude (1997) which allows the correction of the original critical shear stress θ c . The nondimensional critical Shield stress θ co is modified as a function of the bed slope χ, the angle of repose of the sediment φ s and the angle of the current to the upslope direction ψ: and the correction of sediment slide direction is implemented using Talmon et al. [15] formula: with β 2 a coefficient.
2) Secondary currents effects: Three-dimensional effects due to helical flows [17] generated in curves in 2D models can be parameterized with semi-empirical formulation [6], [18]. In this work, the Engelund formula is used. This expression is based on the assumption that the bed shear stress, the roughness and the mean water depth are constant in the crosssection, using the expression tan δ = 7 h r , where δ is the deviation angle between the main flow direction and bed load direction and r the radius of curvature of the bend.

A. Study reach
The Arc River is a mountainous river which springs from its headwaters in the French Alps, flowing until its confluence with the Isère River, through the narrow Maurienne watershed of 1957km 2 (Fig.1a,b). The economic importance of the valley has been shown since the industrial revolution, connecting France to Italy through a diversified transportation network and the implantation of a chain of hydroelectricity generating stations.
The experimental site is a 8km long reach, located next to Sainte-Marie-de-Cuines at downstream. The site is located at approximately 19km of the lower dam of S t Martin-la-Porte (Fig.1b). Constraints due to steep engineered embankments around the channel led to the assumption that active width is equal to the river width (for more details, the reader is referred to [9]). The bed slope varies from 1.1% upstream to 0.5% at the downstream part, showing a transient behavior between suband super-critical flows and presenting a width variation from 35m upstream to 50m downstream. In this configuration, the reach can be divided in two channels connected by the curve located at the KP (Kilometer point upstream the confluence with Isère River) 36.22 (Fig.1d).

B. Hydrology of the reach
The Arc River shows common alpine river characteristics, with a nival hydrological regime which can heavily fluctuate following the seasons (i.e. storage and melting of water supply due to hot/cold seasons). Mean discharges recorded varies from 6-8m 3 /s during autumn/winter and 15-20m 3 /s during spring/summer. In addition to channel deviations, the stream discharge in the study reach is controlled and regulated by the sequence of dams located upstream. According to Hydratec and Cemagref [7], only 5% of the river remained with its natural flow and morphological conditions, highlighting the artificialization of the reach.
Due to the low dams storage capacity, water released from dam flushes do not exceed a one-year return period flood, so extremes floods statistics are hardly ever altered. The 10 and 100-year return-period floods monitored at S t Michel-de-Maurienne (1b) were estimated equal to 300m 3 /s and 660m 3 /s respectively [7]. Every year, a dam flushing event is operated in the beginning of the spring season, in order to release fine sediments trapped by the dams. This artificial flood show high levels of concentrations of fine suspended sediments with a discharge equivalent to the 1-year return-period flood equal to 130m 3 /s. gathered respectively before and after the dam-flushing event according to Jaballah's surveys identification (2012) (Fig1c). M02 (resp. M03) is described using 664 points (resp. 839) with a mean point density of 0.03 (resp. 0.037), showing a diminution of the bar height from 2.40m to 2.25m, and a bar mean surface level decreasing of 19cm. Topography reconstitution from terrain measurements is obtained using the breakline method [8], [11], [19]. A set of cross-sections gathered in 2006 [9] with a sparse mean spatial step of 1km is used to describe the whole reach topography.
2) Hydraulic measurements: A full dataset of hydraulic measurements recorded during the dam flushing event is provided including discharges, water stages, and profile velocities using the LS-PIV (Large Scale Particle Image Velocimetry) technology [10], [11]. The discharge was measured at the stream gauging station of Sainte-Marie-de-Cuines monitored by Irstea [7]. The previous discharge associated a gauge measurement at the same location was useful to draw a rating curve linking stage to discharge. Three water stages (G1-G3) were measured using stream gauges and two LS-PIV cameras were set-up to survey water surface velocities over C and D areas (Fig.1c) at respectively two and four different stages of the event.

BAR-SCALE REACH
In this section, the validation of the model is presented selected reach of the Arc en Maurienne river. The choice of the study zone was motivated by previous studies and the access to measurement data, see [8], [9], [11]. Experimental measurements [11] showed the particularity of the reach to be subjected to subcritical and supercritical flows, depending on several parameters such as the bed slope and the input discharge. Therefore, it has been considered essential to investigate the reach's flow regime by building both subcritical and supercritical models and by imposing different formulations for the boundaries. The imposed upstream discharge was originally measured at the downstream part of the reach, at the stream gauging station of Sainte-Marie-de-Cuines monitored by Irstea [11]. Upstream and downstream discharges were supposed equal due to the small area of the domain and the negligence of source terms such as precipitation or infiltration. A rating curve linking stage to discharge is introduced as the downstream boundary condition of the subcritical model, allowing more flexibility by canceling the time-lag effect between the imposed upstream discharge and the downstream water elevation originally measured at the same location. Water stage G3 is chosen as the second upstream boundary condition of the supercritical model. For the subcritical one the boundary condition was extended of an approximate distance of 400m upstream, allowing more flexibility by flow stabilization.
The supercritical model was unable to generate a relevant stationary flow, showing the importance of the downstream boundary by the predominance of subcritical flow. The hotstart generation (uniform steady flow) using a constant discharge of 5m 3 /s has been well reproduced by the subcritical model using a 20,000s simulated time. At this given discharge, only two areas of supercritical flows located close to the bar head show Froude values in the range [1-1.2].
2) Mesh generation: Bar-scale topographic datasets M02 and M03 (cf III.D.1) were used for this study. Several computational mesh grid were set-up to perform the numerical tests (cf IV.A.3). Previous studies showed that the mesh can efficiently capture the information with an average resolution of 2m, with a minimum cell size approx. 0.5m in steepest slopes areas regions. Therefore, the mesh size has been directly linked to the bed slope S 0 calculated with ArcMap. The operation was only led in the area of interest. Beyond this area, the cell size was kept constant at 2.5m. The mesh density dens was determined by using the following linear-thresholded function: 3) Mesh convergence analysis: In order to assess the ability of the model to reproduce the field data observations, a mesh convergence analysis was proposed. This operation consist on creating several meshes, by using a coarse reference mesh template and increasing its cell density by a factor of 2 at each step. This pre-processing was performed by using the module Stbtel of the TMS. A single triangle is divided into four small ones, where the centered corners are defined by the middle points of each the initial triangle sides (Fig. 3). Three meshes were created, where the coarsest mesh shows a maximal density of 10m (C), 5m for the medium one (B) and 2.5m for the finest one (A). Results of this convergence analysis were compared using the water stages G1-3 and the profile LS-PIV velocities C1,2 and D1-4. The simulations were launched using parallelism MPI library, with a 8-core processors Intel(R) Xeon(R) CPU E3-1240 v3@3.40GHz 3.39GHz.
According to the numerical results (Fig.4a,b), the coarsest mesh (C) showed overestimated water depths, whereas meshes (A) and (B) showed water depths in the same range. We can notice the general tendency of a general water stage decreasing when the mesh resolution increases. The cumulated absolute error was calculated for each gauge where err = t=86400 t=0 (y t −x t ) 2 using x experimental data and y numerical results with time-step of 250s, showing the same conclusions as stated before.  Comparison between computed and observed velocity profiles (Fig.4c) was done for the medium (B) and the fine (C) meshes. The last comparison was made at the D-section, located at the bar's tail. The velocities distribution across the river's section justified the right bank erosion and the local deposits at the bar's tail. The velocity's curve can showed the difficulties of the model to reproduce the local variation due to the topographies changes during the flood and could arise from the model's limits, but on the whole the physicalmechanic process can be assumed to be well reproduced with both meshes. To conclude, the medium mesh (B) was clearly eligible to carry out sensitivity analysis and qualitative tests, giving a satisfactory compromise between the results obtained by a particular mesh and the computational time. The choice of the fine mesh (A) was advantaged to conduct quantitative processes and improve the accuracy of the (B) model's results.

B. Hydrodynamic model calibration
The horizontal spatialization process was handled by the Telemac-2D subroutine corfon.f, where Strickler's coefficient zones are defined using the polygon Fortran function inpoly(). Outcomes together with aerial photographs showing vegetated areas and granulometry maps were useful for the model calibration. Results of the model calibration are shown in Fig.4.

C. Pre-investigation of morphodynamics behavior
Numerical computations of the bed shear stress were performed for a preliminary analysis of the system behavior and its morphological response. Figure 2 shows the differences between the dimensional critical Shields parameter τ c and the dimensional bed shear stress τ . The analysis was carried out for constant grain sizes of 2.5cm, 5cm and 10cm and for a grain size-distribution (Fig.2). For the high discharge values, the difference between τ c −τ can reach values up to 40Pa. The importance of sediment sorting along the domain, which is a characteristic of fluvial reaches presenting gravel-sanded bars, is illustrated in Figure 5). The finer sediment located over the bar surface are prone to erosion, even for medium discharges (Fig.5c).

D. Morphodynamics of the gravel bar unit
Sensitivity analysis were conducted to evaluate the response of the system in function of several parameters, such as the selected bedload transport formula, the activation of bed   1) Sensitivity analysis: A sensitivity analysis was performed to investigate the model response to the variation of different parameters. This investigation started by focusing on the bedload transport capacity (a), calculated with different sediment transport formulas (Meyer-Peter-Müller, Einstein-Brown, etc.) found in the literature. The granulometry effects (b) over the reach morphodynamics were studied by considering uniform and non-uniform sediment distributions, using various particle diameters. Last operations were done by considering the bed slope effect (c) induced by gravity and the secondary currents (d) activation, which allows the reproduction of helical flows induced by 3D effects in 2D-H models.
a) Bedload transport formula: Analysis of the reach morphodynamics evolution in function of the sediment transport capacity formula was done, using a uniform sediment distribution with a particle diameter equal to [1;3;5;7;10] cm. The choice of bedload transport capacity formula has a great impact on the volume of sediments eroded and deposited in the domain of interest (Fig.6). Einstein-Brown and Van Rijn formulas showed results in the same range in comparison with MPM one which slightly overestimated the phenomena of erosion/deposition (e.g. the computed q b is greater). Engelund-Hansen formula is mostly useful to determine the total load transport.
Therefore, results differ from the previous ones, especially for particle diameters inferior to 5cm, showing more sensibility for the sediment diameter fluctuations. However, the erosion/deposition mechanic process was hardly ever perturbed, where locations of erosion and deposition remains quasiidentical.
A sensitivity analysis was done on the MPM α coefficient (cf II.B) for α = [2; 3; 4; 5; 6; 7; 8], using a uniform grain size distribution of 5cm. As expected, variation of α coefficient had an impact on the solid transport rate, which became greater when α increased (Fig.7c,d) and was also confirmed by the computation of positive and negative volumetric changes. MPM α coefficient had strong impacts on the riverbed evolution process, hence the coefficient is suitable to calibration.
b) Bed composition: The bed has been first described with an uniform sediment distribution. Particle diameters of [1;3;5;7;10] cm were considered using several sediment transport capacity formulas (cf IV.D.1.a). Fine sediments seems to be subjected to more transport (e.g. more erosion and deposition), whereas coarser gravels have the tendency to resist to the flow and stay unmoving. This physical phenomena is depicted by a diminution of erosion/deposition in amplitude, but also in space where differences can be mostly found over the bar close to transverse channels and zones of confluence (Fig.7b,c).
As a natural alpine river, the Arc is composed of a large mix of sediments. The river bed was described using a nonuniform distribution of sediments, where the particle diameters are in the range of [1;5;10] cm. Classes of sediments with their respective fraction were specified in the Sisyphe module. Riverbed evolution is calculated by calculating the solid transport discharge q b,i for each class of sediment i. Hence, the total computed solid discharge is equal to The bed-material mixture, described by the fraction of each sediment class, strongly affects the morphodynamic evolution of the reach both in space and amplitude. A poorsorted granulometry (40% of 1cm -20% of 5cm -40% of 10cm) seems to increase the phenomena of erosion-deposition (Fig.7e). While the variation of the particle diameter has little impact in space on the evolution of an uniform bed, the mixture of sediments in non-uniforms riverbeds seems to alter stronger the morphodynamics in space. This can be mostly observed close to the bar's tail and close to the transverse channels.
c) Slope effect: The actual model seemed to be much more sensible to the deviation effects than the magnitude ones. This has been demonstrated using Soulsby's formulation for magnitude (Eq.4) and Talmon's et al. formulation for deviation (Eq.5). Deviation correction effects are controlled by the calibration of β 2 . Low values of β 2 (Fig.7g) showed overestimated erosion of the bar and deposition in the main channel, whereas higher ones (Fig.7h) showed less sediment Fig. 7: Comparison of bed river evolution after the dam flushing event (erosion/deposition) using : (a) experimental data; uniform sediment distribution of (b) 1cm (c) 5cm (d) 5cm and α=5; (e) non-uniform sediment bimodal distribution; (f) secondary currents; bed-slope effects using β 2 = (g) 1.6 (h) 3; non-uniform spatialized sediment distribution with α=5, β 2 =2. 5 and secondary currents for a main channel composed of (i) big boulders (j) mix of big boulders and sand.
motion. Using high values of β 2 diminished the bed slope effect, and tended to converge to the state without bed slope consideration (Fig.7c). However, values in the range [1;3] seems to provide a physically relevant compromise. A general smoothing of results can be noticed once the bed-slope effect is activated due to the diffusivity effects.

d) Secondary currents:
We can notice the bedload movement deviation from the main flow direction, showing sharp differences over the middle part of the bar, where transverse channels connects the secondary and the main channel together (Fig.7f). Deposits fronts directions in this area were changed, taking into account to the helical flows together with gravitational effects.

2) Morphodynamic model calibration:
Morphodynamic simulations presented here were performed for a time step ∆t = 0.25s, with a MPM α coefficient = 5. At the inflow boundary condition, the equilibrium solid discharge is imposed. At the outflow boundary condition, a free or Neumanntype boundary condition has been imposed.
The non-uniform sediment distribution along the reach (cf III.D.1) is given by Fig.2. Horizontal spatialization process is handled by the Sisyphe subroutine called init compo.f, using the Fortran polygon function inpoly() wherein zones composed of different fractions of each sediment class are defined.
Due to difficulties to conduct measurements in the depthdamped areas, the grain size distribution in the main channel remained a mere estimation [11], leading to assumption that the main channel's bed is composed of large boulders (12.8-25.6cm). The riverbank composition was also roughly estimated, mostly composed of bimodal distribution of boulders and sparse fine sediments deposited after flooding events. Several investigations of the morphodynamic response of the system were carried out, depending of the bed-material composition of the main channel and the riverbanks.
The main channel riverbed composition had a predominant effect on the final morphodynamic evolution of the reach. Indeed, large boulders (20% of 10cm -80% of 18cm) were found to be more stable than a bed composed of finer sediments (50% of 1mm -50% of 18cm). The bar was less eroded for the bed composed of finer sediments (Fig.2i,j), whereas the riverside at the opposite of the bar was strongly eroded and the main channel is filled of deposited sediments. Last investigations were done on the riversides. Although the morphological evolution of the reach was altered by the riversides composition, the amplitude seemed to be lessened compared to the effects induced by changes of the main channel granulometry. This was due to the rare contacts between the flows and the riversides, only triggered during highest flood-stages. This effect is non-negligible, where we observe a stronger erosion of the right riverside and the apparition of a deposit front located in the main channel for finer granulometry distribution. Higher flood-stages could play a major role on the reach planform evolution. This is subject of future investigations.

REACH
In this section, preliminary simulations of alternate bars in a 8 km reach, subject to the influence of constant upstream discharges and constant sediment distribution are presented.
A. Model parameters 1) Mesh and bathymetry: The long-time morphodynamical evolution of the 8km long reach (Fig.1b) was investigated. The simulation started with an almost panform bed, enriched by the M02 bathymetric dataset in the downstream part describing the bar unit (cf III.D.1). The computational mesh-grid is composed of 37.330 nodes and was generated using a mesh size around 5m.