Journal of Japan Society of Civil Engineers, Ser.
B1 (HydraulicEngineering),
Vol.
70, No.
4, I_19-I_24,
2014.
Annual Journalof HydraulicEngineering,
JSCE,
Vol.58,
2014, February

EXPERIMENTAL AND NUMERICAL STUDY OF
FLUSHING CHANNEL FORMATION IN
SHALLOW RESERVOIRSTaymaz ESMAEILI1, Tetsuya SUMI2 and Sameh A. KANTOUSH3
1Member of JSCE, Ph.D. Student, DPRI, Kyoto University (Gokasho, Uji, Kyoto 611-0011, Japan)
2Member of JSCE, Professor, DPRI, Kyoto University (Gokasho, Uji, Kyoto 611-0011, Japan)
3Associate Professor, Department of Civil Engineering, German University in Cairo (New Cairo City, Cairo, Egypt)

Sediment flushing is one of the proposed methods for preserving the storage capacity of dam
reservoirs. In flushing with water level drawdown, the incoming flood erodes a flushing channel in thedeposited sediment. Flow pattern as well as flushing channel formation procedure in shallow reservoirs iscomplex phenomenon due to the dynamic interaction between flow field and bed changes. In the present
study, the flow field and flushing channel formation procedure were investigated in various shallow
reservoir geometries using physical experiments and numerical simulation. A fully 3D numerical model
which applies Finite Volume Method (FVM) was utilized. Reasonable agreement was found between the
numerical and experimental outcomes. The results would be useful to understand the influence ofgeometry on flow pattern and flushing process to conduct more efficient sediment management strategies.
Key Words: Shallow reservoirs, 3D numerical simulation, reservoir sedimentation, sediment flushing.

1. 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. This loss of storage volume represents
a huge economic loss and the reduction of flood
control benefits1). Also, the amount of storage loss
varies dramatically from river basin to river basin
due to the different forest cover and geologicalconditions2). The sedimentation problem is serious
for small and medium sized reservoirs with highsediment inflow.
In order to control the reservoir sedimentation,different approaches such as bypassing, dredging,flushing, sluicing and upstream sediment trapping
have been developed. Among several techniques,
the flushing and sluicing plays an important role inthe sediment removal and reduction, as they are
efficient hydraulic sediment removal technique torestore the reservoir storage capacity1),3).The successof flushing depends on the sediment characteristics,
the hydraulic and hydrological conditions1)~3).

Sediment flushing with drawdown is ongoing at
Dashidaira & Unazuki dam reservoirs in the Kurobe
River as pioneer cases in Japan.
Sediment flushing in reservoirs involves several
complex processes. During the flushing withdrawdown, bottom outlets are opened to generate
and accelerate unsteady flow towards the outlet.
This process will initiate the progressive andretrogressive erosion pattern in tail and delta reachesof the reservoirs respectively4). The mentioned
process finally leads to the flushing channel
formation and flushing out both fine and coarse
sediments through reservoirs. Detailed theoreticalexplanation of flushing channel formation inreservoir delta is scarce. The characteristic of theflushing channel when the water level is drawndown
can be briefly introduced by the location, shape,
width, side and longitudinal slope5). For practical
purposes, the pre-assessment of the flushing channel
development characteristic using a 3D numericalmodel would be beneficial to optimize the sedimentflushing operation as well as the flood risk
management in reservoirs.
Application of numerical models enables thereservoir owners to assess the impacts of upcoming

I_19flushing event in the reservoir and implement theappropriate measures aimed to enhance the flushingefficiency. When it is necessary to assess the
distribution of shallow parts, deep parts and sandbar
areas for reservoir sediment management strategies,
numerical models other than one-dimensional one
should be used6). Many calculation methods basedon 2D models have been developed for practicalproblems in rivers7). However, the one and two
dimensional models are not able to directly simulate
the secondary current influences. Also, quasi-3D
models can be categorized whether the shallow
water flow assumed or not8). Still it is difficult foradvanced quasi-3D models to reproduce complex
3D flows9).
Subsequently, 3D numerical models are necessary
to simulate the complex 3D flow pattern and bed
variation in rivers and reservoirs. In the present
study, a fully 3D numerical model was employed toanalyze the surface flow field in a series ofrectangular shallow reservoirs and then the flushingchannel formation and evolution trend were
simulated in the aforementioned reservoirs. The
next step would be simulating the sediment flushingprocess in a real prototype scale to enhance theflushing efficiency by employing various measures.

