7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Modeling Interfacial Viscous Flows and Colloid Transport in a 2D Microchannel
X. Shi* V. Lazouskayat Y. Jint and L.P. Wang*
Department of Mechanical Engineering, University of Delaware, Newark, DE 19716, USA
t Department of Plant and Soil Science, University of Delaware, Newark, DE 19716, USA
graceshi@udel.edu and lwang@udel.edu
Keywords: multiphase flow, contact line dynamics, lattice Boltzmann equation, colloid transport
Abstract
To better understand colloid and contaminant transport and retention in the vadose zone where porescale multiphase
(i.e., airwater) flows through a soil porous medium play a role, we apply the mesoscopic LBM approach to simulate
an interfacial flow in a microchannel. A Lagrangian colloid tracking is then performed to study the trajectories of
colloids near the moving fluidfluid (airwater) interface. The effects of density ratio and viscosity ratio on the detailed
flow near the interface are examined. It was found that the flow on the water side is not sensitive to the density ratio.
The viscosity ratio is important in determining the stability of the interface. For the air front case when the less
viscous fluid invades the more viscous fluid, the fingering instability was realized. The normalized finger widths at
different Capillary numbers were found to be in good agreement with previous experimental and numerical studies.
For both the airfront and the waterfront cases, we found that colloids follow closely the flow streamlines. Even
though the airwater interface is unfavorable collector by design, the colloids can still be brought close to the interface
and may appear to be temporally retained due to the flow stagnation. Colloids may slide along the interface pri
marily due to hydrodynamic force. These observations agree with our parallel micromodel experimental observations.
Introduction
Colloidal particles are present in large quantities in wa
ter systems that enter the subsurface. Colloids, due
to their large available surface area, may help strongly
sorbing contaminants penetrate deeper into the under
ground environment. The role of mobile colloids in the
subsurface transport of contaminants has been recog
nized for some time (McCarthy & Zachara 1989), and
colloidfacilitated transport of contaminants and micro
organisms has been an important research topic in soil
and environmental sciences (Saiers & Ryan 2006; Sen
& Khilar 2006). The widespread use of nanoparti
cles in industrial and commercial products in recent
years also makes the topic of nanoparticle transport rel
evant as the environmental impact of nanoparticles has
been recognized (Fortner et al. 2005; Wiesner et al.
2006). Nanoparticles and natural colloids share simi
lar characteristics in the subsurface environment. There
fore, understanding the fate and transport of colloids and
nanoparticles is of great importance to pollution con
trol and natural resource preservation. A quantitative
description of these subsurface transport phenomena is
required before we can effectively monitor and manage
our environmental quality.
Compared to the saturated system, where the soil ma
trix is filled with water, the unsaturated system is com
plicated by the presence of airwater interface and mov
ing contact line. Wan & Wilson (1993, 1994) performed
the pioneering mechanistic study on the interaction of
colloids with airwater interface. They found that col
loids could irreversibly attach to stationary and moving
air bubbles. However, studies conducted by Chen &
Flury (2005) found that isolated bubbles and airwater
interfaces were unfavorable site to retention of colloids.
The similar conclusion was reached in Crist et al (2004).
Much remains to be understood regarding the mechanis
tic and quantitative description of colloid transport be
havior in the vicinity of airwater interface.
Several colloid retention mechanisms in unsaturated
porous media have been previously reported in the liter
ature: retention at solidwater interface (SWI), retention
at airwater interface (AWI) and retention at the contact
line. For example, Zhuang et al (2007) found that
the lower water saturation results in less colloid con
centration in the bulk, and the retained colloids could
regain mobility when the system is flushed with water.
The reason for a higher retention rate in an unsatura
tion system is likely the presence of airwater interface
and contact line. There have also been studies to isolate
interfacial forces from hydrodynamic effects so their in
dividual roles can be better investigated. For example,
to separate the effect of solid grain surface from that of
airwater interface, isolated bubbles are used by Chen
& Flury (2005) to avoid further complications of con
tact line where additional mechanisms such as the capil
lary force resulting from airwater interface could have a
component tangential to the solid wall (Gao et al 2008).
This current computational study focuses on colloids
transport near AWI and contact lines. Our study is moti
vated by the recent experimental studies of Lazouskaya
et al (2006) and Lazouskaya (2008). They utilized mi
cro channels containing fixed or moving airwater inter
faces to model different aspects of microscale flows rel
ative to AWI in unsaturated porous media. They found
that colloid retention at AWI was assisted by the stag
nant flow conditions. They also found that the hydrody
namic effects dominate transport of colloids away from
interfaces. Near SWI and AWI interfaces, the hydrody
namic effect is coupled to physicochemical processes.
At AWI or near a contact line region, the local flows of
ten contain stagnation points. This increases the interac
tion time between a colloid and the interface, making it
more likely for other interfacerelated interaction forces
to alter the motion of colloids.
There are two components to our computational
model. The first is the simulation of microscale flow
near an airwater interface and the second is the mod
eling of colloid transport. We consider airwater inter
facial flow in a 2D model microfluidic channel. Both
the capillary number and Reynolds number are small,
namely, the fluid inertia may be neglected and the flow is
mainly governed by viscous force and surface tension ef
fects. The mesoscopic lattice Boltzmann method (LBM)
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
is used to simulate the flow.
As an alternative CFD approach for multiphase
flows, Lattice Boltzmann method is still undergoing
rapid development. For problems involving large den
sity/viscosity ratios, LBM suffers from challenges such
as numerical instability and spurious currents. Several
attempts have recently been made to identify the origin
of spurious currents in order to reduce such currents.
Lee & Fischer (2006) showed that the use of potential
form ofintermolecular force can greatly reduce the mag
nitude of spurious currents in freeenergy based LBM
twophase model (He & Doolen 2002). Shan (2006)
argued that the anisotropic form of pressure in Shan
Chen model causes undesirable spurious velocity. He
proposed to use more nonlocal lattice nodes to bet
ter restore isotropy and thus suppress spurious currents.
However, the downside of this approach is that the in
terface becomes more diffusive when two or more belts
of neighboring lattice points are employed to represent
nonlocal interactions. Particular caution is required to
make sure that the thickness of interface is much less
than the characteristic length of flow.
In this work, we employ the LBM multiphase model
of Shan & Chen (1993) to simulate flow field near the
interfacial region, but will incorporate a different equa
tion of state (Yuan & Schaefer 2006) to enhance the
model capability in dealing with large density contrast
while maintaining a sharp interface. When the flow
field is established, we then simulate the motion of col
loids in the interfacial region and compare the simu
lated colloid trajectories with our parallel experimental
observations reported in Lazouskaya (2008); Shi et al
(2010). To study colloid transport behavior and reten
tion mechanism, we track colloids by numerically inte
grating the colloid's equation of motion with physico
chemical, hydrodynamic, Brownian and body forces. It
will be demonstrated that the simulated trajectories are
in reasonable agreement with our experimental observa
tions.
Multicomponent twophase flow
LB model
In the porescale experimental study of Lazouskaya
(2008), a trapezoidal channel with characteristic width
of 40 pm was used. The typical mean flow speed was
around 105 m/s. The resulting Capillary number Ca is
then of the order of 107, namely, surface tension dom
inates the viscous effects. It is challenging to match ex
actly such a small Capillary number in our simulation
as a very small velocity scale would be needed and as
such spurious currents can easily contaminate the phys
ical flow. The Capillary number in our simulation is
made to be much smaller than one. We assume the flow
pattern near the moving airwater interface is not sensi
tive to Ca as long as Ca is reasonably small.
Due to its physical simplicity and relative ease of im
plementation, the multiphase LBM model of Shan &
Chen (1993) has been used to study a variety of mul
tiphase flow problems such as droplet displacement in
a channel (Kang et al 2002), fingering flow of two im
miscible fluids (Kang et al 2004), and coalesce and
breakup (Shan & Chen 1993, 1994). The ShanChen
model is based on the concept of nonlocal fluidfluid
force interaction and can treat fluidfluid and fluidsolid
interaction using a uniform formulation of interaction
potentials.
Consider a complex fluid containing multiple compo
nents. The distribution function for the kth component
is governed by
ff (x, t) f() (x, t)
where Tk is the relaxation time of the kth component,
fk(eq) (x,t) is the corresponding equilibrium distribu
tion function. We use the standard D2Q9 lattice so i =
0, 1, 2,..., 8 and the particle velocities are: eo = (0, 0),
ei = (1,0),e2 = (0,1), e3 = (1,0), e4 = (0, 1),
e5 = (1,1),e6 = (1,1), e7 = (1, 1), and e8 =
(1, 1). The density and velocity of the kth fluid com
ponent are defined as
Pk =EfI' A = fiei. (2)
i i
The equilibrium velocity ueq is modified to take into ac
count the presence of multicomponent fluids and their
interactions, as
eq Ek t7k TFk, (3)
Uk k Pk/Tk + Fk, (3)P
where Fk = Flk + F2k represents the net force due to
mesoscopic fluidfluid interaction Flk and mesoscopic
fluidsolid interaction F2k,.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Shan and Chen (1993, 1994) represented the fluid
fluid interaction by the following nonlocal interaction
force
Flk = Ek(x) Gkk(x, X ',, (x')(x'x), (4)
x' k
where the interaction coefficients Gkk(x, x') will be dis
cussed below.
The scalar function .'. in Flk is known as the ef
fective mass density and is given as a function of par
ticle local number density, which results in a nonideal
equation of state that permits coexistence of two phases
for each component if the effective temperature is be
low a critical temperature. The exact form of .',. deter
mines the density ratio of two phases. More will be dis
cussed when specific forms of .',. are introduced. In the
D2Q9 lattice structure, each lattice node is surrounded
by 4 nearest neighbor nodes with an interaction distance
of 1 lattice unit and the next 4 nearest neighbor nodes
with an interaction distance of /2 lattice units. The
weights GkT must be carefully designed to satisfy sym
metry and isotropy conditions associated with gradient
operators (Shan 2006). When only the 4 nearest neigh
bors and the 4 next nearest neighbors are used, the opti
mal setting is
x x' =
x X' =
otherwise;
where the sign of gkk determines the nature (attractive or
repulsive) of the interaction between component k and
component k.
Similarly, the fluidsolid interaction is modeled as
F2k = Pk(x) Y kwP,(x')(x' x), (6)
x'
where the summation is over lattice nodes inside the wall
only. The relative magnitude and sign of gkw determines
the wetting property of a solid wall. We follow the im
plementation of Kang et al (2002) and consider only the
4 nearest neighbor lattice nodes in the above summation
for solid wall interaction. If we assume that the channel
wall is located half link away from the fluid lattice node
closest to the wall and the wall is parallel to xaxis,
then only one wall link is considered in the above sum
mation. The fluidwall interactions amounts to an ad
ditional wallnormal force applied to those fluid lattice
fik (Xtei6t, tt6t) fik (X, t)
Gkk (XI X.') nken i
nodes closest to the wall before the local collision op
eration. The usual bounceback conditions still apply to
those same lattice nodes to yield the approximate noslip
boundary condition.
It is important to note that, in the Shan and Chen
scheme, the effect of interaction forces is integrated into
the particle distribution function to make the scheme
formally explicit, while in reality the interaction forces
are added using the implicit trapezoidal time integration.
Namely, the distribution functions being solved is re
ally the sum of the original particle distribution func
tions after a net momentum shift resulting from the net
interaction force. As a consequence, the macroscopic
fluid velocity should be computed as pu = k I p, I, +
Ei Fk, where the factor reflects the merging of the
trapezoidal time integration of the force into the distri
bution function.
Finally, the specific detail on the implementation of
boundary conditions at the inlet and outlet is given in Shi
et al (2010).
Singlecomponent twophase
flow LBM simulations
Besides the difficulties of matching exactly the Reynolds
and Capillary numbers in our simulations to the exper
iment of Lazouskaya (2008), another challenge is the
high density ratio (around 1000) and high viscosity ratio
(around 50) for the airwater flow system due to the nu
merical instability and spurious currents. Our approach
here is to use a different equation of state in LBM, as
suggested by Yuan & Schaefer (2006), to handle high
density and viscosity ratios. Still we are currently un
able to treat the true density and viscosity ratio in the
airwater system. Instead, we shall investigate the sen
sitivity of the resulting flow on the assumed density and
viscosity ratios in order to understand if lower density
and viscosity ratios can be used to effectively simulate
the microscale flows in the airwater system for our pur
pose of studying the transport of colloidal particles.
For simplicity, here we shall consider the application
of ShanChen model at the singlecomponent twophase
flow setting; this is a special case (k = 1) of the general
multiphase LB model discussed in the last section. In
this case, the equation of state can be derived in terms of
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
the effective mass density function as
1 3
P = 3P+ 2ip ))2. (7)
While in their original model, Shan and Chen consid
ered only a specific form of t(p), it was later recog
nized by Yuan & Schaefer (2006) that other equations
of state (EOS) can be incorporated to extend the ca
pability of the model. Yuan & Schaefer (2006) con
sidered PengRobinson EOS, modified RedlichKwong
EOS, and a noncubic form of CarnahanStarling EOS
in the extended ShanChen model and showed that these
EOS can greatly enhance the algorithm stability for large
density ratios. By choosing a proper reduced tempera
ture on the coexistence curves, different density ratios
can be modeled. For example, by assuming the follow
ing form of t(p),
) 2pRT 2a2p 2p
3g 3g(1 bp) 3 Tgn (1 + bp) 9911'
(8)
we can model the RedlichKwong equation of state
given as
pRT ap2
S1 bp i(1 + bp)'
where a = 0.4274'' ', b = 0.08664RTc/pc. By
taking the first and second derivatives of pressure with
respect to the density, we can obtain the critical den
sity, pressure and temperature: pc = 2.7292, pc =
0.1784, Tc = 0.1961. The unique feature of this ap
proach is that the reduced temperature is explicitly built
into the model.
We adopted the above RedlichKwong EOS to sim
ulate a twophase interfacial flow in the microchannel
with a density ratio of the order of 100. The walls of the
2D microchannel is parallel to the xaxis, with y being
the cross channel direction. A parabolic velocity pro
file with a maximum velocity of 0.025 (in lattice units)
is imposed at both the inlet and outlet. For the single
component twophase flow setting (i.e., Case 2), the less
dense side (representing the gas side) is advancing while
the more dense side (the liquid side) is receding. For
comparison, a twophase flow with the same density on
two sides of the interface (i.e., Case 1) is also simulated.
The parameters used in the simulations are listed in
Table 1. The same maximum velocity is assumed for
two cases because we intend to compare the detail of the
simulated flow fields at the interfacial region.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 1: Setup of parameters for interfacial flow simulation of different density.
Symbol
Channel width
Channel length
Fluid densities
Viscosity
Maximum inlet velocity
Average front velocity
Reynolds number
Channel length
Initial transition steps
Total time steps
Surface tension
Capillary number
Fluidfluid interaction strength
Fluidsolid interaction strength
Contact angle
H
L
pl /p2
Vo
Vf
Re = HVf/v
L/H
Nsteps
a
Ca = Vfup/a
9
g
glw
0
Single
component(Case
1)
64
800
1.0/1.0
1 /1
0.025
0.0167
6.41
12.5
1000
20,000
0.0794
0.035
912 = 0.2
0.06/0.06
700
We make sure that the steady flow is achieved when
both the interface shape and the velocity field relative
to the interface are independent of time. The interface
locations for Case 1 and Case 2 are found to be only 5
grid spacings apart after 20,000 lattice time steps, cor
responding to a dimensionless time of tVo/H = 5.21,
where Vo is the mean flow speed and H is the channel
width. We found that the interface shape in Case 2 over
laps with that of Case 1, so we simply shift the interface
in order to compare the velocity distributions of the two
cases.
In Figure 1, we compare the distribution of the
streamwise (xcomponent) velocity along the x direc
tion, at two distances from the wall. At both distances
from the walls, the velocity far away from the interface
are identical. In the interfacial region, the velocity dis
tributions on the gas side are different, due to different
local flow patterns for different density ratios near the
interface on the gas side. Interestingly, however, the ve
locity distributions on the liquid side are very similar.
Therefore, the density difference alters the flow on the
low density side but has a negligible effect on the flow
of the high density side. For our application, the flow
on the liquid side is responsible for the transport of col
loids, therefore, a lower than actual density ratio used in
the simulation could be allowed as long as the flow on
the liquid side is of primary concern.
The ycomponent velocity dominates the relative
transport of colloids towards the wall region when col
loids are close to the interface. It determines whether
colloids would slide along the interface or simply slow
down at the stagnation region. In Figure 2, we compare
the distributions of the ycomponent velocity at a dis
tance 1/8 channel width from the channel wall. Again
it appears that the velocity field on the liquid side is not
sensitive to the density ratio, while flow on the gas side
can be different from different density ratios. This same
observation holds for the water front case when the liq
uid is advancing and the gas is receding.
LBM simulations with different
viscosity ratios
In order to further verify our 2D twophase flow simu
lation in a microchannel, in this section we will con
sider finger instability and compare our simulation re
sults with results from previous experimental and nu
merical studies.
It is well known that the interface tends to be un
stable when a less viscous fluid is moving into a more
viscous fluid, and a fingering instability develops. The
RK(Case 2)
64
800
7.53/0.079
1/1
66
0.025
0.0167
6.41
12.5
1000
20,000
0.5079
0.041
gii=0.2
0.06/0.06
700
1.0
I , I , I I
' I I I I '
0.0 0.2 0.4 0.6 0.8 1.0
x/L
rnterfa.
u@18WCsl.
u@1/8W Case.2
Gas Liquid
Locate at x=0.54
0.00 0.20 0.40 0.60 0.80 1.00
Figure 1: Comparison of xcomponent velocity near
the interface: (a) at the center of channel, (b) at a dis
tance of one eighth channel width. Case 1 has a density
ratio of 1 and Case 2 a density ratio of 94. The location
of interface is marked by a vertical line.
most studied topic of this instability is the the relation
ship between the normalized width Wf /W of finger and
the Capillary number, where Wf is the finger width and
W is typically the spanwise with of the HeleShaw cell.
Reinelt & Saffman (1985) simulated the finger penetra
tion of in a 2D channel by solving the Stokes equation
Interface
u@center Case.l
u@center ofCas2
S= 1 0.417 (1 e1.69Cao 25)
(10)
Experimental studies on finger instability went back
to the classical work of Saffman & Taylor (1958) who
showed, both theoretically and experimentally, that the
interface is unstable when a less viscous fluid is invad
ing a more viscous fluid. This was followed by several
other studies including the work by Tabeling et al (1987)
who, based on their experimental data, proposed to use
one single dimensionless number = 12Ca(b)2 to
better fit all data for different ratios of spanwise chan
nel width W of their HeleShaw cell to the distance b
between the two plates. The interface between air and
oil (viscous fluid used in the experiments) is a meniscus
which lives in a threedimensional space. When the in
terface is advancing, the meniscus could leave a thin film
behind it at the cross section of two closely spaced chan
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
2.0 I I
Interface
v@lW Case.l
1.6 v@ Case.2
1.2
Gas Liquid
0.8
0.4
Locate at x=0.54
0.0
0.4
0.00 0.20 0.40 0.60 0.80 1.00
x/L
Figure 2: The distributions of the ycomponent velocity
at the 1/8 channel width away from the channel wall.
with a normalstress boundary condition at the inter
face through a least squares iteration method. Halpem
& Gaver (1994) used a different method known as the
boundary element analysis to study this problem. The
governing equation can be simplified to Stokes equation
when Ca Re << 1. Halpem & Gaver (1994) obtained
a relationship between 3 = Wf/W and the Capillary
number Ca for Ca = 0.35 1.25. They further pro
posed the following regression formula
Liquid
f
Gas
Locate at x=0.543
nel walls. The finger interface has a concave shape at the
center where the maximum film thickness is measured.
Tabeling et al (1987) fitted their experimental results of
film thickness tmax by a stretched exponential as
tmax = 0.119b (l1 e0038(W/b) (l 8.58 ) ,
(11)
which can be transformed to the ratio of finger width to
the width b as
b tmax
1 0.119 (1 e0.038(W/b))
x(i
e8.58C ) (12)
Their prediction is good for Ca < 0.2, which is the Ca
range simulated on our study.
It must be stressed that, when comparing experimen
tal data with our 2D simulations, we must take into ac
count of the effect of the aspect ratio W/b on the liquid
film thickness left by advancing interface. Our 2D simu
lation could be viewed as the special case of W/b + oo,
the limiting case of very large aspect ratio of the Hele
Shaw cell studied in the experiments. The 2D simulation
data of Reinelt & Saffman (1985) would fall on experi
mentally predicted curve in the limit of W/b + oo.
Here we simulate two general situations of twophase
displacement in a 2D channel. The first is the airfront
case when a less viscous liquid is advancing into a more
viscous fluid, while the second is the water front case
when the flow direction is reversed. In the simulations,
the basic parameter settings for the these are the same,
except that the flow direction is reversed. We consider
a slow flow with the Reynolds number changing from
1.06 to 3.47 on the dense liquid (water) side. Two vis
cosity ratios are considered and their values are shown
in Table 2, along with the flow speed used to define the
velocity profile at the inlet and outlet. Case 3 is for the
viscosity ratio equal to 8.6, and Case 4 for 13.3. Other
parameters are the same as Case 1 shown in Table 1. The
resulting ratio of the finger width to the channel width is
also listed in Table 2.
The typical time evolution of the fingering interface
for the air front case is shown in Figure 3, while the sta
ble interface for the water front case is shown in Figure
4. The detail velocity fields and streamlines for these
cases are shown in Figures 5 and 6, respectively. For the
airfront case, the mean velocity is calculated by defini
tion: the average value of the stream velocity cross the
60.
40.0
0.0
0.0
0.0
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
 I IiI
200.0 400.0 I600.0 80.
CONTOUR FROM 0 TO 0 BY 0
Figure 3: The time evolution of the interface shape for
Case 3 at the airfront setting, the steady interface shape
develops into finger shape due to fingering instability.
Umax = 0.025 and the mean velocity of interface is
0.021. The flow goes from left to right.
CONTOUR FROM TO.2 BY.1
Figure 4: Interface evolution for the waterfront set
ting in Case 3, the steady interface is a stable interface.
Same parabolic inlet profile with Umax = 0.025, the
mean velocity is 0.017 which is consistent with the given
parabolic profile. The flow is moving from right to left.
channel. It is also used as a monitoring parameter to
determine whether the steady state is reached. Take the
case of Umax = 0.015 in Case 4 as example, the mean
velocity falls into 0.2%0.8% percent from time steps
10,000 and 20,000 (total time steps used as Case 3 and
Case 4). Therefore, we believe that the steadystate fin
ger shape and propagation of interface are established.
Capillary number is calculated with the mean velocity
and viscosity from the liquid side.
From the vector plots (Figure 5 and Figure 6), we
find that the interface shape does have an effect on flow
field near the interface and contact line and stagnation
regions. A pair of eddies appear on the air side at the tip
of the airfront finger interface (a little bit away from the
interface). This phenomenon has been visualized in the

=
Figure 5: The resulting flow field relative to the interface
and streamlines in the airfront case shown in Figure 3.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Table 2: Setup of parameters of twophase flow simulations of different viscosity.
Viscosity
trighht=0.43,ie/tl =0.05
Viscosity
/tright=0.66,/tIet =0.05
Surface Tension
0.0794
Surface Tension
0.0794
Umax=0.015
Ca=0.081, 0= 0.833
Umax=0.015
Ca=0.099, O= 0.80
Umax=0.025
Ca=0.119,0=0.815
Umax=0.025
Ca=0.153,0=0.76
Umax=0.035
Ca=0.159,0=0.78
Umax=0.035
Ca=0.322,0=0.732
      
  
Figure 6: The resulting flow field relative to the interface
and streamlines in the waterfront case shown in Figure
4.
experimental work of Tabeling et al (1987) by releasing
small fishskin particles into the interfacial flow. It is
also reported in numerical simulations: for example, in
the 2D simulation by Reinelt & Saffman (1985) and the
3D simulation by Wong et al (1995a,b). For the water
front case, the flow relative to the interface exhibits a
smaller vortex pattern near contact line at the interfacial
region.
In Figure 7, we show the normalized finger width
for the air front case from our simulations as a func
tion of the Capillary number. The fitting, Eq. 12,
from the experimental work of Tabeling et al (1987)
with W/b + oo as well as the regression formula, Eq.
(10), from the numerical simulation by Halpern & Gaver
(1994) are plotted in Figure 7 for comparison. We chose
the viscosity ratio so that our Capillary numbers can ex
tend below the minimum of 0.2 simulated in Kang et
al (2004). There is a reasonable level of agreement be
1.10
1.00
0.90
0.80
0.70
I I I I
0.00 0.10 0.20
Ca
Figure 7: The normalized finger width as a function of
Capillary number.
tween our simulation data and these previous results.
This implies that our simulated interface shape are re
alistic for the parameter range we considered here.
We should point out that our Capillary number and
Reynolds number are much larger than those realized
in the experiment of Lazouskaya (2008). At this point,
it is not completely clear whether the fingering for the
air front case would remain in the limit of very small
Capillary number. While in Lazouskaya (2008), a liquid
film was found in the air front case, it is much thinner
in scale. Whether this is a result of smaller Ca number
with the same physical origin macroscopicc fingering in
stability) or it is of a different origin such as molecular
scale adhesion remains to be studied.
Case 3
rightlyl Itleft
Case 4
/right I / left
Targeted area when Ca goes to zero
Targeted area when Ca goes to zero
Transport of Colloids near the
FluidFluid Interface
The simulated flow field near the fluidfluid interface can
be used to study the motion of colloids near the interface
and contact line. This is relevant to colloid transport and
retention in unsaturated soil. We integrate the equation
of motion for colloids under the influence of drag force,
Brown fluctuations, and physicochemical forces. A de
tailed description of the equation of motion and relevant
parameters used for the colloids and solution conditions
can be found in Shi et al (2010). Keeping the same
parameters and treatment for colloid tracking as in Shi
et al (2010) but using the new flow fields discussed here,
we shall briefly study the trajectories of colloids near the
fluidfluid interface.
In the airfront case, the trajectories are plotted in Fig
ure 8 for several colloids released on the water side near
the walls away from the interface. This is because the
relative motion is towards the interface near the wall re
gion but away from the interface in the center region of
the channel (Shi et al 2010). We carefully chose the re
lease locations so that no colloidswall physicochemical
interaction occurred. The airwater interface is unfavor
able retention site since a repulsive electrostatic force is
applied. We set a cutoff distance for colloidsAWI in
teraction to 106 m, as stated in Shi et al (2010). For
the colloids shown in Figure 8, none of them approaches
close enough to the interface to experience this repul
sive interaction force. Thus the motion of colloids is
governed by hydrodynamic force and Brownian motion,
with the hydrodynamic force dominating the motion.
Therefore, colloids appear to follow the flow streamlines
relative to the interface.
When colloids travel into the interfacial region, we
observed that some of them are sliding along the inter
face towards the center of the channel without contacting
the interface. They move very slowly towards the stag
nation point located at the tip of the interface and may
appear to be temporally retained. The curved stream
lines therefore play a major role in determining the tra
jectories relative to the interface. There is also a pos
sibility for a colloid very close to the channel wall to
continue its path into the liquid films between the finger
and the channel walls. These observations are similar to
what were observed in Lazouskaya (2008).
The dependence of colloid transport on the detail flow
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
pattern implies that it only makes sense to compare dif
ferent micromodel experiments with similar geomet
rical and interfacial configurations. Since the contact
line, the interface shape, and corresponding interfacial
flow field are critical to colloids transport, care should
be taken, for example, when generalizing the conclusion
obtained from singlebubble experiments (Chen & Flury
2005) to airwater interface contacting grain surface in
porescale experiments (Johnson & Elimelech 1995).
Our simple 2D simulations so far provide a reasonable
agreement with our parallel experiments with regard to
the hydrodynamic effects on colloids transport. There
are, however, other complications that are not yet con
sidered, including colloidcolloid interactions when col
loids clustered at the AWI and additional Capillary force
when colloids are trapped at interface or near the con
tact line. The minimum velocity in our simulation is also
limited as we wanted to make sure that the flow velocity
is larger than spurious flow for most flow regions, which
prevents us from simulating a much smaller Capillary
number and associated micronsize thin liquid films at
tached to the channel walls. In real soil matrix, the inter
connected pores and irregular channel formed by soil
particles would make flow field much more complicated
than in our simple microchannel.
Colloids trajectories in the waterfront case are shown
in Figure 9. Here colloids released near the center fol
low the streamlines as shown in Figure 6, approaching
the stagnation region at the tip of the interface. Again,
since the flow is so slow there, colloids may appear to
move with the interface visually in the micromodel ex
periment (Lazouskaya et al 2006; Lazouskaya 2008).
The colloids then follow the curved interface and slowly
migrate towards the channel walls. They may reside
close to the contact line before returning to the bulk liq
uid side along the walls. These are similar to the visual
izations from the micromodel experiment, as shown in
Figure 15 of Shi et al (2010).
Conclusions and Summary
In this paper, we presented results on microscale inter
facial flow simulation and transport of colloids near the
fluidfluid interface, in order to better understand colloid
transport behavior in unsaturated soil porous media. At
this stage, the lattice Boltzmann method for twophase
flow is under rapid development. For the purpose of our
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
480 520 560 600 640 680
Figure 8: Colloid trajectories in the airfront setting from Case 3 with Umax = 0.025.
60.0 5 4B B 1 a
r) ,1 I ,1 11
68.0
 M#3 colloAld
