PRACTICAL NUMERICAL SIMULATION METHOD FOR ESTIMATINGLOCAL SCOURING AROUND A PIER ~ BOTTOM VELOCITYCOMPUTATION METHOD BASED ON DEPTH INTEGRATED MODEL
T. Uchida1 and S. Fukuoka2
1. Research and Development Initiative, Chuo University, 1-13-27 Kasuga, Bunkyo-ku, Tokyo, 1128551, Japan, (utida@tamacc.chuo-u.ac.jp)
2. Research and Development Initiative, Chuo University, 1-13-27 Kasuga, Bunkyo-ku, Tokyo, 1128551, Japan, (sfuku@tamacc.chuo-u.ac.jp)

ABSTRACTIn this paper, for estimating local scouring around a pier, the Bottom Velocity Computation (BVC)
method is developed based on the depth integrated model. The governing equations of the BVC
method are composed of depth-integrated horizontal vorticity and water surface velocity equations inaddition to shallow water equations and a depth averaged turbulence energy transport equation. To
compute bed variation, non-equilibrium bed load equation is proposed. The adequacy of the model isdiscussed through the comparisons with the laboratory experimental results for flows and bed
variations around a pier. The pier scour is induced mainly by the horseshoe vortex. The horseshoe
vortex with radial diverging bed-surface flow from the top of the pier is produced by the extension androtation of the horizontal vortex filament in approaching flow. On the other hand, bed-surface flow
converges and deposits sediment downstream of the pier. This paper shows that these important flow
characteristics around a pier are reproduced by the BVC method. And computed scour hole around apier is similar to that of the measurement. The above indicates the BVC method with the bed variation
analysis by non-equilibrium sediment transport model is useful to simulate the severe local scour infront of structures.

1.

INTRODUCTION

River engineers are faced with urgent task of checking the safety degree against floods, because thesafety level of rivers in Japan is predicted to be reduced by the global warming. A numerical model for
bed variations in rivers during floods is required to assess the risk of river structures due to local
scouring.
Local scour around a pier has attracted a lot of interest from engineers and researchers for a long time.
Recently, the horseshoe vortex and local scour around a pier have been calculated successfully by full3D turbulence models coupled with sediment transport models (Olsen & Melaaen, 1993; Olsen &Kjellesvig, 1998; Nagata et al., 2005; Roulund et al., 2005). However, full 3D turbulence models arestill limited for application to floods and bed variations in natural rivers. So, a practical and reliable
model for flows and bed variations during floods are highly required.
Depth integrated models have been used for computations of flood flows and bed variations in rivers(Fukuoka, 2005; Wu, 2008). A 2D model based on the shallow water equations does not give
sufficient accuracy in calculating local scouring around outer bank in curved channels, contracting
sections, river confluences, and river structures. For such conditions, bottom velocity (velocity near
bed surface) which determines sediment motions is different from that for uniform flow, as shown in
Figure 1. For this reason, many improved depth integrated models have been developed, especiallyfor a helical flow in a curved or meandering channel. Those include full-developed secondary flow
model without considering effects of secondary flow on vertical distribution of stream-wise velocity (e.g.
Engelund, F., 1974; Nishimoto et al., 1992), refined secondary flow models with non-equilibrium (e.g.
Ikeda & Nishimura, 1986; Finnie et al., 1999) and mean flow redistribution effects (e.g. Blanckaert &de Vriend, 2003). As more versatile method for evaluating vertical velocity distributions, a quasi-3D
model, in which equations for vertical velocity distribution are computed, has been developed (e.g.

1Ishikawa et al., 1986; Fukuoka et al., 1992; Jin & Steffler, 1993; Yeh & Kennedy, 1993). The above
mentioned models are categorized in the efficiently-simplified 3D model with the assumption ofhydrostatic pressure distribution for each objective of the computation, as shown in Figure 2. However,
it is known that these models are still-inadequate for local scouring induced by complex flows around apier, because they cannot consider effects of flows due to non-hydrostatic pressure distribution on
vertical velocity distributions.

