Reservoir Bank Deformation Modeling: Application to Grangent Reservoir

Within inshore or fluvial environments, submerged fine matter mud banks are characterized by a high water conte spatial variability, and a strong deformability. The study of their instabilities induced by the variation of hydraulic stress requires a modeling of sliding, erosion, and deposition mechanisms. In order to predict the impact of dam reservoir emptying on the s immersed upstream slopes, the method of approach to the problem proposed here combines theoretical developments, nume ing, site observations, and measurements. First, the theoretically achieved sliding criterion is compared with unstable m measurements. For more accuracy in the representation of the natural events, the sliding criterion is then integrated within a code which couples the computation of hydrodynamic conditions, the erosion, and deposition of mud and the banks sliding. F results of the combination of all these mechanisms are compared with the variations in the bathymetric profiles obtained on mental site. keywords Cohesive soils; Saturated soils; Bank erosion; Sliding; Dams; Reservoirs. ore aracnd a atial h can , ent eserentacally im of mud cture ugh ilities, e resuses ersion rates ter. eserthe osioved nomedied to imen-


Introduction
The impact of the water/soil interaction is significant on inshore and fluvial environments.Both environments, indeed, are characterized by the presence of high water content fine soils and a strong deformability.These complex fine soils with great spatial variability are subjected to changing hydraulic stresses, which can be periodic ͑due to swell and tides͒ or exceptional ͑due to floods, seisms or dam reservoir emptying͒.
In fluvial environments, dam erection slows down sediment transport and thus generates fine soil deposits within the reservoirs, which form submerged mud banks.When the sedimentation of the reservoir is considerable and dredging is economically impracticable, emptying processes are conducted with the aim of releasing ͑flushing͒ the sediments accumulated as upstream mud slopes.A reservoir can also be emptied to investigate the structure or carry out maintenance work.Any emptying process, although purposes are different in the two last cases, produces instabilities, and submerged sediment displacements.The lowering of the reservoir level increases the velocity of the flow, which then causes erosion and moves sediments back in suspension.The emersion of the mud banks and the carving action from the flow, generates instabilities, which result in sediments sliding back into water.Some coupled phenomena then occur within the river bed ͑ero-sion, transport, deposition in the downstream part of the reser-voir͒.Upon completion of the emptying process and close to the dam vault, the flow velocity is too high for any sediment deposition.All previously eroded and settled materials are then moved back in suspension.The environmental impact of these phenomena can affect water quality, flora, and wildlife.To enhance sediment storage controls and develop effective solutions, we need to improve our understanding of the different phases of the sedimentary cycle and the connections between them.
Several authors are concerned with some aspects of the impact of reservoir emptying on mud bank stability.Chang et al. ͑1996͒, for instance, are interested in the soil erosion caused by emptying processes.Bouchard et al. ͑1997͒ have developed a numerical code, called Courlis, that combines erosion and deposition phenomena.Paquier and Khodashenas ͑2002͒ calculate the volumes of eroded sediments, assuming a specific repartition of the bottom shear stress across the river width.
Emptying is also at the origin of the sliding of a large proportion of the mud banks.Zhou and Lin ͑1998͒ arbitrarily set a sediment maximum slope, which however, is not easy to estimate since it depends on the soil characteristics, the soil consolidation, the slope of the bedrock, and the thickness of the sediment mass.A more sophisticated approach of the soil mechanic aspect of the channel widening is described by Darby and Thorne ͑1996͒ and Darby et al. ͑1996͒ as bank failure may be planar or by rotational slip.
For accurate representation of slides and consideration of the triple interactions of flow, erosion, and sliding, a sliding criterion developed by the writers and based on soil mechanics has been integrated within the Courlis code.The sedimentary unit calculates an advection-diffusion equation for the concentration of fine sediment, written as follows:

