京都大学防災研究所年報 第 57 号 B 平成 26 年 6 月
Annuals of Disas. Prev. Res. Inst., Kyoto Univ., No. 57 B, 2014

3D Numerical Modeling of Flow Field and Flushing Channel Formation
in Reservoirs
Taymaz ESMAEILI(1), Tetsuya SUMI, Sameh A. KANTOUSH(2) and Yoji KUBOTA(3)(1) Department of Urban Management, Graduate School of Engineering, Kyoto University(2) Department of Civil Engineering, German University in Cairo (GUC)(3) Hydro-soft Technology Institute Co., Ltd.,
Synopsis
Flushing with water level drawdown is one of the proposed methods for preserving
the storage capacity of dam reservoirs. The accelerated shallow flow pattern during thesediment flushing increases the induced bed shear stress and then the deposited sedimentin the dam reservoir would be eroded. Thus, Shallow flow condition which ischaracterized by its complex mixing process and inherent dynamic nature plays majorrole in draw down flushing process. In the present study a fully three-dimensional
numerical model that employs the Finite Volume Approach (FVM) was utilized toreproduce the 3D flow field in various shallow reservoir geometries with fixed anddeformed equilibrium bed. Afterwards, the model with wetting/drying algorithm andadaptive mesh which moves in accordance to water and bed level changes was employed
to simulate the sediment flushing process of Dashidaira reservoir. The total volume offlushed out sediment and bed deformation after flushing compared with the observations.
Numerical outcomes show that the model could reproduce many aspects of the complex
procedure of flushing process and a reasonable agreement was found between the
numerical outcomes and observations.
Keywords: Shallow flows, Reservoir management, 3D numerical simulation, free-flow
sediment flushing1.

Introduction

Sediment deposition is the principal problem
affecting the useful life of reservoirs. A decreased
storage volume reduces the reservoir function forflood control purpose, electricity production andwater supply. Furthermore, the sediment yields of
the Japanese rivers are high in comparison withother countries due to the topographical, geologicaland hydrological conditions (Sumi, 2006). Among
several techniques, the flushing and sluicing plays
an important role in the sediment removal andreduction, as they are efficient hydraulic sedimentremoval technique to restore the reservoir storage
capacity. During the flushing with drawdown,bottom outlets are opened and water level islowered to generate and accelerate the unsteadyflow towards the outlet (Morris & Fan, 1998). This
process will initiate the progressive andretrogressive erosion pattern in tail and delta
reaches of the reservoirs respectively (Battuca,
2000). Many reservoirs are not flushed because an
accurate prediction of the success of an upcoming
reservoir flushing is time and cost-intensive (Haun― 468 ―

& Olsen, 2012). Loss of the stored water besides thelack of knowledge about the amount of the flushed
out sediment and location of the erosion-depositionarea in the reservoir will be the negative factors forchoosing the free-flow sediment flushing.
Following that, the environmental effects
downstream of the reservoir caused by the flushing
are difficult to estimate because of the complexityof free-flow sediment flushing process. So, it isimportant to have an initial assessment of theefficiency of upcoming free-flow sediment flushing.
Numerical models were developed as an alternative
over the last decades to avoid scaling problems ofphysical models and also because of an expected
time and cost reduction (Chandler et al., 2003).
Three-dimensional numerical models are still underdevelopment for application in this field. Atkinson(1996), Khosronejad (2008), Haun & Olsen (2012),
Esmaeili et al. (2014) used three-dimensionalmodels for simulating the flushing process in thephysicalmodel
scale.
Applicationofthree-dimensional numerical models for simulating
the complex structure of interacting flow andsediment field would be useful for understanding
the different stages of flushing process and utilize
the outcomes for application in the prototype scale.
Free-flow sediment flushing through dam
reservoirs is one of the typical techniques in Japanbecause of sufficient rain fall, economical
advantages, independency from other facilities andpotential for high efficiency. During the free-flowstate, after lowering the water level, the shallow
flow plays major role in sediment transportation
pattern and bed deformation (i.e. flushing channel
formation) in reservoir. Therefore, it is necessary to
be familiar with the behavior of shallow waters
characteristic. 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). This type offlows is prominent in the nature and also emerges invarious engineering applications including sudden
expansions (Shapira et al., 1990), compoundchannels (Ghidaoui & Kolyshkin, 1999 and Chu etal., 1991), storage chambers (Adamson et al., 2003),
settling tanks (Frey et al., 1993) as well as shallow
reservoirs sedimentation (Kantoush et al., 2008a &2010a).