Depth averaged flow path
Secondary flow

Bed surface flow path

 Sediment path

3D modelFigure 1

Helical flow due to the centrifugal force in curved channels

Full 3D model
Hydrostatic pressure distribution[Abbreviated vertical momentum equation]

3D model with the assumption ofhydrostatic pressure distribution

Depth integrated model

Functions of vertical velocity distribution

Quasi-3D model (hydrostatic)

Quasi-3D model (hydrostatic)

Helical flow model due to centrifugal force
in curved flow

Secondary flow model
Abbreviated vertical velocity distributions

2D modelFigure 2

Categorization of depth integrated models by derivations

Ghamry & Steffler (2002) proposed quasi-3D model with non-hydrostatic pressure distribution basedon moment equations. However, the model capability for complex flow fields around structures isremain incompletely understood. Recently, we had developed a new depth integrated model without
the assumption of the vertical pressure distribution, in which shallow water equations and horizontal
vorticity equations are solved to compute vertical velocity distributions. It was demonstrated that themodel can simulate complex flow fields induced by a tributary entry (Uchida & Fukuoka, 2009). In thispaper, for estimating local scouring around a pier, the Bottom Velocity Computation (BVC) method isdeveloped based on the depth integrated model (Uchida & Fukuoka, 2009). To compute bed variation,
non-equilibrium bed load equation is proposed and combined with the BVC method. And theadequacy of the model is discussed through the comparisons with the laboratory experimental results
for flows and bed variations around a pier.

2

2.
2.1

BOTTOM VELOCITY COMPUTATION METHOD
Basic concept

Figure 3 shows Bottom Velocity Computation (BVC) method. We define the bottom zb as the surface
on the shin vortex layer zb on the bed, as shown in Figure 3. The velocity on the bottom ubx iscomputed by the BVC method. With Stokes' theorem, velocity difference ux between water surface usx
and bottom surface ubx is represented by depth-integrated vorticity yh and spatial difference of depthintegrated vertical velocity Wh. The non-dimensional form is :

z

u sx

zs

zs

Ux

Flow
k

ubx

h

y

Psy

Pby

zb
u x   y hFigure 3

zbx

The concept of Bottom Velocity Computation Methodu * x  * y h*  h0 / L0 2 W *h* / x*(1)where, *y=yh/U0, h*=h/h0, x*=x/L0, u*x=ux/U0, W*=W/W0=Wh0/U0L0, h0: representative water depth(vertical scale of flow), L0: horizontal scale of flow, U0: representative horizontal velocity, W0:
representative vertical velocity. For shallow water flow, the second term is neglected. That means, formany flow fields in river, we can evaluate the bottom velocity ubx by water surface velocity usx anddepth-integrated horizontal vorticity yh.
The governing equations of the BVC method are composed of depth-integrated horizontal vorticity andwater surface velocity equations in addition to shallow water equations and a depth averaged
turbulence energy transport equation. The important advantage of the BVC method is the ability tocompute bottom velocity distribution of the complex flow field around the river structure due to nonhydrostatic pressure distribution. The local scouring around the structure is calculated by using thebottom velocity with continuity equation for sediment and non-equilibrium sediment transport
equations.
2.2


Vertical velocity distribution and governing equations
Vertical velocity distribution and bed shear stress

Velocity distributions above very thin vortex layer on bed are computed in the BVC method. The
Polynomial equation of vertical velocity distribution is assumed. Then, by using depth averaged
velocity Ui, velocity difference ui between water surface and bed surface, water surface velocity usi,
and no vertical velocity gradients on water surface, the vertical velocity distribution of a cubic curve isobtained:ui  U i  ui 12 3  12 2  1  ui 4 3  3 2

3



(2)

where, ui usi Ui , ui usi ubi , =zszh, zs: water level, h: water depth. The uniform velocity
distribution is derived by equilibrium water surface velocity usei Uiui ui 2 and equilibrium velocity
difference uei. For the equilibrium condition, shear stress zi is written as: zi   bi  v

