Reservoir Sedimentation – Schleiss et al. (Eds)© 2014 Taylor & Francis Group, London, ISBN 978-1-138-02675-9

Three-dimensional numerical modeling of flow field in rectangularshallow reservoirs
T. Esmaeili & T. Sumi

Department of Urban Management, Graduate School of Engineering, Kyoto University, Kyoto, Japan

S.A. Kantoush

Department of Civil Engineering, German University in Cairo (GUC), Cairo, Egypt

A.J. Schleiss

Ecole Polytechnique de Lausanne (EPFL), Laboratory of Hydraulic Constructions (LCH), Lausanne,
Switzerland

S. HaunNorwegian University of Science and Technology (NTNU), Trondheim, Norway

ABSTRACT: Flow field in shallow waters, which is characterized by its complex mixing process and inherent dynamic nature, is interesting mainly due to its practical importance (e. g. in free
flushing operation and sedimentation in large reservoirs). 3D numerical models make it possible
to track two-dimensional large turbulence coherent structures, which are the dominant phenomenon in shallow reservoirs flow field. In the present study a fully three-dimensional numerical model SSIIM that employs the Finite Volume Approach (FVM) was utilized to reproduce
the 3D flow field. Various shallow reservoir geometries with fixed and deformed equilibrium
bed were considered. The measurements by Large-Scale Particle Image Velocimetry techniques(LSPIV) and Ultrasonic Doppler Velocity Profiler (UVP) over the flow depth were used formodel validation. Outcomes revealed reasonable agreement between the simulated and measured flow velocity field even when an asymmetric flow pattern exists in the reservoir.

1

Introduction

Shallow waters are defined as a flow field in which the vertical dimension of fluid domain is significantly smaller than its horizontal dimensions (Yuce & Chen 2003). Wide rivers, lakes, coastallagoons, estuaries as well as large reservoirs are the examples of shallow waters in the real life. Flowpattern in wide and shallow reservoirs with sudden expansion of inlet section may become unstable,
which produces large-scale transverse motions and recirculation zones due to the transverse disturbance because of the high sensitivity of flow pattern to the initial and boundary condition (Dewals
et al. 2008). This type of flows is prominent in the nature and also emerges in various engineering
applications including sudden expansions (Shapira et al. 1990), compound channels (Ghidaoui &Kolyshkin 1999 and Chu et al. 1991), storage chambers (Adamson et al. 2003), settling tanks (Freyet al. 1993) as well as shallow reservoirs sedimentation (Kantoush et al. 2008a & 2010).
When large-scale transverse motions and turbulent coherent structures emerge in shallow
reservoirs, the sediment transportation pattern would be seriously affected by the flow velocity field. Subsequently, measurement of 2D surface velocity and vertical velocity components
with high spatial resolution is essential to predict the favorable sedimentation zone. Such kindof knowledge leads to more efficient sediment management strategies in reservoirs. Also,
assessment of the flow field is necessary to characterize the domain of main jet flow, reverse
flow and eddies within a shallow reservoir.
11

Kantoush (2007) presented a comprehensive review of experimental tests in shallow reservoirs with transverse flow motions in symmetric channel expansions. The observations
revealed that asymmetric flow pattern can be developed under a certain geometric andhydraulic condition even if the symmetric geometry and hydraulic condition is employed.
Same outcome has been obtained by Adamson et al. (2003) and Stovin & Saul (1996) regarding the storage tank sedimentation and storage chambers, respectively. However, most of thestudies in the literature considered the sudden plane expansions of an infinite length andstudies focusing on various geometric and hydraulic parameters are very limited.
Mizushima & Shiotani (2001) used numerical model for studying the flow instabilities insymmetric channels with sudden expanded section. Dewals et al. (2008) and Dufresne et al.
(2011) used 2D numerical models to investigate the flow pattern distribution in shallow reservoirs. Nonetheless, the one and two dimensional models are not able to directly simulate
the secondary current influences in complex 3D flows, especially on deformed beds, since thecomplexity of flow pattern is further magnified.
In the present study, four different reservoir geometries with different length-to-width(Aspect ratio that is called AR hereafter) and reservoir width to inlet channel width (Expansion ratio that is called ER hereafter) were considered for numerical simulation. The fixed
bed and also equilibrium bed after free-flow flushing were used as the initial bed condition as
well. Numerical results were compared with the experimental measurements of surface velocity for all cases as well as 3D velocity component for one case and a satisfactory agreement
was found between the predicted and observed flow pattern. Furthermore, the numericalmodel could reproduce the steady asymmetric flow pattern when a perturbation quantifier
was introduced.
2
2.1