COURLIS
where C M = suspended matter mean concentration; u f = flow velocity; q = diffusion flux; and S v = volume sources.This equation is integrated in the cross section to get a 1D transport equation where sink and source terms are then the erosion and deposition flux which can be modelized by the empirical formula.
Krone's law is used to model the deposition surface flux ͑q deposition ͒.We write where v s = settling velocity; cd = critical shear stress for deposition; and y = transversal coordinate.The shear stress ͑y͒ is defined by where = density; g = gravity; V = mean fluid velocity; h͑y͒ = local water depth; n = Manning roughness coefficient; and R = hydraulic radius.The erosion flux is represented by Partheniades' law where ce = critical shear stress for erosion; and M ce = erosion rate.Bottom sediments are described more accurately by discretizing the cross sections transversally to the flow axis.Exchange flux with the bottom is then computed on each element along the cross section with the local expression of the shear stress Eq. ͑3͒ and added to get a total sediment flux ͑kilogram/meter second͒ available for a 1D transport equation.
Sediment deposit is described vertically by different layers of matter where sediment characteristics and properties are assumed to be homogeneous At each time step, the geometry of the bottom sediment is updated and hydrodynamic calculation relaunched.This exchange between both codes allows us to adapt the hydraulic to the bottom evolution due to both sediment transport and bank sliding.

Sliding
Under the action of earth's gravity, part of the mud bank initiates a slide, whereas the other part remains steady.In the appearing failure surface, the soil shear strength is fully mobilized.The present study is conducted to examine the stability of sediment masses in limit balance conditions, i.e., when failure is about to happen.In order to address the problem, a safety coefficient ͑F͒ is introduced.This coefficient, originally defined by Fröhlich and generalized by Bishop ͑1954͒, is the ratio of the shear strength to the shear stress at the level of the failure surface.Two different methods can be used: the global method and the slice method.The global method implies that the environment is homogeneous and isotropic ͑Morgenstern 1960͒.The slice method consists of dividing the unstable volume into a certain number of vertical slices and then studying the limit balance of each slice.
According to Lam and Fredlund ͑1993͒, the three-dimensional study of n slices gives ͑4n +2͒ equations and ͑12n +2͒ unknown factors.8n equations are then missing to solve the problem and, therefore, several simplifying hypotheses have to be introduced.It is generally theorized that the normal force component applied to the lower surface of the slice is applied in the middle of the surface.The number of unknown factors is then reduced to ͑9n +2͒.Consequently, 5n assumptions are necessary to determine all the unknown factors.
They usually focus on the shape of the sliding surface, on the origin or expression of the interslice forces.Bishop ͑1954͒, for instance, supposes that the shearplane is circular and the vertical components of each interslice force applied to a considered slice have equal amplitudes and opposite directions.
Using the observations collected on a experimental site, a flat sliding criterion for a rectangular soil slice ͑Alexis et al. 2001͒ has been first developed.This criterion suits simplified cases where both the shearplane and the water/sediment interface are parallel.However, to examine cases closer to reality ͑shearplane and interface are not parallel͒, a sliding criterion for a trapezoid soil slice is developed ͑Marot et al. 2003͒ with the following assumptions: 1.The material is assumed to be: ͑1͒ homogeneous; ͑2͒ saturated; and ͑3͒ normally consolidated ͑overpressure is zero compared to the surrounding pressure͒.
2. The failure occurs in undrained conditions.Justification of assumption 1͑1͒: the profiles of total stress and undrained shear strength obtained by different researchers ͑Mekerta 1995; Alexis et al. 2001͒ for this type of soil allow us to suppose the sediment massif to be homogeneous.
Justification of assumption 1͑2͒: given the low incidence of nonsaturation on the effective stress ͑Marot et al. 1998͒ we can assume that the sediment massif is saturated.
Justification of assumption 1͑3͒: sedimentation and sediment consolidation background storage reservoir is often regular.We can thus suppose that the sediments are normally consolidated prior to emptying.
Justification of assumption 2: to carry out laboratory drained shear tests on fine soils samples of low thickness, the shear velocity must be very low ͑lower than 25 m / min or 3.6 cm/ day ͑AFNOR 1994͒.However in situ shearing involves great thicknesses of sediments, and the speed of the drawdown of the storage reservoir, which heavily affects the shear velocity, is generally about several tens of centimeters per day ͑with a maximum of 10 cm/ h͒.We can therefore suppose that shearing appearing at the time of the slope slide occurs in undrained conditions.
In accordance with the assumptions of Bishop ͑1954͒, we suppose that the vertical components of the interslice forces are equal.
The safety factor is then expressed as Slice i weight is expressed as with ␥ = volumic weight of soil When sliding occurs, the remaining stable underlying soil course becomes overconsolidated.So as to consider this phenomenon, the shear strength profile on the potential slip line defined by dimensions T i and T i+1 in undrained conditions is defined by Coefficient a and c u0 are measured on the experimental site and, in order to differentiate remolded sediments from undisturbed sediments, critical remolded erosion stress and remolded undrained cohesion values are assumed to be half that of the undisturbed sediments ͑Marot et al. 1998͒.
In order to perform a continuous survey of the sediment mass stability, three different cases regarding loading must be considered for each slice ͑fully submerged, partly submerged, and fully emergent͒.A different expression of P w and ⌬f corresponds to each leading to: 1. Submerged massif.The force f i can then be expressed by 2. In part emergent massif.We combine Eq. ͑11͒ and and

