京 都 大 学 防 災 研 究 所 年 報 第 50 号 B 平 成 19 年 4 月
Annuals of Disas. Prev. Res. Inst., Kyoto Univ., No. 50 B, 2007

Morphodynamics of Channels with Groins and its Application in River Restoration

Hao ZHANG, Hajime NAKAGAWA, Yasunori MUTO, Yoshio MURAMOTO*,
Daizaburou TOUCHI**and Yasunori NANBU*** Professor Emeritus of Kyoto University** Graduate School of Engineering, Kyoto University
Synopsis
This study investigates the flow and bed deformation process in a channel with groins in ariver restoration project using both physical model experiment and numerical simulation methods.
The experiment is based on a large-scale physical model and is able to resolve the local flow andsediment transport phenomena with a relatively high accuracy. The numerical model is formulatedon unstructured meshes and allows the exact representation of complex geometries commonly
encountered in the engineering practice. The numerical results have been compared with those of
the experimental measurements. Reasonable agreements have been obtained in terms of both flow
patterns and bed variations.
Keywords: channel morphodynamics, groin, river restoration, large-scale physical model, complex
geometry, unstructured mesh

1.

Introduction

In-stream structures such as groins have been widelyused in the river engineering to restore channel
morphologies and improve aquatic habitats. Comparedwith other restoration structures such as bank covers andweirs, groins have been documented to be the mostdurable solution (Shields, et al., 1995). Consequently,
investigations on the flow and/or bed deformation aroundgroins have been carried out, in the past several decades,by a lot of researchers with field surveys, laboratorial
experiments and numerical simulations, e.g. Kuhnle et al.
(2002), Muto et al. (2003), Sukhodolov et al. (2004) andNagata et al. (2005). These studies provide a substantial
collection of data and information on the flow andsediment transport around this typical kind of in-stream
structures. However, existing research results are usuallyproblem-dependent and have not been developed yet up
to practical knowledge for the design of restoration
projects. Although there do exist some design
suggestions and recommendations based on one orseveral researches for river restoration projects, thecomplexities of the real river conditions have been
seldom taken into account in a cost-effective way.
Approaches based on trial-and-error and limited
experiences still play an important role in thestate-of-the-art for restoration designs and managementin practice. Without a full understanding of the channel
hydraulics and morphodynamics caused by restoration
measures, many projects have failed in the early stage of
their design working life or have been confirmed to beunsustainable. The low success rates usually implicate aloss or a require of a huge amount of investments as wellas potential disasters for the riverine eco-system andpeople living nearby. Therefore, it is necessary todevelop more general design methods based on an
insight into the underlying physics of the problem.
This study investigates the channel morphodynamics
before/after the installation of groins and seeks
practicable design methods for river restoration. In orderto fulfill these tasks, a real restoration project in west
Japan has been selected as a reference. The field site of－ 635 －

the project is in the Yodo River- Lake Biwa basin in west
Japan and located at the Kizu River about 6km upstream
from the confluence where three major rivers, i.e. the Uji
River, the Katsura River and the Kizu River, joint
together (Fig. 1).
Lake Biwa

N

Shiga Prefecture
Kyoto Prefecture
Katsura River
Uji River

Kizu River
Study domain

Mie Prefecture

Yodo River
Osaka Prefecture
Nara Prefecture
Osaka Bay

Fig.1 Yodo River-Lake Biwa basin and the study domain
Impermeable groin

Nagare Bridge
5

Tower foundations
Low flow channel

Flow
4

3

2

1

Old river terrace
Groins destroyed in flood
Permeable groinFriparian restoration. Up to now, the project has beenworking well. However, the knowledge is still very pooron the river morphodynamics caused by these restoration
actions. Intensive studies and long-term monitoring arestill of great necessity.
Both large-scale physical model experiments and 2D
numerical simulations have been carried out. The
physical model is non-distorted and represents the detailsof local hydraulic structures to scale in the area of interest,
which is able to resolve the local flow and sedimenttransport phenomena with a relatively high accuracy.
The numerical model solves the depth-averaged shallow
water equations for the mean flow velocities considering
the effect of secondary flows. The bed deformation is
obtained by solving the sediment continuity equation in
the bedload layer and the sediment transport rate isevaluated with a deterministic empirical formula. The
influence of the local bed slope on bedload transport has
been taken into account, as well as the deviation of thenear-bed flow direction from that of the mean flow.
Furthermore, the 2D model is formulated on an
unstructured mesh, which allows the exact representationof complex geometries such as the irregular bank-linesand groin structures. Selected numerical results arecompared with those of the experimental measurements.
Discussions and conclusions are made based on typicalvelocity fields and bed deformations under various
conditions such as the changing of discharges andinstalling/removing groins.

River terrace
River levee2.

Fig.2 Sketch of the restoration project

The project was built around an ancient wooden
bridge: the Nagare Bridge, which is one of the most
attractive sightseeing and recreation spots in this river
basin. The main objectives are to restore beachlandscapes and improve the riparian ecotones around thebridge. Restoration measures mainly include cutting
down part of the river terrace and constructing 3 groins(Ԙ, ԙ and Ԛ) as shown in Fig. 2. The project was not
successful. It was completed in 2004 but soon failed in a
typhoon-induced flood. The 3 groins (one permeable andtwo impermeable) were almost washed away. As aremedy, a permeable groin (ԛ) and an impermeable
groin (Ԝ) were introduced in the next year after a quick
investigation of the effects of the two groins on the

Physical model experiments

The physical model was constructed in the Ujigawa
Open Laboratory, Kyoto University (Fig. 3). The model
had a physical scale of 1:65 and represented the river
reach from section 4.4km to section 7.4km (here km
stands for the distance in kilometer from the riverconfluence as shown in Fig. 1). Local hydraulicstructures in the neighborhood of the Nagare Bridge havebeen exactly represented, including the piers of thebridge, different kinds of groin structures, etc.
Along the transverse section, the river consisted of alow flow channel with a floodplain on either side. The
floodplain area was treated as fixed bed and artificial
grass was made in some area to model the vegetated areasuch as the tea fields. The low flow channel was movable
and was covered with powdered anthracite after acomprehensive comparison with some other model－ 636 －