material was 50 µm with σg of 2.4 and density of1500 kg/m³. The Froude number range set to (0.05 ≤Fr ≤ 0.43) whereas the Reynolds number range was(14000 ≤ Re ≤ 28000) to ensure the subcritical andfully developed turbulent flow. The flow discharge
rate (Q), suspended sediment concentration (C) andwater depth (h) were constant for all experiments as
0.007m³/s, 3gr/L and 0.2m respectively.
The test procedure had three different stages. In
the first stage, the shallow reservoirs was filled withclear water and after reaching to the stable state, inthe second phase, the mixture of water-sediment
was drained by the gravity into the water-filled
reservoir for a total period of 4.5 hours. In the third
phase, two types of flushing with and without
drawdown were performed. The final bed
topography from the second phase was consideredas the initial bed for two types of flushing.
Afterwards, the clear water inflow introduced into
the reservoir to evaluate the surface velocity pattern,
flushing channel evolution, channel location as
well as the flushing efficiency. As for the flushing
without drawdown in run T1, the clear water
without sediment was injected into the reservoir
with the constant hydraulic condition (Q=0.007m³/s
and h=0.2m). The drawdown flushing was
conducted by opening the outlet gate and lowering
the water level to half of the initial level for runs T8,
T11 and T13. Then, pump was turned on by keeping
a constant discharge (Q=0.007m³/s and h=0.1m).
Table 1 shows geometrical attributes of various
geometries employed in the present study.
It was found that the flushing channel formation
is a very rapid procedure that takes less than 3
hours.

2. MATERIALS AND METHODS(1) Experimental setup
The experimental tests were carried out in arectangular basin with the maximum inner length of
6m and width of 4m. Also, the inlet and outlet
rectangular channel width was 0.25m. The different
shallow reservoir geometry achieved experimentallyby moving the PVC plate walls. A moveable framewith 4m long was mounted on the side walls ofreservoir for installing the measurement devices.
Ultrasonic probes were utilized for measuring thewater level and Large Scale Particle Image
Velocimetry technique (LSPIV) was used formeasuring the surface velocity field. Ultra Sonic
Velocity Profiler device (UVP) was employed forproviding the 3D flow velocity measurements aswell. Bed topography was measured by a miniatureecho sounder which was mounted on the moveable
frame and could scan the whole geometry domain.
Moreover, two SOLITAX sc sensors which utilize
the ultrasonic approach were installed for measuring
the suspended sediment concentration.

(3) Numerical model
A fully 3D numerical model was employed in thisstudy. The numerical model solves the
Reynolds-averaged Navier-Stokes equation together
with mass and momentum conservation in three
dimensions to compute the water motion for
turbulent flow as follows.
U i
0
xi

(1)

U j 
U
U i 1  U i P ij   T ( i U j)  (2)
x j
xi 
x j
 x j tin which i=1, 2, 3 is the representative of three
directions; whereis the Reynolds-averaged
velocity over time t,
is the spatial geometricalscale,
is the water density,
is theReynolds-averaged pressure,is the Kronecker
delta and ν T is the turbulent eddy-viscosity. For
transforming the partial equations into algebraic(2) Experimental conditions
The non-uniform crushed walnut shells wereadded to the mixing tank during the test to represent
the suspended sediment. The median size of thisnon-cohesive light-weight and homogenous grain
2

I_20

Table 1 Experimental configuration. L and B are the length andwidth, Pr is the wetted perimeter, A is the surface area of thebasin, AR equals to (L/B) and SF is the shape factor.
Case

L(m)

B(m)

Pr(m)

T1
T8
T 11
T 13

6
6
5
3

4
2
4
4

19.5
15.5
17.5
13.5

A(m²)
24
12
20
12

AR
1.5
3
1.25
0.75

where , is sediment transportation rate for theis theth fraction of bed load per unit width,
diameter of the ith fraction, is the shear stress,
which was
, is the critical shear stress foris the density
calculated from the Shield's curve,
is the density of the water, g is theof sediment,
gravity acceleration and ν is the kinematic viscosity.

SF=(Pr/√A)AR

5.97
13.42
4.89
2.92

