The simulation of sediment transport by three-dimensional (3D) modelling is linked with the question of how accurate such models are. The current paper provides a test where detailed field measurements of velocities and bed elevation changes over a 3-month period from a prototype reservoir are compared with simulation results. The SSIIM software was used to compute water flow and sediment transport in the Iffezheim reservoir. The numerical model solved the Navier–Stokes equations on a 3D unstructured grid with dominantly hexahedral cells. The k-epsilon turbulence model was used, together with the SIMPLE method to find the pressure. The 3D convection–diffusion equation for suspended sediments was solved for nine sediment fractions. The computed velocity pattern showed good correspondence with the measurements. Grid sensitivity tests showed that the main flow features were computed in different grid sizes, but more accurately so with the finest grid. The sediment deposits were reasonably well computed in location and magnitude. A sensitivity test revealed that the computed bed elevation changes were most sensitive to the fall velocities of the finest cohesive sediment particles and sediment cohesion. The sediment deposition computations were also to some degree sensitive to the sediment discharge formula, bed roughness, discretization scheme and the grid resolution.

INTRODUCTION

Sedimentation is today one of the main problems for the sustainable use of zero-emission hydropower worldwide. A hydropower plant will need a reservoir to store water or to make the intake work properly. In a reservoir, the water velocity and turbulence decrease. The sediments that enter the reservoir will therefore settle and decrease the reservoir volume over time. It is estimated that 1% of the volume of the water reservoirs in the world is lost annually due to sedimentation (Mahmood 1987). Walling & Fang (2003) identified reservoir construction as probably the most important influence on land–ocean sediment fluxes. Suspended sediments may also enter the hydropower intake, causing wear on expensive machinery. The hydropower dam and intake should therefore be designed to reduce the sedimentation problem as much as possible. Traditionally, this design has been aided by the use of physical model studies. However, such studies will not be able to scale suspended load and bed load simultaneously. A promising method for the prediction of the sediment deposits in a reservoir is the simulation by three-dimensional (3D) models using computational fluid mechanics. The question is how accurate such models are. This needs to be tested, and the best data available are from the field. The current paper provides such a test, where detailed measurements of velocities and bed elevation changes from a prototype reservoir are used.

3D numerical models can be differentiated in terms of whether they compute hydrodynamics only or model sediment transport as well. Hydrodynamic models which compute the flow field only have earlier been extensively used in river and hydraulic engineering. Examples include Conway et al. (2013), who used a 3D numerical model to predict the stage–discharge relationship for a two-stage channel, showing that the model could be used in a predictive manner. Applying computational fluid dynamics (CFD) to compute spillway capacity is becoming increasingly popular. Andersson et al. (2013) computed the flow over a spillway for a hydropower plant, and compared this with measurements from a physical model. The velocity profiles were very well reproduced, and the accuracy of the discharge capacity was predicted with an error of between 1 and 8% for different spillway openings. CFD models have also been used to compute more complex flow fields. Baranya et al. (2012) computed the flow around a bridge pier including recirculation zones and secondary currents. The results were verified with measurements from a physical model. Both the calculated velocity pattern and the values of the turbulent kinetic energy compared well with the laboratory study.

Considerable work has been done on computing hydrodynamics combined with sediment transport by using CFD models. Villaret et al. (2013) computed sediment transport in laboratory channels and the river Danube. The results were compared with field measurements of bed load discharge at one location in the river. Fischer-Antze et al. (2008) computed bed elevation changes during a 100 year flood in a reach of the Danube river. The results compared reasonably well with measurements of the bed elevation changes obtained from echosoundings. Sediment transport has also earlier been computed in hydropower reservoirs. Souza et al. (2010) used a depth-averaged model to compute sediment deposition in a shallow reservoir. Harb et al. (2014) computed deposition of plastic particles in a physical model of a hydropower reservoir in Austria. Good agreement was found for the deposition pattern. Jia et al. (2013) computed flow through the reservoir of the Three Gorges Dam, including effects of increased fluid density due to high sediment concentrations. The results were verified with velocity profiles measured in the field. Ruether et al. (2005) computed water velocity and suspended sediment transport in the Kapunga water intake in Tanzania using a 3D Navier–Stokes solver. The computed trap efficiencies of the reservoir were compared with measurements from a physical model. Khosronejad et al. (2008) computed sediment release from a reservoir and compared the bed elevation changes with results from a physical model study.

The current study uses a 3D hydrodynamic model combined with sediment transport computations. The purpose of the model was to simulate the complex flow field and bed elevation changes in the Iffezheim hydropower reservoir, located on the Upper Rhine River in Germany. The Iffezheim reservoir exhibits an annual sediment deposition of about 115,000 m3 (WSA Freiburg 2011). The deposited fine sediments contain contaminants (Pohlert et al. 2011), making flushing difficult from an environmental point of view. The sediments are therefore dredged, which is a fairly costly operation for the hydropower company. Due to the high amount of sediment deposition and their contamination, the Iffezheim reservoir was an interesting case for the prediction of sediment transport by 3D modelling. Also, the amount of available data from Iffezheim was substantial compared to other impoundments. The accuracy of the applied model computations was verified by both velocity measurements and bed changes over a period of 105 days for the same reservoir. Several earlier publications have used data from the Iffezheim reservoir. Hillebrand & Olsen (2011) modelled consolidation of sediments in the reservoir. Hillebrand et al. (2012a) and Zhang et al. (2013, 2015) computed the water and sediment transport in the Iffezheim reservoir with a constant time step that was relatively short. One of the main novelties of the current paper is the application of algorithms using long and varying time steps, enabling 3D simulations of sediment deposition over several months with a much shorter computational time than before. The current study also includes sensitivity tests of more parameters than the earlier studies.

