River Flow 2016 – Constantinescu, Garcia & Hanes (Eds)© 2016 Taylor & Francis Group, London, ISBN 978-1-138-02913-2

A 3D flow calculation around a structure in gravel bed rivers by using
non-hydrostatic depth-integrated model with dynamic rough-wall law
T. Uchida & S. Fukuoka

Research and Development Initiative, Chuo University, Tokyo, Japan

A.N. Papanicolaou

Department of Civil and Environmental Engineering, The University of Tennessee, Knoxville, TN, USA

ABSTRACT: A reliable depth-integrated model is required for practical application to rivers. This paper
presents a non-hydrostatic depth-integrated model of the general Bottom Velocity Computation (BVC )method that employs dynamic wall law (BVC-DWL) to calculate 3D flow calculations around a structurein gravel bed rivers. We applied the BVC-DWL method to a local 3D flow around a structure measured ina field experiment using an enhanced measurement method with the Acoustic Doppler Current Profiler(ADCP). The BVC-DWL method is validated for local 3D flow around structures in gravel bed rivers
though comparisons with measured data, demonstrating the advantages of introducing the DWL insteadof the conventional equilibrium wall law. This paper presents the non-equilibrium flow structures near thegravel bed around a structure in the calculation results, and discusses the limitations of the conventionalrough-wall law and the necessity of introducing a dynamic rough-wall law.
1 introduction
Recent advances in computational technology anddevelopments of numerical schemes allow us tocalculate 3D flows around structures in rivers. The
conventional equilibrium wall law that assumes anequilibrium flow condition near the bed and employs
a function of shear velocity and lowest velocity
is not suitable as a boundary condition for gravel
bed rivers, since the non-equilibrium flow near the
bed causes vertical momentum fluxes near the bed.
Some researchers have developed new methods witha 3D flow calculation model that treats bed material as a permeable resistance, instead of assuming
the impermeable surface bed for the conventional
method, and calculates the non-equilibrium motion
near the bed surface and in the roughness layer
under the bed surface (Olsen and Stokseth, 1995;
Lane et al., 2004; Carney et al., 2006; Nicora et al,
2007a, 2007b; Rameshwaran et al., 2011).
While we recognize the need for a non-hydrostatic three-dimensional method for this purpose, areliable depth-integrated model is also required forpractical applications to rivers. In fact, the applications of a 3D model to river flows have been limited to small-scale laboratory experiments owing
to its large computational cost. However, for a
depth-integrated model, which has been developed for applications to actual river flows, there
remains the challenge of developing an evaluation

method for non-equilibrium fluid motions near thegravel bed, as indicated above for a 3D calculation
method. One reason for this is that most previous
depth-integrated models, including quasi-threedimensional models, have been simplified with theassumption of the hydrostatic pressure distribution (Ishikawa et al., 1986; Fukuoka et al., 1992;
Jin & Steffler 1993; Yeh and Kennedy, 1993).
Recently we developed a new quasi-three-dimensional calculation method without a shallow-water
assumption such as hydrostatic pressure distribution. This new approach is called the general
Bottom Velocity Computation (BVC) method, inwhich equations for velocity and pressure distributions in the vertical direction are solved with depth
averaged velocity equations (Uchida & Fukuoka,
2012, 2014). The BVC method has been validated
for several flows including rapidly varied flows overa structure and three-dimensional flows around a
non-submerged structure where the horseshoe vortex is prominent (Fukuoka & Uchida, 2013). However, it is still challenging to calculate the velocity
distribution near the bed for the 3D flow aroundstructures by using the depth-integrated model
including the BVC method. The key to overcoming
this issue is assumed on the bottom boundary conditions for depth-integrated equations of the BVC
method. Uchida et al. (2014) proposed the dynamic
wall law for rough bed (DWL) by introducing continuity and momentum equations for vortex and1339

