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 m^{3} (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 m^{3}/s (LUBW 2011).

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.

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 [%] |
---|---|---|---|---|---|

1 | 20 | 88 | – | 568 | 0 |

2 | 3 | 8.4 | – | 216 | 0 |

3 | 1.3 | 2.3 | – | 136 | 0 |

4 | 0.4 | 1.3 | – | 56 | 2 |

5 | 0.13 | 0 | – | 11 | 11 |

6 | 0.04 | 0 | – | 1.1 | 19 |

7 | 0.02 | 0 | 0.027 | 0.27 | 35 |

8 | 0.005 | 0 | 0.0069 | 0.017 | 12 |

9 | 0.002 | 0 | 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 [%] |
---|---|---|---|---|---|

1 | 20 | 88 | – | 568 | 0 |

2 | 3 | 8.4 | – | 216 | 0 |

3 | 1.3 | 2.3 | – | 136 | 0 |

4 | 0.4 | 1.3 | – | 56 | 2 |

5 | 0.13 | 0 | – | 11 | 11 |

6 | 0.04 | 0 | – | 1.1 | 19 |

7 | 0.02 | 0 | 0.027 | 0.27 | 35 |

8 | 0.005 | 0 | 0.0069 | 0.017 | 12 |

9 | 0.002 | 0 | 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 m^{3}/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 m^{3}/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 m^{3}/s whereas 450 m^{3}/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.

## METHODS

*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).

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.

*c*, of size

*i*, was computed from the convection–diffusion equation for suspended material: The water velocity is denoted

*U*, the fall velocity of the particle is

*w*, Γ is the turbulent diffusion coefficient and

*F*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

_{R,i}*, c*, at the bed was used (van Rijn 1984a): The particle diameter is denoted

_{R,i}*d*,

_{i}*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*: where τ is the actual shear stress on the bed, computed from the turbulent kinetic energy,

*k*: 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 ρ

*is the sediment density. The acceleration of gravity is denoted*

_{s}*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.

*w*, and

*f*is the fraction of the bed material of size i. The parameter

*f*will vary between 0 and 1. The variable

*c*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.

_{c}*z*, were computed based on the computed sediment deposition and the pick-up rate for all the size fractions: The computed concentration in a bed cell is denoted

*c*, while

_{i}*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.

*t*, according to the following formula: The values used in the current study for the reference case run were Δ

*t*

_{0}*=*20,000 seconds,

*n*= 3 and

*Q*= 1,000 m

_{ref}^{3}/s.

## RESULTS AND DISCUSSION

### Velocity field

The velocity field was measured and computed for a constant water discharge of 1,550 m^{3}/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 m^{3}/s, whereas 450 m^{3}/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.

### 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.

### 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 m^{3} over a period of 3 months. The uncertainty of the field measurements is estimated to be up to 20%, i.e., approximately 10,000 m^{3}. Therefore, the computed sediment volumes are expected to be in the range between 41,000 m^{3} and 61,000 m^{3}. 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.

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) |

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 m^{3} and 61,000 m^{3}, 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 m^{3} 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.