Modeling Bank Deformation
The bottom cross section is divided in a set of finite elements ͑segments͒.When a segment is unstable, sliding will result in a smoothing of the bottom elevation values; for an example if the element ͑i , i +1͒ is unstable, sliding would induce a decrease of Z n i value and an increase of Z i+1 n .Mathematically, such a smoothing can be simulated by a diffusion operator applied to the bottom elevation of the unstable parts of the cross section.It can be written as where ͓M͔ = diffusion matrix, the coefficients of which are fixed according to the stability of the element, and = numerical parameter which controls the deformation speed.This parameter is calibrated to avoid any time lag between the loading condition variations and the bottom shape adaptations.
The diffusion operator induces the conservation of mud volume.Such a description implies that at each time step, under variable external conditions, the cross section tends to come back to the stability of all the segments and that no inertial effect exists in the movement ͑slow deformation͒.
The diffusion matrix is built following the finite elements formalism The contribution of segment ͑i , i +1͒ is defined according to two configurations: • If the segment is stable The matrix is then formed by adding the contributions of each segments.For a clearer definition of the matrix, point i is studied depending on the stability or the instability of segments ͓i −1,i͔ and ͓i , i +1͔.Four configurations are possible: 1.Both segments ͓i −1,i͔ and ͓i , i +1͔ are stable: line i is defined by: ͑ 0 1 0 ͒, thus point i does not move.2. Both segments are unstable: line i is defined by The "diffusion" matrix is obtained, point i moves, according to point i − 1 and i +1.
The "mass" matrix is obtained, point i moves according to point i + 1 only.4. With the inverse configuration ͓͑i −1,i͔ unstable, ͓i , i +1͔ stable͒, the symmetric matrix in relation to the main diagonal is obtained.Whatever the configuration, the total sum of each line of the matrix is equal to one, which guarantees sediment volume conservation.
At each time step, erosion and settling is computed and the sediment transport equation is solved.The cross-section shape is then modified according to the stability criterion to sliding.The final cross-section geometry, is then sent to the hydraulic code for flow characteristic updating.

Application to Grangent Reservoir
In order to validate the coupled model, the theoretical results are compared with the measurements carried out on the experimental site.

Site Features
The Grangent storage reservoir is located near Saint Etienne, France on the longest French river, the Loire river ͑1,020 km͒.The dam structure is a cylindrical concrete vault with a height of 55 m and a crown length of 200 m.It is one of the 2,500 dams found in the world with a height within the range 50-100 m ͑Lemperiere 1994͒.The storage reservoir has a length of 23 km, a volume of 57 Mm 3 , whereas the surface area of the total catchment is 3,850 km 2 .Into this reservoir, the Loire receives three tributaries ͑the Semène, the Ondaine, and the Izeron͒.During partial emptying processes ͑October 1995-May 1996͒, the level of the reservoir was depleted down to an elevation of 404 m, i.e., 16 m below its nominal elevation.Upon completion of the partial emptying, the upstream end of the reduced reservoir was exactly at the confluence of the Loire and the Ondaine.The study zone spreads upstream and downstream from the junction over ϳ4 km.Four zones are distinguished: upstream, Ondaine, afterconfluence, and downstream ͑Fig.2͒.

Topography and Bathymetry
A topographic survey carried out in 1952 at every 400 m was used to know the elevation of the Loire coarse alluvial deposit before the construction of the dam.
In order to assess the impact of the lowering of the reservoir

