River Flow 2014 – Schleiss et al. (Eds)© 2014 Taylor & Francis Group, London, ISBN 978-1-138-02674-2

A numerical calculation method for flow in the presence of isolatedboulders atop a rough bed by using an enhanced depth integrated
model with a non-equilibrium resistance law
T. Uchida & S. Fukuoka

Research and Development Initiative, Chuo University, Japan

A.N. Papanicolaou & A.G. Tsakiris

IIHR-Hydroscience and Engineering, The University of Iowa, USA

ABSTRACT: This paper presents a new calculation method for flow over rough bed with a superimposed
array of wall-mounted boulders. Since the resistance of roughness depends on the relative water depth toroughness height for the case of Large Scale Roughness (LSR), it is necessary for calculating the flow overgravel bed rivers to separate the effects of LSR from that of the Small Scale Roughness (SSR). A quasi-3D
calculation method (BVC method) is developed to evaluate velocity vertical distributions around boulders
without making the shallow water assumption. To evaluate additional resistance on rough bed resulted
from the increment in momentum exchange due to flow around boulders, a non-equilibrium wall law is
proposed and the bottom vortex layer and roughness layer equations are derived. The advantages of BVC
method with the non-equilibrium wall law over the 2D and BVC method with the conventional equilibrium
wall law are demonstrated through comparisons with benchmark experimental results.
1

Introduction

The presence of isolated large roughness such as
boulders or clusters in gravel bed rivers, which insome instances may even protrude through the freesurface, produces challenging conditions for determining near-bed flow conditions and flow resistance. Many researchers have indicated that theconventional wall law is inadequate when the large
roughness element diameter is comparable to thewater depth (Bathurst, 1985; Aguirre-Pe & Fuentes,
1990; Clifford et al., 1992, Katul et al., 2002).
Recently, Computational Fluid Dynamics (CFD)
has been widely used to evaluate flows in gravel
bed rivers (Olsen & Stokseth 1995; ­Nicholas, 2001;
Carney et al., 2006; Nikora et al., 2007a, 2007b;
Rameshwaran et al., 2011). Since the resistance ofroughness depends on the relative water depth toroughness height for the case of Large Scale Roughness (LSR), the concept of the roughness scaledistinction, or multi-scale roughness, in which theeffects of LSR is separated from that of the Small
Scale Roughness (SSR), has been proposed by many
researchers (Clifford et al., 1992; Nicholas, 2001;
Carney et al., 2006; Rameshwaran et al., 2011). Yet
some important issues still remain unresolved.
The first issue is an evaluation method of LSR.
While we recognize the need of a non-hydrostatic
three dimensional method for this purpose(Olsen & Stokseth 1995; Nicholas, 2001;­Carney et al., 2006; Nikora et al., 2007a, 2007b;
Rameshwaran et al., 2011), a reliable depth­integrated model is also required for practicalapplications to rivers. However, previous depth
integrated models including quasi-three dimensional models cannot be used in this case, becauseof their hydrostatic pressure distribution assumption (Ishikawa et al., 1986; Fukuoka et al., 1992;
Jin & Steffler 1993; Yeh and Kennedy, 1993).
The second issue is the effects of LSR on SSR.
The applicability of the conventional wall law forrough bed (log-law with equivalent roughness) is
uncertain because the near bed flow is affected by
the LSR presence. For flows over a gravel bed, theexistence of boulders and the irregularity of the
bed topography and sediment diameter exacerbate
the deviations from an equilibrium bottom boundary condition. The essential solution is to calculate
the non-equilibrium motion around bed surface
including the roughness layer (Olsen and Stokseth,
1995; Lane et al., 2004; Carney et al., 2006; Nicora
et al., 2007a, 2007b; Rameshwaran et al., 2011).
Recently, we have developed a new quasi-three
dimensional calculation method that avoids making shallow water assumptions, such as hydrostaticpressure distribution, the general Bottom Velocity Computation (BVC) method, in which the­equations for velocity and pressure ­distributions335