The sediment transportation pattern in shallow
reservoirs is seriously affected by large-scaletransverse motions and turbulent coherent structures.
Therefore, measurement of 2D surface velocity andvertical velocity components with high spatial
resolution would be necessary to predict thefavorable sedimentation zone. Such kind ofknowledge would be interesting for more efficient
sediment management strategies in the reservoirs.
In addition, assessments of the flow field within ashallow reservoir is necessary for characterizing theflow domain such as location of main jet flow,
reverse flow and eddies.
Kantoush (2007) presented a comprehensive
review of experimental tests in shallow reservoirs
with transverse flow motions in symmetric channel
expansions. The observations revealed thatasymmetric flow pattern can be developed under acertain geometric and hydraulic condition even if
the symmetric geometry and hydraulic condition is
employed. Same outcome has been obtained byAdamson et al. (2003) and Stovin & Saul (1996)
regarding the storage tank sedimentation andstorage chambers, respectively. However, most of
the studies in the literature considered the sudden
plane expansions of an infinite length and studies
focusing on various geometric and hydraulicparameters are very limited. Dewals et al. (2008)
and Dufresne et al. (2011) used 2D numericalmodels to investigate the flow pattern distribution
in shallow reservoirs. Nevertheless, the one and two
dimensional models are not able to directly simulate
the secondary current influences in complex 3D
flows, especially on deformed beds due to thecomplexity of this type of flow fields.
In the present study, first, flow pattern and thepotential instabilities in diverse shallow reservoir
geometries with different length-to-width (Aspectratio that is called AR) and reservoir width to inlet
channel width (Expansion ratio that is called ER)
were considered for numerical simulations. The
fixed bed and also equilibrium bed after free-flow
flushing were used as the initial bed condition as
well.
Numerical outputs were compared with theexperimental measurements of surface velocity by
Large Scale Particle Image Velocimetry (LSPIV)
for all cases as well as 3D velocity component for― 469 ―

one case. Besides the good agreement between themeasured and simulated flow velocity field,
numerical model could reproduce the steadyasymmetric flow pattern when a perturbation
quantifier was employed.
Afterwards, the free-flow sediment flushing
event which was performed during a flood event inDashidaira reservoir in Japan was modeled to assess
the bed topography changes in different stages ofsediment flushing process. The cross sections in a
specific segment of the reservoir were comparedwith the measurements. Results show that model
could reproduce the main governing aspects of thesediment flushing event in a real prototype scale.
However, +50/-50 percent differences between themeasurement and simulation result in the real
prototype scale is not uncommon especially when
the geometry is complex and uncertainty exists inthe data collection.
2.

Material and methods

2.1 Three dimensional numerical model
Fully 3D numerical model SSIIM was employedin this study. The numerical model solves the massconservation and Reynolds-averaged Navier-Stokes
equation in three dimensions (Equation 1 and 2) tocompute the water motion for turbulent flow asfollows.

U i
0
xi

(1)

U j  (2)U i
U i 1  
U
U j

 P ij  T ( i 
)

t
x j
 x j 
x j
xi 

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.
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 thedownstream 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 & Spalding,1972). The unknown pressure field is calculated
employing
Semi
Implicit
MethodforPressure-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 usedwhile 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 empirical wall 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 distance tothe wall and ks is the equivalent roughness.
The sediment transport computation forsimulating the morphological change is divided into
suspended and bed load transport. Suspended load
iscalculatedbysolvingthetransientconvection-diffusion equation formula (Equation 4)
and the default empirical formula computing the
bed load is Van Rijn formula (Equation 5)(VanRijn,
1984a). In addition, Meyer-Peter-Muller bed loadtransport formula has been implemented innumerical model (Equation 6). This formula is more
appropriate than others for steep rivers in which thesediment load transportation is mainly bed load (i.eKurobe river). The latter formula has been
examined in European rivers and sediment flushingin Bodendurf reservoir, Austria (Haun et al., 2012)

c
c
c

c
U j
w 
(T
)
t
x j
z x jx j

(4)Where U is the water velocity, w is the fallvelocity of sediments, Γ is the turbulent diffusivity
and can be expressed by (Equation 7) andis thesediment concentration over time t within thespatial geometrical scales x and z.
Numerical model uses an unstructured adaptive
grid and only the flow phase (water body) ismodeled. Since wetting/drying algorithms were

― 470 ―

implemented to avoid emerging distorted cells,
number of cells in spatial directions can vary during
the computation. Thus, it is necessary to introduce
the boundary values for wetting/drying algorithm.
In this paper for simulation of sediment flushing
event, z 1= 0.3m (lower boundary) and z 2= 1m wereused. In areas with a water depth shallower than
0.3m no cell would be generated and the area will
be dried up (The area disappears from thecomputational grid). In areas with water depth
between 0.3m and 1m one cell is produced and a
two-dimensional computation is conducted whilethe various numbers of cells along the verticaldirection is generated in areas with the water depth
higher than 1m. In the latter case, the vertical cellnumbers is an exponential function of the ratio of
local water depth to the maximum water depth.

q b ,idi

1.5( s   w )gw

 0.053   c ,i   c ,i 

1.5

0.3  (    w ) g 
d i  s2  w(5)
0.1

where qb,i is sediment transportation rate for the ith
fraction of bed load per unit width, d i is thediameter of the ith fraction, τ is the shear stress, τc,iis the critical shear stress for di which was
calculated from the Shield's curve, ρs is the densityof sediment, ρw is the density of the water, g is thegravity acceleration and ν is the kinematic
viscosity.

1   w grI  0.047 g (  s   w )d 50 
qb ,i  2g
1     3w
0.25 w 3  s s 

(6)

In Meyer-Peter-Muller bed load formula d50 isthe characteristic sediment size, I is the slope of theenergy line and r is the hydraulic radius.

T 

T
Sc

(7)

where Sc is the Schmidt number representing theratio of eddy viscosity coefficient ν T to diffusion
coefficient and set to 1.0 as default.
In order to compute the suspended sedimentconcentration in the cells close to the bed Equation

8 was used as boundary condition (Van Rijn,
1984b).
1.5

   c ,i 


0.3
d
  c ,i 
Cbed  0.015 i
a  (    ) g  0.1
s
w


2


w



(8)