Geotechnical Characterization
To characterize the submerged soils before the beginning of the emptying process, a core sampling campaign is conducted and 11 samples with a diameter of 6 cm and a mean length of 2.5 m are taken according to the following distribution: • −4 in the upstream area, station Numbers 2, 3, and 5, • −1 in the bed of the Ondaine, station Number 6, • −3 in the after-confluence area, station Numbers 8 and 9, and • −3 in the downstream area, station Numbers 11, 12, and 13.These samples are subjected to measurement every 10 cm on undisturbed then remolded sediments.One thousand tests have been performed.The parameters measured are: mineralogy of clays, particle size analysis, organic mat-ter content, settling velocity, initial rigidity, water content, bulk density, degree of saturation, methylene blue value, and undrained shear strength.
The degrees of saturation obtained when measurements are carried out under atmospheric pressure conditions, are low ͑S rmoy =77%͒ and reveal the presence of a considerable amount of decomposition gas.In Marot et al. ͑1998͒, the degree of saturation of in situ soils, which is above 90% for in situ materials, has been used to calculate the effective stress of in situ sediments and demonstrate the low impact of the in situ pressure for the water heights considered here ͑between 3 and 28 m͒.
The risk of error is determined by the number of individuals and by considering that the hyperbolic tangent argument of the correlation coefficient follows a Gaussian law when there is no connection between two variables.A more detailed definition of risk of error could be found in Moreau and Mathieu ͑1979͒.

Observation Campaigns
In order to characterize the soil displacements, which occurred during the partial emptying process, two observation campaigns have been conducted.The first campaign was carried out during the lowering of the reservoir level to determine its impact on mud bank stability.The second one took place 4 months after the depletion of the reservoir, as the level was kept low to examine mud flats.The following observations were made: • No sliding above the initial water level; • Mud sliding all along the Loire banks; and • No change in mud bank geomorphology over the period, during which the reservoir was kept at its lowest level ͑see Figs. 4  and 5͒.
Slides are dependent on the water level and the lowering process ͑Fig.6͒.The mud layer which is remained undisturbed on the river banks, is uniform and flat.The type of sliding observed here is of the flat type and is due to shearing of the mud layer.

Comparison between Modeling and Experimental Measurements Calibration
The model parameter calibration is a necessary step for the validation of any calculation code.In the present code, the parameters are three in number: 1. Krone's law for deposition Eq. ͑2͒ implies to know the critical shear stress for deposition and the settling velocity.The studies conducted by the IWW Institute ͑Westrich et al. 1996͒ using a straight flume have made it possible to measure very slow settling velocities ͑50 m/s͒.2. Partheniades' law for erosion Eq. ͑4͒ is based on the coefficient of erodibility and on the critical shear stress for erosion.Both parameters are used to modify the shape of the erosion channel.The first computations were performed on Profile P32 to try to identify the part played by erosion, deposition, and sliding in the modeling of the cross-section evolution during partial emptying.The results are presented below.

Comparison between Erosion-Deposition and Bathymetry
In Fig. 7, the sedimentary surface is computed using erosion and deposition phenomena only.Sliding is ignored.Fig. 7 also shows that erosion and deposition calculations cause a narrow channel, the bottom of which reaches the bedrock.Moreover, this simulation does not allow for the representation of mud displacements taking place at the level of the bank emergent parts.

Comparison between Sliding and Bathymetry
The sliding criterion has been first validated by comparing measured ͑h m ͒ and modeled ͑h e ͒ displaced heights under both limit loading conditions exclusively, i.e., completely submerged or completely emergent massif ͑Marot et al. 2003͒.
The sliding criterion gives a satisfactory estimate of displaced heights.However, the computation produces a height, which is systematically underestimated by 0.3 m.
In order to perform a continuous survey of mud bank geomorphologic changes and validate the integration of the sliding criterion into Courlis, computations are made by neglecting both erosion and deposition phenomena ͑Fig.8͒.The results obtained

Comparison between Erosion-Deposition-Sliding and Bathymetry
Figs. 9, 10, and 11 present the results achieved when the three phenomena are combined for Profiles P32, P31, and P33, respectively.The upper part of the slided mud layer is not modified, whereas erosion computational results are still characterized by the presence of a too deep and too narrow channel.
This first comparison proves that the sliding criterion satisfactorily characterizes mud bank stability.The results, however, are strongly dependent on the coupled erosion computations.The emergent massif sliding initiated on the right of the profile could develop further if erosion had removed the passive earth pressure formed by the submerged massif.Furthermore, the presence of an erosion channel up to bedrock is probably due to shear stress repartition on the bottom.According to relation ͑3͒, the local shear stress strongly increases with the flow depth.Paquier and Khodashenas ͑2002͒ propose some other expressions for the local stress.It could improve the erosion computations and, consequently, improve the coupled sliding results.

