7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Direct Numerical Simulation of Sediment Transport and Seabed
Morphology at Laboratory Scales
A. M. Penko,*t J. Calantoni* and D. N. Slinnt
Marine Geosciences Division, Naval Research Laboratory, Stennis Space Center, MS 39529, USA
t Department of Civil and Coastal Engineering, University of Florida, Gainesville, FL 32611, USA
allison.penko@nrlssc.navy.mil, joe.calantoni@nrlssc.navy.mil, and slinn@coastal.ufl.edu
Keywords: sediment transport, direct numerical simulation, ripples, bedform morphology, bottom
roughness, settling velocity
Abstract
Sediment transport and seabed morphology in the nearshore coastal zone are highly influenced by the settling rate of
sediment. Here we investigated the bulk settling rate of sediment in a mixture theory model. In our mixture theory
approach, the sediment flux is a function of advection and two empirical formulations describing the diffusion and the
hindered settling rate. We investigated the fidelity of the model to the specified hindered settling formulation through
a set of controlled simulations. We simulated the settling of a concentrated plane of sediment through a still water
column using the mixture theory model and compared the results to the empirical hindered settling formulation. As
expected, the model generally reproduced the hindered settling formula; however, the interpretation of the results
strongly depend on the definition of the local sediment concentration.
Introduction
Sedimentation rates directly affect suspended load trans
port and consequently have a strong influence on ero
sion, deposition, and bed morphology in the coastal re
gion. Net sediment transport rates are dependent on the
sediment settling velocity, especially fine grain sands,
when phase lags between sediment concentrations and
nearbed oscillatory velocities are common (Dohmen
Janssen et al. 2002). Therefore, a robust and accu
rate specification of the sedimentation rate in numerical
models of sediment transport is necessary to simulate
bedform profile changes and migration rates.
We implement a threedimensional bottom boundary
layer model (SedMix3D) using mixture theory to simu
late the coupled interaction between the fluid and sed
iment. SedMix3D treats the fluidsediment mixture as
a single continuum with effective properties that param
eterize the fluidsediment and sedimentsediment inter
actions using a variable mixture viscosity, a hindered
settling velocity, and a shearinduced, empirically cal
ibrated, mixture diffusion term. While mixture theory
has been well studied for many types of particle laden
flows (Nir and Acrivos 1990; Nott and Brady 1994;
Miskin et al. 1996a,b; Hofer and Perktold 1997; Sun
et al. 2009), and yields reasonable comparisons with
laboratory experiments, it has not yet been applied to
coastal sediment transport.
The model framework describing the sediment flux
includes the effects of advection, diffusion, and set
tling. SedMix3D parameterizes the hindered settling of
sediments in the sediment flux equation with a classic
bulk settling formulation (Richardson and Zaki 1954).
Richardson and Zaki (1954) performed laboratory ex
periments in a pipe to measure the bulk settling rate as
a function of local sediment concentration. We investi
gate the fidelity of the model to the specified hindered
settling formulation through a set of controlled simula
tions with idealized settling conditions that also include
affects from advection and diffusion. We compared the
bulk settling rates of planes of sediment settling through
a still water column simulated by SedMix3D to the em
pirical results of Richardson and Zaki (1954). In gen
eral, the model reproduces the hindered settling formula;
however, the degree to which the settling rate predicted
by the model agreed with the hindered settling formula
was strongly dependent on the definition of the local sed
iment concentration.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Methodology
The modeling framework of SedMix3D includes gov
erning equations for sediment flux, mixture continuity,
and mixture momentum. The mixture continuity equa
tion was derived by combining the fluid and sediment
phase continuity equations,
Op
+ V (pu) 0, (1)
at
where u is the mixture velocity and p is the mixture den
sity,
P = Ps + (1 )pf, (2)
where 0 is the sediment volumetric concentration, and
ps and pf are the sediment and fluid densities, respec
tively.
The mixture momentum equation was derived from
the sum of the individual phase momentum equations.
The mixture stresses can be defined by assuming the
mixture behaves as a Newtonian fluid (Bagnold 1954)
resulting in,
a9u
at PU Vu
VP + V (pVu) + F pg, (3)
where P is the mixture pressure, p is the effective vis
cosity, F is the external driving force vector per unit vol
ume, and g is gravitational acceleration (981 cm s 2 k).
SedMix3D employs a modified Eilers equation (Eilers
1941) to represent effective viscosity, /, here scaled by
the pure water viscosity, pf,
1 +I 05[] (4)
where [p] is the intrinsic viscosity, a dimensionless
parameter representing the sediment grain shape, and
0.0 < 0 < 0.63, where the lower bound represents pure
water and upper bound roughly corresponds to the max
imum concentration of unconsolidated sediment. Here,
we fix the maximum value of the effective viscosity by
specifying ~, 0.644. We chose the maximum pack
ing concentration to be 0 = 0.63, roughly equivalent
to random close packing for identical spheres. The in
trinsic viscosity parameter, [p], accounts for the particle
shape (Nawab and Mason 1958) and has been well doc
umented for spherical particles as 2.5 (Einstein 1906).
For irregularly shaped particles, [p] is still largely un
certain (Ferrini et al. 1979). Here we use an intrin
sic viscosity value of 3.0 in the modified Eilers equa
tion. Previous work has shown that an intrinsic viscosity
value within the range of 2.53.5 when employing Eil
ers' equation will not significantly affect the SedMix3D
0 0.1 0.2 0.3 0.4 0.5 0.6
Figure 1: Effective viscosity formulations.
output (Penko et al. 2009). Studies have found that vis
cosities modeled by concentration dependent equations
compare well with viscosity measurements of small
particles (0(0.001 cm)) in suspension (Mooney 1951;
Krieger and Dougherty 1959; Krieger 1972; Ferrini et al.
1979; Leighton and Acrivos 1987b; de Cindio et al.
1987; Sudduth 1993; Hunt et al. 2002) and large par
ticles (0(0.01 cm)) in dense concentrations (up to 0.70
volumetric concentration) (Vand 1948; Huang and Bonn
2007). Figure 1 compares common effective viscosity
equations.
The concentration of sediment is modeled with a sedi
ment flux equation (Nir and Acrivos 1990) that balances
the temporal gradients in sediment concentration with
advection, gravity, and shearinduced diffusion,
00 (5)
at+ u VO DV20 (5)
at az
where Wt is the concentration specific settling velocity
(Richardson and Zaki 1954),
W = Wto(l ), (6)
where Wto is the settling velocity of a single particle in
a clear fluid and q is an empirical constant,
4.35Re 0.03
4.45Rep 0.10
2.39
0.2 < Rep < 1,
1 < Rep < 500,
500 < Rep.
Re, is defined as the particle Reynolds number,
where d is the sediment grain size diameter. We consider
only noncohesive grains with the material properties of
* Eilers, 1941; SedMix3D, 2010
.* Mooney, 1951; Vand, 1948
Krieger & Dougherty, 1959
Huang & Bonn, 2007 data
Huang & Bonn, 2007 fit
x Hunt et al., 2002 data
 Hunt et al., 2002 fit
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
c = 0.567, d=0.2 cm
c = 0.567, d=0.2 cm
E
S35
S3C
30
o
o 25
N
20
15
10
1234 1234 1234
xDomain (cm)
Figure 2: Snapshots of concentration (0) contours during
gravitational settling of a sediment plane (averaged in the y
direction) initialized with a concentration of 0.567 and sedi
ment grain size of 0.2 cm.
quartz in water (p 2.66 g/cm3). The shearinduced
diffusion of sediment, D, is a function of grain size, vol
umetric concentration, and mixture stresses (Leighton
and Acrivos 1986). Assuming isotropic diffusion (i.e.,
D = Dy D,),
D = d2P3(,)Vul, (9)
4
where IVu (u : uT) and 3() is,
3() Ca2 (1 + e8) (10)
where a is an empirical constant found to be approx
imately 0.33 for large values of the Shields parameter
(0.5 < 0 < 30) by combining results from the dilute
limit with measurements in dense concentration suspen
sions.
The model uses a finite volume method with a stag
gered grid to solve the timedependent sediment concen
tration function and the mass and momentum conserva
tion equations for the fluidsediment mixture to second
order accuracy in space and thirdorder accuracy in time.
With grid spacing typically about an order of magnitude
less than the smallest (Kolmogorov) turbulent length
scale and time steps nearly four orders of magnitude less
than the smallest turbulent temporal scale, we consider
SedMix3D to be a direct numerical simulation.
OA
0.4
0.3 E 30
o
0
N
1234
time (s)
Figure 3: Contoured time series of the horizontally averaged
concentration for d = 0.2 cm. The solid black line denotes
the location of the centroid of the sediment plane. The black
dashed lines denote the boundaries of the concentration plane
(the point at which the concentration goes to zero).
Results
We examined the settling rate of a plane of highly con
centrated sediment in a still water column using three
dimensional simulations for two grain sizes (d = 0.2 cm
and d = 0.04 cm). The only external forcing in the sim
ulations was gravity. Figure 2 plots the concentration
(y) contours (averaged in the ydirection) for the d = 0.2
cm simulation. The full threedimensional domain was
5 cm x 5 cm x 60 cm. Initially the plane of sediment
was a block with a uniform concentration of 0.567, but
quickly diffused out as it settled through the water col
umn. Figure 3 plots the x and yaveraged concentra
tion though the water column as a function of time. The
dashed black lines denote the boundaries of the concen
tration plane (the point at which the concentration goes
to zero) and the solid black line denotes the location of
the centroid of the concentration plane. We calculated
the settling rate of the sediment plane by taking the time
derivative of the position of the centroid. The time se
ries of the concentration of the sediment plane (Figure
4a) was calculated using two different methods. First by
horizontally averaging the concentration at the location
of the centroid, and also by horizontally and vertically
averaging the concentration between the dashed lines in
Figure 3. Figure 4a plots the concentration calculated
using the two methods as a function of time. The blue
line represents the concentration horizontally averaged
at the location of the centroid of the plane (the solid
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
black line in Figure 3). The green line represents the
concentration horizontally and vertically averaged over
the entire plane of sediment (between the dashed lines in
Figure 3). Figure 4b compares of the settling rate of the
centroid of the plane to the empirical formula for hin
dered settling used in the model (Richardson and Zaki
1954). The blue line is (6) calculated with the concen
tration horizontally averaged at the centroid. The green
line is (6) calculated with the concentration horizontally
and vertically averaged over the entire plane of sedi
ment. The plot illustrates the dependence of the em
pirical formulation on the definition of local sediment
concentration.
Figures 5, 6, and 7 are analogous to Figures 2, 3, and
4, respectively, with d = 0.04 cm. The domain for the d
= 0.04 cm simulation was 6 cm x 6 cm x 12 cm. The
results also show a dependence on the definition of local
sediment concentration.
c = 0.567, d=0.04 cm
S 6
S.
1 2 3 4 5 1 2 3 4 5
xDomain (cm)
1 2 3 4 5
c = 0.567, d=0.2 cm
0.6
0.5
0.4
e
0.3
0.2
0.1
0
(a) SedMix3D: av
SedMix3D: ) pla
) 0.5 1 1.5
(b)
SedMix3D: Settl
RZ54: ) average
RZ54: 4 plane a
0 0.5 1 1.5
time (s)
Figure 4: (a) Time series of the horizon
tration at the location of the centroid of
Figure 5: Snapshots of concentration (<) contours during
gravitational settling of a sediment plane (averaged in the y
direction) initialized with a concentration of 0.567 and sedi
ment grain size of 0.04 cm.
Discussion
Dependence on Local Sediment Concentration
Our results demonstrate the inherent ambiguity of the
classic bulk settling formulation (Richardson and Zaki
1954) in that the prediction of the bulk settling rate is
strongly dependent on the definition of local sediment
2 2.5 3 concentration (Figures 4b and 7b). There was approxi
mately a 22% difference in the settling rates calculated
from the empirical formula (Richardson and Zaki 1954)
using the two different definitions of concentration
chosen in the simulation results. Calculating the local
sediment concentration using a plane average rather
than a centroid average decreased the value of the local
concentration, therefore, increasing the empirically
ing rate of centroid calculated settling rate (6). The settling rate of the
ed at centroid centroid from SedMix3D (black line in Figures 4b and
averaged 7b) is in better agreement with the quantitative value
of the hindered settling rate (6) calculated from the
2 2.5 3 plane averaged concentration (green line in Figures 4b
and 7b) for both simulations; however, the slope of the
settling rate from SedMix3D matches better with the
illy averaged concen
the plane (blue line) slope of the hindered settling rate calculated from the
theconentration at the centroid in both simulations.
concentration at the centroid in both simulations.
and me norizontally and vertically averaged concentration over
the entire plane (green line) for the d = 0.2 cm simulation.
(b) Time series of the settling rate of the centroid of the plane
(black line) plotted with the bulk settling rate (Richardson and
Zaki 1954) calculated using the two concentrations: horizon
tally averaging at the centroid (blue line) and horizontally and
vertically averaging over the entire plane (green line).
Dependence on Grain Size
The concentrations as functions of time in the
two simulations differ minimally, however, the time
averaged difference between the two settling rates is
approximately three times larger in the d = 0.2 cm
eraged at centroid
ne averaged
c = 0.567, d=0.04 cm
0.6
0.5
0.4
e
0.3
0.2
0.1
0
0.5
time (s)
Figure 6: Contoured time series of the horizontally averaged
concentration for d = 0.04 cm. The solid black line denotes
the location of the centroid of the sediment plane. The black
dashed lines denote the boundaries of the concentration plane
(the point at which the concentration goes to zero).
simulation than in the d = 0.04 simulation. This result
could be because the average settling rate is also three
times larger in the d = 0.2 cm simulation than in the d =
0.04 cm simulation. The sediment plane also diffuses
out slightly differently in the two simulations. The
height of the sediment plane in the d = 0.2 simulation
is initially about 15 grain diameters and diffuses to
approximately 100 grain diameters over 1.4 s cyann line
Figure 8). The d = 0.04 cm simulation is initialized with
a plane height of about 8 grain diameters and diffuses to
about 115 grain diameters over the same time (magenta
line Figure 8). The plane in both simulations diffuses
out at about the same rate initially, but then the d = 0.04
cm simulation begins to diffuse faster. The increased
diffusion in the d = 0.04 cm simulation was expected
due to the decreased size of the grains.
Conclusions
The hindered settling rates predicted by the mixture
theory model (SedMix3D) agree well with the speci
fied hindered settling formulation (Richardson and Zaki
1954) despite the inclusion of diffusion and advection
terms in the sediment flux equation. The Richardson
and Zaki (1954) empirical formulation is dependent only
on the local sediment concentration and the particle
Reynolds number. While the particle Reynolds num
ber is readily calculated, the definition of local sediment
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
c = 0.567, d=0.04 cm
(a) SedMix3D: 4 averaged at centroid
SedMix3D: 4 plane averaged
.
0.5 1
time (s)
Figure 7: (a) Time series of the horizontally averaged concen
tration at the location of the centroid of the plane (blue line)
and the horizontally and vertically averaged concentration over
the entire plane (green line) for the d = 0.04 cm simulation.
(b) Time series of the settling rate of the centroid of the plane
(black line) plotted with the bulk settling rate (Richardson and
Zaki 1954) calculated using the two concentrations: horizon
tally averaging at the centroid (blue line) and horizontally and
vertically averaging over the entire plane (green line).
concentration is somewhat ambiguous and strongly in
fluences the predicted bulk hindered settling rates.
60
40
20 d
d
0
0 0.2 0.4 0.6 0.8 1
time (s)
Figure 8: Height of the sediment plane (h) normalized with
the grain size (d) plotted with time for the two simulations.
Acknowledgements
AMP and JC are supported under base funding to the
Naval Research Laboratory from the Office of Naval Re
search (PE#61153N). Partial support for AMP provided
by the Office of Naval Research Coastal Geosciences
Program Code 322 CG. This work was supported in part
by a grant of computer time from the DoD High Perfor
mance Computing Modernization Program at the NAVY
DSRC and the ERDC DSRC.
References
R. A. Bagnold. Experiments on a gravityfree dispersion
of large solid spheres in a Newtonian fluid under shear.
Proceedings of the Royal Society of London, Series A,
225(1160):4963, 1954.
B. de Cindio, L. Nicodemo, and P. Masi. On the non
newtonian behavior of suspensions. Rheologica Acta,
26(1):100101, 1987.
C. Marjolein DohmenJanssen, David F. Kroekenstoel,
Wael N. Hassan, and Jan S. Ribberink. Phase lags in
oscillatory sheet flow: experiments and bed load mod
elling. Coastal Engineering, 46(1):61 87, 2002.
H. Eilers. The viscosity of the emulsion of highly vis
cous substances as function of concentration. Kolloid
Zeitschrift, 97(3):313321, 1941.
A. Einstein. Eine neue Bestimmung der Molekuildimen
sionen (German) [A new determination of molecular di
mensions]. Annalen der Physik, 19:289306, 1906.
F. Ferrini, D. Ercolani, B. D. Cindio, L. Nicodemo,
L. Nicolais, and S. Ranaudo. Shear viscosity of settling
suspensions. Rheologica Acta, 18(2):289296, 1979.
M. Hofer and K. Perktold. Computer simulation of con
centrated fluidparticle suspension flows in axisymmet
ric geometries. Biorheology, 34(45):261279, 1997.
N. Huang and D. Bonn. Viscosity of a dense suspension
in couette flow. Journal Of Fluid Mechanics, 590:497
507, 2007.
M. Hunt, R. Zenit, C. Campbell, and C. Brennen. Revis
iting the 1954 suspension experiments of R. A. Bagnold.
Journal of Fluid Mechanics, 452:124, 2002.
Irvin M. Krieger. Rheology of monodisperse latices. Ad
vances in Colloid and Interface Science, 3(2):111 136,
1972.
Irvin M. Krieger and Thomas J. Dougherty. A mech
anism for nonnewtonian flow in suspensions of rigid
spheres. Journal of Rheology, 3(1):137152, 1959.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
D. Leighton and A. Acrivos. Viscous resuspension.
Chemical Engineering Science, 41(6):13771384, 1986.
D. Leighton and A. Acrivos. The shearinduced migra
tion of particles in concentrated suspensions. Journal of
Fluid Mechanics, 181:415439, 1987b.
I. Miskin, L. Elliott, D. B. Ingham, and P. S. Hammond.
The viscous resuspension of particles in an inclined rect
angular fracture. International Journal Of Multiphase
Flow, 22(2):403415, 1996a.
I. Miskin, L. Elliott, D. B. Ingham, and P. S. Hammond.
Steady suspension flows into twodimensional horizon
tal and inclined channels. International Journal OfMul
tiphase Flow, 22(6):12231246, 1996b.
M. Mooney. The viscosity of a concentrated suspension
of spherical particles. Journal of Colloid Science, 6(2):
162170, 1951.
M. A. Nawab and S. G. Mason. The viscosity of dilute
suspensions of threadlike particles. Journal Of Physical
Chemistry, 62(10):12481253, 1958.
A. Nir and A. Acrivos. Sedimentation and sediment flow
on inclined surfaces. Journal of Fluid Mechanics, 212:
139153, 1990.
PR Nott and JF Brady. Pressuredriven flow of suspen
sions simulation and theory. Journal Of Fluid Mechan
ics, 275:157199, SEP 25 1994.
A. M. Penko, J. Calantoni, and D. N. Slinn. Mixture
theory model sensitivity to effective viscosity in simula
tions of sandy bedform dynamics. In Proceedings of the
OCEANS 2009 MTS/IEEE BILOXI Conference & Exhi
bition, Biloxi, MS, October 2009.
J. F. Richardson and W. N. Zaki. Sedimentation and
fluidisation: part 1. Transactions of the Institution of
Chemical Engineers, 32:3553, 1954.
R. D. Sudduth. A generalizedmodel to predict the vis
cosity of solutions with suspended particles .1. Journal
Of Applied Polymer Science, 48(1):2536, 1993.
D. W. Sun, S. R. Annapragada, and S. V. Garimella. Ex
perimental and numerical study of melting of particle
laden materials in a cylinder. International Journal Of
Heat And Mass Transfer, 52(1314):29662978, 2009.
V. Vand. Viscosity of solutions and suspensions .1. the
ory. Journal of Physical and Colloid Chemistry, 52(2):
277299, 1948.
