Computation of Flow, Turbulence and Bed Evolution
with Sand Waves
Yasuyuki SHIMIZU1 and Sanjay GIRI2
1

Dr. Eng., Professor, Dept. of Hydraulic Research, Hokkaido University(North 13, West 8, Kita-Ku, Sapporo 060, JAPAN)
2
Postdoc Research Fellow, Dept. of Hydraulic Research, Hokkaido University(North 13, West 8, Kita-Ku, Sapporo 060, JAPAN)

Abstract
A vertical two-dimensional morphodynamic model with non-hydrostatic, free surface flow is reported
herein. This numerical model can simulate flow and turbulence characteristics over sand waves. Likewise,model can reproduce the sand wave formation and migration process as well as free surface oscillation
simultaneously in a moving boundary domain with the implication of a generalized coordinate system. A
nonlinear k-ε model is employed as a turbulence closure. CIP numerical technique is used to resolve
advection term of momentum equations. An Eulerian stochastic formulation of sediment exchange
process in terms of pick up and deposition function is incorporated to simulate non-equilibrium sedimenttransport that explicitly considers the flow variability during morphodynamic computation. The numericalmodel is validated with experimental data on flow and turbulence over fixed dunes as well as laboratory
measurement and visualization of dune geometry and migration.
Key Words: Numerical simulation, free surface flow, sediment transport, morphodynamics

1. Introduction
The microscale sand waves in alluvial streams are usually formed in lower flow regime under roughturbulent condition. Several investigations have been carried out so far that have made invaluable
contributions to improve understanding on morphodynamic features of dunes as well as flow and turbulence
characteristics induced by sand waves. The complex hydraulic and morphodynamic aspects associated withbed form evolution have an important bearing on problems; however still remain poorly understood despite
numerous efforts that have been made since long.
Understanding flow and turbulence induced by sand waves is of importance to explicate their formation
mechanism and further development. Number of physical studies in this regard was made in the past1). For
instance, Nelson et al.2) conducted flow and turbulence measurement over two-dimensional fixed dunes, inwhich along with the coupling of mean flow and turbulence, effects on bed form instability and finite
amplitude stability were examined. Schmeeckle et al.3) conducted laboratory visualization of turbulence andsuspended sediment transport over two-dimensional dune. Furthermore, some notable numerical studies werealso carried out. Shimizu et al.4) performed a three-dimensional direct numerical simulation of flow andturbulence over two-dimensional fixed dunes. The numerical model was able to reproduce the coherent
structure induced by dune crest adequately. Most observations yielded strong correlation between coherent
turbulence structure induced by flow separation, bed resistance and sediment transport mechanism. The flow
field and bed configuration are found to be interdependent having direct and converse effect on each otherthat creates difficulties in appropriate physical interpretation and quantitative determination of the said
problem.
Numerous approaches were developed in order to quantify sand wave induced morphology of alluvial
channels. Most early investigations were conducted using theoretical, semi-empirical and empiricalapproaches5), 6),. Among them, theoretical studies are largely based on linear stability theory which is, in
effect, not appropriate approach for the determination of natural phenomena having finite amplitude

features2). Scant attention has been given in regard to numerical modeling of sand waves’ morphodynamics7),
probably because of numerical complexity and limitations on its widespread application to real-world
problem. Nakayama and Shimizu8) carried out a numerical study on suspended sediment transport over sandwaves. They compared the computation result of time averaged concentration of suspended sediment with
Rouse distribution curve as well as laboratory measurements. One of the significant studies was made by
Onda and Hosoda9), in which a depth-integrated flow model was proposed with an allowance for verticalacceleration. They described bottom shear stress based on potential flow analysis considering
acceleration-deceleration effect near the bed. In order to compute water surface profile, a numerical
technique was used to simulate undular bore and flow over dunes with an attempt to eliminate the deficiency
of hydrostatic assumption by introducing a reduction factor to the vertical acceleration terms in case of highdepth gradient. However, this flow model cannot reproduce separation behind dune crest and, thereby, its
effect on flow field and sediment transport. Furthermore, they used sediment pick up and deposition
function10) to simulate sand dune formation. Finally, they analyzed temporal variation of sand wave spectrum
as well as dune geometry using their numerical results.
Present research has been prompted by the need to develop a more complete and reliable computationalmodel to simulate flow and morphodynamic features of sand waves. An earlier proposed three-dimensionalmodel4) seems to be more complete approach for computation of flow and turbulence over dunes. However,coupling this flow model with a morphodynamic module may inhibit it from being applied efficiently
because of rather high computational cost. Consequently, aforesaid three-dimensional hydrodynamic model
is simplified to be a vertical two-dimensional model and enhanced by imposing non-hydrostatic, free surface
flow condition as well as coupling with morphodynamic module. The model is verified with flow andturbulence experiment over fixed dunes as well as laboratory measurement of sand dune geometry andvisualization of dune migration process.