sediment materials. The model sediment had a mean
diameter of 0.83mm and a specific gravity of 1.41. The
sieve analysis result of the sediment samples taken at
typical cross-sections on site and the model sediment isshown in Fig. 4.
Sandbar

Nagare Bridge
6.0km
Instrument carriage

Low flow channel
㪧㪅㪈
㪧㪅㪉

7.2km
㪧㪅㪊
㪧㪅㪋

4.8km

㪧㪅㪌

Groin area
Sediment supplierFig.4 Sieve analysis result of sediment particles(a) General view
Floodplain
Impermeable groin

Permeable groin

Flow

Nagare Bridge

4 piles
6.0km

S2

S1 6.2km(b) Groin area

Side view
46 piles72 piers

Longitudinal profile

Longitudinal profile(c) Permeable groin(d) Piers of the Nagare Bridge

Typical side viewFig.3 Physical model experimental setup

Continuous sediment supply was guaranteed during
the experiments with an automatic sediment supplierinstalled at the inlet boundary. At the outlet boundary, atailgate was used to control the water level. At section
4.6km, a point gauge was set to record the downstream
water stage. Water level and bed elevation were
measured at typical transverse sections during theexperiments. In the proximity of the groins, thelongitudinal distance between two measured sections
was around 10cm (about 6.5m in prototype). The waterlevel was measured with an ultrasonic level system(10s-sampling for a point at a frequency of 100Hz) andthe bed elevation was obtained through a sand-surface
profiler. The surface flow structure in the groin stretch
was studied with PIV (Particle Image Velocimetry)
method (30s-sampling for a point at a frequency of 6Hz).

The experimental studies last some 2 years and wereconducted in 4 phases one by one:
Phase 1: Preliminary experiments to quantify some
important hydraulic parameters such as sedimenttransport rate and bed roughness in the model channel.
Detailed reports were given in Touchi (2006) and Zhang
et al. (2006). It has to be mentioned that the sedimenttransport rates (mainly in bedload mode) to maintain the
dynamic equilibrium bed conditions under representative
water discharges have been investigated. According to
the experimental results, the sediment transport rates
under equilibrium conditions were about 70% of thoseestimated with the widely used Ashida-Michiue formula(Ashida and Michiue, 1972) for Japanese rivers.
Phase 2: Experiments to test the performance of themodel channel, i.e. to what extent, the experimental
channel might reproduce the flow and sediment transportphenomena in the real river. A 7-year fluvial process anda one-year bed deformation process (considering theeffects of the detailed flood process) have been
investigated, respectively. The experimental results havebeen compared with the field data and reasonableagreements have been obtained. Detailed information
was reported in Nanbu (2007).
Phase 3: Experiments to investigate the effects ofthe existing restoration project on channel morphologies.
Under the same scenarios for flow and sediment
discharges, experiments with groins and without groins
were carried out, respectively. Based on the experimental
data, together with corresponding numerical simulation
results, the performance of the existing groin structures
was assessed. It was found that the current restoration
measures seemed not enough to fulfill the ultimate goals.
Therefore, some actions should be taken to improve the－ 637 －

effectiveness of existing groins. Possible improvement
alternatives have been proposed after carefully analyses
of both experimental and numerical results.
Phase 4: Experiments to investigate and compare
different morphological consequences due to different
improvement methods. A strong emphasis has been laidon the local bed deformation in the neighborhood of the
Nagare Bridge. Numerical simulations were conducted at
the same time with the models described below. The bestsolution was finally selected.
3.

Wbx Ucf u u2 v2(4)

Wby Ucf v u2 v2(5)

in which the bottom friction coefficient cf is related to the
Manning’s roughness coefficient n by
cf

n2g

(6)

h1/3

From Eq. 4 and Eq. 5, the near-bed shear velocity isimmediately obtained as

Numerical models

A 2D morphological model has been developed tosimulate the flow and sediment transport in this study.
The model consists of a depth-averaged flow module, abedload transport module and a bed deformation module.
The flow field, sediment transport and bed deformationin the model are coupled and solved in a quasi-steadyway, i.e. when one of them is under calculation, theothers are considered as unchanged.
3.1 Hydrodynamic module
Integrating the 3D (Three-dimensional) RANS
(Reynolds-averaged Navier-Stokes) equations from theriverbed to the free surface, the following equations
governing the 2D open channel flow can be derived.
wH w
w

hu 
hv
wt wx
wy

omitted so far when Eq.2 and Eq.3 are solved, which
results in the following expressions:

0

(1)

u*

º Wby
wv wv wv
wH 1 ª ww
u v g  « hWyx  hWyy » Dyx Dyy (3)wt wx wy
wy Uh¬wx
wy¼ Uh
where t= time; H = water stage (i.e. H = h + zb); h =water depth; zb = bed elevation; x, y= Cartesian
coordinate components; u, v = depth-averaged flow
velocity components in the x and y directions,
respectively; g= gravitational acceleration; ȡ = fluid
density; Ĳxx, Ĳxy, Ĳyx, Ĳyy = depth-averaged turbulent
stresses; Dxx, Dxy, Dyx, Dyy = dispersion terms due tovertical non-uniformities of the mean flow and Ĳbx, Ĳby =bottom shear stresses. Strictly speaking, the bottom shear
stresses are related to the near-bed flow field rather thanthe mean flow field. Nevertheless, the influence due tothe deviation of the near-bed flow from the mean flow is(7)U

It has to be mentioned that the shear stress and shear
velocity near the bed are corrected by replacing the mean
velocities in Eq. 4 and Eq. 5 with the near-bed velocities
when sediment transport rate and direction are estimated.
Approximating the turbulent stresses with
Boussinesq’s eddy viscosity concept, one has§ wui wuj · 2 ¸  kG¨ wxj wxi ¸ 3 ij©
¹

Wij UXt ¨i, j x, y(8)

where ȣt = eddy viscosity, k = depth-averaged turbulencekinetic energy and įij = Kronecker delta.
In a zero-equation model, the eddy viscosity is
estimated throughXtº W
wu wu wu wH 1 ª w
w
u v g  « hWxx  hWxy » bx Dxx Dxy (2)wt wx wy
wx Uh¬wx
wy¼ Uh

Wbx2 Wby2

N