equations, the finite volume method is applied as
discretization method. The change in water-levels
was based on calculated pressure field. The pressure
was extrapolated to the water surface and thepressure difference between a surface node and thedownstream node was used to estimate the water
elevation difference 10). The turbulence is modeledby the standard k-ε model, using the constantempirical values11). The unknown pressure field isalso calculated employing Semi Implicit Method forPressure-Linked Equations, (SIMPLE) method12).
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 assumed. As
for the boundary condition at the walls, where there
is no water flux, the empirical wall laws were used
as follows:
U 1  30 y  lnu *   ks T (1) Flow pattern simulation
The complexity of velocity measurement withhigh spatial resolution in different shallow reservoir
geometries pronounces the necessity to conduct thenumerical
modelingalongwithphysical
experiments. As for the real cases, modeling theflow field will provide us useful information about
the areas with the potential erosion and deposition
during the anticipated floods. This information
would be useful for flood risk assessment inreservoirs near the urban areas.
The mesh cell size for case T1, T8, T11 and T13
in X and Y direction was 5cm  2.5cm, 5cm 
1.5cm, 5cm  2cm and 2.5cm  1cm respectively.
Considering the 11 cells for vertical grid distribution,
the total number of cells were 218240, 174460,
220000 and 528000 respectively. For model
validation, the simulated surface flow velocity field
was compared to that of measured experimentally.
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. Timestep was calibrated as 2 seconds for run T1 and T8
whereas it was 1 second for T11 and T13. Bedroughness was fixed as 0.00015m which equals to 3
times of median sediment size. Fig. 1 illustrates thetransversal flow velocity (Uy) contours over thedepth along with the secondary flow velocity vectors
for case T8 and T13. Fig. 2 shows the simulated andmeasured surface velocity magnitudes (V) in m/s
and distribution pattern with uniform scaled vectors.
Fig. 3 illustrates the simulated streamwise andtransversal surface velocity distribution versus

(3)

di

w

 0.053di0.3

 (  s   w ) g 


2

  w(6)

3. RESULTS AND DISCUSSIONS

Where
is the water velocity,is the fallvelocity of sediments, Γ is the turbulent diffusivity
and can be expressed by (Eq. 6) andis thesediment concentration over time t within the spatial
geometrical scales x and z.
1. 5
   c ,i q b ,i  c ,i (5)( s   w )gSc

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, aspecified concentration was used as boundary
condition14).

,  is thewhere the shear velocity is denoted
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 is
calculatedbysolvingthetransientconvection-diffusion equation formula (Eq. 4) andbed load is simulated by Van Rijn formula13) (Eq. 5).
c
c
c

c
(4)
U j
w 
(T
)
t
x j
z x j
x j

1. 5

T

0. 1

3

I_21

4

Velocity magnitudes (m/s)

3

3.5
3

Y(m)

Z

X=2m_simulated
stremwise velocity

2.5

X=2m_measured
streamwise velocity

2
1.5

X=2m_simulatedtransversal velocity

1

Y

0
-0.10

-0.05

(a)

0.00

0.05

0.10

0.15

1
-0.25

0.20

-0.15

-0.05

X=3m_measured
streamwise velocity

2
1
0.5
0
-0.15

(a2)

-0.05

0.05

0.15

X=3m_measured
surface velocity

2

X=3m_simulatedtransversal velocity

1.5

X=3m_measured
transversal velocity

0.5

0.25

X=3m_simulated
surface velocity

2.5

X=3m_simulatedtransversal velocity

1

X=3m_measured
transversal velocity

0
-0.15

0.35

-0.05

(b2)

(c1)

(c2)

Y

X

(d1)

0.05

0.15

0.25

0.35

Surface velocity (m/s)(2) Simulation of flushing with drawdown
Numerical outcomes as well as the observations
in experimental runs revealed that if the water levelis drawndown significantly, the flow starts to erode
the bed progressively propagating to thedownstream and retrogressively from the outlet
towards upstream. The progressive trend was fasterthan retrogressive one. In the meantime, the initialflushing channel deepened and widened rapidly due
to the strong jet flow and subsequent erosion. Then,after formation of initial flushing channel along thereservoir length, the rate of channel widening
reduced noticeably until reaching a dynamic stablecondition over the whole channel length. The
mentioned process is very quick up to the slow
widening stage so that the measurement of the bed
evolution was difficult. Subsequently, numericalmodel outcomes along with the final measured
results in the equilibrium stage utilized to further
description of the process. Fig. 4 demonstrates the