roughness layers. The BVC method employing
the dynamic wall law (BVC-DWL) demonstrated
the ability to calculate velocity distribution over arough bed with isolated boulders in an experimental channel. However, the BVC-DWL method has
not been validated for flows in gravel bed rivers.
The objective of this paper is to discuss thevalidity of the BVC-DWL method for a 3D flow
around a structure in a gravel bed river. We applied
the BVC-DWL method to the local 3D flow arounda structure measured in a field experiment with theenhanced measurement method using the Acoustic
Doppler Current Profiler (ADCP) (Koshiishi et al.,2012; 2013). The calculation results of the BVCDWL method were compared with the measured
data, demonstrating the advantages of introducing
the DWL instead of the conventional equilibrium
wall law (BVC-EWL).

Figure 1. Framework of the present flow calculation
method with dynamic rough wall law (BVC-DWL).

2 calculation methods
2.1 Framework of BVC-DWL method
In this method, the computational domain is separated in the vertical direction into the main computational domain and the wall-law layer near the bed,
as shown in Figure 1. For the main computational
domain between the water surface and the top of thewall-law layer, we refined the BVC method, which isan advanced non-hydrostatic depth-integrated modelcomposed of horizontal momentum equations andseveral equations for vertical velocity and pressure
distributions (Figure 2), to take into account the
momentum transfer at the top of the wall law layer
owing to non-equilibrium flow near the bed. The
wall law layer is composed of a vortex layer and aroughness layer. As a dynamic boundary conditionof the rough bed, this method employs continuity
and momentum equations for these layers (Table 1).
The momentum equations have been derived so that
this method reduces to the conventional equilibrium
rough-wall law for the uniform flow condition, as
indicated in the following section.
2.2

Equations for main computational domain(general bottom velocity computation
method)The governing equations for the main computational domain of the BVC method are derived
from the continuity equation and the Reynoldsaveraged Navier–Stokes equations assuming avortex layer δzb on the bed surface and a bottom
velocity on the layer, as indicated in Figure 2. The
depth-integrated continuity equation and horizontal momentum equations with vertical flux on thebottom are given as follows:

Figure 2. Unknown variables in the main calculation
domain of BVC-DWL method.

∂h ∂U j h+
− wσ b = 0
∂t
∂x j

(1)

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

(2)

i

j

i

i

where i, j = 1(x),2(y); z: vertical direction; τbi:bottom shear stress; Tij: horizontal shear stress;and dp0 = dpb/2. The horizontal shear stresses Tij

1340

SBB5_FP196_57.indd 1340

5/23/2016 10:59:12 PM

Table 1. Unknown variables and equations.

Main computational
domainVortex layer
Roughness layer

Unknown variables/equations

Equationsh: water depth

DI continuity eq. (1)usi: water surface velocity
Ui: depth averaged horizontal velocity
W: depth averaged vertical velocity
Ωi: depth averaged horizontal vorticity
ubi: bottom velocity
dpb: bottom pressure deviation from thebottom hydrostatic pressure
wσb: bottom vertical flux
uvi: horizontal velocity for vortex layer
wσt: bed surface vertical flux
uri: horizontal velocity for roughness layer

Horizontal momentum eq. on water surface (9)DI horizontal momentum eq. (2)Double integral of the continuity eq. of z (10)DI horizontal vorticity eq. (7)DI of definition eq. of horizontal vorticity (8)DI vertical momentum eq. (11)Continuity eq. for vortex layer (12)Momentum eq. for vortex layer (13)Continuity eq. for roughness layer (12)Momentum eq. for roughness layer (14)

DI: depth integrated