in the vertical direction are solved with depth
averaged velocity equations (Uchida & Fukuoka,2012). The BVC method has been validated forseveral flows including rapidly varied flows overa structure and three dimensional flows arounda non-submerged structure where the horseshoe
vortex is prominent (Uchida & Fukuoka, 2012,2013; Fukuoka & ­Uchida, 2013). However, theapplicability of the BVC method for evaluating theeffects of submerged boulders as LSR has not been
­discussed. In addition, the previous BVC method
has employed the conventional equilibrium wall
law without considering the effect of LSR on SSR.
This paper proposes the improved BVC method to
employ a non-equilibrium wall law by introducing
momentum equations for the bottom vortex layer
and roughness layer. The validity of the present
method is discussed through comparisons withexperimentally measured velocity distributions inwithin an array of isolated boulders atop a roughbed and with calculated results from the 2D and
BVC with the conventional equilibrium wall law
methods.
2 calculation method
2.1

Representation of multi-scale roughness
and the framework of the BVC method

In the present calculation method, the roughness
is divided into Large Scale Roughness (LSR),
which includes boulders with comparable heightto the water depth that may even protrude throughthe free surface, and Small Size Roughness (SSR),
the contributions of which are evaluated using thewall law. The LSR is explicitly taken into accountin the flow calculation as variation in bed topography as indicated in Figure 1. The BVC method(Uchida and Fukuoka, 2011) is refined to ­calculate

Figure 1. Representation of multi-scale roughness.

fluid flow and evaluate the effects of LSR directly,
newly introducing momentum and vorticity
exchange owing to those fluxes on bed. The BVC
method is a quasi-three dimensional calculation
method without the need to make key shallow
water assumptions, such as the hydrostatic pressure distribution assumption, in which equations
for velocity and pressure distributions in verticaldirection are solved with depth averaged velocity
equations as described below.
The BVC method considers the thin vortex layerδzb on bed surface and calculates the bottom velocity acting on the layer as shown in Figure 2. The
bottom equations are derived by depth-integrating
the horizontal vorticity:
 ∂Wh
∂z
∂z 
ubi = usi − ε ij 3Ω j h − − ws s + wb b 
∂
x
∂
x
∂xi 

i
i

(1)

where, i, j = 1(x), 2(y) (x,y: horizontal direction,z: vertical direction), ubi: bottom velocity, usi:
water surface velocity, εijk: Levi-Civita symbol, Ωj:
Depth Averaged (DA) horizontal vorticity, h: water
depth, W: DA vertical velocity, zs: water level, zb:
bed level, ws, wb: vertical velocity on water surface
and ­bottom. Equation (1) implies that the bottom velocity on the thin vortex layer on the bed
can be evaluated once the water surface velocity,­horizontal vorticity integrated over water depth,and spatial variation in vertical velocity are known.
The equations for these unknown variables arederived in the BVC method considering the cubic
vertical velocity equation (2).
u ′ i = ∆ui (12η 3 − 12η 2 + 1) + δ ui ( −4η 3 + 3η 2 )(2)where, u′i = ui − Ui, Ui: depth averaged velocity, δui:usi − ubi, ∆ui: usi − Ui, η = (zs − z)/h. The quadric curveof velocity distributions are used for the ­uniform

Figure 2.

336

BVC method.

flow, substituting ∆ui = δui/3 (usei = Ui + δui /3) in theequation (2) (Uchida & Fukuoka, 2011). On theother hand, the bottom pressure deviation fromthe hydrostatic pressure distribution is given bydepth-integrating the vertical momentum equation
after neglecting the temporal variation term forcomputational stability and simplicity.
∂hτ zj∂z
dpb ∂hWU j+ τ bj b −
=
∂x j
∂x j
∂xj
ρ

(3)

where, dp: pressure deviation from hydrostaticpressure distribution (p = ρg(zs − z) + dp), dpb: dp on
the bottom, τij: shear stress tensor due to molecularand turbulence motions.
To calculate Eqs. (1) and (3), the governing
equations for the following unknown quantities
are solved: water depth h (Eq. (4)), DA horizontal velocity Ui (Eq. (5)), DA turbulence energyk (Eq. (6)), DA horizontal vorticity Ωi, (Eq. (7))­horizontal velocity on water surface usi (Eq. (8)),
DA vertical velocity W (Eq. (9)).
∂h ∂Uj h− wσ b = 0
+
∂xj
∂t

(4)