MATERIAL AND METHODS
Experimental setup and conditions

The experimental tests were carried out at the Laboratory of Hydraulic Constructions of
Swiss Federal Institute of Technology (EPFL) in a rectangular reservoir with the maximum
inner length (L) of 6 m and width (B) of 4 m. Also, the inlet and outlet rectangular channel
width (b) and length (l) were 0.25 m and 1 m respectively. Both channels were installed at thecenter of upstream and downstream side wall of the reservoir. The different shallow reservoir
geometry achieved experimentally by moving the PVC plate walls. The reservoir depth is
0.3 m and the both side walls and bottom is hydraulically smooth and flat. The water levelin the reservoir was controlled by a 0.25 m width and 0.3 m height flap gate set up at the endof outlet channel. A moveable frame with 4 m length was mounted on the side walls of thereservoir for installing the measurement devices. Table 1 shows geometrical attributes of fourgeometries employed in the present study. As for the geometrical parameters, ER = B/b shows
the influence of change in the reservoir width while AR = L/B is appropriate for describing
the effect of variations in the reservoir length (Dewals et al. 2008).
Large Scale Particle Image Velocimetry technique (LSPIV) was used for measuring thesurface velocity field. Ultra Sonic Velocity Profiler device (UVP) was employed for providing the 3D flow velocity measurements as well. Also, within the frame work of experimental
study, sedimentation and sediment flushing from the shallow reservoirs were investigated(Kantoush et al. 2008b). Non-uniform crushed walnut shells were used. The median sizeof this non-cohesive light-weight and homogenous grain material was 50 µm with σg of 2.4
and a density of 1500 kg/m³. The flow discharge rate (Q) and water depth (h) were constantfor all experiments as 0.007 m³/s and 0.2 m respectively except for the cases with deformed
equilibrium bed. Thus, in all tested configurations with fixed bed, the measured Froude andReynolds number at the inlet channel were kept constant as Fr = 0.1 and Re = 1.75 × 104. Incase of measurements of flow field on deformed bed after flushing the water level and discharge were 0.1 m and 0.007 m3/s respectively.
Figure 1 illustrates the observed flow streamlines for case T8 and T13. As can be observed,
the issuing flow jet deviated to the right hand side which forms asymmetric flow pattern in
12Table 1. Experimental configuration and corresponding geometrical parameters.
Case

L(m)

B(m)

Pr(m)*

ER(-)

AR(-)

SF(-)**

T8
T9
T11
T13

6.0
6.0
5.0
3.0

2.0
1.0
4.0
4.0

15.5
13.5
17.5
13.5

8.0
4.0
16.0
16.0

3.0
6.0
1.25
0.75

0.375
0.122
0.99
0.97

* Pr is the wetted perimeter.
** A is the total reservoir area and SF is the shape factor which defined as (A/Pr2) × ER.

Figure 1.

Observed streamlines of the surface flow field for a) case T8 and b) case T13.

case T8. The main eddy rotating anticlockwise in the center part of the reservoir and twoothers rotate clockwise in the upstream corners. Also, symmetric flow pattern with one main
jet trajectory in the centerline and two side eddies has been developed for case T13. Kantoush(2007) concluded that the deviation to the right hand side is due to the random disturbanceof the initial boundary condition and a mirror situation would be easily established by disturbing the initial condition slightly. The flow deviation to one side of the reservoir corresponds to the increase of flow velocity in one side of the jet and consequent reduction of
the pressure. This process will lead to flow deviation to one side of the reservoir and called
Coanda effect (Chiang et al. 2000).
2.2 Numerical model & features
Fully 3D numerical model SSIIM was employed in this study. The numerical model solves
the mass conservation and Reynolds-averaged Navier-Stokes equation in three dimensions(Equation 1 and 2) to compute the water motion for turbulent flow as follows.
∂U i
=0
∂xi