u*h
(9)
6
where ț = van Karman constant (=0.41). The result of
the zero-equation model is used as an input for moreadvanced models such as a k-İ model. In the k-İ model,
the eddy viscosity is a function of the turbulent kinetic
energy k and its dissipation rate İ, i.e.

Xt cP

k2

H

(10)

where cȝ = coefficient and is usually set to be a constant
0.09. The two quantities k and İ are solved from theirtransport equations as follows, respectively.

－ 638 －

wk wk wk w § Xt wk · w § Xt wk ·
u v
¨
¸ ¨
¸  Ph  Pkv H (11)wt wx wy wx ©Vk wx ¹ wy © Vk wy ¹
wH wH wH w § Xt wH · w § Xt wH · H
H2
u v
¨
¸ ¨
¸c1H Ph PHv c2H (12)wt wx wy wx©VH wx ¹ wy©VH wy ¹ kkthe secondary flow (=7.0 for fully developed bend flows,
Engelund, 1974) and r = radius of the curvature of thestreamline. According to Shimizu and Itakura (1991), it
is estimated from

where
2

2º

ª § wu ·2 § wv · § wu wv ·Ph Xt «2¨ ¸  2¨ ¸ ¨  ¸ »« © wx ¹ © wy ¹ © wy wx ¹ »
¬
¼
Pkv

u*3
hc1/2
f

PHv

c2H c1/2P

u*4
1/2 3/4 2
e*Vt cf h

in which ıt = Schmidt number expressing the relation
between the eddy viscosity and the diffusivity for scalartransport (a value of 1.0 is adopted in this study) and e* =adjustable dimensionless diffusivity. The model
parameters generally take the following values assuggested by Rodi (1980).

Vk 1.0 VH 1.3

C1H 1.44

C2H 1.92(15)hV  NU(16)*rwhere aU = coefficient considering the primary velocity
profile in the vertical direction (=1.0 in this study forsimplicity); N* = coefficient representing the intensity of1u2 v2ª § wv wu · § wv wu ·ºu u v ¸ v¨u v ¸»«¬ © wx wx ¹ © wy wy ¹»¼

3/2 « ¨(17)

3.2 Sediment transport module
In the absence of suspended load, the bed variation isdue to the movement of bedload transport in a thin layer
near the riverbed. The sediment continuity equation in
the bedload layer is expressed by(14)

During the depth-averaging of the 3D RANS
equations, the dispersion terms Dxx, Dxy, Dyx, Dyy appear
as shown in Eq.2 and Eq.3 (Rodi, 1980). These terms
exhibit similar physical effect with the turbulent stresses
and can be significant when the secondary flow is
obvious. Minh Duc et al. (2004) suggested to adjust thedimensionless diffusivity coefficient e* in the k-İtransport equations to indirectly account for secondary
flow effects according to their numerical experiments. A
value of 0.15 has been suggested for laboratory flumes.
In natural rivers, a value from 0.6 up to 10 might be used.
Following their idea, the dimensionless diffusivity istreated as an adjustable parameter to include the effect of
the secondary flow on the mean flow in this study. In
practice, this parameter is firstly tuned by conducting
fixed-bed computation. After that, the simulation of
movable bed is started.
It is well accepted that the direction and themagnitude of the near-bed velocities do not coincide withthose of the mean velocities in natural river conditions. In
order to estimate the near-bed velocities, they are
resolved into a primary velocity U and a secondary flow
velocity V. The former is in the same direction as that of
the mean flow, while the latter is normal to the primaryvelocity and has a magnitude related to the water depth,
the primary flow curvature and the primary velocity, i.e.

U aU u2 v2

1
r

(13)

1 O

wzb ª wqbx wqby º«»
wt ¬ wx
wy ¼

0

(18)

where zb= riverbed elevation; O = porosity of sedimenton the riverbed; qbx , qby = bedload transport rate in x andy direction, respectively.
The total bedload transport rate is evaluated with the
Ashida-Michiue formula as belowqbs  1 gd 3§
W ·§ W ·17W *3/e 2 ¨1  *c ¸ ¨1  *c ¸¨
W * ¸¹ © W * ¹©

(19)where qb= bedload discharge per unit width; s= specificgravity of sediment; g= gravitational acceleration; d=diameter of sediment; W * e , W * c , W * = dimensionlesseffective shear stress, critical shear stress and shear stress,
respectively. The dimensionless shear stress is connected2with the squared friction velocity u* and is defined asW*u*2( s  1) gd(20)

In the calculation of bedload transport in a channel
with a relatively gentle slope, the critical shear stress isobtained from the well-accepted Shield’s diagram orother popular criteria and the particle gravity does notaffect the sediment motion significantly. However, thelocal bed slope usually becomes much steeper in thesituation where local scour occurs. The local bed slope
may affect [1] the sediment threshold condition, i.e. thecondition just sufficient to initiate sediment motion, and[2] the sediment transport rate and direction once thesediment is in motion. In the conventional treatment, thelongitudinal velocity is usually assumed predominant
and the local bed slope is resolved into a longitudinal bed
and a transverse bed (e.g. van Rijn, 1993 and Duan et al.,2001). However, this kind of treatment is questionable－ 639 －

when applied to the sediment transport around local
scour holes where the magnitudes of the three velocity
components are almost comparable. In view of these, a
new method is proposed herein.

cosT

cosD cosD1  cos E cos E1  cos J cos J1(24)

Under the threshold condition, the resulted driving
force Fr is balanced by the stabilizing force Fs. The
former reads
Fr¬ªW cos J cos D  F cos D 1  cos T cos D ¼º i
 ª¬W cos J cos E  F cos E1  cos T cos E º¼ j (25)
 ª W sin 2 J  F cos J 1  cos T cos J º k¬¼

The latter acts in the direction opposite to the driving
force and has a magnitude of

Bed plane

Datum plane(26)

W cos J  F cos T tan M

Fswhere ĳ = angle of sediment repose. Using the equationof force balance, i.e. Fr = Fs, one obtains

Fig 5 Sediment particle under threshold condition
sin2 T  cos2 T tan2 M F 2  2WF l  cos J cosT tan2 M