2. Computation of Flow and Turbulence(1) Governing equations
The governing equations for vertical two-dimensional flow in cartesian coordinate system reads as:
∂u ∂v
+
=0
∂x ∂y

(1)

(

)

(

)

(2)

∂
∂v
∂v
∂v
1 ∂p ∂+u +v=−+− u ' v' +− v' v' − gρ ∂y ∂x
∂t
∂x
∂y
∂y

(3)

∂
1 ∂p ∂
∂u
∂u
∂u
+u
+v
=−
+
− u' u' +− u ' v'ρ ∂x ∂x
∂t
∂x
∂y
∂y

(

)

(

)

where x and y = coordinates in horizontal and vertical direction respectively; u and v = components ofvelocity in horizontal and vertical direction respectively, − u 'u ' , − u 'v' and − v'v' = Reynolds stresstensors, ρ = fluid density, g = gravitational acceleration, p = pressure.
Eq. (1)-(3) are transformed from x, y, t cartesian coordinate system to a moving boundary fitted ξ, η and τ
coordinate system4).
Pressure term in momentum equations can be computed considering non-hydrostatic component asfollows:H∫

p = p 0 + p' = g ρ dy + p'(4)yEq.(4) can be rewritten as:

p = ρg ( H − y ) + p'

(5)

in which, p’= non-hydrostatic component of pressure, H = the location of free surface. Calculation is
performed substituting pressure term (p) by Eq.(5) in momentum equations, consequently, computing thenon-hydrostatic component.

(2) Calculation of water surface fluctuation
Computation of time-dependent water surface change is of importance for the realistic reproduction offree surface flow over migrating bed forms. Basically, majority of morphodynamic numerical models
developed earlier have assumed the rigid lid water surface condition to achieve numerically stable solution.
In present study, the kinematic condition is imposed along the free surface (at y = H) to compute thetemporal water surface elevation. The kinematic condition, which constraints fluid particles to remain on thewater surface at any time following the local flow velocity, is expressed as follows:

∂H
∂H
+u
=v
∂t
∂x
H = yb + h
where yb = bed elevation and h = local flow depth.

(6)
(7)

(3) Turbulence closure
In conventional k-ε model, turbulence stress tensors are evaluated using linear relationship. In order toreproduce turbulence characteristics more accurately in shear flow with separation zone, a nonlinear term isadded to the standard k-ε model as follows11):− u i' u 'j =ν t S ij −

3
2
1
k


kδ ij − ν t ∑ C β  S βij − S βαα δ ij ε β =1 33k2

ν t = Cµ

(8)
(9)

ε

where νt= eddy viscosity coefficient, δij= Kronecker delta, k= turbulent kinetic energy, ε= dissipation rate ofk. The detail description and formulation of strain tensors Sij, Sβij, Sβαα, as well as coefficient Cβ can be foundelsewhere11).
k-ε equations can be expressed as follows:∂k∂k∂k ∂  ν t ∂k  ∂  ν t ∂k 
(10)
+ 
 + Ph − ε
+u
+v
= 
∂t
∂x
∂y ∂x  σ k ∂x  ∂y  σ k ∂y εε2∂  ν t ∂ε  ∂  ν t ∂ε 
∂ε
∂ε
∂ε
(11)
 + C1ε Ph − C 2ε

 +

=
+v
+u
κ
k
∂y ∂x  σ ε ∂x  ∂y  σ ε ∂y ∂x∂t
Ph = − u i' u 'j

∂u i
∂x j

(12)

where σk, σε, C1ε and C2ε are standard model constants.
(4) Boundary conditions
The boundary conditions are no slip at the bed and free flow at the surface.
for near-bed region is adopted:

up
u*

=

1

κ

ln

yp
y0

The following expression(13)