(1)

 ∂U i ∂U j  ∂U i∂U i 1 ∂  −Pδ ij + ρυT 
+U j
=
+
 
 ∂x

∂t
∂x j ρ ∂x j 
 j ∂xi  (2)in which i = 1, 2, 3 is the representative of three directions; where Ui is the averaged velocity,
x is the spatial geometrical scale, ρ is the water density, P is the Reynolds-averaged pressure,δij is the Kronecker delta and νT is the turbulent eddy-viscosity. For transforming the partial
equations into algebraic equations, the finite volume method is applied as discretization
method, together with the second order upwind scheme.
13

The change in water-levels was based on calculated pressure field. The pressure was extrapolated to the water surface and the pressure difference between a surface node and the downstream node was used to estimate the water elevation difference (Olsen 2013). The turbulence
is modeled by the standard k-ε model, using constant empirical values (Launder & Spalding1972). The unknown pressure field is calculated employing Semi Implicit Method for Pressure-Linked Equations, (SIMPLE) method (Patankar 1980). The grid is adaptive and moves
with change in the bed and water levels.
The Dirichlet boundary condition for the water inflow (logarithmic velocity distribution)
was used while for the water and sediment outflow zero-gradient boundary condition was
specified. For the boundary condition at the walls, where there is no water flux, the empiricalwall laws introduced by Equation 3 were used:U 1  30 y = ln u* κ  ks 

(3)where the shear velocity is denoted u*, κ is the Karman constant equal to 0.4, y is the distanceto the wall and ks is the equivalent roughness.
The sediment transport computation for simulating the morphological change is divided
into suspended and bed load transport. Suspended load is calculated by solving the transientconvection-diffusion equation formula and the bed load is simulated by Van Rijn formula(Van Rijn 1984a). In order to compute the suspended sediment concentration in the cells
close to the bed, a specified concentration was used as boundary condition (Van Rijn 1984b).
Also, the turbulent diffusivity is calculated by Equation 4.
ΓT =

νT
Sc

(4)

where Sc is the Schmidt number representing the ratio of eddy viscosity coefficient νT to diffusion coefficient and is set to 1.0 as the default value in SSIIM model (Olsen 2013).
Based on the experimental measurements computational mesh for all reservoir geometries
were made. The mesh cell size for case T8, T9, T11 and T13 in X and Y direction was
5 cm × 1.5 cm, 5 cm × 2.5 cm, 5 cm × 2 cm and 2.5 cm × 1 cm respectively. Considering the
11 cells for vertical grid distribution, the total number of cells over the main reservoir geometry were 174460, 52800, 220000 and 528000 respectively.
3 results and discussionS
3.1

Model calibration

As for the real cases in prototype scale, flow field modeling will provide us useful information about the areas with the potential erosion and deposition (sedimentation) during the
anticipated floods. The more accurate the prediction of sedimentation zone, the higher is theefficiency of sediment management strategies.
In the numerical simulation, the time step was calibrated as 2 seconds for run T8 and T9
whereas it was 0.5 second for T11 and T13. In case of flow field modeling on the equilibrium
deformed bed, the bed roughness was fixed as 0.00015 m which equals to 3 times the mediansediment size.
3.2

Flow field modeling on fixed beds