Take a look at a sediment particle resting on a
sloping bed. Under the threshold condition, the forces
acting on the particle are shown in Fig.5. Assuming that
the Reynolds stress is large enough, all forces act throughthe particle center. The forces include the submerged
weight of the particle W and the hydrodynamic force F.
The latter one is usually resolved into a drag force FDtangential to the bed plane and a lift force FL normal tothe bed plane in the analyses of sediment transport. In the
figure, n stands for the normal direction of the local bed
plane and Į, ȕ & Ȗ are the corresponding direction angles(i.e. n = cosĮi + cosȕj + cosȖk). Write the hydrodynamic
force F as the multiplication of its magnitude F and its
direction, i.e.

F F cos D1i  cos E1j  cos J1k(21)

where Į1, ȕ1 & Ȗ1 are the direction angles of thehydrodynamic force F. Consequently, the lift force FL
and the drag force FD have the following forms.
F cos T cos D i  cos E j  cos J k(22)

FD F ª¬ cosD1 cosT cosD i  cos E1 cosT cos E j cosJ1 cosT cosJ kº¼

(23)

FL

 sin2 J  cos2 J tan2 M W 2

(27)

0

in which l is expressed by(28)

l cos J cosT  cos J1

The positive solution of Eq.27 isFWªm l  cosJ cosT tan2 M º¼
sin T  cos2 T tan2 M ¬2(29)

in whichml2  sin2T cos2 J sin2 J cos2T 2l cosJ cosT tan2 Msin2T sin2 J(30)
On the other hand, the submerged particle weight W
and the drag force FD are empirically expressed asfollows.

W

c1 s  1 U gd 3 FD

c2W b d 2(31)

where c1 = coefficient related to the shape of the particle;
c2 = coefficient depending on the local Reynolds number.
Normalizing Eq.29 with Eq.31, one has thedimensionless critical shear stress Ĳ*c for an arbitrary
sloping bed as follows

in which ș is the angle between the flow direction andthe normal direction of the bed plane and is expressed by－ 640 －

W*c

Wbc
c1
sinT
ªm l cosJ cosT tan2M º¼
s 1 Ugd c2 sin2T cos2T tan2M ¬(32)

where Ĳbc = critical near-bed shear stress.
In case of a horizontal bed (Į, ȕ, Ȗ) = (90°, 90°, 0)
and the flow parallel to the bed plane ș = Ȗ1 = 90°, Eq.32
becomes
c1
tan M
c2

W *c0

(33)

where Ĳ*c0 = dimensionless critical shear stress for ahorizontal bed. As has been mentioned before, this canbe easily obtained from some well-tested criteria such as
the Shields’ diagram. If Eq.33 is introduced into Eq.32, a
bed slope factor K is obtained below, expressing thedifference in the critical shear stress between on a
sloping bed and on a horizontal bed.
K

W*c
W*c0

sinT

ªm l cosJ cosT tan2 M º¼
sin T cos T tan M tanM ¬
2

2

2

(34)
Eq.34 indicates that the bed slope factor is not only afunction of the bed slope and the angle of sedimentrepose but also a function of the flow direction. If theflow is parallel to the bed plane, it has a simpler form asK

l 2  sin 2 J  cos 2 J tan 2 M  ltan Mtransverse bed slope although the treatment is physicallyquestionable (e.g. van Rijn, 1993 and Duan et al., 2001).
Comparing these two conventional bed slope factors
with Eq.35, it is found that the new bed slope factor ismore general while still maintains the simplicity for
programming.
According to the above analysis process, it is readily
to assume that the sediment movement follows the
direction of the resulted driving force. Hence the threecomponents of the bedload transport rate are easily
obtained from Eq.25.
The three components of the bedload transport rate
are mapped onto a 2D mesh covering the bed surface andEq. 18 is then solved. Consequently, the bed deformation
at each time step is obtained. Nevertheless, unrealistic
local bed slopes over the angle of sediment repose
probably occur in the simulation. This phenomenon
should be avoided for non-cohesive beds. A sedimentslide process is assumed to take place if the local bed
slope becomes too steep. After the slide process, the bed
level at each mesh is adjusted and all the local bed slopes
are no more than the angle of sediment repose.

(35)
Angle of repose

This bed slope factor is used in this paper for themodel applications. In case of a longitudinal bed slope(denoted as șl) and the longitudinal velocity is
predominant, i.e. (Į, ȕ, Ȗ) = (90°-șl, 90°, șl) and (Į1, ȕ1,Ȗ1) = (șl, 0, 90°-șl), one has the following bed slope factor
Horizontal planeK§ tan T l ·cos T l ¨ 1 ¸tan M ¹©

(36)
Fig 6 Sediment slide process

For a transverse bed slope (denoted as șt) and still thelongitudinal velocity is predominant, i.e. (Į, ȕ, Ȗ) = (90°,90°-șt, șt) and (Į1, ȕ1, Ȗ1) = (0°, 90°, 90°), one goes toanother bed slope factor

K

cos Tt 1 

tan 2 Tt
tan 2 M(37)

It is noted that Eq.36 and Eq.37 or their variants areconventional bed slope factors in the literatures. A simple
multiplication of the two factors has been suggested toaccount for the effect of a combined longitudinal andA plane-view of the CV (Control volume) surface on
the riverbed is shown in Fig. 6 (left). At each time step,
the program checks the nodes on the riverbed and thecorresponding neighbors. If a slope steeper than theangle of sediment repose has been detected connecting
node A (xA, yA, zA) and node B (xB, yB, zB), the elevationsof node A and node B should be adjusted to nodes A’ andB’ as shown in Fig. 6 (right). The change of the elevation
is denoted as įzA and įzB, respectively. Assuming that the
local bed slope connecting nodes A’ and B’ equals theangle of sediment repose, one has－ 641 －

zB G zB  zA G zA

2tanM xB  xA  yB  yA

2

(38)

On the other hand, sediment conservation should beensured during the adjustment, i.e.

¦ S A,neighborG z A¦ SB,neighborG zB= coefficient of the final discretized equations and b =source term excluding the pressure term. Defining
pseudo-variables corresponding to M and N as
¦ anb M nb  b f

Mˆ f¦ anb N nb  b f

Nˆ f

af

af

(42)

(39)
Eq. 41 can be written as

