Introduction

Water is a valuable resource in most parts of the world. It is stored in reservoirs and used for several purposes, such as irrigation or hydropower. Many rivers carry large amounts of sediments, which can deposit in the impounded reservoir behind the dams and seriously reduce the water storage capacity. Mahmood [1] estimated that around 1% of the total water reservoir volume in the world is lost annually due to sedimentation. Several regions in the world show much higher sedimentation values. The most common method to mitigate the problem is to carry out a draw-down flushing by lowering the water levels [2]. This method has shown to be the most cost-effective alternative for many reservoirs. However, the success of reservoir flushing is depending on a variety of parameters [3]. The flushing also has some disadvantages. The environmental impact can be significant in case of an incorrect flushing strategy and/or the absence of sufficient monitoring procedures. Most of the water in the reservoir is lost after the flushing, which can have economic implications. Proper management of a flushing operation must therefore carefully consider the cost and benefits. A crucial question is how much sediments can be flushed out of a reservoir. In recent years, numerical modelling has been applied to a number of reservoirs in order to answer this question [4]. The models can be one-dimensional [5], two-dimensional depth-averaged or three-dimensional [6]. They can also take the large number of processes into account to a varying degree. Some are related to sediment erosion, including cohesion effects and sediment slides. The current study focuses on the process of bank failures during the flushing. The lowering of the water level in the reservoir will cause a lateral groundwater gradient that can force geotechnical slides to occur. The soil movement of the slides can be important in determining the amount of sediments that is eroded during the flushing and to interpret the success of the chosen flushing strategy. A correct modelling of the bank erosion can therefore be significant for the quality of the results from the numerical model in simulating reservoir flushing.

Earlier numerical models of reservoir flushing have mostly been based on approaches using the angle of the repose to take bank failures into account. This approach is fairly simple and assumes that the side slopes of the banks can not be steeper than a prescribed angle of repose. Once the numerical model predicts a steeper bed slope, the bed elevations are changed so that the bed slope is smaller. The physics behind this approach is similar to a small slide occurring in non-cohesive soil without groundwater gradients. The procedure for the movement of the bed elevation can then be thought of as a simple bank slide algorithm. However, there are important shortcomings with this approach. One problem is that cohesion of the bank material is not taken into account. Sediments in a reservoir can often contain a considerable amount of silt and clay, which can stabilize the banks due to cohesion. The banks can in such cases be steeper than the angle of the repose for non-cohesive soil. Another problem is that the main driving force for the bank erosion is most often a horizontal gradient in the groundwater level. As the water level drops during the initial stages of the flushing process, the groundwater causes a hydrostatic pressure gradient in the potential slide direction, which can cause the bank to collapse. A numerical model based on the physics of the bank erosion process must therefore also be able to compute the groundwater gradient.

Geotechnical engineering offer several approaches to analyse bank failures. The main principle is to divide the soil/bank into one or more elements and compute at the forces on the elements. Traditionally, the limit equilibrium analysis was the most used approach for many years. The method was pioneered by Bishop [7] assuming the slide plane would follow a slip circle. The traditional limit equilibrium method investigates several glide planes to determine which of them has the lowest safety factor. The safety factor is defined as the acting forces divided by the stabilizing forces. The bank will fail if the factor is above 1, and it will be stable if the factor is below 1. The number of elements for each slide is usually limited. The elements are also vertical with homogenous properties.

A more modern approach is to divide the bank into a large number of 3D elements and use a finite element analysis to compute the stresses and strains in the soil [8]. The advantage is that more detailed spatial variation of different soil parameters can be used and a more complex situation can be modelled. The finite element approach computes the glide plane directly, without the need for guessing or specifying a glide plane [9,10,11]. It is also possible to use an elasto-plastic model for the soil to include a computation of the soil deformation [12]. This makes it possible to compute where the soil moves during the slide, and where it ends up after the slide has occurred. This is not possible with a traditional limit equilibrium approach. The disadvantage with the complex finite element method is that each slide must be resolved with a large number of grid cells. There can be hundreds of slides in a reservoir, leading to long computational times due to the number of cells.

Several advanced studies based on the limit equilibrium approach have earlier been made on bank erosion at rivers. Darby et al. [13] analyzed bank erosion at a river in Italy including modelling of groundwater levels, lateral erosion and geotechnical bank failures. Bank failure including cohesive material was computed with geotechnical software based on a width-averaged approach. A computer program for groundwater level determination was also used with a 2D width-averaged grid, enabling hydrostatic pressures to be included in the computation. The water in the river was not modelled, and it was assumed that the mass from the failure would instantly be transported away by the water flow. An empirical formula for the lateral erosion (Δy) in meters as a function of the excess shear stress was used:

$$\Delta y = k_{d} (\tau_{b} - \tau_{c} )$$
(1)

