7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Enhanced settling due to particleladen convective instability
Xiao Yu* TianJian Hsu and S. Balachandart
Department of Civil and Environmental Engineering, University of Delaware, Newark, DE 19716, USA
t Mechanical and Aerospace Engineering, University of Florida, Gainesville, FL 32611, USA
yuxiao@udel.edu, thsu@udel.edu and balals@ufl.edu
Keywords: convective instability, ,sc'lli particleladen flow, Equilibrium Eulerian approximation
Abstract
Convective instability may play an important role in determining both the location and timing of initial sediment
deposition off river plume. It has been observed in laboratory experiments for concentration as large as few gl 1, the
convective instability occurs and lead to enhanced apparent settling velocity on the order of few cm/s which could
not be explained by conventional Stokes settling law. In addition to the wellknown doublediffusive convection
mechanism, another mode of convective instability development is by heavier particles settling cross the density
interface,i.e., settlinginduced instability. The EulerianEulerian twophase equations is simplified with equilibrium
Eulerian approximation appropriate for fine sediment of small particle response time. The resulting governing
equations are similar to that of RayleightTaylor problem for densitydriven flow except for the consideration of
particle settling in the sediment conservation equation. A linear stability analysis for settlinginduced convection
is carried out to study the finger formation at the density interface. Relationships between growth rate and initial
sediment concentration, mix layer thickness, grain diameter and water depth are obtained to help us further understand
the dominant mechanisms causing the convective instability in a more realistic setting. Fully nonlinear analysis
is further carried out through numerical solution of the complete equations computed by a 3D pseudospertral
NavierStokes solver. Some discrepancies between results from the linear stability analysis and the fully nonlinear
numerical solutions are discussed.
Introduction
River water intruding into coastal ocean typically results
in a surface layer of sedimentladen water over relatively
dense saline seawater. To predict the location of ini
tial sediment deposition off river plume, it is critical to
estimate the equivalent settling velocity of sediment in
saltstratified sedimentladen river plume. Existing liter
atures in coastal oceanography mostly consider Stokes
settling law for primary particle or floc aggregate, which
usually gives settling velocity of no more than few mm/s
(e.g., Hill et al. 2000). However, limited but valuable
field evidences suggest when sediment concentration in
the river plume exceed 0(10) g/1 during episodic river
flooding, equivalent settling velocity exceed few cm/s,
which cannot be explained by conventional Stokes set
tling law (Warrick et al. 2008).
Convective instability across the density interface may
play an important role in such rapid sedimentation pro
cesses (e.g., Thacker and Lavelle 1978; David et al.
1999; Parsons et al. 2000; McCool and Parsons 2004).
However, the mechanisms causing such instability, such
as gravitational settling and doublediffusion, remains
unclear and hence more detailed studies are necessary.
Using fine sediment approximation (Ferry and Bal
achandar 2001) to the EulerianEulerian twophase
equations for fluidparticle system, we investigate con
vective instability and its implication to the occurrence
of rapid sedimentation of fine sediments off river plume
in the field. The idealized particleladen twolayer
flow problem is first solved by linear stability analy
sis in order to understand critical parameters affecting
the growth rate and characteristic size of the instabil
ity. Fully nonlinear 3D analysis is further investigated
by adopting a NavierStokes solver based on pseudo
spectral scheme previously demonstrated to be capable
of carrying out direct numerical simulations (DNS) for
turbulent flow (Mariano et al. 2008).
Conditions for settlingdriven fingering
Green (1987) presents a criterion to determine whether
doublediffusive convection or gravitational settling
dominates particle transport across the interface based
on the ratio (F) of the flux by double diffusion (Fdd) to
the flux by settling (Fs):
Fdd
F,
For F > 1, double diffusion is the dominant mecha
nism causing the instability, while for F < 1 gravi
tational settling dominates. Green (1987) modified the
doublediffusive flux presented by Schmitt (1079) to ac
commodate particles rather than a dissolved substance in
the upper layer. Assume initially there is no particles in
the lower layer and no dissolved substance in the upper
layer, the doublediffusive convective flux is
1 4
Fdd =pb(gi)3(C .)3
where b is a nondimensional constant determined by ex
periment (1/20), g is the acceleration of gravity, K is the
diffusion coefficient of the fastest diffusing substance
and C, is the sediment concentration in the upper layer
(expressed in units of mass/mass). For a dilute particle
suspension. The volumetric expansion coefficient is de
fined as
0 OP
(p18m) eo
\p09C,^
pS pf
P
where p8 is the density of particles.
The flux by settling could be simply calculated as the
product between settling velocity and concentration
Fs = WC
where C is the sediment mass concentration. For small
particles, the settling velocity can be estimated using
Stokes Law
W, 19(p pf)
18p y
where d is the particle diameter and v is the kinematic
viscosity of fluid.
Although in natural system, both doublediffusive and
settlingdriven fingering may be present at the same
time, we can always estimate the relative importance of
these two different mechanisms using parameter F de
fined above. For typical sediment properties and flow
conditions in natural river, i.e., d = 40m, C
2 1'.. /1, we have Fdd = 3.85 x 10 4kg/(m2 s) and
F, 3.80 x 10 3kg/(m2 s), which gives F z 0.1
and gravitational settling is the dominant mechanism for
instability. Hence, as a preliminary step, we focus on
studying instability induced by gravitational settling in
this paper.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Model formulation
It will be useful and instructive to study the main
characteristics of a complex flow with the help of a
somewhat simple flow configuration that possesses
the real physical and numerical features of the full
problem. For our problem, theoretically the sharp
planar interface might continuously descend until the
sand settles at the bottom and form a dense package.
However, the interface instability grows through time
and develops a finger like patten. Here, the system is
modeled as a mixture with variable density. In the initial
stage, the fluid is horizontally stratified, with denser
watersand mixture on top of clear fluid. Therefore, we
make contact to the classic RayleighTaylor instability
according to which the destabilization of the interface
between the layers is inevitable.
The flow is assumed to be consist of two layers with
heavier sedimentladen water on top of clear water. In
the initial state, the sediment concentration is given by
specifying the mixing thickness 6, which represents the
finite mixing layer thickness typically exists at the bot
tom of the river plume. The density interface is located
at y = 0 at t = 0 as shown in Fig.1. Cartesian coor
dinates (x, y) and velocities (u, v) are equivalently re
ferred to as the streamwise and normal component, re
spectively. The mixture density can be written as a func
tion of sediment concentration, p pf +c(p8 p ) with
c representing volumetric concentration of sediment,
and for typical sand, p8 given as 2.65 x 103kg m3.
Initially the sediment concentration profile is given as
co(z) 0.5Cmax(tanh(ay) + 1) with a defined as
a = 2,where the mixing thickness 6 defined as
I00 (
f"0 4
(c Cmax/2)2>
maCa
which is analogy to the definition of momentum thick
ness of the mixing layer.
Equilibrium Eulerian approximation is adopted to de
scribe the twophase fluidparticle system here. For
small particles, whose particle response time T is suf
ficiently smaller than the flow time scale, the particle
velocity can be expanded in terms of T and local car
rier flow field, see Jim Ferry & S. Balachandar (2000).
Here, only the first order term is kept, which gives
u" = u + Ws, suggesting that the particles are mov
ing with the carrier fluid flow except gravitational set
tling. The sediment concentration is considered to be
low, hence the particle interactions will be ignored. Be
sides, the density variations is considered to be small
such that Boussinesq approximation could be used. The
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
H
Figure 1: Sketch of the flow configuration
governing equation can be written as
Oui 0
=0
dOxi
bui bui 1 9p D2Ui
 +Uj +v
Dt dxj p dxi dxj dxj
dc dc D2c
+ (ui + WA,3) = K
at oxi Oxoxii
form are
Oui
=0
Oxi
Oui Oui
a+ u
Dp 1 D2ui
D + Re dxjdj
Oxi Re Oxj8xj
'R,. ..
W
s s
(s 1)gci3
where c is the concentration of sediment, ,u are fluid
velocities. The initial mixing thickness 6 is chosen as the
characteristic length scale for convenience, although it
should not be the correct scale for the convective fingers.
The velocity scale is chosen as U = ^/g ,where g'
g(s 1)Cma, so that
Xi
Xi*
8'
tU
t* 6
X^=^
ui p
ui* U' pP* U2
c WW
Cmax W U
The resulting governing equations in the dimensionless
Du' Dv'
Dx Dy
Du' Dp' 1 (D2
Dt Dx + ReSc 0e
Dv' Dp' 1 (D2
t D9y ReScc D
Dc' Dc' ,dco
W, + v
9t Dy dy
ui D2u'2
v' 2 v'U,
/+ D ) Ric'
R2 +y2
ReSe d2c' D 'y2 c
ReSc 8x2 Dy)
The incompressibility for 2D flow can be formulated in
terms of a streamfunction q'. Since the perturbations
satisfies V u' 0, it is permissible to introduce a per
turbation streamfunction 7' defined by
y '
The disturbance field 0'(x, y, t) now can be taken of the
following travelingwave form
(x, y, t) = (y)ei(kxwt)
dc dc 1 d2c
+ (ui + W, 6,3)
t D3xi ReSc OxiDxi
where the nondimensional groups are defined as
8U
Re 
(s 1)g6Ac
Ri =
U2
with H representing the half depth of the whole water
column and U is the characteristic velocity scale. For
simplicity, the superscripts have been neglected in the
above equations.
Linear stability analysis
For simplicity, a twodimensional flow system is stud
ies here by stability analysis. For field 0, we take the
decomposition in the form of =< 0 > +9', where
< 4 > is the mean of quantity 0, which could be veloci
ties, pressure or concentration, and 0' is the perturbation
from the mean state. We are interested in the stability of
the system, so introducing small perturbations into the
problem the resulting linear stability equations govern
ing the disturbance now becomes
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
where 0 denotes the complex amplitude, and k is the
real wavenumber of the perturbation in the streamwise
direction. Upon substituting the above ansatz into the
linearized governing equation and eliminating the pres
sure terms in the usual way, we will obtain an Orr
Sommerfeld type fourth order ordinary differential equa
tion.
dco dc 1 d2
iw + +W = k c+
dy dy ReSc y
(k  k k2Ric = 0
(k dy2 i Re( dy2 )
The initial concentration profile is given by co(y)
0.5(tanh(ay) + 1) in the dimensionless form.
The temporal longterm evolution of this type of distur
bance is characterized by w whose real part describes the
frequency w, and whose imaginary part the correspond
ing growth rate wi. The onedimensional eigenvalue
problem A0 w9 with matrix A obtained by Cheby
shev expansion with ChebyshevGuassLobatto colloca
tion points can be easily solved using existing eigenvalue
package.
Linear stability result
In the absence of sediment, the flow becomes clear
fluid at rest, which should be always stable. By setting
Ri 0 excluding the concentration equation and
solving the OrrSommerfeld equation, we obtain wi < 0
for all modes, suggesting that the initial perturbation is
always attenuated and the system is stable.
For Ri / 0, model solution gives wi > 0 for small
wavenumber k and the system is unstable. Fig.2 shows
the growth rate as a function of wavenumber k for
d 40pm,Cmax 1 x 104 and 6/H 1/30.
Evidently, the growth rate is not a monotonic function
of wavenumber k. The growth rate increases sharply
first then decreases smoothly towards large wavenum
ber. The peak corresponds to the most unstable mode,
which grows most rapidly and shall become dominate
mode soon after the onset of instability. Hence, the
most unstable mode can be used to estimate the size of
sediment finger.
Our principal interest is to study the growthrate as
a function of different flow parameters and sediment
properties, such as grain diameter d, mixing layer
thickness 6, initial concentration C,,,, and flow depth
H.
Effect of concentration
Considering the classic RayleightTaylor instability, it is
known that the key controlling parameter is the Atwood
dp=4p m
dp=40 I m
d= 100 p m
05
3 04
S03
5
05 1 15 2 25 3 35
Wavenumber k (cm1)
4 45 5
Figure 2: The growth ratewi as a function of wave num
ber k. The concentration Ac 1 x 10 4,
and mixing thickness 6/H 1/30. The most
unstable mode for different sediment diam
eter differs. It reveals that the finger width
will be a function of sediment diameter and
it decreases when the sediment diameter in
creases.
number defined as
SP P2 (s )Cmax
P1 + P2 2 + (s 1)Cma
with s the specific density of sediment, which is given
as 2.65 for natural sand particles. For small Atwood
number A < 1, the Boussinesq approximation is valid,
which is the case studied here. For increasing concentra
tion, the density differences also increase, which leads
to larger Atwood number and larger growth rate is ex
pected. For very small concentration, we can expand the
above equation using Taylor Series, which gives us
A z 0.5(s 1)Cmax[1 0.5(s 1)Cmx]
S0.5(s 1)Cma + O(C,,x)
The classic RayleighTaylor instability gives the growth
rate proportional to /A. Fig. 3 shows how the growth
rate varies with concentration for two choices of grain
sizes. By fitting the power law we acb using least
square method, we obtain b 0.60 for small particle
settling velocity of grain size dp 4/tm. If we set set
tling velocity to zero, which gives the miscible two fluid
system, the curve overlaps with the one for dp 4pm.
b 0.60 obtained here is different from the classic
RayleighTaylor result of b = 0.5 for twoimmiscible in
viscid fluid system. The presence of viscous effect and
sediment settling through the interface may also intro
duce perturbations to this system. Therefore, the growth
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Scm =001
Figure 3: The growth rate cw as a function of sedi
ment concentration Cax for dp = 4pm and
dp = 40pm. The mixing layer thickness is
6/H = 1/30. In loglog plot, the curve close
to a straight line, which implies a power law
relationship between the growth rate wi and
volumetric concentration C,,x.
rate is deviated from the value of 0.6. Fig.3 also shows
the growth rate as a function of concentration for lager
settling velocity of grain diameter dp 40/m. For
larger sediment settling velocity of dp = 40pm, b in
creases to 0.64.
Effect of initial mixing thickness
We can also change the initial mixing thickness to study
its effect on growth rate. In reality, the mixing layer
thickness at the bottom of a river plume is controlled by
freshsaline water densitydriven flow and the associated
ambient mixing processes by KH waves and Hombole
waves. These complicated processes shall be studied in
details however they are highly idealized here by a mix
ing layer thickness through initial condition.
Fig. 4 shows how the growth rate varies as a function
of 6/H. It can be observed from Fig. 4 that if the ini
tial sediment concentration is strongly mixed, i.e. larger
6/H value, the growth rate decreases. Since the instabil
ity occurs at the interface, which means we should con
sider a local Atwood number instead of using Ac in the
above equation for A. For larger 6/H value, the local
concentration at certain level c is smaller, which gives
a smaller Atwood number, therefore, smaller growth
rate. This may imply that in the field condition, stronger
frontal mixing at the bottom of the river plume may de
lay the development of sediment finger.
Effect of settling velocity and domain height
The settling velocity for fine particle can be given by
Stokes law. According to linear stability analysis the
Figure 4: The growth rate cw as a function of dimen
sionless mixing high 6/H. Sediment diam
eter dp 40pm and Cmax 1 x 10 4,
Cmax 1 x 10 2 respectively.
growth rate is not very sensitive to grain size for small
particle diameter. Fig. 5 shows that for dp < 15pm, the
growth rate does not change with sediment diameter. For
dp > 15pm, the growth rate decreases with increasing
sediment diameter. Our results suggesting the growth
rate is insensitive or decreasing with increasing settling
velocity appears to be inconsistent with limited obser
vation in the laboratory experiment (Hoyal et al. 1999).
This may be due to the limitation of a linear stability
analysis. More future work on fully nonlinear analysis
is necessary.
The domain height also plays minor role in determin
ing the growth rate. The boundary condition writes
u' v' 0 at both the top and bottom boundary, which
means the presence of boundary tends to attenuate the
perturbation and hence suppress instability. For larger
domain height, the growth rate only increases slightly
(comparing dashed curve to solid curve in Fig. 5). For
large enough domain height H, the boundary effect can
be neglected.
Fully nonlinear numerical simulations
Full numerical solution of the nondimensionalized gov
erning equations is computed by extending a 3D Navier
Stokes solver based on pseudospectral scheme (Mari
ano et al. 2008) in order to carry out fully nonlinear
analysis of convective instability and resulting sediment
finger development. An operator splitting method is im
plemented to solve the momentum equation along with
the incompressibility condition. A lowstorage mixed
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
H=O1m
  H=O 2m
Figure 6: Concentration contour at t = 153s
Figure 5: The growth rate cw as a function of sediment
diameter dp. For different domain height, the
growth rate only differs slightly.
thirdorder RungeKutta and CrankNicolson schemes
are used for the temporal discritization of the advection
and diffusion terms respectively. More details on the
implementation of this numerical scheme can be found
in Cortese and Balachandar (1995). Since the interface
descends with the settling of sediment phase, a vertical
coordinate transform is further adopted into the pseudo
spetral code. The coordinate transform is done by set
ting znew = Zold Wst in the vertical direction only.
Fourier expansions for the flow variables in the stream
wise and spanwise direction remain unchanged. In the
inhomogeneous vertical direction, a rational collocation
method with adaptively transformed Chebyshev points
Tee & Trefethen (2006) is also applied.
The computational domain is a box of size L, = 2 x Lz
and L, 0.5 x Lz (with L = 306), which extends from
x 2.0 to x 2.0, from y 0.5 to y 0.5 and
from z = 1.0 to z = 1.0. The flow is initialized from
rest with co = Cmax(tanh(az) + 1)/2 where heavier
watersediment mixture on top of nearly clear fluid (see
Fig. 1). The simulations are performed using a resolution
of N, 256, N = 128 and N, 129.
Periodic boundary conditions are enforced for all the
variables in both streamwise and spanwise directions.
This is due to the characteristics of the spectral method
used. However, the computational domain is taken to
be large enough in these directions to allow the full de
velopment of sediment fingers for a sufficient amount of
time. At the top and bottom walls, noslip and no pene
tration conditions are enforced for the continuous phase
Figure 7: Concentration contour at t = 170s
velocity. For the normalized concentration, we apply
and the mass of sediment should be conserved in the sys
tem.
Fig. 6 and Fig. 7 shows preliminary results of two
snapshots of concentration contour at two different time
t 153sec and t 170sec. When these fingers are
macroscopically observable, as shown in Fig. 6 and Fig.
7, the number of fingers is fewer than that predicted by
the linear stability analysis. In this case, 3D simulation
underpredicts the number of finger exactly by onehalf.
The discrepancy between the 2D linear stability analysis
and full 3D solutions could due to the action of pair
ing mechanism upon the 3D microscopic fingers caus
ing them to merge. Tan and Homsy (1987) explain this
as a secondary modulatory spanwise instability, which
causes the fingers to couple and coalesce. More detailed
numerical study shall be carry out in the future to in
vestigate the 3D fully nonlinear dynamics of sediment
fingers.
Conclusions
The linear stability analysis for the initial evolution of a
watersand interfaces is investigated. Equilibrium Eu
lerian method is used here to describe the twophase
particleladen flow appropriate for fine particles. The
linear stability analysis provides us the growth rate,
which can be used to provide critical condition for the
occurrence of sediment finger during initial deposition
of sediment off river plume.
The results of stability analysis suggest that the growth
rate is a strong function of sediment concentration and
mixing thickness. However, the growth rate is less sen
sitive to variation of grain size and flow depth (i.e., the
height of computational domain). The growth rate is a
power function of sediment concentration. For finer sed
iment, this power law gives ci acb, where b is in the
range of 0.5 0.7 depending on sediment diameter and
mixing thickness 6.
Examination of for the most unstable waves has allowed
for the identification of the characteristic size of the sed
iment fingers, such as finger width. Full numerical so
lutions of the 3D nondimensional governing equations
show pairing mechanism that cannot be captured by lin
ear stability analysis.
Acknowledgements
This study is supported by the U.S. Office of Naval
Research (N000140910134) and National Science
Foundation (OCE0913283) to University of Delaware
and by U.S. Office of Naval Research grant (N00014
0710494) to University of Florida. This work is
also partially supported by the National Center for
Supercomputing Applications under OCE70005N and
OCE80005P utilizing the NCSA Cobalt and PSC Pople.
References
A. Lange, M. Schriter, M.A. Scherer, A. Engel, and I.
Rehberg, Fingering instability in a watersand mixture.
The European Physical Journal B 4, 475484(1998).
Mariano I. Cantero, S. Balachandar and Marcelo H.
Garcia. An EulerianEulerian model for gravity currents
driven by inertial particles. International Journal of Mul
tiphase Flow 34(2008) 484501.
C. VOltz, W. Pesch, and I. Rehberg. RayleighTaylor in
stability in a sedimenting suspension. Physical review E,
Volume 65, 011404.
David C.J.D. Hoyal, Marcus I. Bursik and Joseph F.
Atkinson, Settling driven convection: A mechanism of
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
sedimentation from stratified fluids. Journal of Geophys
ical Research, vol. 104, No. C4, Pages 79537966, April
15, 1999.
Jim Ferry and S. Balachandar, A fast Eulerian method
for disperse twophase flow. International Journal of
Multiphase Flow 27(2001) 11991226.
C.T. Tan and G.M. Homsy, Stability of miscible dis
placements in porous media: Radial source flow. Physics
of Fluid 30(5) May 1987.
Green, T., The importance of double diffusion to the
settling of suspended material, Sedimentology, 34, 319
331,1987.
Donald A. Drew, Stability of a Stokes' layer of a dusty
gas. Physics of Fluids 22(11), November 1979.
W.C. Thacker and J.W. Lavelle, Stability of settling of
suspended sediments, Physics of fluids 21(2), February
1978.
W.B. Zimmerman and G.M. Homsy, Nonlinear viscous
fingering in miscible displacement with anisotropic dis
persion. Physics of fluids A, Vol. 3, No. 8, August 1991.
Tamay M. Ozgokmen and Oleg E. Esenkov, Asymmet
ric salt fingers induced by a nonlinear equation of state.
Physics of fluids, Vol. 10, No. 8, August 1998.
D. Parsons, Jeffrey and Garcia Macelo H. Enhanced sed
iment scavenging due to doublediffusive convection.
Journal of Sediment Research 70(1), 4752, 2000.
T.W. Tee, L.N. Trefethen, A rational spectral colloca
tion method with adaptively transformed chebyshev grid
points, SIAM J. Sci. Comput. 28 (2006) 17981811.
Canuto, C., Hussaini, M., Quarteroni, A., Zang, T.,
1988. Spectral Methods in Fluid Dynamics. Springer
Verlag.
Wayne W. McCool and Jeffrey D. Parsons, Sedi
ment from buoyant finegrained suspensions. Continen
tal Shelf Research 24(2004) 11291142.
JA Warrick, J Xu, MA Noble and HJ Lee, Rapid forma
tion of hyperpycnal sediment gravity currents offshore
of a semiarid California river. Continental Shelf Re
search, Volume 28, Issue 8, 9911009 2008.
Paul S. Hill, Timothy G. Milligan and W. Rockwell
Geyer. Controls on effective settling velocity of sus
pended sediment in the Eel River flood plume, Conti
nental Shelf Research, Volume 20, Issue 16, 20952111,
2000.