consist of shear stresses owing to molecular andturbulent motions and horizontal momentum
transportation by vertical velocity distribution.
Tij = τ ij + ui′ u ′j , τij = vSij, v = vm + vt, vm: kinematic
viscosity coefficient, vt: the eddy viscosity coefficient, Sij: depth-averaged strain velocity, andui = U i + ui′ (ui: horizontal velocity). In this paper,an overline indicates a depth-integrated value. For
the eddy viscosity coefficient, this study employs aone-equation model taking into account the velocity vertical distribution in the production term Pk(Uchida & Fukuoka, 2014):∂k∂k 1 ∂  vh ∂k 
+U j=
+ Pk − ε∂t∂x j h ∂xi  σ k ∂xi (3)where vt = Cµk2/ε, Cµ = 0.09, ε = Cεk3/2/∆, Cε/∆ = 1.7/h(Nadaoka & Yagi, 1998).
The cubic vertical velocity equation (4) isassumed by using the depth-averaged velocity Ui,
the bottom velocity ubi, and the water surface velocity usi with the assumption of zero vertical velocity
gradients at the water surface.

(

)

(

ui′ = ∆ui 12η 3 − 12η 2 + 1 + δ ui −4η 3 + 3η 2)

(4)where δui: usi −ubi, ∆ui: usi −Ui, and η = (zs −z)/h. The
quadric curve of velocity distributions is used forthe uniform flow, substituting ∆ui = δui/3 (usei = Ui +
δui/3) in equation (4) (Uchida & Fukuoka, 2014).
To evaluate the vertical velocity and pressure distributions, the following equations are solved forthe main calculation domain of the BVC method:∂hDω ij∂Ω i h= Rσ i + Pω i ++ wσ bω bi
∂t
∂x j

(7)

 ∂Wh
∂z
∂z 
ubi = usi − ε ij 3Ω j h − − ws s + wb b 
∂xi
∂xi 
 ∂xi

(8)

∂usi
∂u
∂z
+ usj si = − g* s + Psi
∂t
∂x j
∂xi

(9)

 ∂z
∂z  ∂h2 ηui′
Wh = h  m + U i m  −+ hwσ b
∂xi 
∂xi
 ∂t

(10)

∂hτ zjdpb ∂hWU j∂z=+ τ bj b −
ρ
∂x j
∂x j
∂x j

(11)

where εijk: Levi-Civita symbol, and ws, wb: vertical velocity on the water surface and bottom. The
above equations with a time variation term calculate a time-advance scheme for each unknown variable in Table 1. The equations (8), (10), and (11),which do not include time variation terms, are calculated as follows. Before evaluating equation (11)for the bottom pressure difference dpb, the velocity
difference between the water surface, the bottom
δui, and the vertical velocity integrated over thewater depth Wh at the next time step are calculated to satisfy equations (8) and (10) at the n + 1
step, introducing a similar approach to the SMAC
scheme for incompressible fluid. The details of thegoverning equations and numerical scheme for themain calculation domain of the BVC method may
be found in Uchida & Fukuoka (2014).
2.3

Dynamic rough wall law for the boundary
conditions of the main computational domain

Boundary conditions on the bottom for the main
calculation domain should be given in the bed
shear stress terms in momentum equations (2),and the production terms in vorticity equations (7)1341with vertical flux wσb. The BVC method employing the equilibrium wall law (Uchida & Fukuoka,2014) assumes the equilibrium velocity and vortex
vertical distribution for the vortex layer, neglecting
the vertical flux (wσb = 0). Those distributions forthe vortex layer have been decided with the log law
and bottom velocity. For the dynamic rough walllaw, those boundary conditions are given dynamically with continuity and momentum equations
for the vortex and roughness layers (Uchida et al.,2014), as shown in Figure 2.
wσ b = wσ t −∂δ zbuvi∂λδ zt uri, wσ t = −
∂xi
∂xi

(12)

∂uvi
τ
∂u
∂( pb + ρ gzb ) τ bi
+ uvk vi = −+− ti (13)∂tρ∂xiρδ zb ρδ zb
∂xk
∂uri
τ
∂u
∂( pb + ρ gzb )D+ urk ri = −+ ti − i (14)∂tρ∂xiρδ zr ρδ zr∂xk
where λ: porosity in the roughness layer (λ = 0.4);
k = 1,2,3; pb: bottom pressure pb = ρgh + dpb; τbi,
τti: shear stress acting on the bottom and roughness
layer surface; δzb, δzr: the thickness of the vortex
and roughness layers; Di: resistance term owing todrag force acting on roughness. The shear stressacting 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 shear stress and resistance terms in the
momentum equations for the vortex and roughness layers are derived so that the DWL reduces tothe EWL for the equilibrium condition. The shear
stress acting on the vortex layer can be assumed
based on the Boussinesq approximation:

τ bi  ∂ui (u − u )= v= v ⋅ A bi viρ  t ∂z  b tb bh(15)

The coefficients Ab are decided to satisfy τb = ρu*2
(τb2 = τbiτbi), u* = ub/cb = uv/cv, and vtb = αu*h for theequilibrium condition by introducing the conventional logarithmic velocity law:
1 / Ab = α (cb − cv ) = (α / κ ) ln ( zb / zv )(16)

where zb = δz0 + δzb, zv = δz0 + δzb/2, δz0: depth ofthe log law origin from the top of the roughness
height, δz0 = 0.3ks. The eddy viscosities for evaluating shear stress acting on the bottom are expressedby using the distance of the wall and velocity
difference:
vtb = α bδ ub h, δ ub 2 = (ubi − uvi ) (ubi − uvi )(17)

where the coefficient αb is given by αb = α/(cb – cv) = Ab
to reduce vtb = αu*h for the equilibrium condition.
However, using bottom vorticity ωb gives anotherexpression of the bottom eddy viscosity vtb:
vtb =ακω b h22 ln( zs / zb )(18)

This study uses the larger value between equations (17) and (18).
In a similar manner, the shear stress and theeddy viscosity on the bed surface (the top of theroughness layer) are given asτ ti(u − uri )  h + δ zb 
= vtt ⋅ At vi


h
h 
ρ

(19)

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

1 / At = 1 / α t = α (cv − cr )(21)

where cv = (1/κ)ln(zv/ks) + Ar, Ar = 8.5, u* = ur/
cr. For the equilibrium condition, since the righthand side of equation (14) is zero, the magnitudeof the resistance term should be represented by
D = τt + ρgδzrS (S: channel slope). The coefficient
cr is evaluated with the resistance term of equation(14), when we introduce a drag coefficient to represent the resistance term:cr =8  δ zr + δ zb 1+hπ CD 

(22)

Di =

π CD
ρur uri
8(23)

where CD = 0.4.
The production terms Pωi in equation (7) are
given by equation (24) (Uchida & Fukuoka,2014).
Pω i = C pω vtb (ω bei − ω bi ) h + Psω i(24)

where Cpω = κ /α, ωbi: horizontal vorticity on bottom, ωbei: equilibrium ωbi, Psωi: vorticity production
owing to flow separation (Uchida & Fukuoka, 2011,2014). The equilibrium vorticity on the production
term in the vorticity equation (11) is assumed byusing a velocity gradient in vorticity in the layer asω bej = ε ij 3Aω

ubi − uvih(25)

where the coefficient Aω is represented by equation(26) so that equation (25) reduces to equation (27)for the equilibrium condition.

1342

Figure 3. Plan view of the experimental channel (Koshiishi et al., 2012, 2013).

Aω =

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

ω be =

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

3 field measument and 3D local
flow measurement using ADCP
A field experiment in the Jyoganji River in 2011 was
conducted to investigate characteristics of erosion
and bed scouring beside a revetment works (Koikeda
et al., 2012), and to obtain detailed data to validate
the new depth-integrated calculation method for bed
scouring in a flow attacking area (Koshiishi et al.,2012, 2013). The experiments were conducted in acompound meandering channel with total channel
length: 170.0 m, total channel width: 20.0 m, main
channel width: 4.0 m, main channel depth: 1.0 m,
and bed slope: 1/200. The experimental discharge
was 7.6 m3/s. There are vertical revetments made withconcrete walls in two flow attacking zones aroundmeandering sections (Figure 3). The three-dimensional components of velocities by means of ADCP
were measured around the flow attacking zone 1.
The conventional method of ADCP calculates avelocity vertical profile with four different direction
beams (Teledyne RD Instruments, 2006). It cannot
give accurate local velocity vectors for flows withlarge velocity gradients and deep water depth. Thisis because of velocity variations within the focus
spot L surrounded by four beams, which expands
toward the vertical direction by an angle of 20º, asshown in Figure 4.
To obtain accurate local 3D velocity distribution around a structure, individual beam data ineach direction for current profiles of the ADCP
are managed as follows (Koshiishi et al., 2012,2013). Koshiishi et al. (2012, 2013) calculated thelateral and vertical velocity at any point using
superimposing spatially interpolated individual
beam data, which are measured continuously inthe lateral direction and are averaged over time asshown in Figure 5. Variations in the stream-wise