The critical shear stress is denoted τc, the bank shear stress is given as τb, and kd is a proportionality factor which must be determined empirically.

Rinaldi et al. [14] and Lai et al. [15] used a similar approach as Darby et al. [13], but did the analysis in three dimensions. However, the shear stress from the water at the banks (τb in Eq. 1) was still computed from a 2D model. Nardi et al. [16] and Langendoen et al. [17] used a similar approach as Rinaldi et al. [14], but only used Eq. 1 to find the lateral erosion without computing the groundwater levels or geotechnical bank failures directly. It was assumed that the main lateral erosion was computed from Eq. 1. Deng et al. [18] also used the same equation to find bank erosion in a large river in China. Marot et al. [5] computed sediment slides in a water reservoir for a very similar situation as in the current study. However, they only used one-dimensional numerical models for the water flow and sediment slide movements.

The problem with Eq. 1 is to find a correct value for kd. This will depend on the soil properties, groundwater level, geometry of the banks etc. The approach used in the current study is similar to the method used by Rinaldi et al. [14] and Lai et al. [15], but Eq. 1 is not used. Instead, lateral erosion is computed from a combination of vertical erosion at the banks and a bank slide algorithm. The limit equilibrium approach is used to find the location and extension of the slides, similar to the method used by Darby et al. [13] and Rinaldi et al. [14]. However, the actual movement of the slides is computed by solving the Navier–Stokes equation, in all three directions, assuming the moving slide is a fluid with high viscosity. This is a less complex approach than using 3D a plastic-viscous model, since our numerical model for sediment transport already includes a Navier–Stokes solver to compute the water flow.

An earlier study by Olsen and Haun [19] also used a numerical model to compute the sediment slides at the banks of the Bodendorf reservoir. However, they used a much coarser grid, requiring an artificially low angle of repose (16.7°) for the slides to occur. The current study uses a finer grid, where the angle of repose was a more reasonable value of 30°. Another disadvantage of the previous study was the occurrence of stability problems with varying sediment thickness, which are solved in the current study. A much more detailed explanation of the numerical algorithms is included in the current article, in the next chapter, together with improved and more detailed figures. An in depth parameter sensitivity analysis is also reported.

Numerical model

The currently presented numerical model for the soil slides is an extension of the existing computer program SSIIM 2 (Sediment Simulation in Intakes with Multiblock option), designed for simulation of sediment transport in rivers and reservoirs [6, 20,21,22,23,24,25,26]. All algorithms are based on the finite volume method and adaptive 2D and 3D grids. The water flow is computed by solving the Reynolds-averaged Navier–Stokes equations with the k-epsilon turbulence model [27]. The sediment transport is computed by solving the convection–diffusion equation for suspended sediment concentration for each grain size separately. In addition an empirical formula is used for calculating the bed load. The grid in SSIIM is adaptive, meaning that the grid moves with changes in the water and bed level as result of an implemented wetting/drying algorithm. This numerical model is further described by a number of articles [28,29,30]. The extension of the model to simulate sediment slides is also based on the finite volume method. The computations are carried out on three types of grids that have common interfaces: A groundwater grid, a water grid and several slide grids.

The groundwater grid is two-dimensional and depth-averaged. It covers the whole reservoir and the banks, and is used to compute the groundwater levels in the reservoir and the banks. It is also used to compute the water surface elevation changes in the water reservoir and the changes in the bed elevations. A fairly coarse 2D grid of the Bodendorf reservoir, used by Haun et al. [4] is shown in Fig. 1. This grid only has 10 cells in the lateral direction. The computational grid used in the present study has 80 cells in the lateral direction. A detail of this finer grid is shown in Fig. 2.

Fig. 1
figure 1

Plan view of the 2D depth-averaged grid for the Bodendorf reservoir with 10 cells in the lateral direction and 195 cells in the longitudinal direction. The arrows indicate the inflow and outflow regions. The computational grid used in the current study used 80 and 1560 cells in lateral and longitudinal direction, respectively

Full size image
Fig. 2
figure 2

Detail of the 3D grid used for the computation of the water flow and sediment transport (plan view). The uneven edges are due to the wetting/drying procedure

Full size image

The second type of grid is three-dimensional and covers the water in the reservoir. Water velocities, suspended sediment transport and bed load is computed in this grid. The grid is located between the water level and the bed level of the reservoir, on top of the two-dimensional groundwater grid.

The third grid type is also three-dimensional, and covers the moving soil within the slide. These grids do not cover the whole computational domain. Instead, there are several of these grids in the reservoir, depending on where slides occur. Figure 3 shows a cross-sectional profile of the three grids.

Fig. 3
figure 3

Cross-sectional profile of a river bank with the three types of grids, used for the slide movement computation. The groundwater grid is given as a, b is the grid for the soil movement and c is the grid for the water/sediment flow. The light blue line is the water surface and the groundwater level. From Olsen [33]