where up = velocity at near-bed grid point, u* = bed shear velocity, κ = Karman constant, yp= distance frombed to nearest grid, y0= ks/30, ks = roughness height (= 2.5d), d = sediment diameter.
The periodic boundary condition is employed in this computation. Turbulent energy and its dissipation innear-bed region are evaluated as follows:
k
1
(14)
=
2
u*
Cµ

ε=

u*3
ky

(15)

In order to account for the diminution of turbulence length scale at the free surface layer, following
damping function to the eddy viscosity is used:( h − y )ε s 
f s = 1 − exp − B
3/ 2

k
s



(16)

where B=10.
The dissipation rate of turbulent energy at the surface layer is evaluated using following equation:εs =C µ3 /04 k s3 / 2
0.4∆y swhere ∆ys = distance from water surface to first grid and coefficient Cµ0= 0.09
∂k
=0
∂y

(17)

(18)

(5) Numerical algorithm
The transformed equations are numerically solved by splitting them into non-advection and pureadvection phase. The non-advection phase is computed using central difference method. The pressure field isresolved using SOR method. The advection phase is calculated using a high-order Godunov scheme known
as CIP method12). In this scheme, at very small time increment, the change in time of velocity components ata point in space can be split into the time evolution of the inhomogeneous terms and the time evolution at apoint due to the advection of the field. A brief description of CIP numerical scheme to solve advection terms
in generalized coordinate system can be found elsewhere4).

3. Sediment Transport Model
Some earlier investigations13) have revealed that spatially or/and temporally averaged sediment transportmodels are inadequate to describe the extremely inconsistent nature of sediment particle motion particularlyin the region with high near-bed turbulence. In effect, a direct computation of sediment particle motion, viz.
Discrete Element Model (DEM), which responds to the local and temporally variable flow field, is physicallymore complete approach. However, coupling this approach with an advanced hydrodynamic model, atpresent, seems to be rather sophisticated in view of computational complexity and expenses, in particular, forsolving comparatively big scale problem. Computationally efficient and relatively comprehensive approachis believed to be the modeling of pick up and deposition of sediment particles in terms of a stochastic
formulation. Consequently, an Eulerian stochastic formulation of sediment transport proposed by Nakagawa
and Tsujimoto10) is incorporated in our numerical code. This method yielded one of the best predictions ofpick up function for fine sediment grain. Moreover, this approach was effectively used earlier9) as well forbed form development computation.
The dimensionless pick up rate is expressed as follows:p s d ( ρ s ρ − 1) g = 0.03τ * (1 − 0.035 / τ * ) 3(19)
Where ps =sediment pick up rate, ρ and ρs = fluid and sediment density respectively and τ*= dimensionlesslocal bed shear stress.
The sediment deposition rate reads as:
p d = p s f s (s )(20)
where pd = sediment deposition rate and fs(s) = distribution function of step length.
Distribution function of step length is found to be exponential as follows:f s ( s) =1 sexp − 
Λ
 Λ

(21)

where Λ = the mean step length and s = the distance of sediment motion from pick up point.
The mean step length can be calculated as Λ = αd, in which α is an empirical constant and proposed to be100.
In order to consider the instability of bed forms in upper flow regime, a module to compute suspended
sediment transport is also incorporated in proposed numerical model. Considering suspended sedimenttransport, the bed deformation rate can be computed using following sediment continuity equation:∂y b
1  A3+
 ( p d − p s )d + (q sui − w f cb ) = 0∂t 1 − λ  A2
(22)

where yb = bed elevation,λ= porosity of sediment particle, A2 and A3 = shape coefficients of sand grain with
two- and three-dimensional geometrical properties, qsui = upward suspended sediment flux per unit area, wf =falling velocity and cb = reference concentration of suspended sediment.
In case of equilibrium suspended sediment transport or its absence, the third term of Eq. (22) vanishes.
qsui is calculated using equation proposed by Itakura-Kishi: ρ − ρ gd
Ω − w f 
q sui = K  α * s
ρ s u*


(23)

in which,
∞

Ω =

1

∫ ξ π exp(−ξ )dξ
τ
2

* a

B*

∞

1

+

B*η 0

∫ π exp(−ξ )dξ
2

τ*

−1

(24)

a

a=

B*

τ*

−

1

η0

(25)

where B*, η0, α* and K are constants, the values of which are adopted, namely B* = 0.143, η0 = 0.5, α* =0.14 and K = 0.008.
The falling velocity is calculated using well-known Rubey's formula:wf =2
36ν 2 6ν
Sgd +
−
2
3
d
d