∂z τ
∂U i h ∂U iU j h− wσ bubi = − gh s − bi
+
∂xi
∂x j
∂t
ρ
∂hTij  ∂hdp0 dpb ∂zb −++ρ∂x  ρ∂x
ρ ∂x 
j

i

i

(5)

∂k
∂k 1 ∂  vh ∂k 
+ Uj=
+ Pk − ε∂t∂xj h ∂xi  σ k ∂xi (6)∂hDω ij∂Ω i h+ wσ bω bi= Rσ i + Pω i +
∂x j
∂t

(7)


 ∂z
∂u
∂usi
 ∂dp s
+ Psi
+ usj si = −  g + 




∂
z
∂
x
∂x j
∂t
z = zs 
 i


(8)


∂ 
2 ∂φ
P
 (C ∆t ) +φ −φ = 0∂x j 
∂x j 

(9)where, wσb: velocity flux across the bottom layer as
described in the following section (wσb = 0: equilibrium wall law), Tij = τ ij + ui ′ uj ′ (superscript bar
indicates depth averaged value and evaluated byvertical velocity distribution), τij = vSij, v = vm + vt,
vm: kinematic viscosity coefficient, vt: kinematic
eddy viscosity coefficient, Sij: DA strain velocity,
dp0: DA pressure deviation evaluated as dp0 = dpb/2
with the assumption of the linear pressure distribution for simplicity, Pk: the production termof turbulence energy, ε: the dissipation term oft­ urbulence energy, Rσi: the rotation term of verticalvorticity (Rσi = usiωσs − ubiωσb), ωσs, ωσb: rotation ofusi, ubi, Pωi: the production term of horizontal vorticity, Dωij: horizontal vorticity flux due to convection, rotation, dispersion and turbulence diffusion,
Psi: the production term of surface velocity (shear
stress under the water surface layer), C = k1h/∆t,
k1 = 1/20, φ = (Wh)n+1 − (Wh)n, φP = (Wh)P − (Wh)n,(Wh)P: the predicted Wh calculated by the continuity equation using (δui)P, (δui)P: predicted velocity
difference calculated by Eq. (1) with (Wh)n. Theseequations can be solved once a vertical velocity
distribution is assumed. More details on the equations and the numerical computation methods
for the BVC method may be found in Uchida &­Fukuoka (2011, 2012 and 2013).
2.2 Equilibrium bottom boundary conditions
Boundary conditions on the bottom for the bed
shear stress terms in momentum equations (5) andthe production terms in vorticity equations (7) canbe derived from the wall law for the rough bed.
According to the logarithmic law for flat roughbed, the relationship between friction velocity andvelocity near the bed may be expressed as (Schlichting and Gersten, 2000):z ub= cb = 1 ln  b  + Arκ  ks u*(10)

where, ub2 = ubiubi, zb = δzb + δz0, δzb: very thin vortex layer assumed in the BVC method, δz0: depthof the log law origin from the top of the roughness height. The vortex layer thickness is defined
as h/δzb = e3 − 1 (e: base of natural logarithm) tosatisfy the vertical velocity distribution simplifiedof equation (2) for uniform flow conditions withlogarithmic velocity distribution.
The vorticity distribution is given by the differential of velocity equation (2) with respect to z.
The production terms Pωi in equation (7) are given
by equation (11).
Pω i = C pωvtb (ω bei − ω bi ) h + Psω i(11)where, Cpω = κ/α, vtb: vt on the bottom converted into
depth averaged scale, ωbi: horizontal vorticity on bottom, ωbei: equilibrium ωbi, Psωi: vorticity production
due to flow separation (Uchida & Fukuoka, 2011,2013). The equivalent vorticity on the bottom is twiceof that of the depth averaged value ωbe = 2Ωe = 2δue/h.
Equation (12) gives the equilibrium vorticity on thebottom (Uchida & Fukuoka, 2011):

ω be =

337

2u*  zs ln
κ h  zb (12)

