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