Simulations have been conducted by employing the geometry and inflow/outflow condition of
the experimental model. The k-ε turbulence model was utilized and numerical runs were performed until a steady-state flow condition is obtained. Simulations revealed that the model was
not able to reproduce an asymmetric flow pattern when the geometry configuration and hydraulic boundary condition are perfectly symmetric. This is because of the implemented mathematical algorithms that were not aimed to reproduce such kind of artificial numerical outputs.
14Due to the fact that symmetric flow pattern was not observed in the experiments for special
geometries, Dewals et al. (2008) introduced a slight disturbance in the initial boundary condition for 2D numerical simulations. They employed non-uniform cross-sectional discharge in
the inlet boundary to test the stability of the numerical model outputs. In the present study,
the same concept of slight disturbance in the inlet boundary condition was implemented
for all runs and non-uniform cross-sectional velocity distribution was utilized in the inflow
boundary condition. Therefore, the initial velocity magnitude differ 2.5% in one side of theinlet channel comparing to the other side. Such kind of disturbance is inevitable in the experimental set up. Nonetheless, very small perturbation of the inflow condition will impose significant effect on the flow field of case T8 and T9 in the numerical model. The reason can be
attributed to the unstable nature of symmetric flow in such geometries which pronounces thehigh sensitivity of the flow pattern to the inflow boundary condition.
Figure 2 shows the simulated and measured surface velocity magnitudes (V) in m/s and flow
distribution pattern. As can be clearly seen from Figure 2, asymmetric flow pattern has developed in case T8 and T9 whereas symmetric flow pattern is observed for case T11 and T13.

Figure 2. Left: The measured surface velocity field for runs (a1) T8, (b1) T9, (c1) T11 and (d1) T13
respectively and right: corresponding simulated velocity field.

15

Figure 3. The measured streamwise and transversal surface velocity at the middle cross section of thereservoirs versus simulated surface velocity for (a) T8, (b) T9, (c) T11 and (d) T13.

The model could simulate the surface flow velocity pattern almost similar to the measured
one by reproducing the dominant aspects such as the main flow jet trajectory and locationof the reverse flow as well as the main vortices and corner gyres. Nevertheless, numericalmodel outcomes show the straighter and longer reverse flow trajectory than observations
for all cases. Consequently, size of the upstream corner gyres in the numerical outputs issmaller than in the experimental measurements. This situation is predominant for case T11(Figure 2c1 and 2c2).
Figure 3 illustrates the simulated streamwise and transversal surface velocity distribution
versus the measured one in the middle cross-section of the reservoirs. Although the numericalmodel results are globally consistent with the measurements, there is a discrepancy for transversal surface velocity of case T11. The reason could be the concentrated flow pattern andlower diffusion of the main jet as well as the reverse flow in the numerical model outcomes.
The vertical distribution of flow velocity field is important for analyzing the sedimenttransport in reservoirs. Thus, the numerically simulated 3D flow velocity field was comparedwith the measured 3D velocity components provided by the UVP. A set of tree UVP probes
which were inclined at 20° to the vertical axis, allowed measuring the 3D flow field. The firstvalid UVP measurements were located at 12.5 cm away from the side walls and 2.5 cm fromthe free water surface (Kantoush, 2008b). Figure 4a and 4b demonstrates the longitudinalvelocity distribution over the flow depth, at upstream, middle, and downstream area of case
T9. It can be seen that the higher longitudinal velocity has been deflected towards the rightbank side and the reverse flow is reproduced beside the left bank side. Such kind of change
in the flow direction across the reservoir is also qualitatively consistent with the experimental
observations.
Velocity distribution in longitudinal, lateral and vertical directions (U, V, W respectively)
at x = 5.5 m and y = 0.375 m from the right bank of the case T9 has been demonstrated in
Figure 4c in order to quantitative assessment of the numerical model outputs. As for thelongitudinal velocity, the numerical model overestimates the magnitudes compared with theexperimental measurements while the lateral velocity was underestimated. This condition
was also found for other locations over the reservoir depth. Furthermore, the vertical velocity
16

Figure 4. a) The UVP measurements of the longitudinal velocity over the flow depth at x = 1.5 m,
3.5 m and 5.2 m; b) the corresponding simulated velocity; c) measured and simulated velocity profiles inthree dimensions: longitudinal (U), lateral (V), Vertical (W) at x = 5.5 m and y = 0.375 m position.