where SA, neighbor = projected area of the mesh surrounding
node A (i.e. mesh 1 to mesh 6 shown in Fig. 6) and SB,neighbor = projected area of the mesh surrounding node B(i.e. mesh 1, mesh 2 and mesh 7 to mesh10 in Fig. 6).
With Eq. 38 and Eq. 39, the new locations of node A
and node B can be obtained. This process is repeateduntil the bed slopes connecting all the CV nodes andtheir corresponding neighbors are not more than theangle of sediment repose.

MfM

af M faf N fwH
dl f b f
wx
wH
dl f b f¦ anb N nb  ³ gh
wy
¦ anb M nb  ³ gh

Nˆ f 

wH
dl f
wy
af

(43)

M*  M c

N*  N cNH

H *  H c (44)a f M *f

¦ anb M nb  ³ gh*

a f N *f*
¦ anb N nb  ³ gh

wH *dl f b f
wx

wH *dl f b fwy(45)

Subtracting Eq.45 from Eq. 41 and taking Eq.44 into
account, one obtains the velocity correction
M cf

N cf

³
c
¦ anb M nb

af

gh

wH c
dl f
wx
af

(46)

wH c
dl f
³ gh
c
wy
¦ anb N nbafaf

If the first term on the right hand side is omitted, thevelocity correction has a simple relationship with thepressure correction. Introducing this relationship into
Eq.44, one goes to Eq.47.

Mf

Mf * 

³ gh

wHc
dl f
wx
af

³ gh
Nf

Nf * 

wHc
dl f
wy
af

(47)

It is found that Eq. 47 has almost the same form as
Eq.43. In order to assure the mass conservation, thecontinuity equation for each CV should be satisfied, i.e.

(41)

where subscript f = the interface between two
neighboring CVs; subscript nb = the neighboring CV; a

Nf

³ gh

With these definitions and Eq. 41, it is readily to seew § wI · w § wI ·h*
 ¨ h* ¸  sI
wx ¨© wx ¸¹ wy © wy ¹

(40)
where hĭ = transport quantity ( for the momentum
equations, such definitions are used for clarity: M = hu
and N = hv), Ƚ = diffusive coefficient and sĭ = source
term.
Eq. 40 is integrated over a series of CVs (Control
volumes) covering the study domain using the FVM(Finite volume method) with a collocated variable
arrangement and is then submitted to discretization. The
discretization adopts a term-by-term process. For the
temporal term, the Euler scheme is employed. The powerlaw scheme is used for the convection and diffusion
terms. Surface fluxes are interploated following themethod suggested by Rhie and Chow (1983). For thevelocity-pressure coupling, the SIMPLER (Semi-implicit
method for pressure-linked equation, revised version)
procedure proposed by Patankar (1980) is adopted.
After the discretization, momentum equations for M
and N at a CV surface have the form as

wH
dl f
wx
af

In the derivation below, the guessed pressure and thevelocity component resulted from the momentum
equation not satisfying the continuity equation are
denoted with an superscript *, the correction of anvariable is denoted with an index ‘ and the variablewithout any superscript stands for the final value, i.e.

3.3 Numerical schemes and solution procedure
In a Cartesian coordinate system, the governing
equations for the flow have a common form as below.
w
w
w
hI  huI  hvI
wt
wx
wy

Mˆ f 

³ gh³

wH
P  ¦ M f lx  N f l y
wt

0

(48)

where subscript P = the present approximated CV;
Substituting Eq.43 and Eq.47 into Eq.48, one has the－ 642 －

so-called pressure equation and pressure-correction
equation for H and H’, respectively. The calculation
sequences are therefore outlined as follows.
[1] Calculate the coefficients in the momentum
equations. Omit the pressure term and solve for thepseudo-velocities.
[2] With the pseudo-velocities, the pressure equation
is solved. This leads to the water stage and water depth.
[3] Based on the resulted water stage, solve themomentum equations to get the approximated velocity
components.
[4] The resulted velocity components are used to
calculate the mass fluxes through CV faces with theRhie-Chow interpolation method. Solve thepressure-correction equation to obtain the pressure
correction.
[5] From the pressure correction, the velocity
correction is obtained. And correct the velocity field.
[6] Solve the transport equations for the turbulencekinetic energy and its dissipation rate, respectively. And
then update the eddy viscosity.
[7] Repeat the above procedures until the residuallevel becomes sufficiently small or the prescribedmaximum iteration step number is covered.
At the dry-wet boundary, the surface fluxes are notproperly estimated with any linear interpolation method
due to the steep gradient of the pressure. Furthermore,
the zero velocity and water depth at the center of the dry
CV bring numerical problems when calculating thevelocity correction by using the pressure correction.
Hence some additional treatments are needed between
step [4] and step [5] to account for the dry-wet boundary.

around dry CVs (Fig.7), i.e.
1
M*f ,new M*f  gh HAc HPc lx
af(49)
where subscript A = the adjacent CV (i.e. wet CV) and lx,
ly = projected length of the dry-wet boundary. The newfluxes are then introduced into Eq.47. Considering Eq.48
for mass conservation, a new pressure correction for the
dry CV is
ghH Ac ¦H Pcl2
 ¦ M *f , new l x  N *f , new l y
af

(50)

l2
ghH Ac ¦ AP / ' tafwhere Ap = area of present CV (i.e. dry CV) and ǻ t=time step.
Since the dry CV has a zero water depth in theprevious time step, the pressure correction gives the newwater stage directly. Other quantities such as ap at thecenter of the dry CV can be estimated by linear
interpolation from the surface values. Hence the dry CV
can be treated as a common CV. Then program goes tostep [5] with the updated pressure correction. This
treatment is found to provide an effective way to accountfor the water edges shifting without significant increasingof computational cost.
3.4 Boundary conditions
In actual river conditions, the boundaries generally
include the inlet, outlet and solid walls. The hydrographis usually known at the inlet boundary. Hence, thevelocity is prescribed there and is estimated according to
the local conveyance. If the discharge is normal to theinlet, then                  Dry CV

Wet CV1
N*f ,new N*f  gh HAc HPc ly
af

uin K

h2/3
n

where

K

Q

(51)

h5/3
l
¦
n

in which uin = velocity at each inlet mesh; Q = total
discharge at a time step; n =Manning’s roughness
coefficient and l = length of the inlet boundary. The
turbulence kinetic energy k is related to the velocity uin

P

f
A

2