An important question using numerical models for water engineering is how much calibration is required (Bennett et al. 2013). If a new hydropower plant is to be constructed, data for the sediment deposition process would not be available until after the project is built. It would be advantageous for the design of the plant to model the sediment deposition as a function of dam height, gate sizes and operation scenarios. A model that would function without calibration would then be needed. The main calibration parameter for solving the Navier–Stokes equations in a river is the bed roughness. In the current study, the bed roughness is computed by an empirical formula as a function of the grain size distribution on the bed. This formula contains an empirical parameter, just like the methods used to compute the sediment transport. None of these parameters has been changed in the current study. If these empirical formulas were more or less correct, the model should in theory be possible to use without calibration. In practice, there are also other parameters affecting the results. The current study investigates the most important of these parameters and their effect on the final result.

STUDY AREA AND DATA

The Iffezheim reservoir is located on the Rhine River near Karlsruhe in Germany. At this location, the Rhine forms the border between Germany and France, and the Iffezheim hydropower plant is owned jointly by French and German power authorities. The Rhine catchment and an aerial view of the Iffezheim reservoir are given in Figure 1. A top view of the reservoir is shown in Figure 2. The reservoir ends in three branches. The right channel contains locks for ships and has usually only a negligible water discharge. The middle channel leads water to the hydropower plant, while the left channel is used for floods and discharging surplus water. The Iffezheim reservoir is 700 m wide at the dam. The annual sediment load at the reservoir is around 1 million tons/year and the average water discharge is 1,240 m3/s (LUBW 2011).

Figure 1

Rhine catchment and aerial view of the Iffezheim reservoir.

Figure 1

Rhine catchment and aerial view of the Iffezheim reservoir.

Figure 2

Top view of the Iffezheim reservoir.

Figure 2

Top view of the Iffezheim reservoir.

The spillway channel is only operated intermittently during high discharges in the river. At normal river discharges the spillway gates are closed, leading to small velocities in this channel. This causes fine sediments to deposit. A method for prediction of the sediment deposits in the reservoir is of interest. The flow field in the reservoir is quite complex, with secondary currents, recirculation zones and unsteady effects as the spillway gates are opened and closed, affecting the amount of sediments depositing in the channel. This is not possible to compute with a 1D model. A 2D model would be able to model some of the complex flow pattern, but a fully 3D model should give the best results. Therefore, a 3D model was used for the current study.

The grain size distribution was analysed by taking core samples. It was shown that the grain sizes in the reservoir range from clay to gravel. The numerical model computed the sediment transport by nine different fractions. The initial grain size distribution on the bed is given in Table 1. The distribution is based on measurements of the main river bed and was initially given uniformly over the whole geometry. The distribution of the inflowing sediment fractions is also based on measurements (Astor et al. 2014) and is also given in Table 1. This was used as the upstream boundary conditions for the numerical model.

Table 1

Sediment characteristics for the different sizes

Size number Diameter [mm] Initial fraction at bed, reference case [%] Fall velocity Winterwerp & van Kesteren (2004) [mm/s] Fall velocity, Zanke (1977), reference case [mm/s] Inflow sediments [%] 
20 88 – 568 
8.4 – 216 
1.3 2.3 – 136 
0.4 1.3 – 56 
0.13 – 11 11 
0.04 – 1.1 19 
0.02 0.027 0.27 35 
0.005 0.0069 0.017 12 
0.002 0.0028 0.0027 21 
Size number Diameter [mm] Initial fraction at bed, reference case [%] Fall velocity Winterwerp & van Kesteren (2004) [mm/s] Fall velocity, Zanke (1977), reference case [mm/s] Inflow sediments [%] 
20 88 – 568 
8.4 – 216 
1.3 2.3 – 136 
0.4 1.3 – 56 
0.13 – 11 11 
0.04 – 1.1 19 
0.02 0.027 0.27 35 
0.005 0.0069 0.017 12 
0.002 0.0028 0.0027 21 

Acoustic Doppler current profiler (ADCP) measurements of water velocities were taken at three cross-sections of the spillway channel, with locations shown in Figure 2. The measurements were taken on 14 May 2012, when the inflow discharge to the reservoir was 1,550 m3/s. This discharge has an exceedance probability of about 25%. The hydropower plant is designed for mean flow conditions, i.e., an exceedance probability of 50%. This means that for discharges of 1,100 m3/s or less, which statistically occur on 183 days per year, there is no discharge through the weir section. Thus, in the current study, the discharge through the power plant channel was 1,100 m3/s whereas 450 m3/s passed through the spillway channel.