Figure 5.

The final equilibrium deformed bed after flushing for case T8.

Figure 6.

Surface velocity field with vectors after flushing for case T8; a) measured and b) simulated.

over the flow depth remains fixed in the numerical outputs whereas there was a slight change
in the measurements. However, magnitudes of the vertical velocity were very small compare
to the longitudinal and lateral velocities.
3.3 Flow field modeling on equilibrium deformed bed
The final bed morphology, which was obtained after the sediment flushing, was introduced
to the model as the boundary condition and then three-dimensional flow field was calculated.
Such kind of flow field simulation is more complex than flow field on fixed beds due to theshallower flow condition with higher velocity components on the existing friction of varying
bed forms.
Figure 5 reveals the final deformed bed after sediment flushing in case T8. The figure
shows that the flushing channel has shifted to the left side of the reservoir. Figure 6a and 6b
also show the measured surface velocity after flushing by LSPIV technique and simulated
one respectively. As can be clearly observed, the model reproduced the surface velocity pattern very similar to the observed one. Likewise the surface velocity pattern on fixed beds, the17Figure 7. a) Lateral velocity distribution over the depth at x = 3 m for T8 and b) comparison of surface
velocity components from the LSPIV measurements and simulations.

reverse flow trajectory is longer, straighter and also the reattachment length is smaller thanthe observations. However, slight disturbance in the boundary condition was not consideredin the model for the current case.
Figure 7a clearly illustrates the lateral flow velocity (V) contours over the depth along withthe secondary flow velocity vectors for case T8 at the middle cross-section. The streamwise
and transversal simulated velocities were also plotted against the measured one in Figure 7b
and good agreement was found.
4 conclusions
In the present study, numerical model SSIIM was utilized for reproducing the 3D velocity
field in rectangular shallow reservoirs with different geometries on fixed and deformed bed aswell. Measured surface velocity utilizing the LSPIV technique was used for model validation.
In addition, 3D velocity components over the flow depth provided by UVP compared to thenumerical results. The following outcomes have been obtained from the present work:
1. Many hydrodynamic aspects such as jet trajectory, recirculation zones, eddies and theflow distribution pattern, in different shallow reservoir geometries were represented by thenumerical model. The numerical model also reproduced both symmetric and asymmetricflow pattern in symmetric geometry setup, similar to the observations, by introducing a
perturbation in the inflow boundary. The numerical outcomes were in good agreementwith the experimental observations. However, there was some discrepancy for reproducing
the upstream vortex dimension due to the longer concentrated flow and reverse flow jetpattern with lower flow diffusion in the numerical outcomes.
2. The numerical model could successfully show the effect of geometry on the flow pattern.
More specifically, the numerical model showed that the flow pattern in geometries withhigher SF (T11 and T13) is insensitive to the small disturbance in the inflow condition andthe numerical model converged to the steady symmetric flow pattern. On the contrary,calculated flow pattern in geometries with lower SF (T8 and T9) converged to the steadyasymmetric flow pattern due to the slight disturbance in the inflow boundary condition.
3. From the practical point of view, 3D numerical models can be used for reproducing thesymmetric and asymmetric flow pattern under the various hydraulic and geometric conditions. Consequently, it would be possible to model the flow and sedimentation pattern
as well as the effects of diverse measures on the flow field for conducting the sedimentation in preferential zones. As to the reservoirs, the outcome can be the effective prediction
and management of flow field and sediment flushing process during drawdown operation (i.e. Dashidaira reservoir in the Kurobe River, Japan). Therefore, application of
3D numerical modeling beside the experimental studies result in more efficient sedimentmanagement strategies in reservoirs. However, the coupled simulation of flow and sediment field interacting with the movable bed during the drawdown flushing would be achallenging work.
18