Figure 4. Concept of the conventional ADCP method(Koshiishi et al., 2012, 2013).

Figure 5. Concept of the present ADCP method(Koshiishi et al., 2012, 2013).

velocity component are assumed to be negligibly
small compared with those of the lateral and vertical directions. Koshiishi et al. (2012) have shownthat this method could measure 3D local flows infront of a revetment with accuracy, especially inthe vicinity of the bed surface. The details of themeasuring method and experiments may be foundin Koshiishi et al. (2012, 2013).
4 computed results anddiscussions
We applied the BVC-DWL to the real scale 3D
flow in front of a structure measured in the Joganji
River field experiment. Figure 6 shows a comparison of the water surface and bottom velocity
between measured (Koshiishi et al., 2012, 2013)and calculated (BVC-DWL) results. We can see in1343

Figure 7. Comparisons of calculated velocity distributions in the roughness layer between BVC-EWL and
BVC-DWL.

Figure 6. Comparisons of water surface and bottom
velocity between measured (Koshiishi et al., 2012, 2013)and calculated (BVC-DWL) results.

Figure 8. Comparisons of local 3D velocity distributions on flow attacking zone between measured(Koshiishi et al., 2012, 2013) and calculated (BVC-EWL,
BVC-DWL) results.

the measured results that the bottom velocity vectors were bent toward the inner bank, while thewater surface velocity vectors move in the directionof the outer bank. These characteristics measured
in the field experiment are well explained by thecalculation results of the BVC-DWL method.
Figure 7 shows comparisons of calculated velocity distributions in the roughness layer between
BVC-EWL and BVC-DWL. The large difference
that occurs by introducing DWL can be seenaround outer bank, in which the roughness layer
velocity of BVC-DWL is bent toward the innerbank owing to the strong downward flux of the secondary flow, while the velocity of BVC-EWL does
not have as much of an inner bank component.
The above numerical investigation indicates theflow pattern under the bottom varies considerably
for the local 3D flow around a structure in gravel
bed rivers. Figure 8 shows comparisons of local
3D velocity distributions on a flow attacking zone
between measured (Koshiishi et al., 2012, 2013)and computed (BVC-EWL, BVC-DWL) results.
In the calculated results, bed heights are set to zt,
and the velocity distribution between zs and zt isshown. The velocity vectors on the bottom zb are
third vectors from the bed.

The measured results indicate the 3D flow characteristics as follows. We can find weak secondary flows near the bottom on the right bank forthe approaching flow (Section B). For the flow
attacking zone of Sections C–E, clockwise secondary flows are developed in the entire cross section.
For Sections F and G, large stones bend flows near
the bed toward the inner bank. The above 3D flow
characteristics in the experiment are explained by
the calculated results using BVC-EWL. However,the calculated secondary flows are underestimatedby BVC-EWL. The main velocity areas do not shift
to the outer bank for Sections C–F compared withthose of the experiment. The calculated results of
BVC-DWL demonstrate that the above problems of
BVC-EWL are considerably improved by introducing DWL, even though some differences between
BVC-DWL and measurements can be found.
5 conclusions
This paper presents a non-hydrostatic depth-integrated model, the general Bottom Velocity Computation method, which employs the dynamic
rough wall law (BVC-DWL). The local 3D flows1344