The bathymetry of the reservoir was measured with single-beam echosounding equipment at cross-sections with an average distance of 20 m from a boat in April and July 2007. The time series of discharge through the upstream boundary, power plant and weir in this period of 105 days is given in Figure 3. These data were used as boundary conditions for the numerical modelling of bed elevation changes.

Figure 3

Time series (05.04.2007–18.07.2007: 105 days) of water discharges in the reservoir and out of the turbines and spillway, respectively.

Figure 3

Time series (05.04.2007–18.07.2007: 105 days) of water discharges in the reservoir and out of the turbines and spillway, respectively.

METHODS

The numerical model solved the Navier–Stokes equations in three dimensions together with the convection–diffusion equation for sediment concentration. The k-epsilon turbulence model was used (Launder & Sharma 1974). The pressure was computed with the SIMPLE method (Patankar 1980). Wall laws for rough boundaries were used (Schlichting 1979). The wall law predicted the shear velocity and thereby the bed shear stress from the velocity in the bed cell. The bed shear stress was used as a sink term at the bed for the Navier–Stokes equations. A constant wall roughness of 2 cm was used, similar to the maximum observed sediment size. A parameter test was done using the bed roughness formula by van Rijn (1982): 
formula
1
A finite volume method was used for the discretization on a 3D unstructured grid. The numerical model had earlier been validated against measurements in a number of studies, for example, Wilson et al. (2003) and Rüther et al. (2010).

In order to calculate the water flow only, the grids were generated based on the bathymetry from February 2012. The bathymetry from April 2007 was used in generating the grid for the sediment computations. The grids are made up of mostly hexahedral cells, but some pentahedral and tetrahedral cells were used at the side boundaries to facilitate modelling of the complex geometry. A plan view of the downstream part of the coarser grid and a vertical slice through the grid are given in Figure 4. The grid covered the whole of the geometry shown in Figure 2. The grid had 100 × 40 × 10 cells in the streamwise, lateral and vertical direction, respectively. The grid contained approximately 300,000 3D cells with horizontal dimensions between 5 and 10 m. The grid is reasonably uniform throughout, since there are no areas in the model domain which have to be refined in the grid locally. The cell resolution was chosen in a way to achieve a compromise between accuracy and computational efficiency. The finer grid had approximately twice as many cells in each of the two horizontal directions, giving a grid of 1.1 million cells. The grids were generated from bed levels computed from echosoundings of the bathymetry. The water reservoir level then defined the other surface of the grid. A 2D depth-averaged grid was extended to three dimensions by assuming a uniform distribution of the grid in the vertical direction. The number of cells in the vertical direction was varied according to the water depth, with between one and ten cells. The minimum and maximum cell heights were 0.20 and 3.3 m, respectively. The average height for all the cells in the grid was 2.4 m. More details about the grid generation are given by Olsen (2003).

Figure 4

Downstream part of the coarse grid seen in plan view (left). Vertical slice through the grid (right).

Figure 4

Downstream part of the coarse grid seen in plan view (left). Vertical slice through the grid (right).

The current study used an adaptive grid to model the changes in the water and bed levels. The grid was then moved vertically according to the computed free water surface and bed level changes. The free surface was computed using the pressure computed with the SIMPLE method when solving the Navier–Stokes equations. The measured water level at the downstream gates was used as a boundary condition. The method to generate the free surface is further described and verified by Tritthart & Gutknecht (2007) against field measurements in the Danube River.

The sediment concentration, c, of size i, was computed from the convection–diffusion equation for suspended material: 
formula
2
The water velocity is denoted U, the fall velocity of the particle is w, Γ is the turbulent diffusion coefficient and FR,i the pick-up rate. The diffusion coefficient was set equal to the eddy-viscosity from the k-epsilon turbulence model. The reference computations were done with the fall velocities computed by Zanke's formula (Zanke 1977), as given in Table 1. The third term on the left side of the equation is a drift-flux term, taking into account the fall velocity of the particles. An empirical formula for the equilibrium sediment concentration, cR,i, at the bed was used (van Rijn 1984a): 
formula
3
The particle diameter is denoted di, T is a dimensionless shear stress parameter, and D is a dimensionless particle parameter. The parameter a is the distance from the bed to where a reference concentration is computed, which is set equal to half the bed cell height in the current study. T and D are given with the following formulas for size number i: 
formula
4
 
formula
5
where τ is the actual shear stress on the bed, computed from the turbulent kinetic energy, k: 
formula
6
The turbulent kinetic energy is computed from the k-epsilon turbulence model. The cμ parameter is a constant in the k-epsilon model, equal to 0.09.

The critical shear stress, τc, for movement of the bed particles was computed from the Shields (1936) curve, where ρ is the water density and ρs is the sediment density. The acceleration of gravity is denoted g and ν is the viscosity of water. For cohesive particles in the silt and clay fraction, the critical shear stress was increased by adding a constant value.

