Numerical modelling of bank failures during reservoir draw-down

Numerical algorithms are presented for modeling bank failures during reservoir flushing. The algorithms are based on geotechnical theory and the limit equilibrium approach to find the location and the depth of the slides. The actual movements of the slides are based on the solution of the Navier-Stokes equations for laminar flow with high viscosity. The models are implemented in the SSIIM computer program, which also can be used for modelling erosion of sediments from reservoirs. The bank failure algorithms are tested on the Bodendorf hydropower reservoir in Austria. Comparisons with measurements show that the resulting slides were in the same order of magnitude as the observed ones. However, some scatter on the locations were observed. The algorithms were stable for thick sediment layers, but instabilities were observed for thin sediment layers.


Introduction
A flushing operation is often the most inexpensive method to maintain water storage capacity in hydropower reservoirs that fill up with sediments.Erosion of sediments from reservoirs can also be an important process to consider in decommissioning of dams [1].The current article describes bank failure algorithms that have been included in a fully 3D numerical model for modelling the sediment movement during the flushing [2][3][4][5][6].The novel approach is an incorporation of a soil slide algorithm based on geotechnical theory.Lateral soil slides are common during the flushing process, as the eroding water channel cuts down into the sediments.Widening of the channel and possible meandering/braiding often involve soil slides [7][8][9][10].The width of the flushing river system will determine how effective the flushing is and how much sediment can be removed.The bank erosion algorithm is therefore important in assessing the effectiveness of the flushing, and thereby also in sediment management of a reservoir.This question is also relevant for the engineering design of flushing gates in the reservoir dam.
The current article describes progress in developing algorithms for lateral bank failures during reservoir drawdown.The numerical models are presented first, and then they are applied to the Bodendorf reservoir in Austria.Considerable data from this reservoir have been collected before and after flushing operations, making it a valuable test case for the current study.Fig. 1 shows the bank of the reservoir during a flushing operation in 2007.

Numerical model
The numerical algorithms used in the current study are based on a finite volume method, where the domain is divided into water, soil and slide cells.The soil domain uses a twodimensional depth-averaged grid.The water grid and the slide grid are three-dimensional.Fig. 2 shows a profile of the three grids, along a cross-section of the bank of a river.The soil grid only has one cell in the vertical direction.The water grid and the slide grid have multiple cells in the vertical direction, depending on the depth.The grids are adaptive, and move over time with the changes in bed level, water level and slides.The water grid will cover the reservoir water and the channel that is formed during the flushing.The Navier-Stokes equations are solved on this grid, together with the convection-diffusion equation for suspended sediment transport and a bedload equation [11].The main water flow direction in the reservoir is normal to the plane in Fig. 2.
The computation of the soil slide can be divided into two parts: 1.Where is the location of the slide and what is its magnitude?2.Where does the soil move?
The current study solves the first question with geotechnical theory: The Limit Equilibrium approach [12].The slide is then divided into vertical elements, as seen in Fig. 3: The Limit Equilibrium approach computes the forces on each of the elements (8 in the example in Fig. 3).There are four forces: gravity and groundwater pressure are driving forces, while friction and cohesion are stabilizing forces.Summing up all the forces for the elements, a Safety Factor (SF) can be computed as the ratio of stabilizing forces divided by the driving forces.If SF > 1, the bank is stable, and if SF < 1 the bank will slide.The sum of all four forces on one slide, F, is according to [13]: The area of a slice projected in the horizontal plane is A z, rs is the density of the sediments, rw is the density of the water, p is the porosity of the soil, h is the height of the slice and hw is the height of the groundwater level through the slice.The critical shear stress of the cohesive soil is denoted tc, q is the bed angle of the slide plane, f is the angle of repose for the sediments and y is the slope of the groundwater level.
A problem with the Limit Equilibrium Approach is that the location and the magnitude of the slide can not be computed directly.The procedure involves guesses of where the slide is, and computations of the SF for each slide.The slide with the lowest SF will be the critical slide.In a CFD model it is advantageous to compute the slide location and magnitude directly, without the need to make a large number of guesses.In the current study, Eq. 1 was transformed to the differential equation for the slide movement, dp, of one element, p, as a function of the movement of the neighbouring elements (subscript nb).The equation was solved on the depth-averaged 2D grid, so that each slice had four neighbours, two in each of the two horizontal directions.
The source term, Sd, is made up by the forces from Eq. 1.
The solution of Eq. 2 determined which vertical slices moved and which did not move.The next step was to compute the vertical extension of the slide.This was done by solving another differential equation similar to Eq. ( 2), where the vertical elevation of the bottom of the slide, z, was the unknown instead of d.
( ∑ a nb ) z p = ∑ (a nb z nb )+ S z (3) The source term, Sz in Eq. 3 was negative and chosen so that the maximum slide thickness was 1/10 of its length.
By solving Eq. 2 and 3, the vertical and horizontal location of the slides could be determined.This information was used to generate the grid shown in Fig. 2. The actual movement of the slide was computed by assuming the soil in the slide was a very viscous fluid.The velocity field in the slide was computed from the Navier-Stokes equations, together with the resulting bed elevation changes [11].The currently used algorithms were implemented in the SSIIM program, which already had a Navier-Stokes solver.