Full size image

The groundwater grids are made up of 2D quadrilateral cells and have a structured layout. The other two grid types are non-orthogonal, unstructured and made up of dominantly 3D hexahedral cells. The hexahedral cells are chosen to improve stability and accuracy of the numerical solution.

The groundwater levels were computed by solving Eq. 2, which is based on the water continuity equation and Darcy’s equation for the relationship between the water velocity and pressure head in a soil element:

$$\frac{{\partial z_{w} }}{\partial t} = \frac{\partial }{{\partial x_{i} }}\left( {K\frac{{\partial z_{w} }}{{\partial x_{i} }}} \right)\quad {\text{i}} = 1,2$$
(2)

where zw is the ground water level to be computed, t is the time, K is the conductivity of the soil, and x is the horizontal spatial dimension.

The water velocities in the reservoir were computed from the Reynolds averaged Navier–Stokes equations. The k-epsilon model was used to determine the turbulence. The water flow solver had previously been tested on a number of cases, for example Wilson et al. [31]. The sediment transport was computed by a procedure including the solution of the convection–diffusion equation for each grain fraction of the suspended sediments. The Engelund and Hansen [32] formula was used as sediment boundary condition at the bed. The continuity equation for sediment transport was used to compute the bed elevation changes due to the resulting erosion/deposition pattern in the reservoir. This procedure is described in more detail by Haun and Olsen [24]. The water level in the reservoir was computed from algorithms solving the diffusive wave equation [29].

The main numerical challenge in the present study was the computation of the geotechnical sediment slides. The approach used for this purpose was similar to the procedure used by Olsen [33], who computed a 2D slide from a laboratory experiment. In the first step, the horizontal locations of the slides were found. This algorithm was based on the limit equilibrium approach [7]. A slide is divided into elements, as shown in Fig. 4. The force on each side in the direction of movement can be computed from Eq. 3.

$$F_{EX} = (1 - p)A_{z} (\rho_{s}\, - \,\rho_{w} )gh\tan (\theta )\, + \,A_{z} \rho_{w} gh_{w} \tan (\psi )\, - \,(1 - p)A_{z} (\rho_{s}\, - \, \rho_{w} )gh\tan (\phi )\,- \,A_{z} \tau_{c} /\sin (\theta )$$
(3)
Fig. 4
figure 4

2D slip circle with seven vertical elements. The length, L, of the slide is defined on the figure

Full size image

The cross-sectional area of the slide normal to the soil flow direction is denoted Az and the depth of the slide is given as h. The soil porosity is denoted p, and the distance from the glide plane to the groundwater level is given as hw. The buoyancy of the submerged soil particles is found by subtracting the water density, ρw, from the particle density, ρr [34]. The slope of the glide plane is denoted θ and ϕ is the friction slope. The cohesive shear strength of the soil is denoted τc. and the slope of the groundwater level is given as ψ. The gravitational acceleration is denoted g.

There are four terms in Eq. 3. The two driving forces are the gravity and the water pressure from the groundwater. These are the first and second terms on the right side of Eq. 3. The two remaining terms in Eq. 3 are stabilizing forces due to friction and cohesion of the soil.

In addition to the forces given from Eq. 3, there are also forces acting between the grid cells. If we look at one cell in the grid shown in Fig. 1, it will have four neighbour cells (2D plan view). The grid cell and its neighbours are shown in Fig. 5a. The center cell is denoted P and the neighbours N, S, E and W, after the four directions North, South, East and West, respectively. The indexes will follow the grid lines, so that the East–West direction is the main flow direction along the river, and North–South is the direction of a cross-section of a river. Since grid and river may flow in any geographical direction, the definition of the North/South/East/West direction in Fig. 5 does not follow the geographical directions.

Fig. 5
figure 5

a Computational molecule for cell P and the four neighbours, S, W, N and E. b Normal (σ) and shear (τ) stresses on the sides of cell P

Full size image

The banks are on the side of the reservoir, so the soil slides will mainly move in the North–South or in the South-North direction of the grid. We initially denote the x-direction to be following the East–West direction of the grid and y is in the North–South direction. The force, FS, on the South side of the element P in Fig. 5b due to the moving soil can be expressed using Hooks law that relates the force to the horizontal movement, δ:

$$F_{S} = A_{S} \sigma_{S} = A_{S} E\varepsilon_{S} = A_{S} E(\delta_{S} - \delta_{P} )/\Delta y$$
(4)

AS is the area of the south side of the cell, while E is the module of elasticity of the soil, σS is the normal stress on the south side and εS is the strain rate of the soil. The horizontal distance between the cells is denoted Δy. A similar equation can be set up on the West side of the cell.

A shear stress due to the movement of the soil will give a force on the West side of cell P, as show in Fig. 5b. This force can be expressed from the relationship including the shear modulus, G, of the soil:

$$F_{W} = A_{W} \tau = A_{W} G(\delta_{W} - \delta_{P} )/\Delta x$$
(5)

The total force on an element is made up by the result of Eq. 2, together with Eqs. 3 and 4 applied on the four sides of cell P. The resulting equation can be divided by the height of the slide cell, resulting in an equation for the movement of cell P, δP

$$\delta_{P} = \frac{{a_{E} \delta_{E} + a_{W} \delta_{W} + a_{S} \delta_{S} + a_{N} \delta_{N} + S_{P} }}{{a_{P} }}$$
(6)

The weighting coefficients, a, are getting the following values:

$$\begin{array}{*{20}c} {a_{E} = \frac{{E\Delta y_{E} }}{\Delta x}} & {a_{W} = \frac{{E\Delta y_{W} }}{\Delta x}} & {a_{N} = \frac{{G\Delta x_{N} }}{\Delta y}} & {a_{S} = \frac{{G\Delta x_{S} }}{\Delta y}} \\ \end{array}$$
(7)
$$a_{P} = a_{E} + a_{W} + a_{S} + a_{N}$$
(8)

The source term, Sp, is made up of the forces from Eq. 3, divided by the element height. The next approximation is to use a limiter on the source term, according to the following formula:

$$S_{P} = \max (\frac{{\left| {F_{EX} } \right|}}{h},0)$$
(9)

The rationale behind this limiter is that the direction of the soil movement is not considered. The soil moving in the positive and negative direction will get the same positive value of δ. The necessary assumption is that each soil slide moves mainly in one direction and does not interact with other slides. The purpose of computing δ in Eq. 6 is not to determine the movement itself, but to find in which cells a slide will occur.

Dividing Eq. 3 by the slide height eliminates the variable h in three out of the four terms. The last term with the cohesion do not contain the variable h. An assumption for how to evaluate cohesion divided by h, needs to be made for this term. The concept of a cohesive height is introduced in the model as shown in Fig. 6. Looking at a bank where cohesion is the dominant stabilizing force, a linear slide plane is assumed along an element of triangular shape with a height hc. The critical value of hc for when the slide moves is a function of the cohesion of the material. A force equilibrium in the vertical direction then gives a relationship between the cohesive shear stress τc and hc.

$$F_{g} = \frac{1}{2}\rho gh_{c} (1/3h_{c} ) = h_{c} \tau_{c}$$
(10)
Fig. 6
figure 6

Definition sketch for determining the cohesive height, hc. Fg is the weight of the sliding element and τ is the cohesive shear stress

Full size image

This gives the following expression for hc:

$$h_{c} = \frac{{6\tau_{c} }}{\rho g}$$
(11)

It is here assumed that the horizontal length of the edge is 1/3 of hc. With a cohesive shear stress of 1000 Pa, hc becomes 0.6 m. This is used as a cohesive height in Eq. 9.

Equation 6 was solved on the same 2D grid that was used in the computation of the groundwater level (Fig. 1). Figure 7 shows a section of the computational domain with location of slides. The grey areas are cells where slides occur. Note that the figure shows the situation at one time step. The location of the slides will change during the simulation, depending on the solution of Eq. 6.

Fig. 7
figure 7

Section of the computational domain (plan view). Grey shaded cells indicate cells where slides occur. The blue lines represent the sides of the domain/reservoir

Full size image

Classical soil mechanics relates the E and G modules through the following equation:

$$G = \frac{E}{2(1 + \nu )}$$
(12)

The ν parameter is the Poisson ratio, which is around 0.4 for soil. This means that the G module is around 1/3 of the E module. The default parameters in the current study assumes that E = G. Then the limiter in Eq. 9 means that Eq. 6 gives the soil slide direction independent of the direction of the grid. The soil slide does not have to move in the North–South direction of the grid, it can also move East–West or in any other direction in the horizontal plane. The assumption E = G will be tested later.

The next step is to determine the vertical magnitude of the slides, h. In the current study, it was assumed that the maximum slide height, H could be computed from the following equation [33]:

$$H = k_{L} L$$
(13)

where L is the length of the slide and kL is a user-defined parameter set to 0.15. The parameters H and L are also shown in Fig. 4. An equation similar to Eq. 5 was then solved to find the vertical elevation, z, of the slide plane elevation in each cell:

$$z_{P} = \frac{{a_{E} z_{E} + a_{W} z_{W} + a_{S} z_{S} + a_{N} z_{N} + S_{P} }}{{a_{P} }}$$
(14)

The weight factors, a, were set to 1 and the source term, S, was varied until the maximum slide depth was equal to the H value from Eq. 13. The boundary conditions for Eq. 14 were the z values outside the slide area being equal to the river bed elevations. This gave a smooth slide plane, also for complex bed topographies.