The resulting sediment pick-up rate of size i was computed from the following equation: 
formula
7
The fall velocity of the sediment particle is denoted w, and f is the fraction of the bed material of size i. The parameter f will vary between 0 and 1. The variable cc is the concentration in the bed cell computed in the previous time step. The limiter on the sediment pick-up flux in Equation (7) ensures that the deposition of sediments is always larger than or equal to the erosion. This makes the results less sensitive to the initial grain size distribution on the bed and the size of the time step.
Bed load was computed in addition to the suspended load according to van Rijns (1984b) formula: 
formula
8
The bed level changes, Δz, were computed based on the computed sediment deposition and the pick-up rate for all the size fractions: 
formula
9
The computed concentration in a bed cell is denoted ci, while w is the fall velocity of the particles, Δt is the time step and r is a conversion factor between volume of sediment particles and volume of deposits on the bed. Measurements from the field showed that the solid fraction of the bed material was 24% on average (Hillebrand 2014). This value was used in the numerical model.
Modelling sediment movements in a 3D grid over several months can take long computational times. A time series of discharges over a longer period often shows periods with low and high discharges. Usually, the main bed changes occur during the high discharges. It is therefore more important to use a short time step during high discharges than during low discharges. The numerical model therefore used a varying time step, Δt, according to the following formula: 
formula
10
The values used in the current study for the reference case run were Δt0= 20,000 seconds, n = 3 and Qref = 1,000 m3/s.

RESULTS AND DISCUSSION

Velocity field

The velocity field was measured and computed for a constant water discharge of 1,550 m3/s, a discharge with an interesting flow pattern in the reservoir. The fine grid was used, with twice as many cells in each horizontal direction compared with the grid in Figure 4. The discharge through the intake channel was 1,100 m3/s, whereas 450 m3/s passed through the spillway channel. The lock channel on the right side was not operated and had zero discharge.

Figure 5 shows the computed velocity field in the reservoir. The water flow is fairly uniform in the central hydropower intake channel, except for a small recirculation zone on the left side of the hydropower intake channel at the divide wall between the left and right channel. There is a large recirculation zone on the left side of the spillway channel. Both recirculation zones are at the inside of a curve or a curved flow field. Similar recirculation zones at the inside of a bend have been observed by Blanckaert et al. (2013) in a physical model study. The recirculation zone covers almost the whole length of the spillway channel. This flow feature is also found in Figure 6, showing velocities measured with the ADCP instrument in three profiles. The profiles' locations are given in Figure 2. The recirculation pattern is shown in profile 2 and 3 by negative velocity vectors at the left side of the profiles. Profile 1 is located next to the weir, where the computed recirculation zone is very weakly established. This is reflected by the near absence of negative velocity vectors in this profile. The measurements in the recirculation zones show a particularly irregular pattern. This could be due to unsteady vortexes in the recirculation zones. In some parts of the channel, the velocity vectors seem more uniform. This is where the numerical model also predicts such a velocity pattern. The main flow pattern is computed correctly, in that the highest velocity is on the right side and follows the main flow direction. Both measured and computed velocities are lower on the left side of the channel, in the recirculation zone. Profile 1 shows the largest deviation between the measured and computed velocity profiles. In this profile, the numerical model shows a slightly different direction of the main flow and higher velocities at the right bank than in the measurements. This could be because the profile is located close to the weir, and the velocities are thereby affected by which gates were open at the time of the measurements. The numerical model assumed an outflow over the whole cross-section.

Figure 5

Water velocity vectors at the surface computed with the fine grid. The streamtraces in bold black indicate the recirculation zones on the left bank in the weir channel and on the left side in the hydropower intake channel (enlarged view).

Figure 5

Water velocity vectors at the surface computed with the fine grid. The streamtraces in bold black indicate the recirculation zones on the left bank in the weir channel and on the left side in the hydropower intake channel (enlarged view).

Figure 6

Measured (above) and computed (below) velocity vectors at profiles 1, 2 and 3.

Figure 6

Measured (above) and computed (below) velocity vectors at profiles 1, 2 and 3.

Bed elevation changes

The bed elevation changes over a continuous time period of 105 days from April to July 2007 were computed by the numerical model. The time period was selected based on dredging records. The reservoir is often dredged, but in this time period no dredging took place. Dredging is not recorded in detail, and is therefore difficult to take into account when computing 3D changes in the bed geometry. The time series of discharges from these 105 days are given in Figure 3.

The inflow sediment concentrations were based on a rating curve from records of measured suspended sediments. A total of 0.6 million tons of sediments were estimated to flow into the reservoir during the 105 days of the simulation. The spatial variation of the upstream concentrations was determined from measurements using an ADPC instrument. The procedure was calibrated by taking bottle samples (Hillebrand et al. 2012b). The computations were done with both a coarse (Figure 4) and a fine grid. The fine grid had approximately twice as many cells as the coarse grid in each of the two horizontal directions. To save computational time, most of the parameter tests were done on the coarse grid.

Bed elevations were measured by echosounding equipment in April and July 2007. The resulting measured bed elevation changes are shown in Figure 7(left), together with the computed sediment deposits in Figure 7(right), for the reference case with the coarse grid. The majority of sediments deposit in the spillway channel, both in the numerical model and in the measurements. Figure 7(right) also shows some computed deposition in the central hydropower intake channel. This channel was not measured with the echosounding equipment, and this is why Figure 7(left) shows no deposition there.

Figure 7