The Bodendorf reservoir
The numerical algorithms were tested on the Bodendorf reservoir in Austria.This is a runof-the river reservoir with a weir that includes two 12 m wide segments with heights of 8.5 m.The modelled geometry is approximately 100 meters wide and 3 km long.A plan view of the modelled section of the reservoir is given in Fig 4 .The water is flowing from the left to the right.The dam is at the right boundary of the figure.Profiles of ten cross-sections were measured before and after the flushing in 2004 [13].The profiles and their numbers are given in Fig. 4. The flushing of 2004 caused erosion of the bed in some parts of the reservoir.The water level was drawn down 7.1 meters at the dam during the flushing, and sediment erosion only took place once the water velocity was sufficiently high.This only occurred when the water level reached a low value.The bed elevation changes at the higher banks could not have been caused by erosion.Fig. 5 shows cross-section number 60, with bed levels before and after the flushing.Fig. 5. Bed levels of cross-section 60, before (A) and after (B) the flushing.Fig. 9 shows the computed bed elevation changes in the reservoir due to the sediment slides.The changes varied between -4 and 4 meters.Note that continuity of the slide must be fulfilled, so that each slide will give the same volume for lowering of the bed at the upper end of the slide as the increased bed level at the lower end of the slide.This is seen in Fig. 9 as each slide gives corresponding black (rise) and white (lowering) of the bed levels.
Fig. 9 also shows the measured bed elevation changes at the end of each of the crosssections given in Fig. 4. The numbers were obtained by looking at all the cross-sections in figures similar to Fig. 5, subtracting the bed levels before and after the flushing at the ends of the cross-sections.Fig. 9 shows that there is correspondence between computed and measured values at some locations of the reservoir, but not everywhere.A the downstream part of the left bank, the numerical model computed a fairly large slide which was not observed in the field.
The results shown in Fig. 9 depend on several calibration parameters.The most important parameters are the angle of repose for the sand and the viscosity of the slide.The tangens to the angle of repose for the sand was set to 0.3, which is lower than the physical correct value of around 0.6.This parameter affects how many slides are formed and their horizontal magnitude.A higer value would give smaller and fewer slides.Looking at Fig. 8, the grid is relatively coarse compared to the grid in Fig. 2. It is possible that a finer grid would produce better results with a higher angle of repose for the sediments.
The vertical magnitude of the slides were affected by the viscosity of the slide.The slide was assumed to be a fluid with high viscosity.A value of 8x10 7 m 2 /s was used in the results from Fig. 9.A lower value would give more elevation changes.
The main problem for the current study was to get sufficient stability when solving the Navier-Stokes equations for the slide.Looking at Fig. 8, the vertical resolution of each slide is fairly good, as there are several cells over the depth and each cell has a reasonable height/length ratio.This works well as long as the sediment magnitude is fairly large.In the current study, it was set to 5 meters.However, in reality the sediment magnitude varies and can be much lower.If it is too low, than the slide can not be resolved properly with grid cells of a certain hight/length ratio.Then instabilities were observed and the solution of the Navier-Stokes equation would diverge.

Conclusions
A numerical model for computing geotechnical slides reservoirs during flushing was developed and tested against data for the Bodendorf hydropower reservoir.The main calibration parameters are the angle of repose for the sediments and the viscosity of the sediment slides.The computed slides look reasonable in number and magnitude compared to the field data, but the locations are not always correct.The new implemented algorithms work well for thick sediment layers, but instabilities occur for thin layers.Further research is therefore required.

Fig. 2 .
Fig. 2. Bank profile with grid.GL is the groundwater level.

Fig. 3 .
Fig. 3. River bank profile showing slide and eight elements above the slip circle.

Fig. 4 .
Fig. 4. Plan view of the modelled section of the Bodendorf reservoir.