The algorithm to find the length, L, of the slide was to first find the lowest and highest cell in the slide. Then the direction vector between the two cells was computed, and decomposed in the two directions of the grid. This meant the horizontal slide length in the lateral and longitudinal direction was found. The smaller length was then chosen as L. This meant for most of the slides that the lateral length was used.

Another concern was to find the maximum slide height, H. In the current study, a midpoint was selected, between the highest and lowest slide cell. The slide depth in this point was taken to be the maximum value for the slide. Due to an irregular topography and slide geometry, this point may not always be located correctly in the middle of the slide. Therefore, slides might emerge that had lower or higher magnitude than what was given from Eq. 13.

The procedure determined the slide depth and therefore the initial geometry of each slide. A grid was then made around each slide. This was embedded in the grid for the groundwater level simulations and the water/sediment grid for water flow and sediment transport computations, as shown in Fig. 3.

The final algorithm calculated the actual movement of the slides. The computer program used in the current study was originally developed to compute the water flow and sediment transport/erosion in the reservoir. This meant that the Reynolds-averaged Navier–Stokes equations were solved to find the water velocities. In the current study, this Navier–Stokes solver was also used to compute the slide movement, assuming the soil in the slide was a very viscous fluid. Olsen [29] had earlier made algorithms for computing the free surface in combination with the Navier–Stokes equations. The CGA method (Continuity, Gravity and Adaptive grid) was modified so that it could be used on the interface between the top of the slide and the water/air. This procedure then computed the changes in the bed elevation due to the geotechnical slides.

The Bodendorf reservoir

The numerical model was tested on a data set from the Bodendorf reservoir, located at the river Mur in the province of Styria, Austria. The reservoir is part of the Bodendorf hydropower plant, and has a length of about 2.5 km and an average width of around 40 m [35]. This reservoir experience considerable sedimentation problems. The initial volume was 900,000 m3 at the plants commission in 1991, but it was reduced to 300,000 m3 in only 12 operation years. The reservoir has been flushed several times since 1994. The flushing is invoked by drawing down the water level at the dam, causing erosion of the sediments as a result of increasing water velocity and bed shear stresses. A time series of the downstream water level during the flushing in 2004 is given in Fig. 8, together with the occurring outflow discharge. In 2004 a detailed survey programme was carried out in the reservoir by the Verbund AG and the Technical University of Graz. Bed elevation measurements were taken before and after the flushing. A sonar instrument was used to take around 130 cross-sectional profiles before the flushing. These bed elevation points are shown in Fig. 9. The data gave a detailed map of the reservoir bed, which was used as geometry input for the numerical model. Eleven cross-section profiles were also taken after the flushing, so that the bed elevation changes could be computed. During the 2004 draw-down flushing 47,300 m3 of sediments were flushed out [36]. The profiles are shown in Fig. 10. One of the profiles is given in Fig. 11, including the bed levels before and after the draw-down flushing. The bed elevation from 1980, before the construction of the dam, is also shown. This level indicates a fixed level under which erosion or slides can not occur. Slides are seen on both sides of the reservoir. On the right side of the reservoir, the bed level closest to the bank is higher in May 2004 (before the flushing) than in June 2004 (after the flushing). The lowering of the bed level can be due to a slide, as the bed level has risen further to the left after the flushing.

Fig. 8
figure 8

Time series of downstream water level and outflow discharge during the flushing event and used for the simulations [4], modified)

Full size image
Fig. 9
figure 9

Plan view of the Bodendorf reservoir, showing the bathymetric survey before the flushing in 2004. Different colours represent bed elevations

Full size image
Fig. 10
figure 10

Plan view of the Bodendorf reservoir. The numbers indicate the cross-sectional profiles where bed level elevations were measured before and after the flushing

Full size image
Fig. 11
figure 11

Cross-section 374.6 with measured bed elevations before and after the reservoir flushing in 2004, together with the bed level from 1980 before the dam was built

Full size image

Another indication of slides at the banks is observations from the reservoir. Figures 12 and 13 shows photographs of the banks during the flushing. Slides can be seen along most of the bank where there are sediment deposits.

Fig. 12
figure 12

Photograph of a sediment slide at the bank of the Bodendorf reservoir

Full size image
Fig. 13
figure 13

Photograph of a bank of the Bodendorf reservoir during flushing

Full size image

Haun et al. [4] computed bed elevation changes during the 2004 flushing of this reservoir without taking the sediment slides at the bank into account. Reasonable correspondence between measured and computed bed elevation changes were found in the middle of the reservoir. However, the numerical model was not able to replicate fairly large lowering of the bed levels along parts of the banks. This was thought to be due to geotechnical slides, which is the motivation of the current study.

Results from the numerical model