dui
2vui


dz
h

(3)

where, bi i direction bed shear stress defined as: bi  C0 2U iU  cb 2ubiub(4)

where, U2 Ui Ui, ub2 ubi ubi. ui for the equilibrium condition is given by:ui  C0 2U iUhh2 cb ubiub
2v
2v

(5)

With the depth averaged eddy-viscosity vu*h, 6, the relationship between cb and C0 is describedas:cb  C0 / 1  2C0 /  (6)

C0 is described by manning roughness coefficient, C02gn2/h1/3.
Water surface velocity

The equations for the water surface velocities are derived by assuming very thin layer just below water
surface zs (see Figure 3), neglecting water surface curvature:

usi
u
z
 usj si   g s  Psi
t
x j
xi

(7)

where, Psi , shearstress acting on zs, is represented by Equation (8).

Psi 

2vt
12C ps u sei  u si   3ui  6ui h2(8)where, Cps=3.


Depth averaged horizontal velocity and water depth

Equations for water depth h and depth averaged velocity vectors Ui (shallow water equations) withoccupancy ratio of water described as:h U j  h

0
t
x j

(9)

h ij
U i h U iU j hz   g s  0i  swi hthx jxi h Rsw hx j(10)where, swi =shear stress acting on side wall, Rsw=hydraulic radius of side wall. ij is horizontal shear
stress: ij  2 t Sij  ui ' u j '  2 / 3   ij  k

4

(11)

where, k=depth averaged turbulence energy, Sij=strain ratio tensor of depth averaged velocity vectors
Ui. The second term is produced by vertical velocity distribution, which is defined as:

ui ' u j ' 13ui u j  2uiu j  2u j ui  3uiu j

(12)

35

k is computed by Eq.(13), which is known as one-equation model:kk1   vh k  Dsw


U j

 Pk  tx j h xi   k xi  Rsw(13)

where, vCk2, Ck3/2, C0.09, C=1.7/h (Nadaoka & Yagi, 1998). With using vertical velocity
distribution (2), the production term Pk is derived as:



Pk
13
4
3

 4C
 2 S 2  S 2  S ijS ij  S 2   2h 8u 2  7uiui  2u 2
vt
35
35
35

 5h



(14)

where, S2=SijSij, S2=SijSij, S2=SijSij, u2=uiui, u2=ui ui , Sijstrain ratio tensor of ui , Sijstrain
ratio tensor of ui . Ch9C4C3, vte/uhvte: equilibrium depth averaged kinematic eddy
viscosity. The model is equivalent of zero-equation model for the equilibrium condition of k. And in thismodel, the last term of equation (11) is eliminated, because the normal stress is assumed ashydrostatic pressure distribution.


Depth integrated horizontal vorticity equations

The velocity difference ui is computed by using depth integrated vorticity for shallow water conditions,as discussed above:ui   ij 3 j h(15)

where, ij3permutation symbol (1232131,113223=0), idepth averaged vorticity vectors. The
equation for i is described by:

hDiji h ERi  P i 
t
x j

(16)

where,ERi  usis  ubib , Dij  U j i  U i  j   ' j u 'i   'i u ' j 
65

vt i ,  x j u j u j  , i   ij 3u j / h

2
5



 'i u ' j  iu j  i 
2.3

Bottom vorticity and production of vorticity from vortex layer

The production term Pi of the vorticity equation (16) represents net vorticity flux from the thin vortex
layer on bed. For flat bed condition, the flux is induced by the turbulent mixing (Pbi ). In addition, theflux is supplied by the separation downstream of the abrupt change of bed slope (Psi ), as shown in
Figure 3. The production term is defined as the sum:P i  Pb i  Ps i

5

(17)

Pbi is approximated by the product of the kinematic eddy viscosity and the vertical gradient ofvorticity near the bed:Pbi  C p vtb

bei  bih(18)