Fig. 7 Treatment of the dry-wet boundary

The pressure correction in Step [4] is firstly used tocorrect surface fluxes (denoted with subscript new)and the turbulence intensity I via k 1.5uin I . The
dissipation rate İ is estimated from Eq.10 by specifying aviscosity ratio (Ȟt/Ȟ) around 10-20.
At the outlet, a zero gradient boundary is assumed
there. Moreover, an ad hoc treatment is applied to correct
the velocity components in order to assure the global
mass conservation, i.e.

－ 643 －

uout

Qin old
u
Qout

vout

Qin old
v
Qout

(52)

in which the discharge Q with subscripts in and out is thetotal discharge computed at the inlet and outlet
boundaries, respectively. The velocity component with asuperscript old stands for the computed value before
correction. It is noted that Eq.52 provides a direct way tocorrect the velocity at the outlet. Therefore, there is noneed to correct the outlet velocity again in the SIMPLER
procedure.
At solid walls such as channel banks and emerged
hydraulic structures, the wall function approach is used
for the computational convenience. However, theboundaries between wet and dry CVs are not treated assolid walls as has been mentioned in the previous
contexts.
4.

Results and discussions

sandbar has increased the sandbar height obviously.
Moreover, the stream way on the right side of thechannel is almost colonized by the head of the sandbar.
The local scour at the head of the impermeable groin isfilled with sediment transported from upstream. Behind
the impermeable groin, the bed is degraded significantly.
This is due to the flow separation as well as overtopping.
These observations suggest that flood conditions should
be prudently considered in groin designs.
Water discharge (m3/s)(m3/s)

Accumulated sediment amount (m3)

(m3)

Peak

1000

7200

700

Post-peakPre-peak

5400

500

3600

250

1800

0

0
0.0h

Selected results are presented in this section to further
the understanding of channel morphodynamics aroundgroins and to seek design methods for groins in riverrestoration. Firstly, the flow field and bed deformation
process around the two existing groins during a typicalflood season in 2005 is investigated. After that, the bed
deformation characteristics in two different cases (withgroins and without groins) are compared under the samewater discharges and sediment conditions. The
simulation is based on the experimental scale, but all the
results have been transformed to prototype scale forclarity.

12.2h

18.7h

24.5h

Fig 8 Modeled flood process

Impermeable groin
Permeable groin
Nagare Bridge
Fig 9 Hybrid mesh system

4.1 Bed deformation in a flood process
The bed evolution process during a flood in
September 2005 is investigated here. The flood is
simplified to a 3-stage process: pre-peak (0~12.2h), peak
(12.2~18.7h) and post-peak (18.7~24.5h). At each stage,
water discharge and sediment discharge are assumed
constant. According to the hydrograph and calibrated
sediment discharge by preliminary experiments, themodelled flood process is shown in Fig.8. Hybrid
unstructured mesh is used in the computation for a better
resolution. Part of the mesh system is shown in Fig.9.
From Fig. 10 to Fig. 13, the simulated bed evolution
process has been plotted. The flood discharges almost
trigger sediment movement in the whole channel.
Deposition of sediment in the upstream part of the big

Fig.13 and Fig.14 show the final bed contour after
the flood process according to computation andexperiment, respectively. Compared with the bed contour
before flood as given in Fig.10, the bed morphology
exhibits significant changes after the flood. Both thecomputation and experimental results indicate that thescour-deposition morphology around the impermeable
groin is destroyed. In some area, the scour and bed
degradation seem to be under-estimated. It is probably
caused by the under-estimation of the effect of thevertical velocity component, which contributes a lot tothe sediment transport inside of scour holes and aroundgroins.

－ 644 －

Bed level (Unit: m)
y (m)

y (m)

1250

1250

Sandbar

1150

1150

1050

1050

t=24.5ht=0.0h
950

950

850

850

1350

1450

1550

1650

1750

1850

1950

2050

1350

1450

1550

1650

1850

1950

2050

       x (m)x (m)Fig. 13 Final bed contour (Simulation)

Fig. 10 Initial bed contour
y (m)

y (m)

1250

1250

1150

1150

1050

1050

t=12.2ht=24.5h

950

850

1750

950

1350

1450

1550

1650

1750

1850

1950

2050

850

1350

1450

1550

1650

1750

1850

1950

2050

       x (m)x (m)Fig.11 Bed contour before peak (Simulation)

Fig.14 Final bed contour (Experiment)y (m)Water level (m)Experiment (18.7h)

1250

Simulation (18.7h)
1150

1050

t=18.7h

Experiment (12.2h)   Experiment (24.5h)

950

Simulation (12.2h)   Simulation (24.5h)
850

1350

1450

1550

1650

1750

1850

1950

2050

x (m)



Fig.12 Bed contour after peak (Simulation)

Distance from the confluence (km)Fig.15 Water level in the flood process－ 645 －

Bed level (m)Before flood
After flood (Experiment)
After flood (Simulation)

Distance from the left bank (m)Fig. 16 Bed level before permeable groin (6.2km)
Bed level (m)Before flood
After flood (Experiment)
After flood (Simulation)

The computed water levels along the centerline of the
low flow channel are compared with those of theexperiments as shown in Fig. 15. In general, water levels
exhibit no significant differences at the end of the
pre-peak stage (12.2h) and at the end of the post-peakstage (24.5h). Higher water level is observed at the endof the peak (18.7h) due to the higher flood discharge.
Simulation results are in reasonable agreement withthose of the experiments. Nevertheless, slightover-estimation is found upstream, in particular under the
peak discharge. It may be attributed to the use of thesame set of Manning’s roughness coefficients. The
development of ripples and its effect on the bed
resistance are not accurately accounted for in the model.
The roughness, in the upstream area and in case of higher
discharges, has probably been over-estimated.
y (m)1.0m/s1300
1250

Distance from the left bank (m)1200

Fig. 17 Bed change between the two groins (S1)1150

Bed level (m)

1100
1050

Before flood
After flood (Experiment)
After flood (Simulation)

1000
950
1600 1650 1700 1750

1800 1850 1900 1950 2000 2050

x (m)Fig. 19 Surface velocity (Experiment, t=18.7h)

Distance from the left bank (m)Fig. 18 Bed change after impermeable groin (S2)y (m)1.0m/s1300