(26)

where S = relative density of sediment in water, namely 1.65, ν = kinematic viscosity of water.

4. Dune Migration Experiment
An experimental investigation on sand dune formation and migration process was carried out under the
department of hydraulic research of Hokkaido University. Experiments were performed in laboratory flumeof CERI of Hokkaido. The total length of the flume was 30m, width was 10cm and height was 30cm having
1/500 slope, regulated by automatic control system. The middle part of flume (10m) was covered with 5cm
thick sediment layer with 0.28mm median diameter, which was trapped by using sediment stoppers at thepart of upstream and downstream having length of 10m each. The sand feeding was done simply by fillingup manually the space between movable bed and sand stopper in the upstream boundary, formed because oftransported sand. The migration process of sand dunes was visualized using high resolution video camera.
Maximum likelihood evaluation of sand dune geometry was made at the end of each experiment. The
experimental cases and condition are listed in Table 1. It is to be noted that dune migration feature(particularly, length of dune) seems to be irregular for run A-4 and A-5.

5. Model Verification(1) Mean flow and turbulence characteristics
In order to evaluate the performance of hydrodynamic model, an experiment2) (along with data obtainedby personal correspondence with the author of said work) on flow and turbulence over fixed dunes isreproduced by using proposed numerical model. A typical example of instantaneous flow and vorticity field,
simulated by numerical model, is depicted in Fig.1. From this, it can be seen that the numerical model is able
to reproduce the realistic feature of flow separation and vortices generation behind the dune crest.
The time-averaged streamwise and vertical velocity profiles in all regions over dunes are reproduced
quite well by numerical model (Fig.2). The turbulence characteristics, particularly in the separation region
behind the dune crest, are found to be under-predicted in some points (Fig.3). It seems that proposed modelis still unable to resolve turbulence in small scales with compare to three- dimensional model4). However, on
the whole, result can be regarded as satisfactory.

Table 1 Experimental cases and condition
Run

Discharge[cm3/sx103]

A-1
A-2

Time
[min]
250
105

2.0
4.0

Flow depth
[cm]
7.0
9.9

A-3

90

6.0

12.1

A-4

75

8.0

15.2

A-5

60

10.0

16.7

(2) Water surface fluctuation
Proposed numerical model is found to be able to reproduce realistically the water surface fluctuation due
to form drag induced by migrating dunes. The increase in depth after dune formation is found to be in sameorder as observed in our experiments. However, for run A-1 and A-2, average water depth seems to beunder-predicted (Table 2). The reason is not so clear. Nonetheless, it is thought to be associated with
under-prediction of dune geometry in numerical computation as well as some experimental discrepancies.
(3) Dune migration simulation
An attempt has been made to simulate sand wave formation and migration by adding an initial
perturbation on sand bed. The compared shape and size of dunes after analogous time increment forexperiment and computation has been listed in Table2. The dune geometry is found to be poorly predictedby numerical model, particularly for low flow rate, i.e. run A-1 and A-2. On the other hand, the dunecharacteristics are not so regular, especially for run A-4 and A-5, which is reproduced by numerical model inthe same manner as in experiments. Overall, computed results show acceptable prediction capability even forthe migration speed of dunes. Furthermore, simulated dune migration process is compared with theexperimental visualization for the qualitative assessment of model performance. Some qualitative
comparison between simulated and visualized dune geometry is depicted in Pic.1. The numerical simulationof dune migration reveals identical qualitative characteristic as in laboratory observation.

Dots indicate particle simulation of flow field

Fig.1 Numerical simulation of flow separation and vorticityu20(dot-experiment, line-numerical)

47.8cm/s

20

4.1cm/s15y (cm)y (cm)

15

v

10

5

10

5

0

0
70

80

90

100

110

120

130

140

150

70

80

90

100

x (cm)

110x (cm)

Fig.2 Mean flow comparison (Nelson’s experiment, 1993)

120

130

140

150

20

u'u'

213.5cm**2/sec**2

y (cm)

15

10

5

0

70

80

90

100

110

120

130

140

150

120

130

140

150

120

130

140

150

x (cm)
20v'v'

63.7cm**2/sec**2

y (cm)

15

10

5

0

70

80

90

100

110

x (cm)
20u'v'

50.6cm**2/sec**2

y (cm)

15

10

5

0

70

80

90

100

110

x (cm)

Fig.3 Turbulence characteristics comparison (dot-measured by Nelson et al., line-computed)