where, Cp3, vtb: equivalent depth averaged value of bottom kinematic eddy viscosity,
vtbmax2h2b, 2h2be, vt), b2bi bi, be2bei , vt: depth averaged kinematic eddy viscosity, bi:bottom vorticity, bei : equilibrium bottom vorticity for bottom velocity ubi, bei ij3cbubih, b: bottom
vorticity. The bottom vorticity is defined by Equation (19), minimizing the production term between
bi bci andbi bqi :bci bqi  0 :
bci bqi  0 :bi  bei bci  bqi : bi  bci bci  bqi : bi  bqi

(19)

where,  bcibeibci, bqi bei bqi, bci ui ui)h, bqii.
The flux supplied by the separation Psi is represented by the product of the bottom layer integrated
vorticity and the bottom layer averaged velocity: ubx / 2  (ubx ) (z x  tan  ) uby / 2  (uby ) (z y  tan  )Psx x  , Psy y  (z x  tan  )(z y  tan  )

0
0

(20)

where, x, y: computational grid size of x, y direction, zx=zxxx, zy=zyyy, zxx, zyy: second orderdifferential of bed level computational by x, y.
2.4

Numerical schemes

The explicit conservative CIP scheme for shallow water flows (Uchida, 2006) is adopted forcomputations of water depth and depth averaged velocity vectors. In the model, to directly capture theeffect of distributed parameters in the Cartesian coordinate system, each control volume ij has threekinds of variables, i.e., the value at the intersection of the grid (Point Value), the averaged value along
the side of the grid (Line-averaged Value) and the averaged value over the grid (Area-averaged
Value). Those values are computed all together based on CIP-CSL scheme (Nakamura et al., 2001).
All the variables in the governing equations are set on the same location. The utilization of multivaluables on a computational cell enables us to capture complex geometry even on the Cartesian
coordinate system. The above scheme is not adopted for computations of the equations for water
surface velocity, the depth averaged turbulence energy and vorticity. Those averaged value of thecomputational cell are computed by traditional scheme.

3.
3.1

COMPUTATIONAL METHOD FOR BED VARIATION
Continuity equation and momentum equation for sediment

The time variation of bed topography is computed by the continuity equation for sediment:(1  B )z qBi0t xi(21)

where, B: sediment porosity, qBi: bedload sediment transport rate. The equation for the bedload rate isderived by momentum equation for sediment:6q Bi u Bj q Bi Pu BPi  DuBDi   m* hB  ei   i 
t
x j

(22)

where, uBi: sediment particle velocity of bedload, PuBPi: gain momentum from bed material, DuBDi: loss
momentum by particle deposition, mksgcoss+1+CM, CM : coefficient of added mass, k: dynamic
friction coefficient of sediment, s: specific gravity of sediment in water,  : maximum grade of bed,
hBqBuB, qBqBiqBi, uBuBiuBi ,i, ei: i component of unit vector of sediment movement and its equilibrium
value. uB, uBehB are given by reference to Ashida & Michiue (1974).
The first term in right side for momentum exchange between bed load and bed materials is describedas:

Pu BPi  DuBDi  qBe uBei / Le  qBuBi / L(23)

where, qBe: equilibrium sediment transport rate of bed load, L, Le: average particle step length and thatfor equilibrium condition. The step length formulae are assumed by reference to Phillips & Sutherland(1989). The second term in right side is derived by forces acting on sediment particles.
3.2

Bed tractive force and equilibrium bed load

The total resistance of bed is evaluated as bed shear stress of equation (3) in the present model. The
total resistance of bed with sand waves is divided into form drag acting on sand waves and surface
resistance acting on sediment particles. The latter is considered as bed tractive force which induces
bed load (Ashida & Michiue, 1974). Generally, the surface resistance is evaluated by: bse  (cbsub )2(23)

where, 1/cbsr1lnzbksr8.5, ksequivalent roughness. However, the above is assumed as a
uniform flow condition and inadvisable in flow conditions around structures with local scouring. In thisstudy, the surface resistance is evaluated by the Boussinesq approximation near bed to evaluate nonequilibrium flow conditions:v  u  u   v c u  tb cv vte      tb bs bh
 z b vte  z b  bs  cv vtb (24)