.aa1 #2 cilod
##1 collkid
20.0
0.0
360 400 440 480 520
Figure 9 Colloids trajectory in the waterfront setting from Case 3 with U 0.025.
Figure 9: Colloids trajectory in the waterfront setting from Case 3 with Umax = 0.025.
application, the ability to treat large density and viscos
ity ratio is crucial. We have demonstrated that the use
of the RedlichKWong EOS in LBM can handle a high
density ratio. An important observation is that the den
sity difference has a small effect on flow field at dense
fluid (water) side, although the flow field on the low den
sity (gas) side is greatly affected with larger spurious
currents. On the water side, not only do the interface
shapes from different density ratio simulations overlap
with each other, the velocity profiles are also insensitive
to the density ratio. As far as the waterside flow field
in the vicinity of the interface is concerned, we may ne
glect the influence of density difference and just focus
on the effect of the viscosity ratio.
We then studied the effect of viscosity ratio. First, the
interface shape evolutions for the airfront and the water
front cases are simulated. Consistent with previous ex
perimental observations (Saffman & Taylor 1958), we
found that the finger instability is associated with the
case of a less viscous fluid driving a more viscous fluid
(the air front case), while a stable interface is maintained
for the water front case. Quantitatively, we compared
our simulated finger width and its dependence on the
Capillary number to previous experimental and numeri
cal results (Tabeling et al 1987; Halpem & Gaver 1994)
and concluded that a good agreement was achieved. The
flow pattern near the interface is also found to be quali
tatively similar to previous results (Shikhmurzaev 1997,
2008).
Colloid trajectories relative to the interface were ob
tained by integrating the equation of motion. For both
the airfront and the waterfront cases, we found that col
60.0
40.0
20.0
0.0
loids follow closely the flow streamlines. Even though
the airwater interface is unfavorable collector by design,
the colloids can still be brought close to the interface and
may appear to be temporally retained due to flow stag
nation. They could slide along the interface primarily
due to the hydrodynamic forces. These observations are
similar to experimental observations in Lazouskaya et
al (2006); Lazouskaya (2008). These results must be
viewed as preliminary as the parameter conditions are
not exactly the same as the experimental conditions and
there are a number of simplifications in our current nu
merical models.
Acknowledgements
This study is supported by the US Department of Agri
culture (NRI200602551, NRI200802803), US Na
tional Science Foundation (NSF CBET0932686), and
National Natural Science Foundation of China (Project
No. 10628206). We thank Dr. Qinjun Kang for provid
ing his LBM codes and for his comments on LB multi
phase flow models.
References
Chen, G., and M. Flury Retention of mineral colloids
in unsaturated porous media as related to their surface
properties, Colloids Surf. A, vol 256, pp207216, 2005
Crist J.T., McCarthy J.R, Zevi Y., Baveye P., Throop J.A.
and Steenhuis T.S., Porescale visualization of colloid
transport and retention in partly saturated porous media,
Vadoze Zone J., vol 3, pp444450, 2004
Fortner, J. D.; Lyon, D. Y.; Sayes, C. M.; Boyd, A. M.;
Falkner, J. C.; Hotze, E. M.; Alemany, L. B.; Tao, Y.
J.; Guo, W.; Ausman, K. D.; Colvin, V. L.; Hughes,
J.B., C60 in water: Nanocrystal formation and microbial
response. Environ. Sci. Technol., Vol 39, pp43074316,
2005
Gao H, Qiu CQ, Fan D, Jin Y.. Threedimensional mi
croscale flow simulation and colloid transport model
ing in saturated soil porous media. Comput. Math. appl.
2008
Halpem D. and Gaver D.P, Boundary Element Analysis
of the time dependent motion of a semiinfnite bubble
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
in a channel, J. Computational Phys. vol 115, pp366
375, 1994
He X., Doolen G.D., Thermodynamic foundations of
kinetic theory and lattice Boltzmann models for multi
phase flows, J. Stat. Phys.vol 107,pp309328,2002
Johnson, W.P, and Elimelech M. Dynamics of colloid
depositionin porous media: Blocking based on ran
dom sequential adsorption. Langmuir vol 11, pp801
812, 1995 transport and retention in partly saturated
porous media, Vadoze Zone J., vol 3, pp444450, 2004
Johnson W.P., Li X. and Yal G., Colloids retention
in porous media: Mechanistic confirmation of wedg
ing and retention in zones of flow stagnation. Environ.
Sci.Technol. Vol 41, ppl2791287, 2007
Kang Q, Zhang D. et al. Immiscible displacement in a
channel: simulations of fingering in two dimensions. Ad
vances in Water Resources, vol 27,ppl322, 2004
Kang Q, Zhang D et al, Displacement of a two
dimensional immiscible droplet in a channel, Physics of
FLuids, vol 14, pp32033214,2002
Lazouskaya V.I., Jin Y. and Or D., Interfacial interac
tions and colloids retention under steady flows in a cap
illary channel, J.Colloids and Interface Scie, vol 303, pp
171184, 2006
Lazouskaya V.,Ph. D dissertation, Porescale investiga
tion on mechanisms of colloid retention in unsaturated
porous media, University of Delaware, 2008.
Lee T and Fischer P, Eliminating parasitic currents in the
lattice Boltzmann equation method for nonideal gases,
Physical review E, vol 74, 046709, 2006
McCarthy J.F, Zachara J.M. Subsurface transport of
contaminants mobile colloids in the subsurface envi
ronment may alter the transport of contaminants. En
rivon, Sci. Technol., Vol 23 pp496502, 1989
Reinelt D.A. and Saffman P.G., The penetration of a fin
ger into a viscous fluid in a channel and tube, J. Sci. Stat.
Comput., vol 6 No.3, 1985
Shi G. X., Gao H.,Lazouskaya V.I. et al, Visous flow
and colloids transport near airwater itnerface in a mi
crochannel, Computers and Mathematics with Applica
tions, in press, 2010
Shan X., Chen H., Simulation of nonideal gases and
liquidgas phase transitions by the lattice Boltzmann
equations. Phys. Rev. E. Vol 49, pp29412948,1994
Shan X, Analysis and reduction of the spurious current
in a class of multiphase lattice Boltzmann models. Phys
ical Review E, vol 73, 047701, 2006
Shikhmurzaev Y D, Moving contact lines in liq
uid/liquid/solid systems, J. Fluid Mech. vol. 334, pp.
211249,1997
Shikhmurzaev Y D., Capillary flows with forming inter
face, by Taylor & Francis Group, 2008
Saffman P.G and Taylor G, The penetration of a fluid
into a porous medium or HeleShaw cell containing a
more viscous liquid. c. Roy. Soc. A., Vol 245, pp312
329,1958
Shan X., Chen H., Lattice Boltzmann model of simulat
ing flows with multiple phases and components, Phys.
Rev. E, vol47, ppl8151819, 1993
Saiers J.E. and Ryan J.N., Introduction to special section
on colloid transport in subsurface environments. Water
Resources Res., Vol 42, W12S01, 2006.
Sen T.K. and Khilar KC, Review on subsurface colloids
and colloidassociated contaminant transport in satu
rated porous media. Adv. Colloid and Interf. Sci. Vol
119, pp7196, 2006
Tabeling P, Zocchi G. and Libchaber A., An experi
mental study of the SaffmanTaylor instability, J. Fluid
Mech. vol 177, pp6782, 1987
Saffman P.G and Taylor G, The penetration of a fluid
into a porous medium or HeleShaw cell containing a
more viscous liquid. c. Roy. Soc. A., Vol 245, pp312
329,1958
Wong H, Morris S, Radke CJ.. The motion of long bub
bles in polygonal capillaries. Part I Thin films. J. Fluid
Mech. vol 292, pp7194, 1995a
Wong H, Morris S, Radke CJ.. The motion of long bub
bles in polygonal capillaries, Part II Drag, fluid pressure,
and fluid flow. J. Fluid Mech. 292, pp:95 110,1995b
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Wan, J., and J. L. Wilson, Vizualization of the role of the
gaswater interface on the fate and transport of colloids
in porous media, Water Resour. Res., vol 30, ppll1123.
1994
Wan, J., and J. L. Wilson, Visualization of colloid trans
port during single and twofluid flow in porous media,
in: J. F. McCarthy, F. J. Wobber (Eds.), Manipulation
of Groundwater Colloids for Environmental Restoration,
Proceedings of a Workshop in Manteo, North Carolina,
October 1518, 1990, pp. 335339. 1993
Wiesner, M. R.; Lowry, G. V.; Alvarez, P.; Dionys
iou, D.; Biswas, P., Assessing the risks ofmanufactured
nanomaterials. Environ. Sci. Technol., VO140, pp4338
4345, 2006
Peng Y.,Schaefer L., Equations of state in a lattice Boltz
mann model, Phys. Fluids., Vol. 18, 042101, 2006 trans
port and retention in partly saturated porous media,
Vadoze Zone J., vol 3, pp444450, 2004
Zhuang, J., J. F. McCarthy, J. S. Tyner, E. Perfect, and
M. Flury In situ colloid mobilization in Hanford sedi
ments under unsaturated transient flow conditions: Ef
fect of irrigation pattern, Environ. Sci. Technol., v41,
pp31993204.2007