The previously described numerical model was used to compute the reservoir flushing of the Bodendorf reservoir. The simulations were similar to the work conducted by Haun et al. [4], except that newly developed algorithms were used, taking geotechnical bank slides into account. A time series of outflow discharges and downstream water levels, shown in Fig. 8, was used as input for the numerical model. A time step of 10 s and 13,300 time steps were used to model the flushing event which lasted 37 h. Bed elevations taken by sonar before the flushing were used to make the 2D and 3D grid for the computation Fig. 9. Haun et al. [4] originally used a grid with 10 cells in the lateral direction, which was too coarse to resolve the slides. Therefore, a new grid with 80 cells in the lateral direction was made. The grid cell size was then between 0.7 and 1 m. The aspect ratio of the 2D cells was not changed, so the number of cells in the longitudinal direction was 1,560. This 2D grid was used to compute the water level in the reservoir, the groundwater level and the location of the slides. The water velocities, pressure, turbulence and sediment concentrations in the reservoir were computed on a 3D unstructured grid, made on the basis of the 2D grid. The Reynolds-Averaged Navier–Stokes equations were solved, together with the k-epsilon turbulence model for each time step. Convection–diffusion equations were used to compute the sediment concentrations of 10 size fractions, together with empirical formulae for the sediment transport at the bed. The results were erosion pattern in the reservoir for each time step of the flushing. The new geotechnical slide module was then invoked, starting with computing the groundwater level. As the water level in the reservoir was drawn down, the banks of the reservoir dried up. A lateral groundwater gradient thereby occurred, which was an important factor for the slide movement. The numerical algorithms for finding the location of the slides and computing the slides itself were then applied. The sediment slides were computed in each time step, following the previously described procedures. The slides also changed the bed elevations, similarly to the erosion due to the water velocities. This resulted in lower bed levels at the higher levels of the slides and an increase of the bed elevation where the movement of the slides stopped. The bed evolution due to slides was then computed over time. The bed levels were adjusted according to the new geometry before the grids were regenerated and a new time step was computed. A wetting/drying algorithm was used to change the 3D grid over time, due to changes in the bed level and the water surface level. The wetting/drying algorithms are further described by Olsen [26].

A question was how the parameters in the geotechnical model could be determined: E-modulus, G-modulus, cohesion, angle of repose for the soil and its viscosity. Data for geotechnical parameters were not available for the Bodendorf site, as investigations of slides was not considered during the monitoring in 2004. Therefore it was decided to use data from another study computing a geotechnical slide in the laboratory. Jia et al. [37] carried out a very large scale laboratory experiment of a bank slide in a water tank. The tank was 15 m long, 5 m wide and 6 m high. An idealized geometry was built in the tank, with a 3.8 m high bank, similar to what can be seen in Fig. 4. In the experiment, the water level was initially at the top of the bank and then it was gradually drawn down over time. The experiment is therefore very similar to the situation described in the current study. Since the current study lacked measurements of basic geotechnical parameters, the parameters from Jia et al. [37] were used. They found the angle of repose, ϕ, (Eq. 3) for the soil material to be 30°. The cohesion of the soil was measured to 1,000 Pa (Eq. 3). The ground water conductivity of the soil, K, (Eq. 2) was determined to be 4.2 × 10–6 m/s.

Another parameter for the slide movement was the viscosity of the moving soil. This value was found numerically by modelling the experiment by Jia et al. [37] and varying the viscosity until the measured slide geometry emerged. This case had previously been modelled by Olsen [33], so all input data was available for this test. The optimum viscosity was found to be 150 kPa s.

The current study focused on the geotechnical slides, so only these results are shown in the following. Movement of the bed due to slides and erosion were stored in two separate arrays, so that it was possible to present just slide movements. Figure 14 shows a plan view of the bed elevation changes due to the slides. Since the reservoir is fairly long compared to the width, it is difficult to see details of the slide in this figure. An enlargement of Fig. 14 is therefore shown in Fig. 15, where the reservoir is divided into four parts (compare Fig. 14), to see more details. The slides are displayed in red and blue in Figs. 14 and 15. The higher part of a bank will be lowered due to the slide (blue colour), and the lower part of the bank will rise (red colour). Figures 14 and 15 show that the slides occur along the banks and the lowering (blue) takes place closest to the side of the reservoir. The corresponding bed level rise (red) occurs further away from the sides. This is a result of the slip circle movement, meaning that bank failures at higher banks will result in a rise of bed elevations closer to the centre of the reservoir. The bed level changes due to the slides are partly connected, forming lines along the reservoir banks. This corresponds with observation from the field (Fig. 13).

Fig. 14
figure 14

Plan view of the reservoir, showing computed vertical bed elevation changes in meters due to the slide movements. Water is flowing in from the left and out to the right. Figure 15 shows more details

Full size image
Fig. 15
figure 15

Detailed plan view of four parts of the reservoir, where the colours show the computed bed elevation changes in meters due to the slide movements. The four sections are numbered in the flow direction, where 1 is the most upstream located part and 4 is the most downstream