where, cvvtb: bottom kinematic eddy viscosity, vte: equilibrium depth averaged kinematic eddy
viscosityThe bed tractive force and the critical shear stress of sediment particle on bed slope are
given by Fukuoka & Yamasaka (1983). The equilibrium bed load is given by Ashida & Michiue (1974).
And where the local bed slope exceeds the repose angle, the sand slide occurs in the model.

4.
4.1

APPLICATION TO LOCAL SCOURING AROUND A PIER
Experimental and computational conditions

The present model is applied to the experiment for local scour around a pier by Fukuoka et al. (1997)
to validate the model performance. Experiments were conducted in the channel with 27.5 m long, 1.5
m width, 1/600 initial bed slope and uniform sand material d=0.80 mm. A pier with 0.20 m diameter
was installed at the center of the cross section in the channel. The experimental discharge was 0.066
m3/s, and sediment supply rate at the upstream end was 8.8 cm 3/s. The average water depth h0 and
Fr number of the experiment were h0=11.6 cm and Fr=0.4, respectively.
The computational domain includes the overall channel width from 5m upstream to 6m downstream ofthe pier. From about 2.5m upstream to about 4m downstream of the pier, uniform grid ofdx=dy=0.025m is used, as shown in Figure 4. The experimental discharge and sediment supply rate
are given at the upstream end of the domain and the fixed water level is given at the downstream end.

7

The downstream end water level and manning roughness n of channel bed are decided to adjust
longitudinal water surface profiles between computation and measurement. And the equivalentroughness for bed tractive force is given by ksd. The initial flow condition of the bed variation
analysis is given by computation results for fixed flat bed with z=0 at the pier center section with initial
bed slope 1/600.
1.5
1

0.3

Flow

Water surface velocity distribution
 29grid
Mar 2010

Uniformof dx=dy=0.025
m

0.2

0.5

0.1

0

0
-0.1

-0.5

-0.2

-1
-1.5

pier

Grid around a pier
-0.3
-0.4

-3

-2

-1

0

1

Figure 4

2

3

4

-0.2

0

0.2

0.4

Grid around a pier5

The computational domain and grid
0.5 m/s

0.2
Water surface velocity distribution  09 Mar 2010 
Water surface velocity distribution  09 Mar 2010 

0

-0.2
-0.4

-0.2

0

0.2

0.4

0.6

0 min: Flat bed

Bed variation(cm)

0.5 m/s

0.2

Figure 5

0.2

0

0

-0.2

-0.2

-0.4

18

0.5 m/s 16-0.4 -0.2

-0.2

0

0

0.2

0.2

0.4

0.4

0.6 0.6

14
12
10
8
6
4
2
0
-2
-4
-6
-8
-10
-12
-14
-16
-18

18
16
14
12
10
8
6
4
2
0
-2
-4
-6
-8
-10
-12
-14
-16
-18

30 minThe computed bottom velocity fields on flat bed and scour bed after 30min.
0.2

0

-0.2
-0.4

Figure 6

-0.2

0

0.2

0.4

0.6

The measured bottom velocity fields for equilibrium stage (65 hours)(Fukuoka et al., 1997)

8

4.2

Computed results and discussion

Figure 5 shows computed temporally averaged bottom velocity fields on initial flat bed and computed
local scour after 30 minutes. Figure 6 shows the measured bottom velocity field around the pier for theequilibrium condition (after 65 hours). In the experimental results, it is fined that radial bottom flows
from the top of the pier induces local scour and the bottom flow runs around to the back of the pierwith sediment deposition. These important characteristics of the bottom velocity around the pier are
shown in the present model, especially for the radial bottom flow induced by the horseshoe vortex.
The radial bottom flow after 30 minutes by the present model is developed from that on initial flat bed
due to bed scouring. This phenomenon is similar as examined in previous experimental results by
Melville & Raudkivi (1977) and computed results by the full 3D turbulence model (Nagata et al., 2005).
As indicated above, the present model is reasonable to evaluate the bottom velocity fields around thepier, because the horseshoe vortex is mainly developed by the extension and rotation of thetransverse vorticity in approaching flow.