where a is the distance reference level and set equal
to roughness height.
2.2 Experimental setup for flow field modeling inshallow reservoirs
The experimental tests were carried out at theLaboratory of Hydraulic Constructions of Swiss
Federal Institute of Technology (EPFL) in arectangular reservoir with the maximum innerlength (L) of 6m and width (B) of 4m (Kantoush,2007). Also, the inlet and outlet rectangular channel
width (b) and length (l) were 0.25m and 1m
respectively. Both channels were installed at thecenter of upstream and downstream side wall of thereservoir. The different shallow reservoir geometry
achieved experimentally by moving the PVC plate
walls. The reservoir depth is 0.3m and the both sidewalls and bottom is hydraulically smooth and flat.
The water level in the reservoir was controlled by a
0.25m width and 0.3m height flap gate set up at theend of outlet channel. A moveable frame with 4m
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 forthe geometrical parameters, ER=B/b shows theinfluence of change in the reservoir width while
AR=L/B is appropriate for describing the effect ofvariations 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 theshallow reservoirs were investigated (Kantoush etal., 2008b). Non-uniform crushed walnut shells
were used as suspended load in sedimentation andflushing tests. The median size of this non-cohesive― 471 ―

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) wereconstant for all experiments as 0.007m³/s and 0.2m
respectively except for cases with deformed
equilibrium bed. Thus, in all tested configurations
with fixed bed, the measured Froude and Reynolds
number at the inlet channel were kept constant as
Fr=0.1 and Re=1.75×10 4 respectively. In case ofmeasurements of flow field on deformed bed after
flushing, the water level and discharge were 0.1m
and 0.007m3/s respectively.
Figure 1 illustrates the observed flow streamlines
for case T8 and T13. As can be observed, theissuing flow jet deviated to the right hand sidewhich forms asymmetric flow pattern in case T8.
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 isdue to the random disturbance of the initialboundary condition and a mirror situation would be
easily established by slight disturbance in theboundary condition. The flow deviation to one sideof the reservoir corresponds to the increase of flow
velocity in one side of the jet and consequent
reduction of the pressure. This process will produce
the flow deviation to one side of the reservoir andcalled Coanda effect (Chiang et al., 2000).
Table 1. Geometrical parameters of the reservoirs
Case

L(m)

B(m)

Pr(m)*

ER

AR

SF**

T8

6

2

15.5

8

3

0.375

T9

6

1

13.5

4

6

0.122

T 11

5

4

17.5

16

1.25

0.99

T 13

3

4

13.5

16

0.75

0.97

*Pr is the wetted perimeter.
** A is the total reservoir area and SF is the shapefactor defined as (A/Pr 2)×ER.

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

2.3 Case study description: Dashidaira reservoir,
Kurobe River, Japan
Kurobe River located in the eastern part of the
Toyama prefecture. The river originates from the
Washibadake Mountain and then extends to theplain area at the Aimoto point which is 13km
upstream of the coastal area of Japan Sea. The
Catchment area of the Kurobe River is 682 km 2 andthe length of the river is 85km. River bed is steep
and the bed slope varies between 1% and 20%. The
dominant sediment type is weathered granite andowing to the steep valleys and consequent sandslides the sediment production potential is very highin this river basin (Minami et al., 2012). The annual
rainfall in the mountain area of the basin is 4000mm(PWRI, 1993). Dashidaira dam with 76.7m height
has been constructed in 1985 by Kansai ElectricCompany over the Kurobe River that has thehighest sediment concentration among the rivers inJapan. This dam provides 124MW hydroelectric
power. Also, this dam has been equipped withbottom outlet facilities to flush out the deposited
sediment in the reservoir. The gross and effective
storage capacity of the dam reservoir is 9.01 and1.66 MCM respectively. Until year 2000 50% of thegross storage capacity and 29% of the effective
storage capacity was filled due to the sedimentation.
The sediment flushing operation has been
performed since 1991. From June 2001 thecoordinated sediment flushing in both Dashidairaand Unazuki reservoir, which is 7 km downstreamof Dashidaira dam, is ongoing (Liue et al., 2004).
Coordinated sediment flushing is executed at thefirst major flood event every year and the sedimentsluicing is done at the successive bigger ones by thesame reservoir operation to prevent additional
sediment deposit in the reservoir (Sumi, 2006).
Figure 2 illustrates the location of Dashidairareservoir in Kurobe River Basin. Although
Dashidaira reservoir is a narrow shaped, there are
local wide segments. Figure 3 shows the
topographic map of Dashidaira reservoir which was
extracted from the surveying data before thesediment flushing operation in 2012.
As for the coordinated sediment flushing, if theflood satisfies a stipulated standard, both reservoirs
water level is lowered based on the instructions toincrease the induced bed shear stress. As a result the― 472 ―

incoming sediment and deposits in the reservoirs
are flushed out through the sediment flushing gatesto downstream channel (Minami et al., 2012).
Figure 4 demonstrates the sediment flushing
through the Dashidaira reservoir during thefree-flow state.

contributed in preserving the storage capacity ofDashidaira and Unazuki reservoir noticeably.

Fig. 4 Aerial photo from free-flow sedimentflushing at Dashidaira Reservoir.

Nowadays, sediment flushing operation in the
Kurobe River is aimed to keep a constant bed form
without trapping sediment in the reservoir as muchas possible by executing the flushing within the

280

B

Discharge (m³/s)

A

C
A

B
C

335
Dischrage

260

Water level

330

240

325