2.3 Non-equilibrium bottom boundary conditions
An issue with the equilibrium bottom boundary condition is the assumption of uniform flow
condition for a bottom velocity regardless ofthe velocity and vorticity distributions above the­bottom. For example, the same velocity distribution for vortex and roughness layer is assumed for
different velocity distributions above the bottom(solid and dotted lines of velocity vertical distributions in BVC calculation layer in Fig. 3). This
study defines the unknown variables in the vortex
and roughness layers as shown in Figure 3 to calculate non-­equilibrium flow under the bottom of
the BVC calculation layer. The assumption of verythin layer gives the following momentum equations
for the vortex and roughness layers:∂( pb + ρ gzb ) τ bi
∂u
∂uvi
τ
+
− ti
+ uvk vi = −
∂xk
∂t
ρ∂xi
ρδ zb ρδ zb(13)

D∂( pb + ρ gzb ) τ ti
∂u
∂uri
+
− i
+ urk ri = −
∂xk
∂t
ρ∂xi
ρδ zr ρδ zr(14)

where, k = 1,2,3, uvi, uri: i direction velocity in the­vortex and roughness layers, pb: bottom pressure pb =
ρgh + dpb, τbi, τti: shear stress acting on the bottom
and roughness layer surface, δzb, δzr: the thicknessof vortex and roughness layers, Di: resistance term
due to drag force acting on ­roughness. The shear
stress acting on the lower surface of the roughness
layer is neglected in equation (14). The roughness
layer thickness is assumed as δzr = ks in this paper.
The fluxes wσb, wσt crossing the bottom and roughness surface are calculated with the following continuity equations for vortex and roughness layers:∂δ zbuvi∂λδ zt uri
wσ b = wσ t −, wσ t = −
∂xi
∂xi

(15)

where, λ: porosity in the roughness layer (λ = 0.4).
The shear stress and resistance terms in the
momentum equations for vortex and roughness

layers are derived by using the equilibrium wall
law as follows. The shear stress acting on thoselayers following the Boussinesq approximation­(Shlichting & Gersten, 2000) can be assumed as:(ubi − uvi )τ bi  ∂ui =  vt = vtb ⋅ Abhρ  ∂z  b(uvi − uri )τ ti  ∂ui =  vt = vtt ⋅ At

∂z t
h
ρ

(16)
(17)

The coefficients Ab and At are decided to satisfy
τb = τth/(h + δzb) = ρu*2 (τb2 = τbiτbi, τt2 = τtiτti) and
vtb = vtt = αu*h with the conventional logarithmic
velocity vertical distribution for the equilibrium
condition:Ab =1
1  δ z0 + δ zb / 2 , c = ln + Arα (cb − cv ) v κ 
ks

(18)

At =

h + δ zb
1
1 δz , ct = ln  0  + Ar
α (ct − cr ) hκ  ks (19)

where, cr: coefficient for the ratio between thevelocity in the roughness layer and shear velocity
for uniform flow condition.
The kinematic eddy viscosities on the bottom andthe top of the roughness layer are expressed by using
the distance of the wall and velocity difference:
vtb = α bδ ub h, δ ub 2 = (ubi − uvi ) (ubi − uvi )(20)

vtt = αtδ ut h, δ ut = (uvi − uri ) (uvi − uri )

(21)

2

where, the coefficient ab and at are given by ab = α/(cb - cv) and at = α/(cv - cr) to reduce vtb = vtt = αu*h forthe equilibrium condition. However, using bottom
vorticity ωb gives another expression of the bottom
kinematic eddy viscosity vtb as presented in previous
research (Uchida & Fukuoka, 2011, 2013):
vtb =ακω b h22 ln( zs / zb )(22)This study uses the average value of equation (20)and (22) for vtb.
The resistance term Di is expressed with theassumption of the roughness layer as the one layerof spheres array and drag coefficient CD:
Di =Figure 3. Calculation domains for non-equilibrium
wall law.

πCDρur uri
8(23)

where, ur2 = uriuri, CD = 0.4. By considering the equilibrium condition in equation (14), the coefficient
cr is given by equation (24).

338

cr =

8  δ zr + δ zb 
1 +

πCD 
h

(24)

The equilibrium vorticity on the production
term in the vorticity equation (11) is assumed byusing velocity gradient in vorticity in the layer as:ω bej = 2Ωej = 2ε ij 3Aω

ubi − uvih(25)

where, the coefficient Aω is obtained to reduce
equation (25) to equation (12) for the equilibrium
condition.
Aω =

z 
1
ln  s κ (cb − cr )  zb (26)