Measured (left) and computed (right) sediment deposition pattern. The greyscale shows the bed elevation change. Only region A was measured in the field. Region B was used in the comparison of measured vs. computed sediment deposition volume.

Figure 7

Measured (left) and computed (right) sediment deposition pattern. The greyscale shows the bed elevation change. Only region A was measured in the field. Region B was used in the comparison of measured vs. computed sediment deposition volume.

Sensitivity analyses

The numerical model contains a number of input parameters that will affect the final result. In order to evaluate the uncertainty of these parameters on the model results, sensitivity analyses are essential and a prerequisite for further sediment transport computations in river impoundments. In the current study, the simulation sensitivity to each chosen parameter was analysed by comparison of the computed total amount of deposited sediments in the spillway channel. The deposition amount from the field measurements was estimated at 51,000 m3 over a period of 3 months. The uncertainty of the field measurements is estimated to be up to 20%, i.e., approximately 10,000 m3. Therefore, the computed sediment volumes are expected to be in the range between 41,000 m3 and 61,000 m3. The chosen sensitivity parameters and the resulting volumes from the numerical model are given in Table 3, where also the percentage of deviation between computations and measurements are shown.

The numerous numbers of input parameters which will affect the sediment transport simulation can be distinguished in three groups. One group is numerical parameters which influence the stability and convergence of the computation. This includes, for example, the time step, grid resolution, discretization scheme and the thickness of the active sediment layer. Another group is physical parameters which are used to analyse their effect on sediment transport. This would, for example, be the initial grain size distribution, roughness and the applied sediment transport formula. When calculating bed elevation changes caused by deposition of fine, cohesive sediments, sediment-specific parameters also have to be taken into account. In this context, flocculation is an important sediment process, i.e., the aggregation and disaggregation of single particles and larger aggregates. Input parameters which are implicitly correlated to flocculation are the fall velocities of the finest particle sizes. The cohesion of the sediments is also an essential factor which will affect the simulation results in terms of its influence on the critical shear stress of the particles.

On the basis of these considerations, the above-mentioned sensitivity parameters were chosen in this study. As a result, a total of nine parameters were investigated by applying physical feasible ranges for each parameter. The calculated deposited amount of each sensitivity simulation was compared to the simulation results of a reference case computation and to the amount of deposits calculated from measurements (given in Table 3). The reference case computations were done on the coarser gird to save computational time. The definitions of the input parameters for the reference case are given in Table 2.

Table 2

Settings of the input parameters for the reference case

Numerical parameters 
 Grid resolution Coarser grid (100 × 40 × 10) 
 Time step Varying time step with an average time step of 10,000 seconds 
 Discretization scheme Second-order scheme 
 Active layer thickness 10 cm 
Physical parameters 
 Initial grain size distribution See Table 1 
 Roughness 2 cm 
 Bed load transport formula van Rijn (1984b)  
Parameters correlated to cohesive sediments 
 Fall velocity Zanke (1977), see Table 1 
 Critical shear stress for bed particles Shields (1936)  
Numerical parameters 
 Grid resolution Coarser grid (100 × 40 × 10) 
 Time step Varying time step with an average time step of 10,000 seconds 
 Discretization scheme Second-order scheme 
 Active layer thickness 10 cm 
Physical parameters 
 Initial grain size distribution See Table 1 
 Roughness 2 cm 
 Bed load transport formula van Rijn (1984b)  
Parameters correlated to cohesive sediments 
 Fall velocity Zanke (1977), see Table 1 
 Critical shear stress for bed particles Shields (1936)  
Table 3

Sediment deposits in intake channel

Method Deposits [m³] Deviation, reference [%] Deviation, measured [%] 
Measured 51,000 −10 0.0 
Computed, reference case 56,850 0.0 11 
Computed, 300 seconds fixed time step (30,000 time steps) 54,721 −3.7 7.3 
Computed, fine grid and 300 seconds fixed time step 50,031 −14 −1.9 
Computed, with initial bed grain size distribution from reference case 56,239 −1.1 10 
Computed, fall velocities from Winterwerp formula for 3 finest fractions 45,305 −20 −13 
Computed, 5 × larger fall velocities for 3 finest fractions 83,961 48 65 
Computed, roughness from 2 to 5 cm 47,769 –16 −6.7 
Computed, roughness from 2 cm to a value calculated using the van Rijn formula in equation (1) 51,947 –8.6 1.8 
Computed, active layer thickness from 0.1 to 0.2 m 57,635 1.4 13 
Computed, active layer thickness from 0.1 to 0.02 m 57,144 0.5 12 
Computed, Engelund–Hansen formula 46,140 −19 −11 
Computed, van Rijn suspended load formula only 60,458 6.3 19 
Computed, first-order scheme instead of second-order scheme 59,616 4.9 17 
Computed, cohesion equivalent to 1 Pa on 4 finest fractions 70,395 24 38 
Computed, cohesion equivalent to 0.1 Pa on 4 finest fractions 67,022 18 31 
Computed, cohesion equivalent to 0.01 Pa on 4 finest fractions 58,264 2.5 14 
Method Deposits [m³] Deviation, reference [%] Deviation, measured [%] 
Measured 51,000 −10 0.0 
Computed, reference case 56,850 0.0 11 
Computed, 300 seconds fixed time step (30,000 time steps) 54,721 −3.7 7.3 
Computed, fine grid and 300 seconds fixed time step 50,031 −14 −1.9 
Computed, with initial bed grain size distribution from reference case 56,239 −1.1 10 
Computed, fall velocities from Winterwerp formula for 3 finest fractions 45,305 −20 −13 
Computed, 5 × larger fall velocities for 3 finest fractions 83,961 48 65 
Computed, roughness from 2 to 5 cm 47,769 –16 −6.7 
Computed, roughness from 2 cm to a value calculated using the van Rijn formula in equation (1) 51,947 –8.6 1.8 
Computed, active layer thickness from 0.1 to 0.2 m 57,635 1.4 13 
Computed, active layer thickness from 0.1 to 0.02 m 57,144 0.5 12 
Computed, Engelund–Hansen formula 46,140 −19 −11 
Computed, van Rijn suspended load formula only 60,458 6.3 19 
Computed, first-order scheme instead of second-order scheme 59,616 4.9 17 
Computed, cohesion equivalent to 1 Pa on 4 finest fractions 70,395 24 38 
Computed, cohesion equivalent to 0.1 Pa on 4 finest fractions 67,022 18 31 
Computed, cohesion equivalent to 0.01 Pa on 4 finest fractions 58,264 2.5 14 