Table 2 Comparison of measured quantities with computed
Flow depthDune length

Dune height

[cm]

[cm]

[cm]

Migration
speed[cm/sec]

A-1

7.0/4.9

23/15

3.0/1.0

0.2/0.5

A-2

9.9/7.0

25/20

3.5/2.0

2.8/3.0

A-3

12.1/13.6

40/25

2.4/2.0

4.7/3.5

A-4

15.2/15.8

irregular

2.2/2.5

6.3/6.5

A-5

16.7/16.0

irregular

3.6/3.0

9.4/7.0

Run

Note: values in numerator and denominator denote measured and computed respectively.

Initial bed and water level (Run A-3)

Run A-2

Run A-3

Run A-5

Pic.1 Comparison between numerical simulation and experimental visualization of sand dunes

6. General Conclusion
Apparently, a significant attempt has been made here to embody a reliable and physically basedmodeling approach relating to the flow, turbulence and non-equilibrium sediment transport. The proposed
approach is believed to be more appropriate, for computing flow-field and non-equilibrium sedimenttransport mechanism with bed forms, than generally used morphologic modeling approach. Additionally,
module to compute suspended sediment flux could also be relevant to consider the bed form instability
during upper flow regime, though still remains to be verified. In general, computation results seem to berealistic and promising, which demonstrate scope and usability of numerical simulation technique that could
be a supplementary tool to comprehend a sophisticated river engineering phenomenon.

Acknowledgement: We gratefully acknowledge Kazutake Asahi who carried out experiments as a part ofhis graduate study. Likewise, we acknowledge Jonathan Nelson for providing experimental data.
References
1) Bennet, J.S. and Best, J.L.：Mean flow and turbulence structure of fixed, two-dimensional bed forms;
implication for sediment transport and bed form stability, Sedimentology, Vol.42, pp.491-513, 1995.
2) Nelson, J.M., McLean, S.R. and Wolfe, S.R.：Mean flow and turbulence over two-dimensional bed
forms, Water Resour. Res., Vol.29 (12), pp.3935-3953, 1993.
3) Schmeeckle, M.W., Shimizu, Y., Hoshi, K. and Tateya, K.：Turbulent structures and suspended sedimentover two-dimensional dunes, Proc. Int. Conf. Riv. Coast. Morph. Dyn., Genova, pp.261-270, 1999.

4) Shimizu, Y., Schmeeckle, M.W. and Nelson, J.M.：Direct numerical simulation of turbulence overtwo-dimensional dunes using CIP methods, J. Hydrosci. Hydr. Eng., Vol. 19 (2), pp.85-92, 2001.
5) Engelund, F.：Instability of erodible beds, J. Fluid Mech.,Vol.42, pp.225-244, 1970.
6) Engelund, F. and Fredsoe, J.：Sediment ripples and dunes, Ann. Rev. Fluid Mech., Vol.14, pp.13-37,
1982.
7) Fredsoe, J. and Tjerry, S.：Morphological computation of dunes, Proc. Int. Conf. Riv. Coast. Morph.
Dyn., Obihiro, Japan, pp.225-231, 2001.
8) Nakayama S., and Shimizu, Y.：Numerical calculation of suspended sediment over sand waves, Proc.
Int. Conf. Riv. Coast. Morph. Dyn., Obihiro, Japan, pp.217-224, 2001.
9) Onda, S. and Hosoda, T.：Numerical simulation on development process of micro scale sand waves andflow resistance, J. Hydrosc. Hydr. Eng, pp.13-26, 2005.
10) Nakagawa, H. and Tsujimoto, T.：Sand bed instability due to bed-load motion, J. Hyd.Div., ASCE, 106,
pp.2029-2051, 1980.
11) Kimura, I. and Hosoda, T.：A nonlinear k-e model with realizability for prediction of flows around bluff
bodies, Int. J. Num. Meth. Fluids, Vol.42, pp.813-837, 2003.
12) Yabe, T., Ishikawa, T., Kadota, Y. and Ikeda, F.：A multidimensional cubic-interpolated pseudoparticle(CIP) method without time splitting technique for hyperbolic equations, PSJ, Vol.59, pp.2301-2304,
1990.
13) Schmeeckle, M.W. and Nelson, J. M.：Direct numerical solution of bedload transport using local,
dynamic boundary condition, Sedimentology, Vol.50, pp.279-301, 2003.