0.6
-2

0.4

-2

0

2

-6
-8

0

-14
-10
-8

-4-2

0.2
0

2

2
-2
-1
Water surface velocity distribution  28 Mar 2010 
-0.2
4
-6 -

2

0

-2

0

-0.4

0

-2

0

-0.6
-0.5

0
0.5
The measured bed topography

1

1.5

0.6
0.4

0

-2

0

-2
2

-4
-2

0

-4
-2

2

2
0

-6

-0.2

0

2

-12

-10
-8

-6

0.2

-4

0

-0.4
-0.6
-0.5

0

0.5

1

1.5

The computed bed topography
Figure 7

The comparison of between measured (Fukuoka et al., 1997) and computed
bed topographies around the pier after 30 minutes

Figure 7 shows the comparison of bed topographies around the pier between measured (Fukuoka etal., 1997) and computed results after 30minutes. And Figure 8 shows the comparisons of longitudinal

9and transverse bed profiles along the center-line of the pier. At the back of the pier, computed scour is
not developed but computed sediment deposition is developed compared with those of themeasurement. Roulund et al. (2005) also indicated that the scour depth at the back of the pier by thefull 3D turbulence model with the bed variation model by the bed load sediment transport was smaller
than that of the experiment. To improve this issue, a suspended load sediment transport model isrequired, as indicated by Roulund et al. (2005). However, in the upstream part of the pier, computed
scour hole is similar to that of the measurement. And the upstream and lateral bed profiles of scour
hole of the computation are in good agreement with those of the experiment. As indicated above, thepresent model is useful to simulate local scour in front of structures due to horseshoe vortex.

(cm)
3
0
-3
-6

-9
-12

Exp.

-15

Cal.

-18
-0.60

-0.40

-0.20

0.00

0.20

0.60 (m)0.40

The longitudinal bed profiles

(cm)
3
0
-3
-6

-9
-12

Exp.

-15

Cal.

-18
-0.6

Figure 8

5.

-0.4

-0.2
0
0.2
The transverse bed profiles

0.4

0.6 (m)The comparison of bed between measured (Fukuoka et al., 1997) andcomputed bed profiles

CONCLUSIONSThis paper indicates that the bottom velocity is evaluated by water surface velocity and depthintegrated horizontal vorticity for shallow water flows. Based on this consideration, this study proposesa bottom velocity computation (BVC) method, in which depth-integrated horizontal vorticity and water
surface velocity equations are solved in addition to shallow water equations and a depth averaged

10

turbulence energy transport equation. To develop a bed variation model for local scour aroundstructures, non-equilibrium bed load equation is derived based on the momentum equation of theparticle in the bed load sediment. It is demonstrated that the model can reproduce bottom velocity
fields around the pier and the scour hole induced by the horseshoe vortex.

6.

REFERENCES