3 calculation condition
The conditions for the benchmark experimental
results are described in the following. The experiments were conducted in a 21.0 m long, 0.91 m
wide and 0.53 m deep channel with slope S = 0.003
by Papanicolaou et al. (2012). Three layers of glass
beads with diameter d0 = 0.0191 m were placed ina fixed, well-packed arrangement on the bottom.
48 immobile, spherical boulders with diameterdc = 0.055 m were fixed in a staggered arrangement atop the flat porous bed of the glass beads.
The diagonal spacing between boulders was l = 6dc
to establish isolated roughness conditions often
found in natural mountainous streams (Fig. 4).
The experiment was conducted for fully submerged
condition h0/dc = 3.5 (h0: water depth). Velocity vertical distributions were measured by ADV along
longitudinal sections across the center of boulders
and between the boulders.
The calculation domain includes the completeboulder arrangement (Fig. 4). At the upstream
boundary, discharge Q = 0.127 m3/s is applied tomatch the calculated average velocity with averaged velocity measured by ADV (bulk velocity,

Figure 4.

U0 = 0.72 m/s). The water depth at the downstream
end is given to reproduce the experimental water
depth (h0 = 0.193 m) at the center of the boulder
arrangement. The bed height is set to the top of thebeads. The origin of the axis x = y = z = 0 is definedon the bed surface at the center of the boulder
arrangement (x: stream-wise direction, y: leftward
direction toward downstream, z: upward direction). While boulders are taken into the calculation
directly as bed height undulations, the equivalentroughness height ks and the log law origin heightδz0 are used to evaluate the resistance of the beads.
The ks and δz0 are determined to be ks = 0.9 d andδz0 = 0.3 d, using measured velocity vertical distributions on the flat rough bed without the boulders
under identical flow conditions. The results of
the BVC method with Non-equilibrium Wall Law(BVC-NWL) are compared with the results fromthe 2D calculation method and the BVC methodwith Equilibrium Wall Law (BVC-EWL) for flow
over rough bed with boulders. In the 2D calculation
method, the quadric velocity vertical distribution
and hydrostatic pressure distribution are assumed
with calculating equations (4)–(6). The BVC-EWL
employs equations (10) and (12) without calculating momentum equations (13) and (14), assuming
vertical flux across the bottom wσb = 0. The parameters, ks = 0.9d and δz0 = 0.3d, are the same for all
three calculation methods. The drag coefficient
CD = 0.4 and the roughness layer thickness δzr = ks
are additionally given for BVC-NWL.
4 calculation results
and discussion
Figure 5 shows comparisons of the velocity within
the water surface and bottom horizontal planes
between the ADV measurements (Papanicolaou
et al., 2012) and the calculations using the 2D,
BVC-EWL and BVC-NWL methods. The water
surface velocity and bottom velocity were measured at η = 0.1 (highest measurement points)
and 0.95, respectively. The experimental results

Calculation domains for non-equilibrium wall law.

339

Figure 5. Water surface velocity and bottom velocity distributions (SI unit).

in F­ igure 5 illustrate the complex bottom velocity field with flow separations behind boulders,
while the water surface flow shows a little spatial­variation. All calculated water surface velocity distributions are in good agreement with the experimental results, showing almost uniform flow at
the water surface, even though the 2D calculates
a little lower water surface velocity than others.
However, there is large difference in the calculation results for the bottom velocity distribution.
Since the 2D method cannot take into account­ eformation in the ­velocity ­vertical distributionsddue to the ­submerged boulders, the calculated bottom ­velocity is almost uniform with considerable
high velocity magnitude compared to the measured results. The bottom velocity by BVC-EWL
captures the velocity reduction upstream of theboulders and the reverse flows in the separation
zone behind them. ­However, the bottom velocity
in-between boulders is still overestimated compared to the measured results. The bottom velocity distribution by BVC-NWL method can capture