The reference case computation estimated the deposition to be 11% too high. This computation was done with a varying time step, according to Equation (10), with an average time step of 10,000 seconds. The numerical model used an implicit solver, enabling long time steps without causing stability problems. However, the accuracy of the results was a function of the magnitude of the time step. Therefore, this parameter was included in the parameter sensitivity test. In addition to using Equation (10), a constant time step was tested. Reducing the time step to 300 seconds gave 3.7% less deposited volume compared with the reference run. The results were therefore not very sensitive to the time step. Using a constant time step of 300 seconds, the computations were also done on the finer grid. The finer grid had twice as many cells in each of the horizontal directions. This gave a 14% underprediction of the deposited volume compared with the reference run, or 10% lower deposited volume compared with the coarse grid. This showed that the results were not completely grid-independent. However, the difference is not very large, considering that there are around four times more cells in the fine grid compared with the coarse grid.

The next numerical parameter to test was the discretization scheme. Using a first-order upwind scheme instead of the second-order scheme increased the deposited volume by 4.9%. The first-order scheme introduces false diffusion, which is especially pronounced for a coarse grid. Another parameter is the thickness of the active sediment layer. The current numerical algorithms use only two layers: an inactive and an active layer. The thickness of the active layer was initially set to 10 cm, a factor five times larger than the largest sediment size fraction modelled. Table 3 shows that varying the value to 2 cm or 20 cm only changes the deposited volume by under 1.4%. In other words, the size of the active layer was not very important in the current study.

The physical parameters were analysed by varying the spatial variation of the initial grain size distribution. This sensitivity was tested by saving the grain size distribution at the end of the reference run and using it as a start-up value for a new run. The result was only 1.1% less deposition than the reference run. The result was therefore not very sensitive to the initial grain size distribution. This also makes sense from a physical point of view, where deposition dominates over erosion processes, as there were no large flood events in the chosen time period between April and July 2007. The dominance of deposition is confirmed by the numerical model.

A physical parameter that is often used to calibrate numerical models is the roughness. The initial value used in the current study was 2 cm, similar to the largest size of the sediment particles. Table 3 shows that using a roughness of 5 cm or using Equation (1) decreased the deposited sediment volume by 16% and 8.6%, respectively. Equation (1) also increases the roughness on the bed, causing a higher bed shear stress. The increased bed shear stress causes a higher pick-up rate from Equation (3), reducing the sediment deposition. The roughness is an important parameter for the results.

The sediment transport formula used in the current study was developed by van Rijn for CFD modelling. An alternative formula was tested in the current study: the Engelund–Hansen formula (Engelund & Hansen 1967). Table 3 shows that the resulting deposited sediment volume decreased by 19% using the Engelund–Hansen formula. Removing van Rijn's bedload formula (Equation (8)) in the computation increased the deposited volume by 6.3%.

An important sediment process is flocculation. The finer sediments stick together to form coarser sizes. This mainly affects the fall velocity of the particles. There exist very complex models for flocculation, but in the current study the process was only investigated by varying the fall velocity of the three finest sizes. Table 3 shows two variations. Instead of using the Zanke (1977) formula, fall velocities of the three finest sizes were computed using the Winterwerp & van Kesteren (2004) formula. This gave 20% less deposited sediment volume. The other test was to increase the fall velocities of the three finest fractions by a factor of 5. This had a very large effect on the results, causing 48% increase in the deposited volume. Of all the parameters tested in this study, the fall velocity of the particles had the largest effects on the results.

The cohesion of the sediments was also tested by increasing the critical shear stress of the particles. Table 3 shows an increase in the deposited volume from 2.5 to 24% for the highest increased critical shear stress of 1 Pa. Field measurements show that the increase in critical shear stress due to cohesion in the Iffezheim reservoir generally is in the order of 1 Pa, but can be up to 5 Pa (Hillebrand 2014). The cohesion of the sediments is therefore also a very important parameter for the results.