220

320

200

315

180

310

160
140

Fig. 3 Topography of Dashidira reservoir before
flushing event in May 2012 and considered crosssections for further quantitative assessment.
natural flood events between June and August. As aresult, the water quality deterioration such as
happened in the 1991 sediment flushing event, has
not occurred recently. The mentioned approach also

Water level (m)Fig. 2 Location of Dashidaira and Unazuki
reservoir in Kurbe River Basin (Kantoush et al.,2010b).

2.3.1. Data collection
Due to the importance of bed level assessment inDashidaira reservoir for sustainable sediement
management, three times annual bed changes
measurement are done in the reservoir (i.e. before
and after flushing). The bed level measurements
before flushing were employed to prepare the initial
bed topography of the reservoir.
The water level at the dam site as well as theInflow and outflow discharge is recorded during thefree-flow sediment flushing operation with the 60
minutes intervals. These data were utilized as thehydrodynamic boundary condition to introduce thewater level draw-dawn procedure as well as
discharge changes to the numerical model. Figure 5
shows the water level and discharge fluctuations
during the free-flow sediment flushing operation
event in June 2012.

305

0

4

8

12

16

20

24

28

32

36

40

44

48

Time (hrs)Fig. 5 Water level and discharge rates during theflushing operation in June 2012.
The grain size distribution of the bed materials in
the reservoir was determined based on the taken
onsite samples. The sorting process of the bed
material sizes in the reservoir can be clearly― 473 ―

observed. While d50 of the sediment at the entrance
area of the reservoir was in the range of
32mm-76mm (very coarse gravel and cobble), thed50 decreased to less than 0.5mm in downstream
areas. Bed material size distribution within thereservoir area was obtained based on the sizeinterpolation of the available bed material samples.
7 sediment size fractions were considered that is
ranged between 0.374mm and 316mm. The densityof the sediment determined 2.65gr/cm3 while theporosity assumed to be in the range of 50%-70%.
The free-flow sediment flushing event in June2012 was a special case in which the major
sediment inflow from the upstream boundary waswash load. Wash load is basically assumed to betransported without deposition in the reservoir
because of their small size that is mainly in therange of clay and silt. Moreover, a reliable washload measurement when the water level anddischarge fluctuates during the free-flow sedimentflushing is scarce. Thus, in the sedimenttransportation process (i.e. erosion & deposition) in
the reservoir, wash load effect was neglected.
3.

Results & Discussions

3.1 Numerical model ability to represent
instabilities in flow field
Simulations have been carried out by employing
the geometry and inflow/outflow condition of theexperimental model. The k-ε turbulence model was
employed and numerical runs were done until asteady-state flow condition is achieved.
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 5cm  1.5cm,
5cm  2.5cm, 5cm  2cm and 2.5cm  1cm
respectively. Considering the 11 cells for verticalgrid distribution, the total number of cells over themain reservoir geometry were 174460, 52800,
220000 and 528000 respectively.
Simulations showed that the model was not able
to reproduce an asymmetric flow pattern when thegeometry configuration and hydraulic boundary
condition are completely symmetric. This is owing
to the implemented mathematical algorithms whichwere not aimed to reproduce such kind of artificial

numerical results.
In order to have such kind of asymmetric flow
pattern, slight disturbance in the inlet boundary
condition was implemented for all runs andnon-uniform cross-sectional velocity distribution
was utilized in the inflow boundary condition. Thus,
the initial velocity magnitude in one side of the inlet
channel differs 2.5% in comparison to the other sidein all cases. Such kind of disturbance is anunavoidable phenomenon in the physical modeling.
This small perturbation of the inflow condition will
impose significant effect on the flow field of narrow
reservoirs (i.e. case T9) in the numerical model.
This would be because of the unstable nature ofsymmetric flow in such geometries which reveals
the high sensitivity of the flow pattern to the inflow
boundary condition.
Time step was calibrated as 2 seconds for run T8
and T9 whereas it was 0.5 second for T11 and T13
in the numerical simulation. In case of flow field
simulation on the equilibrium deformed bed, thebed roughness was fixed as 0.00015m which equals
to 3 times the median sediment size.
Figure 6 shows the simulated and measured
surface velocity magnitudes (V) in m/s and flow
distribution pattern. As can be clearly observed
from Figure 6, asymmetric flow pattern has
developed in case T9 whereas symmetric flow
pattern is seen for case T13. The model could
simulate the surface flow velocity pattern almost
similar to the measured one by reproducing thedominant aspects such as the main flow jet
trajectory and location of the reverse flow as wellas the main vortices and corner gyres. Nonetheless,numerical model outcomes show the straighter andlonger reverse flow trajectory than observations inboth cases. Consequently, size of the upstream
corner gyres in the numerical outputs is smaller
than the experimental observations.
Figure 7 illustrates the simulated streamwise andtransversal surface velocity distribution versus themeasured one in the middle cross-section of case T9& T13. Outputs reveal that the numerical model
results are globally consistent with the observations.
However, small discrepancy between the simulation
results and measurements is due to the concentrated
flow pattern and lower diffusion of the main jet aswell as the reverse flow in the numerical model― 474 ―

overestimates the longitudinal velocity magnitudes.
Furthermore, numerical model results reveal a more
uniform distribution of longitudinal velocity invertical and lateral direction than measurements.

outcomes.

X=3m_simulated
streamwise velocity
X=3m_measured
streamwise velocity

1