around a structure in a gravel bed river measured
with the enhanced measurement method using
the Acoustic Doppler Current Profiler (ADCP)
are compared with the BVC-DWL method. The
water surface and bottom velocity distributions
calculated by BVC-DWL are in good agreementwith those of the measurement. The roughness
layer velocity distribution with BVC-DWL differsfrom the calculation result that employs the conventional equilibrium rough wall law (BVC-EWL)
around a structure owing to the downward fluxof the secondary flow. This paper presents local
3D velocity distributions in several cross sections
around the structure in the gravel bed river measured by enhanced ADCP. The detailed comparisons with measured results indicate the problemsof BVC-EWL in which secondary flows are underestimated and the main velocity areas do not shift
to the outer bank. Although there are some differences between BVC-DWL and the measurements,
the BVC-DWL is able to improve the above problems of BVC-EWL considerably.
references
Fukuoka, S., Watanabe, A., and Nishimura, T. (1992)“On the groin arrangement in meandering rivers”
Journal of Hydraulic Engineering, JSCE, No. 443/
II-18, pp. 27–36, in Japanese.
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.
Koikeda, S., Ishii, A., Iwai, H., Ishikawa, T., and Fukuoka, S. (2012) “Nature friendly bank protection works
using sandbars with boulders in the Jyoganji river,”
Advances in River Engineering 18: 233–238, JSCE, inJapanese.
Koshiishi, M., Uchida, T. and Fukuoka, S. (2012) “A
New method for measuring local Flows around riverbank using ADCP” Proceedings of the 10th International Conference on Hydroscience & Engineering,
Orlando, FL, USA.
Koshiishi, M., Uchida, T. and Fukuoka, S. (2013) “A new
method for measuring and calculating flows and bed
forms around river banks” Advances in River SedimentResearch, Proceedings of 12th International Symposium
on River Sedimentation, ISRS-2016, Kyoto, Japan.
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.
Lane, S.N., Hardy, R.J., Elliott, L., and Ingham, D. B.
(2004) “Numerical modeling of flow processes overgravelly surfaces using structured grids and a numerical porosity treatment,” Water Recourses Research 40:
W01302.
Nadaoka, K. and Yagi, H. (1998) “Shallow-water turbulence modeling and horizontal large-eddy computation of river flow” Journal of Hydraulic Engineering
124(5): 493–500.
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 nu-merical modelling of water flow in a riverwith large bed roughness” Journal of HydraulicResearch 33(4): 571–581.
Teledyne RD Instruments (2006) “Acoustic Doppler
current profiler principles of operation a practicalprimer.”Uchida, T. and Fukuoka, S. (2011) “A bottom velocity
computation method for estimating bed variation ina channel with submerged groins” Journal of JSCE,
Ser. B1 (Hydraulic Engineering), 67(1): 16–29, inJapanese.
Uchida, T. and Fukuoka, S. (2012) “Bottom velocity
computation method by depth-integrated model without shallow water assumption” Journal of JSCE, Ser.
B1 (Hydraulic Engineering) 68(4): I_1225–I_1230, inJapanese.
Uchida, T. and Fukuoka, S. (2014) “Numerical calculation for bed variation in compound-meandering channel using depth in-tegrated model without assumptionof shallow water flow” Advances in Water Resources
72: 45–56.
Uchida, T., Fukuoka, S., Papanicolaou, A.N., andTsakiris, A.G.(2014) “A numerical calculation
method for flow in the presence of isolated boulders
atop a rough bed by using an enhanced depth-integrated model with a non-equilibrium resistance law”
Proceedings of International Conference on Fluvial
Hydraulics, River Flow, Lausanne, Switzerland, pp.
335–343, 2014.
Yeh, K.-C. and Kennedy, J.F. (1993) “Moment modelof nonuniform channel-bend flow. I: Fixed beds” J.
Hydraul. Eng. 119(7): 776–795.

1345