The bed deformation at typical transverse sections isshown from Fig.16 to Fig.18. The locations of thesections are depicted in Fig.3. Due to the wavy bed ofripples, the measured bed profiles exhibit great
irregularity. Irrespective of this, the computational result
is quite encouraging. Fig.18 deserves special attention.
Again, the scour-deposition pattern (i.e. pool-riffle
morphology) caused by the impermeable groin before
the flood is found to be destroyed during the flood. It
gives quantitative evidence to previous observations in
the change of the bed contour.

1250
1200
1150
1100
1050
1000
950
1600 1650 1700 1750 1800 1850 1900 1950 2000 2050

－ 646 －x (m)Fig.20 Mean velocity (Simulation, t=18.7h)

particular the impermeable one. At the end of the flood,
one may observe some area where the flow is almost
stagnant in the experimental result (Fig.21). However,this phenomenon is not well predicted with the model(Fig.22). The over-estimation of the water level as
depicted in Fig.15 is probably responsible for thisdiscrepancy. The bed forms change under different flow,
sediment and bed conditions, whose effects on the bed
roughness necessitate quantitative estimation. Thisshould be included in the future research.

1.0m/s

y (m)

1250

1200

1150

1100

1050

1000

950
1600

1650

1700

1750

1800

1850

4.2 Bed deformation in cases with/without groins
In order to investigate the impacts of groins on thechannel morphodynamics, experiments with and without
groins have been carried out. Considering the importanceof flood discharges in the design of groins, cases underflood discharges are selected and presented in thissection, which are detailed in Table 1.

1900

x (m)
Fig.21 Surface velocity (Experiment, t=24.5h)

1.0m/s

y (m)

1250

Table 1 Selected experimental conditions
1200

Cases

Discharge (m3/s)

GroinsTime (hour)

1
2
3
4

2000
2000
500
500

Yes
No
Yes
No

16
16
80
80

1150

1100

1050

1000

950
1600

1650

1700

1750

1800

1850

1900

x (m)
Fig.22 Mean velocity (Simulation, t=24.5h)

The surface velocities were measured in theexperiment considering the shallowness of the river flow.
The surface velocities around groins at the end of theflood peak (18.7h) and at the end of the flood (24.5h) are
shown in Fig.19 and Fig.21, respectively. The computed
mean velocities at corresponding time are also plotted in
Fig.20 and Fig.22. It was found that the measured
surface flow patterns are quite similar to the computed
mean flow patterns. The flow fields are closely related tothe bed morphologies and water depths. The highvelocity zone basically follows the river thalweg. Since
the groins are submerged during the flood, they do nothave significant influence on the flow directions.
However, the flow was accelerated around the groins, in

The flood discharges 2000m3/s and 500m3/s
correspond to the maximum annual discharge and thedominant discharge of the river, respectively. The bed
contours before the flood and after the flood are plottedfrom Fig. 23 to Fig. 34 according to the experimental andsimulation results.
In case of flood discharge of 2000m3/s, the numericalresults give very obvious evidence that the groin
structures are effective for the bank protection.
Comparing Fig. 24 with Fig. 27, one may observe
several scour in the area where the groins are removed.
Unfortunately, the groins do not show advantages increation of deposition area in desired locations. This isprobably caused by the over-topping flow during theflood.
Similar conclusions may be drawn from the results incase of discharge of 500m3/s (Fig. 29 to Fig. 34).
However, the movement of sandbars and development of
bed forms under bigger discharge (Case1 and Case2) aremuch more significant compared with those undersmaller discharge (Case3 and Case4).

－ 647 －

8

8.5

9

9.5 10 10.5 11 11.5 12 12.5 13 13.5 14 14.5 15 15.5 16

Bed level (Unit: m)

Fig.23 Bed contour (Case1, t=0h)

Fig.26 Bed contour (Case2, t=0h)

Fig.24 Bed contour (Case1, t=16h, Simulation)

Fig.27 Bed contour (Case2, t=16h, Simulation)

Fig.25 Bed contour (Case1, t=16h, Experiment)

Fig.28 Bed contour (Case2, t=16h, Experiment)－ 648 －

8

8.5

9

9.5 10 10.5 11 11.5 12 12.5 13 13.5 14 14.5 15 15.5 16

Bed level (Unit: m)

Fig.29 Bed contour (Case3, t=0h)

Fig.32 Bed contour (Case4, t=0h)

Fig.30 Bed contour (Case3, t=80h, Simulation)

Fig.33 Bed contour (Case4, t=80h, Simulation)

Fig.31 Bed contour (Case3, t=80h, Experiment)

Fig.34 Bed contour (Case4, t=80h, Experiment)－ 649 －

5.

Conclusions

This paper presents some results of a research workon the channel morphodynamics due to the constructionof groins with both experimental and numerical methods.
The experiment is based on a large-scale physical modeland is able to resolve the local flow and bed deformation
process with a high accuracy. The numerical model takes
into account the effects of secondary flows and local bed
slopes on the flow and sediment transport. Unstructured
meshes are employed in the discretization process of thegoverning equations, making possible the accurate
resolution of commonly encountered complex
geometries in actual river conditions. The numericalresults have been compared with those of theexperimental measurements and reasonable agreements
have been obtained.
Compared with most of the existing researches for
simplified channels and groin structures under idealized
flow and sediment conditions, the results of this research
are of practical value. These results may be expected toprovide important information for the sustainable designof groin structures. However, comparisons between theexperimental and numerical results are not enough yet.
Numerical simulations based on the prototype scale aswell as the long-term monitoring of the project are a
crucial consideration in the next step. Development of a
3D model is also of great need to resolve the detailed
flow structure and sediment movement around groins.
Acknowledgements
The authors are grateful to Dr. Yasuyuki Baba and
Mr. Akira Nakanishi, who provide a lot of advices andconveniences for the numerical simulations. Mr.
Michinari Tani and Mr. Hirokazu Ikeda are also
sincerely acknowledged. Without their expertise andkindness, the experimental data and analyses would not
be possible.
References
Ashida, K. and Michiue, M. (1972): Studies on bed loadtransportation for nonuniform sediment and river bed
variation. Disaster Prevention Research Institute
Annuals, Kyoto Univ., No.14B, pp. 259-273. (inJapanese)
Duan, J.G., Wang, S.S.Y. & Jia, Y. (2001): The