The sensitivity studies above show that the computed deposition volume, which was expected to range between 41,000 m3 and 61,000 m3, was not underestimated in any of the tests. Only in the case of larger fall velocities for the three finest fractions, a significant sediment deposition of 83,961 m3 occurred. This sediment amount corresponds to an overestimation of 65% compared with the measured deposition volume.

CONCLUSIONS

The numerical model is able to compute the velocity field in the Iffezheim reservoir fairly well. The main recirculation zone on the left side of the spillway channel has a reasonable size compared with measurements. The computed bed elevation changes occur in the spillway channel, which is also observed in the measurements. A parameter sensitivity test shows that the resulting bed elevation changes are most sensitive to the flocculation process and fall velocity of the sediments (48%). It is also fairly sensitive to the cohesion of the bed material (24%) and the sediment transport formula (19%). A fair sensitivity to the roughness (16%) and the grid size (10%) was also observed. The results were less sensitive to the discretization scheme (5%) and the time step (4%).

It can be concluded that the 3D model is sufficiently accurate to compute the complex flow field and sediment deposition in the Iffezheim reservoir. However, the sensitivity analyses also showed that there are uncertainties in terms of the flocculation process. Hence, further investigations include detailed in situ measurements of the fine, cohesive sediments. Measured data on floc settling velocities or floc sizes can thereby serve as input data for numerical models applied for future studies on sediment transport in the Iffezheim reservoir.

REFERENCES