Conclusion
When dam reservoir mud banks are subjected to stress variations under the action of exploitation processes ͑emptying, for in-stance͒, instabilities occur.For a more accurate representation of the natural phenomena resulting from the combination of erosion, sliding, and sediment deposition, an alternative approach, which makes it possible to integrate a sliding, criterion into an erosiondeposition code, is developed.
The erosion-deposition code is based on the external coupling of a hydraulic unit and a sedimentary unit.The hydraulic unit solves the 1D shallow-water equations in cross sections along the flow axis and the sedimentary unit calculates the erosion and deposition flux with each element along the cross section transversally to the flow axis.The flat sliding criterion proposed here takes the specific features of high water content soils and is based on the theoretical formulation of the safety factor when the massif is submerged, partly and fully emergent.When an element of the bottom cross section is unstable, sliding will result in a smoothing of the bottom elevation values.This smoothing is simulated by a diffusion operator, which guarantees sediment volume conservation.At each time step, the geometry of the bottom sediment is updated and hydrodynamic calculation relaunched.Both the sliding criterion and the erosion-deposition laws used in the numerical code require a number of parameters, which have to be calibrated.Two bathymetric surveys were carried out on the whole Grangent reservoir ͑before partial emptying of the reservoir and upon completion͒.The submerged soils are characterized by a geotechnical campaign.The computations are shown to agree with the measurements carried out on the selected site, whereas mud heights implicated in slides are satisfactorily estimated.The results, however, are strongly dependent on the coupled erosion computations and the presence of an erosion channel up to bedrock is probably due to shear stress repartition on the bottom.Another expression for the local stress could improve the modeling.The approach presented in this paper for reservoir sediment control will be useful for other reservoirs.
has been developed at the Laboratoire National d'Hydraulique et d'Environnement, a Department of Electricité de France.It is based on the external coupling of a hydraulic unit and a sedimentary unit described in Bertier et al. ͑2002͒.The hydraulic unit is the software MASCARET which solves the onedimensional ͑1D͒ shallow-water equations, in cross sections along the flow axis.It allows the computation of unsteady supercritical flows.A more detailed description of MASCARET could be found in Goutal and Maurel ͑2002͒.

͑6͒
where ⌬f = horizontal components of the interslice forces; f i = horizontal component of the force applied to the slice i by the slice i −1; c u = undrained shear strength of the soil on the potential slip line; L = length of the slice; b = width of the slice; W = weight of the slice; P w = force resulting from the pressure of the water applied to the upper surface of the slice; ␣ = inclination of the lower surface of the slice; and = inclination of the upper surface of the slice.In Fig.1Z i 0 and Z i+1 0 = elevations of the top slice in the initial configuration; Z i n and Z i+1 n = elevations of the top slice in configuration n; T i and T i+1 = elevations of the potential line of sliding, i.e., the topographic dimensions; and H i n and H i+1 n = heights of water above the top slice in configuration n.

Fig. 3 .
Fig. 3. Total stress profile; effective stress profile; and undrained shear strength profile

Fig. 11 .
Fig. 11.P33 profile modeling using Courlis ͑erosion-deposition-sliding coupling͒ ϭ width of slice; C M ϭ suspended matter mean concentration; c u ϭ undrained shear strength on potential slip line; c u0 ϭ undrained shear strength on zero point of profile; F ϭ safety factor; f i ϭ horizontal component of force applied to slice i by slice i −1; ͓M͔ ϭ diffusion matrix; M ce ϭ erosion rate;N ϭ number of individuals; n ϭ Manning roughness coefficient; P w ϭ force resulting from pressure of water applied to upper surface of slice; q ϭ diffusion flux; q deposition ϭ deposition flux; q erosion ϭ erosion flux; R ϭ hydraulic radius; r ϭ correlation coefficient; S v ϭ volume sources; T i ,T i+1 ϭ elevations of potential line of sliding, i.e.,