Full size image

The slides can also be visualized by showing cross-sections. Figure 16 shows cross-section 373.6 (location shown in Fig. 10) at different times during the flushing. The figure shows the different grids, with the groundwater grid at the bottom and the water/sediment grid at the top. The slide grid is located in between these grids. Slides form initially at both sides of the cross-section, and they move the soil towards the centre of the channel with different speed. The slides emerge and disappear over time. Figure 16 only shows the slides at some selected times. The intermittent nature of the slides is probably also taking place in the prototype. It was observed in the laboratory experiment of Jia et al. [37].

Fig. 16
figure 16

Cross section 373.6 Fig. 10 showing the soil and water grid at different times (in seconds) after the start of the flushing. The upper legend is the water velocity in m/s. The lower legend is the soil velocity in mm/s

Full size image

In order to numerically compare measured and computed bed elevation changes due to slides, a search algorithm was used to find this computed value in the grid at the end of the cross-section profiles shown in Fig. 10. The minimum value of the 10 cells closest to the edge was used. The average computed value of the two edges of the 11 cross-section was found to be 0.4 m. The measured average lowering of the end cross-sections of Fig. 10 was 0.48 m, which represents a 20% deviation.

Discussion

A number of input parameters will affect the quality of the results from the numerical model. Unknown parameters governing the slides (spatial distribution and the magnitude) are the cohesion, soil E and G modules, angle of repose and the soil viscosity. Due to the absence of measured values for these parameters, a sensitivity analysis was carried out. The results are shown in Table 1 and Fig. 17. The figure shows a plan view of the downstream part of the reservoir with the bed elevation changes due to the slides. The downstream part of the reservoir was selected because the slides were larger there.

Table 1 Computed average bed lowering from the parameter sensitivity tests
Full size table
Fig. 17
figure 17

Plan view of downstream part of the reservoir, where the colours show the vertical slide movement [meters]. The letters indicate results with various parameter changes, as given in Table 1

Full size image

One of the parameters varied in the test was the soil viscosity. The result is given in Table 1 and Fig. 17 as test C. The value of the soil viscosity will affect the velocity of the slide. The soil viscosity was increased by 67%, which gave 3% smaller slides. Comparing the Default computation and test C in Fig. 17 shows that the slides are slightly smaller at some locations. However, the main pattern in the bed elevation changes due to the slides remains the same. The algorithms of the limit equilibrium model will check the stability of the slide and stop it once the safety factor is above 1. The change in the viscosity will therefore not affect the final slide pattern very much.

Parameter test F is a modification of the kL value in Eq. 14 from 0.2 to 0.15. Table 1 shows a 20% decrease in the bed level movements due to the slides. Figure 17 shows that the main pattern of the slides remains the same. However, some slides have been reduced in size and a few are missing. A 20% change means that the value of the parameter has a significant effect on the result. In theory, sediments without cohesion should have a fairly flat slide plane. An increasing cohesion would mean a slide plane closer to the shape of a circular segment. There is room for improvement of the current numerical algorithm computing the depth of the slide.

Parameter test D shows that the values of the E and G modulus of the soil were not very sensitive to the slide magnitude. Parameter test D is a 40% reduction in the E and G modulus from 6.5 × 107 to 4 × 107 Pa. This gave an increase in slide volume of only 5%. Figure 17 shows that the main pattern of the slides remained the same. Looking in detail at Fig. 17, it is possible to see that at some locations the slides are slightly larger, but the difference is not substantial.

Another concern in the current study is the choice that the E and G modulus have the same value. Geotechnical analysis shows that this is not the case. A normal Poisson ratio from measurements of soils can be around 0.4. Using this value in Eq. 12 means that the E modulus is 2.8 times larger than the G modulus. A parameter test was therefore carried out where the E modulus was set to 9.6 × 107 Pa and the G modulus to 3.4 × 107 Pa. This gave a reasonable Poisson ratio of 0.4. The sum of these E and G modulus was then the same as for the Default computation. The test was done by assuming the slide would move in the lateral direction, so that the values could be used as given in Eq. 7. Looking at the Figs. 14, 15 and 17, this seems like a reasonable assumption. This parameter test was called E, and the results are shown in Table 1 and Fig. 17. The bed lowering from this test was 0.41 m, compared to the Default value of 0.4. This is only 3% difference. Looking at Fig. 17, the slides are slightly larger at some locations, but the main pattern in bed elevation changes remains the same. The reason that the results are not so sensitive to the Poisson ratio or the E and G modulus is that the shear/normal forces between the elements in the grid is smaller than the other forces: water pressure gradients, gravity and cohesive shear forces under the cells. The assumption that the G module is the same as the E module is therefore not very critical for the current computation. This also means that Eq. 6 can be used independently of the direction of the slide without major errors being introduced. This again means that the slide does not have to be in the lateral direction of the grid.