340the velocity reduction upstream and downstreamof the boulders, reproducing the measured averagebottom velocity. A detailed comparison of bottom
velocity between the measurement and BVC-NWL
shows that the calculation overestimates effects ofboulders such as flow convergences beside boulders and flow separations behind them. This is
attributed to the gaps between the lower surface ofthe boulders and the beads surface, which were nottaken into account in the calculation; as the calculation bed surface is defined at the boulder top­surface. Notwithstanding these differences between
the experiment and calculation conditions, the calculation results provide good explanations of theexperimental results.
Figure 6 shows stream-wise velocity vertical distributions along the longitudinal sections
across boulders (y/dc = 0) by measurement andcalculations. The differences in calculated velocity vertical distributions between BVC-EWL and
BVC-NWL are small. The complex velocity distribution of the experiment around the boulders
could not be reproduced perfectly by the methods
assuming a cubic velocity vertical distribution,
especially for separation zone behind boulders(x/dc = −4). The resolution of this issue requires
a high-resolution vertical distribution function
which needs additional momentum equations.
A high-resolution three dimensional calcula-

Figure 6.

tion method with advanced turbulence modelis required for the precise calculation of flow
around separation zone. However, calculation
velocity vertical distributions capture the wholecharacteristics of those of the measurements
except for local velocity variations due to the flow
separation.
Figure 7 shows the velocity distributions
between boulders (y/dc = 1). We can see BVCNWL provides a better agreement with the measurements than the BVC-EWL. In the experimental
velocity distribution, the flow in the vicinity of the
bed surface is disturbed due to the boulder presence which increase the mass and momentum
exchange across the vortex and roughness layers.
In other words, the presence of isolated boulders
on rough bed increases the total resistance not
only by drag force acting on the boulders but also
by additional resistance acting on the bed due tothe enhanced momentum exchange on bed (secondary resistance). Since BVC-EWL assumes
equilibrium velocity distribution and no verticalflux, it underestimates deformations in velocity
vertical distribution but overestimates the bottom
velocity.
The reduction in the computational error
by combining the BVC method with the non­equilibrium wall law is shown in Figure 8 for thelongitudinal cross sections across (y/dc = 0, 3) andStream-wise velocity vertical distribution along longitudinal section across boulders (y/dc = 0).

341

Figure 7. Stream-wise velocity vertical distribution along longitudinal section between boulders (y/dc = 1).

5 conclusion

Figure 8. Calculation errors of bottom velocity by several calculation methods.

between boulders (y/dc = 1, 2). The computational
error is estimated as the root mean square of­velocity difference between measured and computed results normalized by the bulk velocity U0.
The improvement by calculating vertical velocity
distribution (2D—BVC-EWL) reflects the effectof the drag forces on velocity distributions, which
are dominant for the line across-boulders. The
error reduction for BVC-NWL compared to BVCEWL indicates the effects of secondary resistancethat the isolated boulders exert. As discussed earlier, for flow between boulders, it is important toevaluate not only the boulders drag force but also
the boulders secondary effect due to the enhanced
vertical momentum exchanges.

The non-hydrostatic quasi-3D calculation method(BVC method) with non-equilibrium wall law is
proposed to calculate flow over rough bed withboulders. The present method has three calculational domains in the vertical direction in order
to separate the effects of Large Scale Roughness(LSR) from that of Small Scale Roughness (SSR),
namely: a main domain that takes into accountthe form effect of boulders, a thin vortex layer on
the bed surface and a roughness layer under the
bed surface. The main domain between the water
surface and the vortex layer is calculated via the
BVC method. For the non-equilibrium wall law,
momentum and continuity equations in vortex
and roughness layers are derived and solved withequations of the BVC method. Since the present
method is physically based and reduces to the conventional wall law for rough bed for the uniform
flow condition, we can use same parameters for theconventional equivalent wall law, such as equivalent roughness height.
The presence of isolated boulders atop a roughbed increases flow resistance not only by the dragforce acting on the boulders but also by producing
secondary resistance, which acts on the rough bed
due to enhancing momentum exchange across thevortex and roughness layers. The BVC method can342