V:

(b1)

0.35

(c)
(d)
Fig. 3 The measured streamwise and transversal surface
velocity at the middle area of the channel versus simulated
surface velocity for (a) T1, (b) T8, (c) T11 and (d) T13.

V:

V:

0.25

3
Y(m)

2.5

Surface velocity (m/s)

(a1)

0.15

4
X=3m_simulated
streamwise velocity

1.5

V:

0.05

Surface velocity (m/s)

3.5

3
Y(m)

X=3m_simulatedtransversal velocity
X=3m_measured
transversal velocity(b)Surface velocity (m/s)

4Fig. 1 Transversal velocity distribution over the depth at the
middle cross section of (a) T8 and (b) T13.

X=3m_measured
streamwise velocity

1.5

3.5

(b)

X=3m_simulated
streamwise velocity

2

X=2m_measured
transversal velocity

0.5

(a)

2.5
Y(m)

Uy:

(d2)

Fig. 2 Left: The measured surface velocity magnitudes (V) andvelocity vectors for runs (a1) T1, (b1) T8, (c1) T11 and (d1) T13
respectively and right: corresponding simulated velocity
magnitudes and velocity vectors.

Bed elevation (m)0.05
t=Initial
t=15 mint=30 min
t=40 min
t=80 min
t=120 mint=180 min
t=480 min

0.04
0.03
0.02
0.01
0

measured one in the middle part of the channel
length. As can be clearly observed form the Fig. 2
and Fig. 3, the model could simulate the surface
flow velocity pattern almost similar to the measured
one by reproducing the dominant aspects such as themain flow jet trajectory and location of the reverse
flow as well as the main vortices and corner gyres.
However, simulation results show the higher flow
velocity magnitudes than observations over the main
jet trajectory for all cases and over reverse flow incase of T1 and T8. Also, the size of the upstream
corner gyres in case T11 is smaller than themeasured one because of the straighter simulated
surface velocity pattern than hydraulic model test.

0

(a)

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

6

Reservoir length (m)Bed elevation (m)

0.05
0.04
0.03
0.02
0.01
0

(b)

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

5.5

Reservoir length (m)Bed elevation (m)0.05

Initial

0.04

t=5 min

0.03

t=10 min0.02

t=15 min0.01

t=20 min0(c)t=Initial
t=15 mint=30 min
t=40 min
t=60 min
t=80 min
t=120 min
t=480 mint=25 min
1.5

2

2.5

3

3.5

4

4.5

Reservoir length (m)Fig. 4 Longitudinal development of flushing channel along thecenterline for (a) T8, (b) T11 and (c) T13 provided by numericalmodel.

4

I_22

longitudinal progressive and retrogressive pattern,
reproduced numerically, along the centerline of T8,
T11 and T13 reservoirs. Numerical computation
was performed until starting the slow channel
widening stage. Thus, as soon as cumulative
sediment pass variation within 30 minutes interval
became smaller than 1%, simulations were stopped.
Described flushing channel formation process wasconsistent with the real observations during thedrawdown flushing in Dashidaira and Unazuki dam
reservoirs. The side bank erosion and lateraldevelopment of flushing channel in the middle
length of the reservoirs has been shown in Fig. 5.
Fig. 6 illustrates the plan view of final measured
bed level contours (Z) after 48 hours and simulated
one after initiating the slow channel wideningstage.The size, shape, location and evolution patternof the flushing channel was simulated up to theearly stage of slow widening phase. It takes longtime to simulate whole 48 hours by 3D numericalmodel. The characteristics of flushing channel
components, more or less, have been reproduced by
the numerical model except for run T8 in which thelocation of the channel was different than measured
one. Although the geometry was symmetric, thechannel did not develop along the shortest path frominlet to outlet during the experiment. The reason canbe attributed to the accidental small disturbance in
the inflow discharge distribution along the inlet
channel width. In such a condition, location of theflushing channel is very likely to change from thecenterline of the reservoir. As the numerical model
assumes the uniform inflow discharge distribution
and symmetric initial condition, the flow direction
would be straight and subsequently model is notable to break the symmetry of input data.
As to the run T11 and T13, both experimental
measurements and numerical outputs show that thechannel width increased in the downstream direction