The parameter having most effect on the sediment slide magnitude was the soil cohesion. Table 1 shows the effect of reducing or increasing the cohesion in test A and B. Test A is a reduction of the default value of 1000–500 Pa, which gives an increase in slide magnitude of 50%. Figure 17 shows that the main slide pattern is the same, but the magnitudes of the slides are significantly larger. The choice of using 1000 Pa as the default value came from the measurements of Jia et al. [37]. The soil used in this experiment was chosen from typical values along a river bank in China, consisting of 12% sand, 80% silt and less than 5% clay size particles. The average sediment particle diameter, d50, was 0.035 mm. The grain size distributions measured at Bodendorf [38] show a large variation along the reservoir. The d50 at the entrance of the reservoir was 32 mm and in front of the weir the d50 showed a value of 1.5 mm. Even though the sediments close to the banks are most probably finer, it is likely that the Bodendorf sediments are coarser than the material used by Jia et al. [37]. This means that the cohesion value of 1000 Pa reported by Jia most likely is lower at Bodendorf. Using a cohesion of 1000 Pa gave a computed average lowering of the end cross-sections of 0.4 m. Selecting a value of 500 Pa gave 0.56 m. Considering that the measured value was 0.48 m, it is possible that the actual cohesion in the soil at Bodendorf could be somewhere between 500 and 1000 Pa.

Deng et al. [18] reported measured cohesive values in bank sediments of up to 10 000 Pa along the Upper Jingjiang river in China. Parameter test B uses this value, which reduces the vertical slide movement to 3% of the default value, as given in Table 1. Figure 17 also shows that almost all the slides have disappeared with such a high value of the soil cohesion. The lateral transport of sediments from the bank into the river would be reduced considerably if this value would have been present for the Bodendorf reservoir.

Looking at Figs. 14 and 15, the computed slides in the downstream part of the reservoir are larger than in the upstream region. The reason is that the water level is lowered more in the downstream location as the water depth in the reservoir increases in the downstream direction. The groundwater gradient will then cover a larger area of the bank. The measured bed elevation changes due to the slides do not show a similar trend. One reason is that there could be systematic variations in some of the soil parameters. The downstream part of the reservoir has lower velocities under normal operation. The downstream sediment deposition will therefore contain more fine sediments than the upstream deposits. The cohesion of the material is a function of the content of the fine material. It is therefore possible that there is a naturally occurring systematic increase in bank cohesion in the streamwise direction of the reservoir. That could counteract the effect of the lateral groundwater gradients. However, in the current model a constant value for the cohesiveness of the sediment was chosen as no measurements of this parameter were available.

Ideally, a study of soil slides in a reservoir should include more detailed field measurements. The geometry of the current study was given from 130 cross-sections, giving 23 m between each section. Given the currently used grid resolution of around 1 m, the vertical elevation of each grid line was linearly interpolated from a fairly long distance. A side-scan sonar or other remote measuring devices could be used to get a much more detailed reservoir topography before and after the flushing. Not only the reservoir bed itself could be investigated, but also the reservoir banks some meters above the highest water level. Soil samples could be taken for measurements of the geotechnical parameters, especially cohesion. Spatial variation of these values along the reservoir could be obtained using a large number of samples. It would also be useful to know the thickness of the sediments and variation of soil parameters over the sediment depth. The actual sediment magnitude could ideally be obtained by scans of the river/valley bed before the dam was built. An alternative would be to use ground-penetrating sonar to find the thickness of the sediments.

Conclusions and outlook

The suggested numerical model for bank erosion and sediment slide movement gives reasonable results for many of the slides observed in the Bodendorf reservoir during a draw-down flushing, both for the location of the slides and the magnitude of the bed elevation changes. The predicted slide movements coincide with steep lateral banks where the largest groundwater gradients occur. Including a groundwater model in the prediction of the slides was therefore found to be an important factor for obtaining reliable results.

The parameter tests show that the sediment cohesion is the most important geotechnical parameter in the determination of the slide magnitude. This value is not so easy to obtain, but there are field techniques with which it can be measured. The current study indicates that the E and G modules also have an effect on the results, but they are not as important as the cohesion. The reason is probably that the friction and cohesion forces on the bottom of the slide are more important than the inter-element forces in the depth-averaged grid.

The numerical algorithms have been developed for sediment slides in water reservoirs. However, it is possible that the algorithms would also work on general land slides. The algorithm used to predict the location of the slides could be used separately to identify areas of slide risks, without computing the actual slide movements itself.

Further work should include testing the numerical model on data sets where the geotechnical soil parameters have been obtained from the actual reservoir where the model is applied, instead of using data from soils of other sites. More detailed geometrical data could be obtained by using a side-scan sonar. A ground penetrating sonar could be used to measure the sediment thickness in more detail, where the values of the original bed levels are not well known.