Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Direct, Dynamic Measurement of Interfacial Area within Porous Media
D. Crandall1, G. Ahmadi2, D.H. Smith3, and G. Bromhal4
1URS Washington Division, National Energy Technology Laboratory
Morgantown, WV 26507 USA
Dustin. Crandall@ur.netl.doe.gov
2Mechanical and Aeronautical Engineering Department, Clarkson University
Potsdam, NY 13699, USA
gahmadi@tclarkson.edu
National Energy Technology Laboratory
Morgantown, WV 26507 USA
Duane. Smith @netl.doe.gov
4National Energy Technology Laboratory
Morgantown, WV 26507 USA
Grant.Bromhal@netl.doe.gov
Keywords: Twophase flow, microscale experiments
Abstract
Standard models of twophase flow in porous media have been shown to exhibit several shortcomings that might be partially
overcome with a recently developed model based on thermodynamic principles (Hassanizadeh and Gray, 1990). This
alternative twophase flow model contains a set of new and nonstandard parameters, including specific interfacial area. By
incorporating interfacial area production, destruction, and propagation into functional relationships that describe the capillary
pressure and saturation, a more physical model has been developed. Niessner and Hassanizadeh (2008) have examined this
model numerically and have shown that the model captures saturation hysteresis with drainage/imbibition cycles. Several
static experimental studies have been performed to examine the validity of this new thermodynamically based approach; these
allow the determination of static parameters of the model. To date, no experimental studies have obtained information about
the dynamic parameters required for the model.
A new experimental porous flow cell has been constructed using stereolithography to study twophase flow phenomena
(Crandall et al. 2008). A novel image analysis tool was developed for an examination of the evolution of flow patterns during
displacement experiments (Crandall et al. 2009). This analysis tool enables the direct quantification of interfacial area between
fluids by matching known geometrical properties of the constructed flow cell with locations identified as interfaces from
images of flowing fluids. Numerous images were obtained from twophase experiments within the flow cell. The dynamic
evolution of the fluid distribution and the fluidfluid interface locations were determined by analyzing these images. In this
paper, we give a brief introduction to the thermodynamically based twophase flow model, review the properties of the
stereolithography flow cell, and show how the image analysis procedure has been used to obtain dynamic parameters for the
numerical model. These parameters include production/destruction of interfacial area as a function of saturation and capillary
pressure. Our preliminary results for primary drainage in porous media show that the specific interfacial area increased
linearly with increasing gas saturation until breakthrough of the displacing gas into the exit manifold occurred.
Introduction applications range from the movement of fluids and
particulates in the subsurface, to brain and liver cancer
A number of engineering problems require knowledge of treatment, to processes occurring during paper
flow in porous media, including applications in the manufacturing and within fuel cells. Modeling of these
environment, biology, industry, oil & gas recovery, and systems is usually performed with a standard set of
geologic sequestration of carbon dioxide. These equations based upon Darcy's Law (Darcy, 1856); a century
Paper No
and a half old phenomenological relationship describing the
flow of a single fluid through a homogenous porous domain.
For systems consisting of two fluids moving within a
porous medium, an extended form of Darcy's law to
determine fluid phase velocities and saturations has been
used for decades (Helmig, 1997). From these velocities and
saturations functional relationships to determine the
capillary pressure (pc) and relative permeabilities (kr) of the
system are typically used. This method is deficient with
respect to physics in that it completely neglects the role and
presence of fluidfluid interfaces and leads to hysteretic
relationships. Interfaces between fluids are not only crucial
quantities for interfacedependent processes, such as mass
transfer between phases, but there is a growing amount of
experimental evidence that shows a relationship between pc,
interfacial area, and saturation (S) accounts for the observed
hysteresis in the krS relationship by collapsing the curves
upon a single surface (Culligan et al. 2004, Brusseau et al.
2006, Chen and Kibbey 2006, PyrakNolte et al. 2008).
Alternative theories of multiphase flow and transport in
porous media have been developed. Among these are
theories based on rational thermodynamics by Hassanizadeh
and Gray (1980, 1990) and theories based on
thermodynamically constrained averaging theory by Gray
and Miller (2005, Miller and Gray 2005). In addition, there
is an approach which separates balance equations for
percolating and nonpercolating part of each phase (Hilfer
and Doster, 2009). In this work we use the approach of
Hassanizadeh and Gray (1980, 1990), which is based on
thermodynamic principles and is more physicallybased
than the standard modified Darcy's Law model. This new
approach not only includes balance equations for the bulk
phases, but also an interfacial balance equation.
Furthermore, in this theory pc is not a function of S only, but
also depends on specific interfacial area (an). A modeling
approach using this theory with idealized parameters was
put into practice by Niessner and Hassanizadeh (2008) that
includes balance equations for two bulk fluid phases and the
fluidfluid interface, a functional relationship between an,
Pc, S, and a production rate term of interfacial area (Ewn).
In order to properly use this new model the necessary
parameters must be obtained. While many static tests have
been performed to calculate relations between awn, Pc, and S,
e.g. (Culligan et al. 2004, Brusseau et al. 2006, Chen and
Kibbey 2006, PyrakNolte et al. 2008), there are some
parameters that cannot be determined by standard, static
experimental procedures. In this work, we present a
combination of new techniques that overcome this
discrepancy and allow for an estimation of the needed
dynamic model parameters.
A micromodel of an idealized porous medium, i.e. a flow
cell, was created from a computer aided drafting (CAD)
model using stereolithography (Crandall et al. 2008). The
flow cell geometry was created with an inhouse computer
code. The fabricated physical model geometry was known
from the CAD model and thus the size of each individual
interface within the flow cell could be determined by
combining image analysis with the original computer model.
Experiments of air invading an initially waterfilled flow
cell were performed. By recording images periodically
during flow experiments, data on dynamic evolution of
interfaces, awn, Pc, S, and velocity of interfaces were
determined.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Nomenclature
an Specific interfacial area between wetting and
nonwetting fluids (m2)
Ca Capillary number ()
Ew Production rate of awn (m2s1)
g Gravitational force (N)
h, Height of throat (m)
j, Number of interfaces
K Permeability (m2)
kr Relative permeability
p Pressure (Nm2)
pc Capillary pressure (Nm2)
q Mean volumetric flow rate (m2s')
S Fluid saturation ()
t Time (s)
v Fluid velocity (ms1)
Wt Width of throat (m)
Greek
p
p
letters
Porosity ()
Fluid density (kgims.
Absolute viscosity (Pas)
Subscripts
n Nonwetting fluid
w Wetting fluid
TwoPhase Flow Model
Hassanizadeh and Gray (1990) have derived general
balance equations for twophase flow based on
thermodynamic principles. However, it is not possible to
numerically model this general system as it would require
an incredible amount of computational power and because a
number of the involved constitutive relationships have not
been determined. Niessner and Hassanizadeh (2008)
modeled a simplified version of the general equation system,
which represents a relatively simple extension of the
standard twophase flow model, but still includes an.
Their procedure resulted in a system of three partial
differential equations, one for each of the two phases and
one for the fluidfluid interface. Additional constraints
included summing the wettingphase saturation (S,) and
nonwetting phase saturation (Sn,) to unity, setting the pc as
the difference between nonwetting pressure (pnw) and
wettingphase pressure (p,), and a functional relationship
between awn, Pc, and S was provided. This yielded the
following system of equations (Niessner and Hassanizadeh
2008):
a0 + V v.v =0
at
k
v =Kk (Vp.
vww
PSw
n +V. =0
at
k 
v, K (Vp,
fun
pg)
pg)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
anw + V av }
at
vw = K wn Vawn
S, + S, = 1
P, Pw = Pc
aw, = aw,(Sw, Pc)
resistances of other flow devices in the literature is listed in
(5) Table 1 of Crandall et al. (2008). Two manifolds were
created on opposing edges of the porous matrix, arbitrarily
(6) labeled #1 and #2 in Figure 1. Air was injected into the
centered port of one of these manifolds and through the
(7) initially watersaturated flow cell to perform the
(8) experiments.
where v is the Darcy velocity, p is the dynamic viscosity, t is
time, 0 is the porous medium porosity, K is the intrinsic
permeability, v. is the interfacial velocity, and K,, is
interfacial permeability. The an and the E,, are dependent
on the primary variables S, and pc. Relative permeabilities
are only dependent on S, (Helmig 1997)
To model the system given in Eqn. (1) to Eqn. (9) the S,, pc,
and aw, must be determined. A functional relationship
an, (S,, pc) was used by Niessner and Hassanizadeh (2008).
Ideally, a number of subsequent static drainage and
imbibition experiments would be needed to obtain a reliable
relationship. The variation of the a,, (S,, pc) relationship
with different porous media will become more apparent as
various researchers study different applications (Brusseau et
al. 2006, PyrakNolte et al. 2008), but currently a generic
relationship from porethroat modeling suffices
(JoekarNiasar et al. 2008)
Dynamic properties that have not been determined in
previous studies had to be estimated in the recent model of
Niessner and Hassanizadeh (2008). These include the
interfacial velocity v.n, interfacial permeability K., and the
production rate of specific interfacial area Ewn. The
remainder of this paper will review to experimental devices
and techniques we used to determine these parameters and
present preliminary results.
Experimental Setup
The porous medium used for our experiments was
constructed using stereolithography (SL) rapid prototyping,
a production technique that has been used to make complex
parts from computer aided design models. SL models are
constructed by curing successive layers of photosensitive
resin with a laser to form 3D objects. Our SL flow cell was
constructed using a 3D Systems Viper Si2
stereolithography apparatus out of DSM 11120 Watershed
(DSM Somos, New Castle DE, USA), a waterresistant
resin. The flow cell production procedure is reviewed in
detail by Crandall et al. (2008).
The porous matrix of the flow cell, labeled in Figure 1, was
constructed as a squarelattice of over 5000 throats within a
10.16 cm by 10.16 cm region. The individual throat widths
varied from 0.35 mm to 1.0 mm. These throat widths were
randomly distributed throughout the porous matrix. To
introduce a greater range of porelevel resistances into the
flow cell, seven different throat heights were assigned to
these throats, varying from 0.2 mm to 0.8 mm. These throat
heights were assigned in such a manner so as to reduce the
aspect ratio of the throats. That is, the narrowest throats
were assigned to be the shortest (smallest throat width =
smallest throat height), the widest throats were made the
tallest, and so forth. A detailed comparison of the throat
capillary and inertial resistances of this SL flow cell to
10.16 cm
Figure 1: Photograph of stereolithography flow cell.
A relatively simple experimental setup was used to control
the fluid flow into the flow cell and to capture the images. A
schematic of the experimental setup is shown in Figure 2
and consists of the SL flow cell, a constantrate syringe
pump (KD Scientific KDS 200), a CCD camera (NTSC
COHU 4915400/000), lighting, a collection vessel at
atmospheric pressure, and a data acquisition computer. A
LabVIEWTM (version 7.1) module was used to control the
imagecapture rate. The same experimental setup was used
by Crandall et al. (2009) to capture images of fluid motion
in the flow cell. A postprocessing routine was written to
determine the size and distribution of fluid bursts associated
with Haines jumps and these were shown to be well
described by the model of SelfOrganized Criticality
(Crandall et al. 2009). For the current study this
postprocessing code was changed to identify the location
and size of individual interfaces between fluids in the flow
cell. By linking the known height and width of the flow cell
at each location an interface was identified the fluidfluid
interface was determined. For the results presented here the
individual interfaces were estimated as the crosssectional
area of the throat they were identified in. While the
curvature of the interface at each location would increase
the total interfacial area, for this preliminary study this
approximation was deemed satisfactory.
Camera
FI ll
Flowcell
IPump
Computer
Figure 2: Schematic of experiment.
Paper No
Paper No
Injection of air at different flow rates was conducted. The
ratio of inertial to capillary forces for the different flow rates
was quantified using the following definition of the
capillary number,
Ca = (10)
where q is the mean volumetric flow across the flow cell, p.
is the wetting fluid viscosity and yis the interfacial tension
between the injected air and the defending water (y = 72
mi/m). Here q is defined as the injected volumetric flow rate
divided by the mean crosssectional area of the flow cell
perpendicular to the flow direction. For this study the
volumetric flow rate was varied from a maximum of
20 m"'l/ to a minimum of 0.002 ml/mm, which corresponds to
4.63(109) < Ca < 4.63(1013). All experiments were
conducted with the flow cell initially saturated with distilled
water and air injected into either Manifold #1 or Manifold
#2 at a constant rate. This fluid pairing has a viscosity ratio
(M = i,,r/Iater) of 0.0178, which indicates the drainage case
studied is in the unstable regime (Lenormand et al. 1988).
In order to determine the macroscale parameters that are
needed for the numerical model the size of an averaging
volume (Representative Elemental Volume, REV) needed to
be determined for the flow cell. We subdivided the flow cell
in n x n subvolumes and calculated the 0 of each REV as a
test parameter. As the flow cell was constructed using a
random distribution of throats and pores the flow cell as a
whole is considered a homogeneous porous medium with a
porosity of 35.79% (Crandall et al. 2008). The subdivision
of the flow cell to determine the REV for which 0 varies
little using the largest possible value for n allowed us to
consider multiple regions within the porous medium as
homogeneous. Table 1 shows average, maximum, minimum
values of 0, as well as the standard deviation for n = 2 to 6
sized REVs. As is shown the standard deviation between the
0 values were below 1% for n = 2 and 3, i.e. 4 and 9
separate REVs. These lower values of variation between
separate volumes were thought to be acceptable for
reporting of a,, values, and thus the results of E,, are
presented for each. Two images of air invasion into the flow
cell are shown in Figure 3 with the n = 2 and n = 3 REVs
identified.
Avg. 4 (%)
Min. ( (%)
Max. 4 (%)
Std. Dev. (%)
2 35.89 35.41 36.71 0.57
3 35.61 34.23 37.21 0.92
4 35.99 34.16 38.82 1.28
5 35.77 32.08 38.87 1.66
6 35.63 29.67 39.14 2.01
Table 1: Average, minimum and maximum porosity (4) of
potential REVs within the flow cell, and the standard
deviation of the values.
Results and Discussion
The number of identified interfaces within the entire flow
cell for an 11 hour, 0.002 ml/mm flow rate experiment is
shown in Figure 4 as a function of time. Also the interfaces
identified in each of the four n = 2 REVs are shown. The
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
overall number of interfaces increases approximately
linearly while the number of interfaces in each REV is
observed to increase sporadically as the invading air moves
in only one REV at a time at this low flow rate. This pattern
is similar to the measured S,, of air shown in Figure 5 as a
function of time. These S values were determined by
summing the volume of invaded flow cell regions and
dividing by the total volume of the flow cell, or the volumes
within the REVs.
To determine the a,, within the cell the interfacial area was
summed and divided by the total volume of the flow cell.
These values are shown in Figure 6 for the same 11 hour,
0.002 ml/mm flow rate experiment in both the entire cell and
within the n = 2 REVs as a function of the S,,. As can be
seen the a,, increases approximately linearly for all the
measured regions, with similar slopes both in the entire cell
and in each REV
n=2 n=3
Figure 3: REVs for n = 2 and n = 3 shown on top of
experimental images of air in the water filled flow cell.
1500
1000 Entire Cell
Number
of
n = 2
S REV
10
Time (hr)
Figure 4: Identified interfaces in the entire flow cell and in
each of the n = 2 REVs.
Entire
50 cell
40 n=Z
io
of 30
Time (hr)
Figure 5: Saturation of air in the entire flow cell and in
each of the n = 2 REVs.
Paper No
1.5
Spedfic
Interfacial
Area 1.0
(l/n)
O.S
Entire
Cell
n..
=I REV
*O 1 2 0 0 5
0 10 20 30 40 s0
saturation of Invadhg Air (%)
Figure 6: a., in each of the n = 2 REVs and the flow cell.
The pressure across the entire flow cell was measured with a
pressure transducer, as shown in the experimental schematic
in Figure 2. These pressure values are shown in Figure 7 for
the same 11 hour, 0.002 m"l/m flow rate experiment shown in
Figures 4 6. The pressure increased gradually to
approximately 675 Pa, with fluctuations of about 100 Pa
every 15 minutes. To estimate the average Pc of the
interfaces the individual pc at each throat was determined,
2Y
Pc=
h, + w,
where h, is the height of the throat the interface was
measured in and w, is the width. The average p, was
determined by summing the individual pc and dividing by
the total number of interfaces measured. This average Pc is
shown in Figure 8 for the entire cell and within each of the n
= 3 REVs. The average Pc value also fluctuates significantly
over the experiment and is within the same range of values
as the macroscopically measured pressure difference across
the cell.
7 600
500
s.oo
400
0 2 4 6 8 10
Time (hours)
Figure 7: Pressure measured across the flow cell.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
S REV
Time (min)
Figure 9: v,, in each of the n = 2 REVs the x direction.
2.0 n = 2
REV
(cm/s)
0 2 4 6 8
Time (min)
Figure 10: vw in each of the n = 2 REVs the y direction.
Our preliminary measurements of the interfacial velocity
(vw,) reveal a complex and sporadic pattern of interfacial
motion. The vw, was determined for both the x and the y
directions, where y is in the bulk flow direction. By
averaging the number of interfaces that moved a distance in
the x (Ix) and y (I,) directions in between successive images
the v,, was determined for both the x and the y direction.
,C 
I ly
(v ~) =412
kny I At
IY 1 1I1 Atd
As shown in Figures 9 and 10 this resulted in a chaotic
pattern of velocity values in both positive and negative
directions. This jumpy pattern is due to the occurrence of
Haines jumps as the nonwetting air quickly invades water
filled pores (Crandall et al. 2009). By averaging over a
reasonable amount of time we are hoping to see a smoother
description of the v,, emerge from these initial results. The
overall average of the v,, in the x direction is approximately
zero, as would be expected since the bulk fluid motion is in
the y direction. Once we have a reasonable set of velocity
measurements we will determine the K,, using the
following relation,
Kjw wn / oawn
/ aj
Time (hr)
Figure 8: an, in each of the n = 2 REVs and the entire cell.
and E,, will be determined via Eqn. 5. Both of these
dynamic quantities rely on an accurate description of the
Vw, on a macroscopic scale.
Paper No
Conclusions
We have performed experiments of air flowing within a
stereolithography flow cell. Images of air invasion were
captured during the experiments. With knowledge of the
porous medium geometry the location and size of the
interfaces within the flowing air structure was determined.
The number of interfaces was observed to increase in a
linear fashion during the primary drainage flows evaluated
for this paper. The capillary pressure determined from the
image analysis technique was shown to be in good
agreement with the pressure measured across the flow cell
via a pressure transducer. Initial measurements of the
interfacial velocity show a great deal of scatter due to the
sporadic motion of the interface at the pore level. This is
because of the small capillary resistances being overcome
as the nonwetting fluid displaces the water. By averaging
these velocity values we hope to be able to determine
interfacial production and permeability values for our
experiments.
Acknowledgements
The authors extend their gratitude to Jen Niessner and Majid
Hassanizadeh for the impetus to study this topic and
conversations on the matter. Also we wish to thank Martin
Ferer for his knowledgeable talks on flow in porous media.
References
Brusseau, M.L., S. Peng, G Schnaar, and M.S.
CostanzaRobinson. 2006. Relationships among airwater
interfacial area, capillary pressure, and water saturation for
a sandy porous medium. Water Resour. Res. 42:W03501.
Chen, L., and T.C.G. Kibbey. 2006. Measurement of
airwater interfacial area for multiple hysteretic drainage
curves in an unsaturated fine sand. Langmuir 22:66746880.
Crandall, D., G. Ahmadi, M. Ferer, and D.H. Smith. 2009.
Distribution and occurrence of localizedbursts in
twophase flow through porous media. Physica A
38:574584.
Crandall, D., G. Ahmadi, D. Leonard, M. Ferer, and D.H.
Smith. 2008. A new stereolithography experimental porous
flow device. Rev. Sci. Instrum. 79(4):044501.
Culligan, K.A., D. Wildenschild, B.S.B. Christensen, W.
Gray, M.L. Rivers and A.F.B. Tompson. 2004. Interfacial
area measurements for unsaturated flow through a porous
medium. Water Resour. Res. 40(12):112.
Darcy, H. 1856. Les fontaines publiques de la ville de Dijon.
Dalmont, Paris.
Hassanizadeh, S.M. and W.G Gray. 1980. General
conservation equations for multiphase systems: 3.
Constitutive theory for porous media flow. Adv. Water
Resour. 13(4):169186.
Hassanizadeh, S.M. and W.G Gray. 1990. Mechanics and
thermodynamics of multiphase flow in porous media
including interphase boundaries. Adv. Water Res.
13(4):169186.
Gray, W.G and C.T. Miller. 2005.Thermodynamically
constrained averaging theory approach for modeling of flow
in porous media: 1 Motivation and overview. Adv. Water
Resour. 28(2):161180.
Helmig, R. 1997. Multiphase Flow and Transport Processes
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
in the Subsurface. Springer, Berlin.
Hilfer, R. and F. Doster. 2009. Percolation as a basic
concept for macroscopic capillarity. Trans. Porous Media.
JoekarNiasar, V, S.M. Hassanizadeh, and A. Leijnse. 2008.
Insights into the relationship among capillary pressure,
saturation, interfacial area and relative permeability using
porescale network modeling. Trans. Porous Media
74(2):201219.
Lenormand, R., E. Touboul, and C. Zarcone. 1988.
Numerical models and experiments on immiscible
displacements in porous media. J. Fluid Mech.
189:165187.
Miller, C.T. and W.G. Gray. 2005.Thermodynamically
constrained averaging theory approach for modeling of flow
in porous media: 2. Foundation. Adv. Water Resour.
28(2):181202.
Niessner, J., and S.M. Hassanizadeh. 2008. A model for
twophase flow in porous media including fluidfluid
interfacial area. Water Resour. Res. 44, W08439.
PyrakNolte, L.J., D.D. Nolte, D.Q. Chen, and N.J.
Giordano. 2008. Relating capillary pressure to interfacial
areas. Water Resour. Res. 44:W06408.