similar to a T shape head (Fig. 6(b2) & (c2). The
longer simulation period, the similar channel width
as well as T head shape will be reproduced close to
the outlet.
(3) Flushing efficiency
In this study the Flushing Efficiency (FE) was
defined as the volume ratio of flushed out sedimentto cumulative deposited sediment after thedeposition (second) phase. Fig. 7 illustrates the
temporal sediment discharge change in milliliter per
second and cumulative sediment pass through theoutlet for run T8, T11 and T13. The erosion rate isvery high at early stage owing to the rapid
deepening and widening of flushing channel. Then,the sediment discharge rate decreases remarkably as
the initial flushing channel was formed along thewhole length of the reservoir. Afterwards, during the
Z:

(c1)

Bed elevation (m)

1

1.25

1.5

1.75

2

2.25

2.5

2.75

3

Reservoir width (m)Fig. 6 Left: The measured flushing channel shape and locationfor (a1) T8, (b1) T11 and (c1) T13 respectively and right: thecorresponding simulated results.

t=Initial
t=15 mint=30 min
t=40 min
t=80 min
t=120 mint=180 min
t=480 minEquilibrium

25

0.03
0.02
0.01
0

(b)

0

0.5

1

1.5

2

2.5

3

3.5

4

Reservoir width (m)t=Initial
t=15 mint=30 min
t=40 min
t=60 min
t=80 min
t=120 min
t=480 minEquilibrium

0.02
0.01
0

(c)

0

0.5

1

1.5

2

2.5

3

3.5

4

Sediment
discharge

10

Comulative
sediment pass

5
0
0

30

60

90

120 150 180 210 240 270 300 330 360 390 420 450 480

Time (min)30

Sediment
discharge

20

Comulative
sediment pass

100(b) 0 30 60 90 120 150 180 210Time240(min)270 300 330 360 390 420 450 480

Initialt=15 mint=30 min
t=60 mint=90 min
t=120 mint=180 minEquilibrium

50

Sediment discharge(ml/s)

Bed elevation (m)

0.03

15

40

0.04

Reservoir width (m)

140
120
100
80
60
40
20
0

20

(a)

0.05
0.04

X

40
Sediment
discharge

30
20

Cummulative
sediment pass

10

(c)

Cumulative sedimentpass (kg)

0

(a)

(c2)

0
0

20

40

60

80

100

120

140

160

140
120
100
80
60
40
20
0

Cumulative sedimentpass (kg)

0.01

(b2)

70
60
50
40
30
20
10
0

Cumulative sedimentpass (kg)

0.02

(b1)

Sediment discharge(ml/s)

0.03(a2)

Sediment discharge(ml/s)

Bed elevation (m)

0.04

(a1)

Y

0.06
0.05

Z:

180

Time (min)Fig. 7 Calculated flushed out sediment for the run (a) T8, (b)T11 and (c) T13.

Fig. 5 Bank erosion and lateral development of flushing channel
for (a) T8, (b) T11 and (c) T13.

5

I_23

formation and evolution. In terms of the flushingchannel location and Flushing Efficiency (FE), theinfluence of reservoir width is higher than reservoir
length. Thus, higher Shape Factor (SF) leads tohigher efficiency of flushing with drawdown. In thecase of numerical simulation of flushing withdrawdown, the efficiency can be further increased if
the non-symmetric inflow pattern, as a kind of
perturbation, can be included in the numericalmodel to initiate wide meandering channel form.

slow channel widening stage, the sediment
discharge rate remains almost stable. The measured
and calculated FE has been plotted versus reservoir
shape factor in Fig. 8.
In case of drawdown flushing, FE increases
significantly with the higher reservoir shape factor
in experimental runs. In other word, the flushingefficiency with drawdown will be high for narrow
reservoir geometries. On the contrary, the FE influshing without drawdown is low because of very
local erosion pattern in the inflow and outflow
area.Right illustration in Fig.8 clearly shows local
erosion during the flushing without drawdown.
The numerical model outputs showed slightly
lower FE rate than measured ones in Case T11 andT13. The main reason is due to the computation stop
as soon as slow widening stage starts. Nevertheless,in case of T8, measured and numerically modeled
FE has great discrepancy. This is because of thewide meandering flushing channel formation in theexperiments whereas the width and length of thesimulated straight flushing channel is smaller thanthat of the measured one (Fig. 6(a1) & Fig.
6(a2)).Subsequently, the volume of the eroded
sediment as well as the flushing efficiency would behigher in the experimental model of T8.

REFERENCES
1)