Ashida, K. and Michiue, M.(1972). Study on hydraulic resistance and bed-load transport rate in alluvial
stream. Proceedings of the Japan Society of Civil Engineers, Vol.206，pp.59-69, in Japanese.
Blanckaert, K. and de Vriend, H. J.(2003). Nonlinear modeling of mean flow redistribution in curved
open channels. Water Resources Research, Vol.39, No.12, 1375, doi:10.1029/2003WR002068.
Engelund, F.(1974). Flow and bed topography in channel bends. Journal of hydraulics division, Proc.
of ASCE, Vol.100, HY11, pp.1631-1648.
Finnie, J., Donnell, B., Letter, J., and Bernard, R.S.(1999). Secondary flow correction for depthaveraged flow calculations. Journal of Engineering Mechanics, ASCE, Vol.125, No.7, pp.848-863.
Fukuoka, S.(2005). Flood hydraulics and channel design. Morikita, Japan, in Japanese.
Fukuoka, S. and Yamasaka, M.(1983). Alternating bars in a straight channel. Proceedings of theJapanese Conference on Hydraulics, JSCE, pp.703-708, in Japanese.
Fukuoka, S., Miyagawa, T. and Toboishi, M.(1997). Measurements of flow and bed geometry aroundA cylindrical pier and calculation of its fluid forces. Annual journal of hydraulic engineering, JSCE,
pp.729-734, in Japanese.
Ghamry, H. K. and Steffler, P.M.(2002). Two dimensional vertically averaged and moment equations
for rapidly varied flows. Journal of hydraulic research, Vol.40, No.5, pp.579-587.
Ikeda, S. and Nishimura, T. (1986). Three-dimensional flow and bed topography in sand-silt
meandering rivers. Proceedings of the Japan Society of Civil Engineers, Vol.369/II-5，pp.99-108, inJapanese.
Ishikawa, T., Suzuki, K. and Tanaka, M. (1986). Efficient numerical analysis of an open channel flow
with secondary circulations. Proceedings of the Japan Society of Civil Engineers, Vol.375/II-6 ，
pp.181-189, in Japanese.
Jin Y.-C. and Steffler, P.M.(1993). Predicting flow in curved open channels by depth-averaged method.
Journal of Hydraulic Engineering, ASCE, Vol.119, No.1, pp.109-124.
Melville, B. W. and Raudkivi, A. J. (1977) Flow characteristics in local scour at bridge piers, Journal of
Hydraulic Research, 15(4), 373-380.
Nagata, N., Hosoda, T., Nakato, T. and Muramoto, Y.(2005). Three-dimensional numerical model forflow and bed deformation around river structures. Journal of Hydraulic Engineering, ASCE, No.131,
No.12, pp.1074-1087.
Nadaoka, K. and Yagi, H.(1998). Shallow-water turbulence modelling and horizontal large-eddy
computation of river flow. Journal of Hydraulic Engineering, Vol.124, No.5, pp.493-500.
Nakamura, T., Tanaka, R., Yabe, T., and Takizawa, K. (2001). Exactly conservative semi-Lagrangian
scheme for multi-dimensional hyperbolic equations with directional splitting technique. Journal ofComputational Physics 174, pp. 171–207.
Nishimoto, N., Shimizu, Y. and Aoki, K. (1992). Numerical simulation of bed variation considering thecurvature of stream line in a meandering channel. Proceedings of the Japan Society of Civil
Engineers, Vol.456/II-21，pp.11-20, in Japanese.

11Olsen, N. R. B. and Melaaen, M. C. (1993). Three-dimensional numerical calculation of scour aroundcylinders. Journal of Hydraulic Engineering, ASCE, No.119, No.9, pp.1048-1054.
Olsen, N. R. B. and Kjellesvig, H. M. (1998). Three dimensional numerical flow modeling for estimationof maximum local scour depth. Journal of Hydraulic Research, IAHR, 36(4), pp.579-590.
Phillips, B.C. and Sutherland, A.J. (1989). Spatial lag effects in bed load sediment transport, Journal of
Hydraulic Research, Vol.27, No.1, pp.115-133.
Roulund, A., Sumer, B.M., Fredsøe, J. and Michelsen, J. (2005). Numerical and exparimantal
investigation of flow and scour around a circular pile. Journal of Fluid Mechanics, 534, pp.351-401.
Uchida, T. (2006). A CIP-based method for shallow water flows in complex geometries using
Cartesian grids. Proceedings of 7th International Conference on Hydroscience and Engineering,
ICHE-2006.
Uchida, T. and Fukuoka, S. (2009). A depth integrated model for 3D turbulence flows using shallow
water equations and horizontal vorticity equations. 33rd IAHR Congress: Water Engineering for a
Sustainable Environment, pp.1428-1435.
Yeh, K.-C. and Kennedy, J..F. (1993). Moment model of nonuniform channel-bend flow. I: fixed beds.
Journal of Hydraulic Engineering, ASCE, Vol.119, No.7, pp.776-795.
Wu, W. (2008). Computational river dynamics, Taylor & Francis, London.

12