(a1)
Y(m)

0.75

X=3m_simulatedtransversal velocity
X=3m_measured
transversal velocity

0.5
0.25

(a2)

0
-0.15

-0.05

0.05

0.15

0.25

Surface velocity (m/s)(a)
X=3m_simulated
surface velocity

4

X=3m_measured
surface velocity

Y(m)

3

X=3m_simulatedtransversal velocity
X=3m_measured
transversal velocity

2

1

(b1)

0
-0.10

0.00

0.10

0.20

0.30

Surface velocity (m/s)(b)
Fig. 7 The measured streamwise and transversal
surface velocity along the middle cross section of
the reservoirs versus the simulated surface velocityfor (a) T9 and (b) T13.

(b2)
Fig. 6 Measured surface velocity field withvelocity vectors for runs (a1) T9 and (b1) T13;
Corresponding simulated velocity field for runs(a2) T9 and (b2) T13.
(a)
The vertical distribution of flow velocity field isalso important for assessing the sediment transport
in reservoirs. Thus, the numerically simulated 3D
flow velocity field was compared with the measured
3D velocity components provided by the UVP.
Figure 8a and 8b 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 right bank side andthe reverse flow is reproduced beside the left bankside. Such kind of change in the flow direction
across the reservoir is also qualitatively consistent
with the experimental observations. Compared withthe experimental measurements, numerical model(b)
Fig. 8 (a) The UVP measurements of thelongitudinal velocity over the flow depth at x=1.5m,
3.5m and 5.2m; (b) the corresponding simulated
velocity.
3.2 Flow field simulation on deformed beds
The final bed morphology, which was obtained
after the sediment flushing, was introduced to themodel as the boundary condition and thenthree-dimensional flow field was calculated. This
information would be useful for flood risk

― 475 ―

flushing event in June 2012.
Since Dashidaira reservoir has been located in aassessment especially in reservoirs near the urban
areas.
Flow field calculation on a deformed bed ismuch more complex than that of on the fixed bed
owing to the shallower flow condition with higher
velocity components on the existing friction of bed
forms. The final deformed bed after sedimentflushing in case T8 and T11 has been shown in
Figure 9. Figure shows that the flushing channel has
shifted to the left side of the reservoir in case of T8.
Figure 10 also illustrates the measured surface
velocity after flushing by LSPIV technique versus
the simulated one for case T8 and T11. As can beclearly observed, model reproduced the surface
velocity pattern similar to the observed one in case
T8 while it was different in case T11. In case T11,
the main and reverse jet trajectory is narrower andstraighter than that of the measured one. Moreover,slight disturbance in the boundary condition was notconsidered in the model for the current cases.

(a1)

(a2)

(b1)
(a)

(b2)
Fig. 10 Measured surface velocity field with
vectors on deformed bed for runs (a1) T8 and (b1)
T9. Corresponding simulated velocity field forruns (a2) T8 and (b2) T11.

Figure 11 demonstrates the quantitative
comparison between the measured surface velocity
by LSPIV and the simulated one for case T8 in the
middle cross-section.

Surface velocity (m/s)(b)
Fig. 9 Equilibrium deformed bed after flushing for(a) case T8 and (b) case T11.

0.25

X=3m_simulated streamwise velocity

X=3m_measured streamwise velocity

0.15

X=3m_simulated transversal velocity

X=3m_measured transversal velocity

0.05

-0.05
-0.15
0.00

0.50

1.00

1.50

2.00

Y( coordinate across the reservoir (m))

3.3 Numerical modeling of flushing channel
formation in Dashidira reservoir
The modeled bed level changes after the flushingwere validated with the measured data after the

Fig. 11 Measured surface velocity field versus thesimulated one for case T8.
steep mountainous region with high sediment― 476 ―

transportation in the bed load form, theMeyer-Peter-Mueller
sediment
transportation
formula was employed. Due to the implicit
free-water surface algorithm in which the computedpressure field is used for updating the water surface,large time steps can be utilized for numericalsimulations. Consequently, time step of 10 seconds
was selected and the active sediment layer thickness
adjusted to 1 m in the model. The angle of repose
for the sediment ranged between 30 °-32° and valueof 32° was applied in the sand slide algorithm. As tothe hydrodynamic boundary condition, water leveland discharge rates in Figure 5 were introduced tothe model and numerical results were presented.
Second order upwind scheme (SOU) was
employed in numerical modeling since it can handle
more complex flow field than first order power law
scheme (POW), especially in the channel bending,where the secondary currents play major role inerosion and deposition pattern (Haun & Olsen,2012). Figure 12(a) illustrates the measured bed
levels before flushing. In Figure 12(b) and 12(c),the measured and simulated bed topography after
flushing has been shown respectively. Formation of
the flushing channel during the free-flow state,similar to the prototype, was observed in the
numerical model outcomes as shown in Figure 13.
Regarding the surface flow velocities during thefree-flow condition, super critical flows emergedwith the Froude number between unity and 1.45 inseveral sections. Meanwhile, the horizontal flow
velocities in the flushing channel increased up to
5.5 (m/s) during the draw down stage (Figure 14).
Because of the presence of the air in supercriticalflows, the free-water surface is disturbed andsubsequently it is very likely to have uncertainty inthe outputs of this stage.
As can be seen from Figure 12, in the area close
to the dam the measured bed levels are different
than measured one (Area Ӏ). The simulation results
show narrower flushing channel in comparison tothe observations. This can be attributed to thecomplex flow field due to the secondary currents
and reverse flow pattern in the wide bending area
close to the dam site. Moreover, bed forms (i.e.
dunes and anti-dunes) which are developed in thisarea will increase the complexity of the flow andsubsequently the sediment transportation pattern.