2)
3)

4)
5)

4. CONCLUSION

6)

In the present study, a fully three-dimensional
numerical model was utilized for reproducing thesurface velocity distribution pattern and flushingchannel formation in shallow reservoirs. Main
outcomes of the study were as follows:(a) Many hydrodynamic aspects such as jet
trajectory, recirculation zones, eddies and the flow
distribution pattern, in different shallow reservoir
geometries were presented quantitatively by thenumerical model. Application of the hydrodynamic
module of the numerical model would be beneficial
for flood risk management in reservoir areas.
(b) As for the flushing channel formation, acomplex dynamic interaction exists between various
effective parameters. However, three-dimensional
numerical model could successfully reproduce themain governing procedures in flushing channel

Flushing efficiency (%)

100

7)

8)

9)

10)

11)

Calculated FE with darwdown
Measured FE with drawdown
Measured FE without drawdown

80

12)

60

13)

40
20
0
2

14)

6
10
14
Reservoir shape factor (SF)

Fig. 8 Left: Measured and calculated FE with drawdown forvarious reservoir geometries. Right: Illustration of local erosion
obtained during the flushing without drawdown.

Morris, G. L., and Fan, J.: Reservoir Sedimentation
Handbook: Design and Management of Dams, Reservoirs
and Watersheds for Sustainable Use, McGraw-Hill, New
York, 1998.
White, R.: 2. Review of sedimentation in reservoirs,
Thomas Telford, London, 2001.
Liu, J., Minami, S., Otsuki, H., Liu, B., and Ashida, K.:
Prediction of Concerted Sediment Flushing, Journal of
Hydraulic Engineering, ASCE, Vol. 130, No. 11, pp.
1089-1096, 2004.
Batuca, D. G. and Jordaan, J.M.: Stilling and desilting ofreservoirs, A. A. Balkema, Rotterdam, 2000.
Kantoush, S. A., Sumi, T., Suzuki, T., and Murasaki, M.:
Impacts of sediment flushing on channel evolution andmorphological processes: Case study of the Kurobe River,Japan, Proc.5th River Flow Conference, Germany,
Braunschweig, 2010.
Fukuoka, S., Sumi, T., and Horiuchi, S.: SedimentManagement on the Arase Dam Removal Project,
Proc.12th
International
SymposiumonRiver
Sedimentation, Japan, Kyoto, 2013.
Shimizu, Y., Itakura, T., and Yamaguchi, H.: Numericalsimulation of bed topography of river channels using
two-dimensional model, Proceedings of the Japanese
Conference on Hydraulics, Vol. 31, pp. 689-694, 1987.
Uchida, T. and Fukuoka, S.: A computation method forflow over structures, Advances in River Engineering, Vol.
18, pp. 351-356, 2012.
Fukuoka, S. and Uchida, T.: Toward Integrated
Multi-Scale Simulations of Flow and Sediment Transportin Rivers, Journal of JSCE, Ser. B1 (HydraulicEngineering), Vol. 69, No. 4, pp. II_1-II_10, 2013
Olsen, N. R. B.: A Three-Dimensional Numerical Model
For Simulation of Sediment Movements In Water Intakes
With Multi-block Option, Department of Hydraulic andEnvironmental Engineering, The Norwegian University of
Science and Technology, Trondheim, 2012.
Launder, B.E. and Spalding D. B.: Lectures inMathematical Models of Turbulence, Academic Press,
London, 1972.
Patankar, S.V.: Numerical Heat Transfer and Fluid Flow,
McGraw-Hill, New York, 1980.
Van Rijn, L. C.: Sediment Transport. Part II: Suspended
load Transport, Journal of Hydraulic Engineering, ASCE,
Vol. 110, No. 11, pp. 1431-1456.
Van Rijn, L. C.: Sediment Transport. Part I: Bed load
Transport, Journal of Hydraulic Engineering, ASCE, Vol.
110, No. 10, pp. 1733-1754.
(Received September 30, 2013)6

I_24