in general take into account the drag force by calculating flows around submerged boulders. However,the BVC method with Equilibrium Wall Law (BVCEWL) underestimates the velocity profile deformations especially in-between boulders, which is not
the case when the Non-equilibrium BVC method(BVC-NWL) is used. The BVC method with Nonequilibrium Wall Law (BVC-NWL) demonstrates
the ability to calculate velocity distribution overrough bed reach with isolated boulders.
references
Aguirre-Pe, J. & Fuentes, R. 1990. Resistance to flow insteep rough streams, Journal of Hydraulic Engineering
116(11): 1374–1387.
Bathurst, J.C. 1985. Flow resistance estimation in mountain rivers, Journal of Hydraulic Engineering 111(4):625–643.
Carney, S.K., Bledsoe, B.P. and Gessler, D. 2006. Representing the bed roughness of coarse-grained streams
in computational fluid dynamics, Earth Surface Processes and Landforms 31(6): 736–749.
Clifford, N.J., Robert, A. and Richards, K.S. 1992.
­Estimation of flow resistance in gravel-bedded rivers:
A physical explanation of the multiplier of roughness
length, Earth Surface Processes and Landforms 17(2):111–126.
Fukuoka, S. and Uchida, T. 2013. Toward integrated
multi-scale simulations of flow and sediment transport in rivers, Journal of JSCE Ser.B1 (HydraulicEngineering) 69(4): II-1–II-10.
Fukuoka, S., Watanabe, A. and Nishimura, T. 1992. On
the groin arrangement in meandering rivers, Journalof Hydraulic Engineering JSCE 443/II-18: 27–36.
Ishikawa, T., Suzuki, K. and Tanaka, M. 1986. Efficient
numerical analysis of an open channel flow with secondary circulations, Proc. of JSCE 375/II-6: 181–189.
Jin Y.-C. and Steffler, P.M. 1993. Predicting flow incurved open channels by depth-averaged method,
J. Hydraul. Eng. 119(1): 109–124.
Katul, G., Wiberg, P., Albertson, J. and Hornberger G.
2002. A mixing layer theory for flow resistance in shallow
streams, Water Resources Research 38(11): 321–328.
Nicholas, A.P. 2001. Computational fluid dynamics
modelling of boundary roughness in gravel-bed rivers:
an investigation of the effects of random variability
in bed elevation, Earth Surf. Process. Landforms 26:
345–362.

Nikora, V., McEwan, I., McLean, S., Coleman, S.,­Pokrajac, D. and Walters, R. 2007a. Double-­Averaging
Concept for Rough-Bed Open-Channel and Overland
Flows: Theoretical Background, J. Hydraul. Eng 133:
873–883.
Nikora, V., McLean, S., Coleman, S., Pokrajac, D.,­McEwan, I., Campbell, L., Aberle, J., Clunie, D.
and Koll, K. 2007b. Double-Averaging Concept for
Rough-Bed Open-Channel and Overland Flows:
Applications, J. Hydraul. Eng 133: 884–895.
Olsen, N.R.B. and Stokseth, S. 1995. Three-dimensional
numerical modelling of water flow in a river with large
bed roughness, Journal of Hydraulic Research, 33(4):571–581.
Papanicolaou, A.N., C.M. Kramer, A.G. Tsakiris,
T. Stoesser, S. Bomminnayuni, and Z. Chen. 2012.
Effects of a fully submerged boulder within a boulder array on the mean and turbulent flow fields:­Implications to bedload transport, Acta Geophysica.
60(6): 1502–1546.
Rameshwaran, P., Naden, P. and Lawless, M. 2011. Flowmodeling in gravel-bed rivers: rethinking the bottom
boundary condition, Earth Surf. Process. Landforms
36: 1350–1366.
Shlichting, H. and Gersten, K. 2000. Boundary Layer
Theory, 8th Revised and Enlarged Edition, Springer,
Verlag Berlin Heidelberg.
Uchida, T. and Fukuoka, S. 2011. A bottom velocity
computation method for estimating bed variation in achannel with submerged groins, Journal of JSCE, Ser.
B1 (Hydraulic Engineering) 67(1): 16–29.
Uchida, T. and Fukuoka, S. 2012. Bottom velocity computation method by depth integrated model without
shallow water assumption, Journal of Japan Society of
Civil Engineers, Ser. B1 (Hydraulic Engineering) 68(4):I_1225–I_1230.
Uchida, T. and Fukuoka, S. 2013. Quasi 3D numericalsimulation for flow and bed variation with various
sand waves, Advances in River Sediment Research,
Proceedings of 12th International Symposium on River
Sedimentation, ISRS, Kyoto, Japan: 221–229.
Yeh, K.-C. and Kennedy, J.F. 1993. Moment modelof nonuniform channel-bend flow. I: Fixed beds,
J. Hydraul. Eng. 119(7): 776–795.

343