(a)

(b)

III

II

I

(c)
Fig. 12 (a) Measured bed topography of Dashidairareservoir before flushing; (b) Measured bed levels
after flushing, and (c) simulated bed levels.

(a)

6 0 0 .0 m

Le ve l 2(b)

6 0 0 .0 m

Le ve l 2

Fig. 13 (a) Computational grid at the beginning offlushing simulation. (b) computational grid during
the free-flow condition which shows the flushingchannel shape and location.

― 477 ―

It should be noted that, numerical model utilizes
an empirical formula by Van Rijn (1984c) that isbasically developed for small Froude numbers.
Thus, it is very likely to simulate different bed
forms than that is formed during the sedimentflushing process. Moreover, in the areas close to thedam the sediment transport capacity is reached to its
highest level that would be smaller than its originalvalue in the prototype. Thus, enough sediment
could not be removed from this area and bed levels
are different than measured ones. It should be also
noted that the complex drawdown procedure isconducted by contribution of spillways and bottom
outlets in the prototype. In the numerical model,
this procedure has been simplified.

Fig. 14 Simulated flow velocity field in the flushingchannel during the free-flow condition.
In Area ӀӀ, which is the wide middle section of
the reservoir, the flushing channel shape andlocation is more or less consistent with theobservations (Fig. 12a-12c). In this area, there is alarge amount of deposited sediment. Red hatchedrectangles in Figure 12c reveal the regions withhigh amount of sediment deposition. The flushingchannel geometry in Area II has a noticeable
contribution in preserving the reservoir storage
capacity as well as enhancing the flushingefficiency. Three cross-sections (A-C; Fig. 3) fromthe outcomes of the numerical model have beenpresented in Figure 15. In addition, Figure 16 shows
the longitudinal profile at the centerline of thereservoir before and after flushing.
During the free-flow flushing operation the
narrow flushing channel location is deflected to theright bank region of Area II and consequently there
is no chance for erosion of deposited sediment in
the left bank region. A large amount of deposited
sediment in the form of a big sand bar will remain

after flushing operation. Subsequently, applicationof appropriate measures is likely to increase thevolume of sediment evacuation from this segment
which may enhance the flushing efficiency. Thus,
this study mainly focuses on flushing channel shape,
location and geometry within Area II.
Figure 15 clearly shows the bed degradation as
well as the flushing channel geometry in Area ӀӀ.
Similar to the measurement, numerical results showthat thalweg of the flushing channel has been
formed close to the right bank of wide middle
section (Area II). This is expected since bed
degradation usually occurs along the outer part ofthe channel bending. In all cross sections, bed
degradation has been over estimated in right half ofthe cross section. Concerning the bed levels in left
half, it has been underestimated in case ofcross-section A whereas it was more or lessconsistent with the measurements in cross-section B
and C.
In this study, the volume of the eroded sedimentfrom the reservoir was also computed. The total
computed volume of the eroded sediment was
308,900 m3 while it is 408,700 m3 based on thesurveying data. The lower volume of the eroded
sediment in numerical simulation is due to theunderestimation of the eroded sediment quantity in
Area I and II.
As can be observed from Figure 12, thesimulated flushing channel width in Area I isnarrower than measurements. As for the simulated
flushing channel in Area III, the channel could not
be propagated enough towards the far upstream
region. The reason can be attributed to the type ofthe bed materials in this area which are mainlycoarse gravels. Subsequently model could not
simulate the erosion pattern properly in thementioned region due to the complex velocity
profile and effective hydrodynamic forces that is
emerged especially during the free-flow state. The
other reason is elimination of the cells with shallow
water head, from the computational domain, along
the left bank of lower part of Area III in where there
is a sever bed degradation. It also should be taken
into account that the bed measurements is notconducted right after flushing event and also thesediment flushing is conducted in the rainy season
inDashidairareservoir.
Therefore,some― 478 ―

uncertainties would be possible in the bed levelmeasurements because of the dynamic changes inhydrodynamic boundary condition right after the
flushing event until the measurements. The otheruncertainty is related to the bed material
distribution that can have noticeable effect on theaccuracy of the numerical model outcomes.
However, the model results are more or less inagreement with the measurements in the main study
segment (Area II). An extensive sensitivity analysis
on the numerical algorithms and empiricalparameters (i.e. in bed load transport formula) in
the model will offset the uncertainties to someextent and also result in higher accuracy.
350

Measured bed level_before flushing
Measured bed level_after flushing
Simulated bed level

Elevation (m)

345
340
335
330
325

320
315
0

20

40

60

80

100

120

140

160

180

Distance (m)(a)
350
Measured bed level_before flushing
Measured bed level_after flushing
Simulated bed level

Elevation (m)

345
340
335

330
325
320

315
0

20

40

60

80

100

120

140

160

180

200

180

200

220

260

280

Distance (m)

(b)
350

Measured bed level_before flushing
Measured bed level_after flushing
Simulated bed level

Elevation (m)

345
340

335
330

325
320
315
0

20

40

60

80

100

120

140

160

220

240

Distance (m)(c)
Fig. 15 Bed level changes in (a) cross-section A; (b)cross-section B and (c) cross-section C.