application of the enhanced CCHE2D model to study
the alluvial channel migration processes, Journal of
Hydraulic Research, IAHR, Vol. 39, No. 5, pp. 1-12.
Engelund, F. (1974): Flow and bed topography inchannel bends, Journal of the Hydraulics Division,
ASCE, Vol. 100, No. HY11, pp. 1631-1648.
Kuhnle, R.A., Alonso, C.V. and Shields, F.D.Jr. (2002):
Local scour associated with angled spur dikes. Journalof Hydraulic Engineering, ASCE, Vol.128, No. 12,
pp.1087-1093.
Minh Duc, B., Wenka, T. and Rodi, W. (2004):
Numerical modeling of bed deformation in laboratory
channels, Journal of Hydraulic Engineering, ASCE,
Vol.130, No.9, pp.894-904.
Muto, Y., Kitamura, K., Khaleduzzaman, A.T.M. andNakagawa, H. (2003): Flow and bed topographyaround impermeable spur dykes, Advances in River
Engineering, JSCE, Vol.9. (in Japanese)
Nagata, N., Hosoda, T., Nakato, T. and Muramoto, Y.
(2005): Three-dimensional numerical model for flow
and bed deformation around river hydraulic structures,
Journal of Hydraulic Engineering, ASCE, Vol. 131, No.
12, pp. 1074-1087.
Nanbu, Y. (2007): Experimental Study on the
Restoration of Channel Morphology with Groinsaround the Kouzuya Bridge on the Kizu River,Bachelor’s Thesis, Kyoto University, Kyoto, Japan. (inJapanese)
Patankar, S.V. (1980): Numerical Heat Transfer andFluid Flow, McGraw-Hill, New York.
Rhie, C.M. and Chow, W.L. (1983): Numerical study ofthe turbulent flow past an airfoil with trailing edge
separation, Journal of AIAA, Vol. 21, pp. 1525-1532.
Rodi, W. (1980): Turbulence Models and Their
Application in Hydraulics - a State of the Art Review,
University of Karlsruhe, Karlsruhe, Germany.
Shields, F. D. Jr., Cooper, C.M. and Knight, S.S. (1995):
Experiment in stream restoration, Journal of HydraulicEng., ASCE, Vol.121 (6), pp.494-502.
Shimizu, Y. and Itakura, T. (1991): Calculation of flow
and bed deformation with a general non-orthogonal
coordinate system, Proc. 24 IAHR Congress, Madrid,
C-2, pp. 241-248.
Sukhodolov, A., Engelhardt, C., Kruger, A. andBungartz, H. (2004): Turbulent flow and sediment
distributions in a groyne field. Journal of HydraulicEngineering, ASCE, Vol.130, No. 1, pp.1-9.
Touchi, D. (2006): Investigation of Bed Evolution－ 650 －

around Spur Dykes with Large Scale Physical Model
Experiments, Bachelor’s Thesis, Kyoto University,
Kyoto, Japan. (in Japanese)van Rijn, L.C. (1993): Principles of Sediment Transportin Rivers, Estuaries and Coastal Seas, AQUA
Publications, Amsterdam, The Netherlands.

Zhang, H., Nakagawa, H., Muto, Y., Touchi, D. and
Muramoto, Y. (2006): 2D numerical model for river
flow and bed deformation based on unstructured mesh,
Journal of Applied Mechanics, JSCE, Vol. 9, pp.
783-794.

᳓೙ࠍ᦭ߔࠆᴡ㆏ߩᴡᐥᄌേߦ㑐ߔࠆ⎇ⓥߣߘߩᴡᎹౣ↢߳ߩᔕ↪ᒛ ᶈ࡮ਛᎹ ৻࡮ᱞ⮮ ⵨ೣ࡮᧛ᧄ ཅ㓶*࡮⇧ౝ ᄢਃ㇢**࡮ධㇱ ᵏ▸***੩ㇺᄢቇฬ⹷ᢎ᝼

**੩ㇺᄢቇᄢቇ㒮Ꮏቇ⎇ⓥ⑼

ⷐ ᣦᧄ⎇ⓥߢߪ㧘ᴡᎹౣ↢੐ᬺߢ⸳⟎ߐࠇߚ᳓೙Ꮏ๟ㄝߩᵹࠇߣᴡᐥᄌേߦ㑐ߒߡ㧘⒖േᐥታ㛎ߣᢙ୯ࠪࡒࡘ࡟࡯࡚ࠪ
ࡦࠍⴕߞߚ‫ߩ࡞࡯ࠤࠬ ޕ‬ᄢⷙᮨᮨဳ᳓〝ࠍ↪޿ߚ᳓ℂታ㛎ߢߪ㧘ዪᚲ⊛ߥᵹࠇߣ࿯⍾⒖േ⃻⽎߇㧘චಽߥ♖ᐲߢౣ⃻ߢ߈ߚ‫ޕ‬ᢙ୯ࡕ࠺࡞ߢߪ㧘㕖᭴ㅧᩰሶࠍ↪޿ࠆߎߣߦࠃߞߡ㧘ታᴡᎹߩࠃ߁ߥⶄ㔀ߥ࿾ᒻ߿᭴ㅧ‛ߩᒻ⁁ࠍㆡಾߦ⴫⃻ߔࠆߎߣ߇น⢻ߦߥߞߚ‫⎇ᧄޕ‬ⓥߢ㐿⊒ߒߚᢙ୯ࡕ࠺࡞ࠍ᳓ℂታ㛎ߩᵹࠇߦኻߒߡㆡ↪ߒߚߣߎࠈ㧘߶߷ḩ⿷ߥ♖ᐲߢᵹᴫߣᴡᐥᄌേࠍౣ⃻ߔࠆߎߣ߇␜ߐࠇߚ‫ޕ‬ࠠ࡯ࡢ࡯࠼: ࡕ࡞ࡈࠜ࠳ࠗ࠽ࡒࠢࠬ㧘᳓೙㧘ᴡᎹౣ↢㧘ᄢⷙᮨ᳓ℂታ㛎㧘ⶄ㔀Ⴚ⇇㧘㕖᭴ㅧᩰሶ

－ 651 －