REFERENCES
Andersson
A. G.
Andreasson
P.
Lundström
T. S.
2013
CFD-modelling and validation of free surface flow during spilling of reservoir in down-scale model
.
Engineering Applications of Computational Fluid Mechanics
7
(
1
),
159
167
.
Astor
B.
Gehres
N.
Hillebrand
G.
2014
From Source to Mouth, a Sediment Budget of the Rhine River: Grain Size Distribution of Suspended Sediment Samples in the Rhine and its Tributaries
.
BfG-1798. Bundesanstalt für Gewässerkunde
,
Koblenz
(In German).
Baranya
S.
Olsen
N. R. B.
Stoesser
T.
Sturm
T.
2012
Three-dimensional RANS modeling of flow around circular piers using nested grids
.
Engineering Applications of Computational Fluid Mechanics
6
(
4
),
648
662
.
Bennett
N. D.
Croke
B. F. W.
Guariso
G.
Guillaume
J. H. A.
Hamilton
S. H.
Jakeman
A. J.
Marsili-Libelli
S.
Newham
L. T. H.
Norton
J. P.
Perrin
C.
Pierce
S. A.
Robson
B.
Seppelt
R.
Voinov
A. A.
Fath
B. D.
Andreassian
V.
2013
Characterising performance of environmental models
.
Environmental Modelling & Software
40
,
1
20
.
Blanckaert
K.
Kleinhans
M. G.
McLelland
S. J.
Uijttewaal
W. S. J.
Murphy
B. J.
de Kruijs
A.
Parsons
D. R.
Chen
Q.
2013
Flow separation at the inner (convex) and outer (concave) banks of constant-width and widening open-channel bends
.
Earth Surface Processes and Landforms
38
(
7
),
696
716
.
Conway
P.
O'Sullivan
J. J.
Lambert
M. F.
2013
Stage-discharge prediction in straight compound channels using 3D numerical models
.
Proceedings of the Institution of Civil Engineers – Water Management
166
(
1
),
3
15
.
Engelund
F.
Hansen
E.
1967
A Monograph on Sediment Transport in Alluvial Streams
.
Teknisk Forlag
,
Copenhagen, Denmark
.
Fischer-Antze
T.
Olsen
N. R. B.
Gutknecht
D.
2008
Three-dimensional CFD modeling of morphological bed changes in the Danube River
.
Water Resources Research
44
(
9
),
W09422
.
Harb
G.
Haun
S.
Schneider
J.
Olsen
N. R. B.
2014
Numerical analysis of synthetic granulate deposition in a physical model study
.
International Journal of Sediment Research
29
(
1
),
110
117
.
Hillebrand
G.
2014
Sediment Analysis of Drill Core Samples from the IKSR Investigation of the Reservoirs in the Upper Rhine from the Years 2000 to 2002
.
KLIWAS-Report series, KLIWAS-56/2014
,
Bundesanstalt für Gewässerkunde
,
Koblenz
,
63
pp.
(In German).
Hillebrand
G.
Olsen
N. R. B.
2011
Towards Modeling Consolidation of fine Sediments Upstream of the Iffezheim Barrage, Upper Rhine River, Germany
.
River, Coastal and Estuarine Morphodynamics (RCEM)
,
Hamburg
,
Germany
.
Hillebrand
G.
Klassen
I.
Olsen
N. R. B.
Vollmer
S.
2012a
Modelling fractionated sediment transport and deposition in the Iffezheim reservoir
. In:
10th International Conference on Hydroinformatics (HIC)
,
Hamburg
,
Germany
.
Hillebrand
G.
Otto
W.
Vollmer
S.
2012b
Findings from ADCP-measured flow velocities and suspended sediment concentrations at the Upper Rhine
. In:
2nd IAHR Europe Conference
,
Munich
,
Germany
.
Khosronejad
A.
Rennie
C. D.
Neyshabouri
A. A. S.
Gholami
I.
2008
Three-dimensional numerical modeling of reservoir sediment release
.
Journal of Hydraulic Research
46
(
2
),
209
223
.
LUBW Landesanstalt für Umwelt, Messungen und Naturschutz Baden-Württemberg
(ed.).
2011
German Hydrological Yearbook, Rhine Area, Part 1 – High and Upper Rhine (In German).
Mahmood
K.
1987
Reservoir Sedimentation: Impact, Extent and Mitigation
.
World Bank Technical Paper 71
,
Washington, DC
,
USA
.
Olsen
N. R. B.
2003
3D CFD modeling of a self-forming meandering channel
.
Journal of Hydraulic Engineering
129
(
5
),
366
372
.
Patankar
S. V.
1980
Numerical Heat Transfer and Fluid Flow
.
McGraw-Hill Book Company
,
New York
,
USA
.
Pohlert
T.
Hillebrand
G.
Breitung
V.
2011
Trends of persistent organic pollutants in the suspended matter of the River Rhine
.
Hydrological Processes
25
,
3803
3817
.
Ruether
N.
Singh
J. M.
Olsen
N. R. B.
Atkinson
E.
2005
Three-dimensional modelling of sediment transport at water intakes
.
Proceedings of the Institution of Civil Engineers, UK, Water Management
158
,
1
7
.
Schlichting
H.
1979
Boundary Layer Theory
.
McGraw-Hill Book Company
,
New York
,
USA
.
Shields
A.
1936
Use of Dimensional Analysis and Turbulence Research for Sediment Transport, Preussen Research Laboratory for Water and Marine Constructions, Publication no. 26
,
Berlin
,
Germany
(In German).
Souza
L. B. S.
Schulz
H. E.
Villela
S. M.
Gulliver
J. S.
Souza
L. B. S.
2010
Experimental study and numerical simulation of sediment transport in a shallow reservoir
.
Journal of Applied Fluid Mechanics
3
(
2
),
9
21
.
Tritthart
M.
Gutknecht
D.
2007
3-D Computation of flood processes in sharp river bends
.
Proceedings of the Institution of Civil Engineers – Water Management
160
(
4
),
233
247
.
van Rijn
L. C.
1982
Equivalent roughness of alluvial bed
.
Journal of Hydraulic Engineering
108
(
10
),
1215
1218
.
van Rijn
L. C.
1984a
Sediment transport. Part I: bed load transport
.
Journal of Hydraulic Engineering
110
(
10
),
1431
1456
.
van Rijn
L. C.
1984b
Sediment transport. Part II: suspended load transport
.
Journal of Hydraulic Engineering
110
(
11
),
1613
1641
.
Villaret
C.
Hervouet
J.-M.
Kopmann
R.
Merkel
U.
Davies
A. G.
2013
Morphodynamic modelling using the Telemac finite-element system
.
Computers and Geosciences
53
,
105
113
.
Walling
D. E.
Fang
D.
2003
Recent trends in the suspended sediment loads of the world's rivers
.
Global and Planetary Change
39
(
1–2
),
111
126
.
Wasser- und Schifffahrtsamt (WSA) Freiburg
2011
Sachstandsbericht oberer Wehrkanal Staustufe Iffezheim
[Technical report – Upper weir channel of the Iffezheim hydropower reservoir].
Wilson
C. A. M. E.
Boxall
J. B.
Guymer
I.
Olsen
N. R. B.
2003
Validation of a 3D numerical code in the simulation of pseudo-natural meandering flows
.
Journal of Hydraulic Engineering
129
(
10
),
758
768
.
Winterwerp
J. C.
van Kesteren
W. G. M.
2004
Introduction to the Physics of Cohesive Sediment in the Marine Environment. Developments in Sedimentology 56
.
Elsevier
,
Amsterdam
,
The Netherlands
.
Zanke
U.
1977
Berechnung der Sinkgeschwindigkeiten von Sedimenten. Mitteilungen des Franzius-Instituts für Wasserbau und Küsteningenieurwesen der TU Hannover, Vol. 46
[Calculation of settling velocities of sediments. Announcements of the Franzius-Institute for Hydraulic, Estuarine and Coastal Engineering, University Hannover].
Zhang
Q.
Hillebrand
G.
Klassen
I.
Vollmer
S.
Olsen
N. R. B.
Moser
H.
Hinkelmann
R.
2013
Sensitivity analysis of flow field simulation in the Iffezheim reservoir in Germany with the 3D SSIIM model
. In:
Proceedings of 2013 IAHR World Congress
,
Chengdu
,
China
.
Zhang
Q.
Hillebrand
G.
Moser
H.
Hinkelmann
R.
2015
Simulation of non-uniform sediment transport in a German reservoir with the SSIIM model and sensitivity analysis
. In:
Proceedings of 2015 IAHR World Congress
,
Delft
,
The Netherlands
.