350
Measured bed level_before flushing
Measured bed level_after flushing
Simulated bed level

Elevation (m)

340
330

320
310
300
0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

1.3

1.4

1.5

1.6

1.7

Distance from upstream (km)Fig. 16 Longitudinal bed profile at the centerline of
the reservoir.

4.

Conclusions

Due to the major role of shallow flow conditionin flushing channel formation and development
during the free-flow state, it is necessary to be
familiar with the behavior of such kind of flows.
In this paper, first, numerical model SSIIM wasutilized for reproducing the 3D velocity field inrectangular shallow reservoirs with different
geometries on fixed and deformed bed as well.
Measured surface velocity utilizing the LSPIV
technique was used for model validation. In
addition, 3D velocity components along the flow
depth provided by UVP were compared with thenumerical results. Afterwards, the numerical model
along with wetting/drying algorithm was utilized tomodel a coupled computation of flow and sedimentfield in a real prototype scale during a free-flow
flushing event. The following outcomes wereobtained from this study:
1Many hydrodynamic aspects such as jet
trajectory, recirculation zones, eddies and the flow
distribution pattern, in different shallow reservoir
geometries were represented by the numericalmodel. The numerical model also reproduced both
symmetric and asymmetric flow pattern insymmetric geometry setup, similar to theobservations, by introducing a perturbation in the
inflow boundary condition. The numericaloutcomes were in good agreement with theexperimental observations. However, there wassome discrepancy for reproducing the upstream
vortex dimension due to the longer concentrated
flow and reverse flow jet pattern with lower flow
diffusion in the numerical outcomes.
2The numerical model could successfully
show the effect of geometry on the flow pattern.
More specifically, the numerical model showedthat the flow pattern in geometries with higher SF(i.e. T13) is insensitive to the small disturbance in
the inflow condition and the numerical model
converged to the steady symmetric flow pattern.
On the contrary, calculated flow pattern ingeometries with lower SF (i.e. T9) converged tothe steady asymmetric flow pattern due to theslight disturbance in the inflow boundary
condition. In prototype scale reservoirs withnarrow shape (i.e. Dashidaira reservoir, Japan),― 479 ―

small disturbance in the inflow condition is very
likely to change the flow pattern and subsequently
different erosion and deposition pattern would
occur. Due to the high probability of disturbancein the inflow condition within a flushing event,accurate modeling of flushing channel size, shapeand location in prototype scale would be difficult.
3Meyer-Peter-Muller bed load transportation
formula was employed for flushing simulation in
Dashidaira reservoir. The numerical model could
simulate the volume of the flushed out sedimentwith an appropriate accuracy. The simulated
volumeofflushedout
sediment
wasunderestimated about 25%. The reason is very low
erosion rate in areas close to the dam and also far
upstream of the reservoir. In areas close to thedam, complex flow field emerges because of thereservoir geometry. Besides, water leveldrawdown procedure is simplified in the
numerical model compared to the reality. Thus,
transport capacity was small in this region.
Regarding the far upstream area, bed material type
is coarse sediment and model could not represent
the erosion pattern properly since the flow field
simulation over such kind of material is complex.
The complexity further magnifies when theaccelerated shallow flow condition emerges. As an
overall conclusion, numerical results disclosedthat secondary currents have been overestimatedwhich result in deeper erosion pattern along theouter part of the channel bending whereas smaller
erosion is obtained for the inner part. Inaccuracies
also could be attributed to the bed forms (i.e.
dunes and anti-dunes) that have a significant
effect on erosion pattern in numerical modeling.
An empirical formula that is basically developed
for small Froude numbers is used to represent the
bed-forms in the model. Nevertheless, numericalmodel can predict an appropriate range for total
volume of eroded sediment. Therefore, it would bepossible to examine the effect of some measureson flushing efficiency for conducting more
efficient sediment management strategies in thisreservoir.
Acknowledgement
The authors would like to give their sincere thanks

to Mr. Shuhei Minami from NEWJEC company forproviding the data of Dashidaira reservoir in thisstudy.
References
Adamson, A., Stovin, V. & Bergdahl, L. (2003):
Bed shear stress boundary condition for storagetank sedimentation. Journal of Environmental
Engineering, ASCE, Vol. (129), No. 7, pp.
651-658.
Atkinson, E., (1996): The feasibility of flushingsediment from reservoirs, Tech. report, HR
Wallingford, UK.
Batuca, D. G. and Jordaan, J.M. (2000): Stilling anddesilting of reservoirs, A. A. Balkema, Rotterdam.
Chandler, K., Gill, D., Maher, B., Macnish, S. andRoads, G. (2003): Cooping with maximum
probable flood-an alliance project delivery for
Wivenhoe Dam, Proceedings of 43 rd conference,
Hobart, Tasmania.
Chiang, T. P., Sheu, T.W.H. & Wang, S. K. (2000):
Side wall effects on the structure of laminar flow
over plane symmetric sudden expansion.
Computers & Fluids, Vol. (29), No. 5, pp.
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, Vol. (117), No. 10, pp. 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, Vol. (8), No. 1,
pp. 31-54.
Dufresne, M., Dewals, B., Erpicum, S.,
Archambeau, P. & Pirotton, M. (2011):
Numerical investigation of flow patterns inrectangular shallow reservoirs, Engineering
Application of Computational Fluid Mechanics
Vol. (5), No. 2, pp. 247-258.
Esmaeili, T., Sumi, T. & Kantoush, S. A. (2014):
Experimental and numerical study of flushingchannel formation in shallow reservoirs, Annual
Journal of Hydraulic Engineering, JSCE, Vol.
(70), No. 4, pp. I19-I24.
Frey, P., Champagne, J. Y., Morel, R. & Gay, B.