references
Adamson, A. Stovin, V. & Bergdahl, L. 2003. Bed shear stress boundary condition for storage tanksedimentation. Journal of Environmental Engineering, ASCE 129(7): 651–658.
Chiang, T.P. Sheu, T.W.H. & Wang, S.K. 2000. Side wall effects on the structure of laminar flow overplane symmetric sudden expansion. Computers & Fluids 29(5): 467–492.
Chu, V.H. Wu, J-H. & Khayat, R.E. 1991. Stability of transverse shear flows in shallow open channels.
Journal of Hydraulic Engineering, ASCE 117(10): 1370–1388.
Dewals, B.J. Kantoush, S.A. Erpicum, S. Pirotton M. & Schleiss, A.J. 2008. Experimental and numerical analysis of flow instabilities in rectangular shallow basins. Environmental fluid mechanics 8(1):31–54.
Dufresne, M. Dewals, B. Erpicum, S. Archambeau, P & Pirotton M. 2011. Numerical investigationof flow patterns in rectangular shallow reservoirs. Engineering Application of Computational Fluid
Mechanics 5(2): 247–258.
Frey, P. Champagne, J.Y. Morel, R. & Gay, B. 1993. Hydrodynamics fields and solid particle transport
in a settling tank. Journal of Hydraulic Research, 31(6): 736–776.
Ghidaoui, M.S. & Kolyshkin, A.A. 1999. Linear stability of lateral motions in compound open channels. Journal of Hydraulic Engineering, ASCE 125(8): 871–880.
Kantoush, S.A. 2007. Symmetric or asymmetric flow patterns in shallow rectangular basins with sediment transport. JFK paper competition, In International Association of Hydraulic Engineering andResearch (ed.), 32rd IAHR Congress: Harmonizing the Demands of Art and Nature in Hydraulics;
Proc. intern. conf., Italy, 1–6 July 2007. Venice.
Kantoush, S.A. Bollaert, E. & Scleiss, A.J. 2008a. Experimental and numerical modelling of sedimentation in a rectangular shallow basin. International Journal of Sediment Research, 23(3): 212–232.
Kantoush, S.A. De Cesare, G. Boillat, J.L. & Schleiss, A.J. 2008b. Flow field investigation in a rectangular shallow reservoir using UVP, LSPIV and numerical modelling. Flow Measurement & Instrumentation, 19: 139–144.
Kantoush, S.A. Sumi, T. & Schleiss, A.J. 2010. Geometry effect on flow and sediment deposition patternin shallow basins. Annual Journal of Hydraulic Engineering, JSCE 54: 212–232.
Launder, B.E. & Spalding, D.B. 1972. Lectures in mathematical models of turbulence. London: Academic press.
Mizushima, J. & Shiotani, Y. 2001. Transitions & instabilities of flow in a symmetric channel with a suddenly expanded and contracted part. Journal of Fluid Mechanics 434: 355–369.
Olsen, N.R.B. 2013. A three dimensional numerical model for simulation of sediment movement in water
intakes with multiblock option. Department of Hydraulic and Environmental Engineering, The Norwegian University of Science and Technology, Trondheim.
Patankar, S.V. 1980. Numerical heat transfer and fluid flow. New York: McGraw-Hill.
Shapira, M. Degani, D. & Weihs, D. 1990. Stability and existence of multiple solutions for viscous flow
in suddenly enlarged channels. Computers & Fluids 18(3): 239–258.
Stovin, V.R. & Saul, A.J. 1996. Efficiency prediction for storage chambers using computational fluid
mechanics. Water Sciences & Technology 33(9): 163–170.
Van Rijn, L.C. 1984a. Sediment transport. Part I: Bed load Transport. Journal of Hydraulic Engineering,
ASCE 110(10): 1733–1754.
Van Rijn, L.C. 1984b. Sediment transport. Part II: Suspended load Transport. Journal of HydraulicEngineering, ASCE 110(11): 1431–1456.
Yuce, M.I. & Chen, D. 2003. An experimental investigation of pollutant mixing and trapping in shallow
costal re-circulating flows. In Uijttewaal, W.S.J. & Gerhard, H. J (ed.), Shallow flows; Proc. intern.
symp., Netherlands, 16–18 June 2003. Delft: Balkema.

19