― 480 ―

(1993): Hydrodynamics fields and solid particle
transport in a settling tank. Journal of HydraulicResearch, Vol. (31), No. 6, pp. 736-776.
Ghidaoui, M. S. & Kolyshkin, A. A. (1999): Linear
stability of lateral motions in compound open
channels. Journal of Hydraulic Engineering,
ASCE, Vol., (125), No. 8, pp. 871-880.
Haun, S. and Olsen, N. R. B. (2012):
Three-dimensional numerical modeling ofreservoir flushing on a prototype scale, Int. J.
River Basin Management, Vol., (10), No. 4, pp.
341-349.
Haun, S., Dorfman, C., Harb, G. and Olsen, N. R. B.
(2012): 3D numerical modeling of the reservoir
flushing of the Bodendorf reservoir, Austeria,
Proceedings of IAHR European conference,
Munich, Germany.
Kantoush, S. A., Sumi, T. & Schleiss, A. J.
(2010a): Geometry effect on flow and sedimentdeposition pattern in shallow basins. Annual
Journal of Hydraulic Engineering, JSCE, Vol.
(54), pp. 212-232.
Kantoush, S. A., Sumi, T., Suzuki, T., andMurasaki, M. (2010b): Impacts of sedimentflushing on channel evolution and morphologicalprocesses: Case study of the Kurobe River, Japan,
Proceedings of 5 th River Flow Conference,
Germany, Braunschweig.
Kantoush, S. A., Bollaert, E. & Scleiss, A. J.
(2008a): Experimental and numerical modellingof sedimentation in a rectangular shallow basin.
International Journal of Sediment Research, Vol.
23, No. 3, pp. 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 andnumerical
modeling.
Flow
Measurement & Instrumentation, Vol. (19), pp.
139-144.
Kantoush, S. A. (2007): Symmetric or asymmetricflow patterns in shallow rectangular basins withsediment transport. Proceedings of 32 rd IAHR
Congress, Venice, Italy.
Khosronejad, A., Rennie, C.D., Salehi Neyshabouri,
A.A. and Gholami, I. (2008): Three dimensional
numerical modeling of sediment release, Journalof Hydr. Research, Vol. (46), No. 2, pp. 209-223.
Launder, B.E. and Spalding D. B. (1972): Lecturesin Mathematical Models of Turbulence, Academic
Press, London.
Liu, J., Minami, S., Otsuki, H., Liu, B. and Ashida,K. (2004): Prediction of concerted sedimentflushing, Journal of Hydr., Eng., ASCE, Vol.
(130), No. 11, pp. 1089-1096.
Minami, S., Noguchi, K., Otsuki, H., Fukuroi, H.,
Shimahara, N., Mizuta, J. and Takeuchi, M.
(2012): Coordinated sediment flushing and effectverification of fine sediment discharge operationin Kurobe River, Proceedings of the International
symposium on dams for changing world, ICOLD,
Kyoto, Japan.
Morris, G. L. and Fan, J. (1998): Reservoir
Sedimentation
Handbook:
DesignandManagement of Dams, Reservoirs and Watersheds
for Sustainable Use, McGraw-Hill, New York.
Olsen, N. R. B., (2011): A three-dimensional
numerical model for simulation of sediment
movements in water intakes with multiblock
option, www.ntnu.no, Online User's manual.
Patankar, S.V. (1980): Numerical Heat Transfer andFluid Flow, McGraw-Hill, New York.
Public Works Research Institute of Ministry ofConstruction (PWRI) (1993): The property andplanning of the channel in Kurobe river (inJapanese).
Shapira, M., Degani, D. & Weihs, D. (1990):
Stability and existence of multiple solutions forviscous flow in suddenly enlarged channels,
Computers & Fluids, Vol., (18), No. 3, pp.
239-258.
Stovin, V. R. & Saul, A.J. (1996): Efficiency
predictionforstoragechambers
using
computational fluid mechanics. Water Sciences& Technology, Vol. (33), No. 9, pp. 163-170.
Sumi, T., 2006: Reservoir sediment management
measureandnecessary
instrumentation
technologies to support them, The 6 th
Japan-Taiwan Joint Seminar on Natural Hazard
Mitigation.
Van Rijn, L. C. (1984a): Sediment transport. Part I:
Bed load Transport, Journal of HydraulicEngineering, ASCE, Vol. (110), No. 10, pp.
1733-1754.
Van Rijn, L. C. (1984b): Sediment transport. Part
II: Suspended load Transport, Journal of
Hydraulic Engineering, ASCE, Vol. (110), No.

― 481 ―

11, pp. 1431-1456.
Van Rijn, L. C. (1984c): Sediment transport. Part
III: Bed forms and alluvial roughness, Journal of
Hydraulic Engineering, ASCE, Vol. (110), No.
12, pp. 1733-1754.
Yuce, M. I. & Chen, D. (2003): An experimental
investigation of pollutant mixing and trapping inshallow costal re-circulating flows. Proceedingsof International Symposium on Shallow flows,
Balkema, Delft, Netherlands.
(Received June 10, 2014)― 482 ―