vvEPA
United States
Environmental Protection
Agency
Office of Research and
Development
Washington DC 20460
EPA/600/R-99/022
March 1999
Evaluation and Analysis of
Microscale Flow and
Transport During
Remediation
-------
-------
EPA/600/R-99/022
March 1999
EVALUATION AND ANALYSIS OF MICROSCALE FLOW AND
TRANSPORT DURING REMEDIATION
by
Yanis C. Yortsos
and
Katherine Shing
University of Southern California
Department of Chemical Engineering
Los Angeles, CA 90089-1211
CR-824592
Project Officer
Jong Soo Cho
Subsurface Protection and Remediation Division
National Risk Management Research Laboratory
Ada, Oklahoma 74820
NATIONAL RISK MANAGEMENT RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
CINCINNATI, OHIO 45268
Printed on Recycled Paper
-------
NOTICE
The U.S. Environmental Protection Agency through its Office of Research and
Development funded and managed the research described here under Cooperative Agreement
Number CR-824592 to the University of Southern California, Los Angeles, CA. It has been
subjected to the Agency's peer and administrative review, and has been approved for publication
as an EPA document. Mention of trade names or commercial products does not constitute
endorsement or recommendation for use.
All research projects making conclusions or recommendations based on environmentally
related measurements and funded by the Environmental Protection Agency are required to
participate in the Agency Quality Assurance Program. This project did not involve
environmentally related measurements and did not involve a Quality Assurance Project Plan.
-------
FOREWORD
The U.S. Environmental Protection Agency is charged by Congress with
protecting the Nation's land, air, and water resources. Under a mandate of national
environmental laws, the Agency strives to formulate and implement actions leading to a
compatible balance between human activities and the ability of natural systems to support
and nurture life. To meet these mandates, EPA's research program is providing data and
technical support for solving environmental problems today and building a science
knowledge base necessary to manage our ecological resources wisely, understand how
pollutants affect our health, and prevent or reduce environmental risks in the future.
The National Risk Management Research Laboratory is the Agency's center for
investigation of technological and management approaches for reducing risks from
threats to human health and the environment. The focus of the Laboratory's research
program is on methods for the prevention and control of pollution to air, land, water, and
subsurface resources; protection of water quality in public water systems; remediation of
contaminated sites and ground water; and prevention and control of indoor air pollution.
The goal of this research effort is to catalyze development and implementation of
innovative, cost-effective environmental technologies; develop scientific and engineering
information needed by EPA to support regulatory and policy decisions; and provide
technical support and information transfer to ensure effective implementation of
environmental regulations and strategies.
The design of in-situ remediation is currently based on a description at the macro-
scale. Phenomena at the pore and pore-network scales are typically lumped in terms of
averaged quantities, using empirical or ad hoc expressions. These models cannot address
fundamental remediation issues at the pore and pore network scales, including: The
emplacement in-situ of the contaminant NAPL, and the displacement patterns; the mass
transfer of the contaminant from the NAPL to the groundwater or from the groundwater
to a sparging fluid, and of the remedial agents to the NAPL; and possible microscale flow
instabilities during injection of a remedial fluid. The objective of this work is to obtain a
fundamental understanding by conducting theoretical, experimental and computational
pore-scale studies. Emphasis is placed at the pore network scale. Use of this information
can be incorporated in macroscopic simulators to provide fundamentally correct
expressions for the various coefficients or parameters, currently treated empirically.
Clinton W. Hall, Director
Subsurface Protection and Remediation Division
National Risk Management Research Laboratory
in
-------
ABSTRACT
The design of in-situ remediation is currently based on a description at the
macroscopic scale. Phenomena at the pore and pore-network scales are typically lumped
in terms of averaged quantities, using empirical or ad hoc expressions. These models
cannot address fundamental remediation issues at the pore and pore network scales,
including: The emplacement in-situ of the contaminant NAPL, and the displacement
patterns; the mass transfer of the contaminant from the NAPL to the groundwater or from
the groundwater to a sparging fluid, and of the remedial agents to the NAPL; and possible
microscale flow instabilities during injection of a remedial fluid.
The objective of this work is to obtain a fundamental understanding by
conducting theoretical, experimental and computational pore-scale studies. Emphasis is
placed at the pore network scale. Use of this information can be incorporated in
macroscopic simulators to provide fundamentally correct expressions for the various
coefficients or parameters, currently treated empirically. The theoretical findings are
compared with findings from experiments in glass micromodels and Hele-Shaw cells.
Specific tasks of the work described include: (1) The determination of immiscible
displacement patterns that develop during seepage ofNAPLs in the groundwater or
during air sparging of contaminated groundwater. (2) The pore-scale study of the mass
transfer from a trapped NAPL to an injected remedial fluid. (3) A sensitivity analysis of
the effective mass transfer coefficient on parameters such as the flow rate and the patern
geometry. (4) Studies of pore-scale instabilities developing during miscible
displacements hi porous media, particularly when non-monotonic viscosity profiles are
involved. (5) The up-scaling of the information obtained to the macroscopic scale. (6)
Finally, the comparison of theoretical predictions with experimental results obtained in
etched-glass micromodels and Hele-Shaw cells in order to validate and improve the
theoretical models.
The research described here is expected to improve the fundamental
understanding of remediation processes, to help in the more accurate design of
remediation projects, to assess in the screening of new or proposed remediation
techniques and to help in the design of new ones. More generally, the research will
contribute to the advancement of the state of the science in pollution abatement.
This report was submitted hi fulfillment of CR-824592 by the University of
Southern California under the partial sponsorship of the United States Environmental
Protection Agency. This report covers a period from October 1, 1995 to September 30,
1997, and work was completed as of March 30,1998.
IV
-------
ACKNOWLEDGEMENTS
This research was supported by U.S. EPA Cooperative Agreement No. CR-
824592. The encouragement, help, and support of the Project Officer, Dr. Jong Soo Cho,
is gratefully acknowledged. We also acknowledge useful interactions with Dr. Carl
Enfield, Dr. Candida West, and Dr. Lynn Wood of the EPA's Subsurface Protection and
Remediation Division, Ada, Oklahoma.
Parts of this report are in various stages of publication in scientific journals.
Certain chapters of this report contain material from the PhD thesis of Chunsan Jia, from
current doctoral research by Maryam Shariati and from research to satisfy MS degree
requirements by Yuyong Zhang, all in the Chemical Engineering Department at the
University of Southern California. The various co-authors of each chapter are identified
at the beginnning of each chapter. We also note that the research reported in chapters 2
and 6 was also supported by DOE Contract No. DE-FG22-96BC14994/SUB.
-------
-------
TABLE OF CONTENTS
NOTICE ii
FOREWORD iii
ABSTRACT iv
ACKNOWLEDGEMENTS v
TABLE OF CONTENTS vii
LIST OF FIGURES viii
LIST OF TABLES xiii
Chapter 1. INTRODUCTION 1
Chapter 2. DISPLACEMENT PATTERNS USING PORE-NETWORK MODELS 8
Chapter 3. VISUALIZATION AND SIMULATION OF NAPL SOLUBILIZATION
IN PORE NETWORKS 44
Chapter 4. CONVECTIVE MASS TRANSFER FROM STATIONARY SOURCES
IN POROUS MEDIA 75
Chapter 5. PORE-SCALE INSTABILITIES DURING SOLVENT INJECTION 111
Chapter 6. A NOTE ON THE EFFECT OF DISPERSION ON THE STABILITY OF
NON-MONOTONIC MOBILITY PROFILES IN POROUS MEDIA 150
Chapter?. CONCLUSIONS 159
vii
-------
LIST OF FIGURES
2.1 Displacement patterns from numerical simulation of drainage in a 100 x 100 lattice,
with M = 0.1 and different values of the capillary number: (a) Ca = 3.2 x 10"8 (IP),
(b) Ca = 3.2x 10"6,(c) Ca = 3.2x 10'5and(d) Ca = 3.2x 10"4 23
2.2 Schematic of the front region and the notation used 24
2.3 Numerical solution of (2.5) vs. the modified capillary number for three different
values of M 25
2.4 Percolation probability (dashed fine) and saturation (solid line) profiles vs.
normalized length, for 3-D VGP: (a) B = 10"4 in a lattice of 200 x 200 x 65, (b) B
= 10'6 in a lattice of 200 x 200 x 500 26
2.5 Single finger in GP in a destabilizing gradient 27
2.6 Pore network simulations (saturation profile and displacement patterns) for
uniform permeability and Ca = 1.5 x 10"4 28
2.7 Pore network simulations (saturation profile and displacement patterns) for
permeability increase (1:2) and for Ca = 1.5 x 10"5 29
2.8 Pore network simulations (saturation profile and displacement patterns) for
permeability increase (1:2) and for Ca = 1.5 x 10"4 30
2.9 Displacement pattern with parallel gradient in a 2-D 800x 400 lattice: (1) Bx =
-0.1, (2)Bx = -0.01, (3)Bx = -0.001, (4) BJC = -0.0001 31
2.10 Displacement pattern with parallel gradient in a 3-D 128 x 128 x 128 lattice: (1)
Bx = -0.03, (2)5* = -0.003, (3) 5.x = -0.001, (4) Bx = -0.0003 32
2.11 The dependence of finger width on the Bond Dumber for (a) 2-D and (b) 3-D
lattices 33
2.12 Displacement pattern with transverse gradient in a 2-D 800 x 400 lattice: (1) By =
-0.1, (2) By = -0.01, (3) By = -0.001, (4) By = -0.0001 34
2.13 Displacement pattern with gradient in a 3-D 128 x 128 x 128 lattice: (1) By = 0.1,
(2) .By = -0.01, (3) By = -0.001, (4) By = -0.0005 35
2.14 The dependence of finger width oy on By for (a) 2-D and (b) 3-D lattices 36
2.15 Displacement pattern with both transverse and parallel gradient in a 2-D 800 x 400
lattice: (1) Bx = 0.01, By = -0.001, (2) Bx = - 0.01, By = -0.005, (3) Bx = 0.01, By
= -0.01, (4) Bx = 0.005, By = -0.01, (5) Bx = 0.001, By = -0.01, (6) Bx = 0.001,
By = -0.001 37
2.16 Displacement pattern with both transverse and parallel gradient in a 2-D 400 x 200
lattice: (1) Bx = 0.1, By = -0.01, (2) Bx = 0.1, By = -0.05, (3) Bx = 0.1, By = -0.1,
(4) Bx = 0.05, By = -0.1, (5) Bx = 0.01, By = -0.1 38
2.17 Displacement pattern with both transverse and parallel gradient in a 2-D 400 x 200
lattice: (1) Bx = -0.01, By = -0.01, (2) Bx = -0.001, By = -0.01, (3) Bx = -0.0005,
By = -0.01, (4) Bx = -0.01, By = -0.001, (5) Bx = -0.001, By = -0.001, (6) Bx =
-0.0005, By = -0.001 39
2.18 Displacement pattern with both transverse and parallel gradient in a 3-D 128 x 128
x 128 lattice: (1) Bx = -0.01, By = -0.1, (2) Bx = -0.01, By = -0.005, (3) Bx =
-0.01, By = -0.001, (4) Bx = -0.0005, By = -0.1, (5) Bx = -0.0005, By = -0.005, (6)
vin
-------
Bx = -0.0005, By = -0.001 40
2.19 The dependence of finger width oz on Bx for a 3-D lattice and various values of
By 41
2.20 The scaling function/(z) from simulations in: (a) 2-D lattice, (b) 3-D lattice 42
2.21 Modified destabilizing IPG patterns to model air sparging for B = -0.1 and P = 0,1
and 2, respectively 43
3.1 Schematic of the configuration of trapped NAPL in a porous medium 59
3.2 Typical pattern used for constructing a 2-D pore network 60
3.3 Schematic of the experimental apparatus 61
3.4 Photograph showing the steady-state NAPL distribution in the glass micromodel
for typical experimental conditions 62
3.5 A typical effluent curve for the 29 x 15 micromodel corresponding to the
solubilization of benzonitrile by water 63
3.6 Typical NAPL interface configurations in the porespace of the micromodel. Note
the cavity-like characteristics 64
3.7 Overall pressure drop vs. injection rate in single-phase flow in the micromodel for
two different cases: (a) in the absence of trapped NAPL (top part of the Figure),
and (b) in the presence of trapped NAPL (bottom part of the Figure) (for example
as shown in Figure 3.4). The solid line is the theoretical prediction and
experimental data are bracketed by error bars 65
3.8 Comparison between experimental (top) and numerical (bottom) concentration
profiles during miscible displacement of 50% glycerol solution for 0.091 pore
volumes injected 66
3.9 Comparison between experimental (top) and numerical (bottom) concentration
profiles during miscible displacement of 50% glycerol solution for 0.21 pore
volumes injected 67
3.10 Comparison between experiments and simulations for mass transfer from trapped
benzonitrile for pattern 1. The concentration at the effluent is normalized with its
theoretical solubility value. The solid line is the theoretical prediction,
experimental data are bracketed by error bars 68
3.11 Comparison between experiments and simulations for mass transfer from trapped
benzonitrile for pattern 2. The concentration at the effluent is normalized with its
theoretical solubility value. The solid line is the theoretical prediction and
experimental data are bracketed by error bars 69
3.12 Comparison between experiments and simulations for mass transfer from trapped
benzonitrile for pattern 1. Only diffusion was used for the simulation of local mass
transfer. The concentration at the effluent is normalized with its theoretical
solubility value. The solid line is the theoretical prediction, experimental data are
bracketed by error bars 70
3.13 Sensitivity analysis of the effect of parameter a in the local mass transfer
coefficient equation (3.13) on the effluent concentration curve. The effect of a is
minimal as all curves practically coincide (particularly at large values of the Peclet
IX
-------
number) 71
3.14 Sensitivity analysis of the effect of parameter b in the local mass transfer
coefficient equation (3.13) on the effluent concentration curve. Four different
curves corresponding to b = 0.2, 0.4, 0.664 and 1 (from left-to-right) are shown (a
= 1, c = 0.5. The effect of & is to shift the curve parallel to itself towards the right
(essentially rescaling the Peclet number) 72
3.15 Sensitivity analysis of the effect of parameter c in the local mass transfer coefficient
equation (3.13) on the effluent concentration curve. Seven different curves
corresponding to c = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, and 0.7 (from left-to-right) are
shown (a = 1, b = 0.664). Note the asymptotic approach to a straight line at large
values of the Peclet number 73
3.16 Numerical simulation results for the dependence of the slope s of the effluent
concentration curve on the parameter c of equation (3.13). The curve essentially
coincides with the theoretical result (3.14) 74
4.1 Photograph showing the steady-state NAPL distribution in the glass micromodel
for typical experimental conditions (from Jia et al., 1997) 94
4.2 Schematic of the unit cell used by Quintard and Whitaker (1994) 95
4.3 Typical NAPL interface configurations in the porespace of the micromodel. Note
the cavity-like characteristics (from Jia et al., 1997) 96
4.4 The normalized Sherwood number as a function of parameter a 97
4.5 Numerical results of Sh^ vs Pe^, for mass transfer over a flat plate and for three
local exponents (c = 0, c = 1/3 and c = Vz, from bottom to top, respectively) 98
4.6 Numerical results of Sh^g vs Pe^ for mass transfer over a sphere and for three local
exponents (c = 0, c = 1/3 and c = Vz, from bottom to top, respectively) 99
4.7 A Koch curve of 3 generations used for the simulation of mass transfer over a self
similar surface 100
4.8 Velocity profile for flow over the Koch curve of Figure 4.7 101
4.9 The local mass flux profile over the perimeter of the Koch curve of Figure 4.7 for
various values of the Peclet number 102
4.10 The local concentration profile over (a) the streamwise coordinate, (b) the
perimeter of the Koch curve of Figure 4.7, for various values of the Peclet number .. 103
4.11 Numerical results of ShOTg vs PeL for mass transfer over the Koch curve of Figure
4.7 for three different generations and for c - 0 104
4.12 Numerical results of Sh^ vs PeL for mass transfer over the Koch curve of Figure
4.7 for three different generations and for c = 1/3 105
4.13 Invasion percolation cluster used for the simulation of flow over a source with
percolation characteristics 106
4.14 Numerical results of Sh^g vs Pe^ for mass transfer over the percolation cluster of
Figure 4.13 and for three local exponents (c = 0, c = 1/3 and c = 1/2, from bottom
to top, respectively) 107
4.15 Numerical results of ShOTg vs PeL for mass transfer over the percolation cluster of
Figure 4.13 and for three local exponents (c = 0, c = 1/3 and c = Vz, from bottom
x
-------
to top, respectively) 108
4.16 Uniformly-distributed sources to simulate flow over a distributed source 109
4.17 Numerical results of Sh^ vs J>eL for mass transfer over a uniformly-distributed
source and for three local exponents (c = 0, c = 1/3 and c = 1/2, from bottom to
top, respectively) 110
5.1 Typical displacement geometries (from Yang and Yortsos, 1997) 124
5.2-a Concentration profiles for miscible displacement with a monotonic viscosity-
concentration profiles between parallel plates at time step = 0.3 and M = 100 for
various values pf the dipsersion coefficient N^. (from Yang and Yortsos, 1997) ... 125
5.2-b Concentration profiles for miscible displacement with a monotonic viscosity-
concentration profiles between parallel plates at time step = 0.3 and M = 0.01 for
various values pf the dipsersion coefficient NTO. (from Yang and Yortsos, 1997) ... 125
5.3 Monotonic viscosity-unfavorable mobility ratio (a= 5, NTO = 0, Time = 0.5) 126
5.4 Non-monotonic viscosity- unfavorable mobility ratio (a= 5, nm = 10, cm= 0.5, NTO
= 0, Time = 0.5) 127
5.5 Monotonic viscosity. Effect of dispersion (a= 5, NTO=0.01, Time=0.5) 128
5.6 Non-monotonic viscosity. Effect of dispersion (a= 5, nm = 10, cm = 0.5, NTO =
0.01, Time = 0.5) 129
5.7 Monotonic viscosity- favorable mobility ratio (a= 0.2, NTO = 0, Time = 0.5) 130
5.8 Non-monotonic viscosity- favorable mobility ratio (a = 0.2, jum= 10, cm = 0.5,
NTO=0.01, Time=0.5) 131
5.9 Monotonic viscosity- favorable mobility ratio. Effect of dispersion (a = 0.2, NTO =
0.01, Time=0.5) 132
5.10 Non-monotonic viscosity- favorable mobility ratio. Effect of dispersion (a = 0.2,
pm= 10, cm = 0.5, Nro=0.01, Time = 0.5) 133
5.11 Monotonic viscosity- Poiseuille flow (a= 1, NTO = 0, Time=0.4) 134
5.12 Non-monotonic viscosity at Poiseuille flow conditions (a= 1, jum = 10, cm = 0.5,
NTO = 0, Time = 0.4) 135
5.13 Effect oft on non-monotonic viscosity for a= 1 (jum = 7, cm = 0.5, NTO = 0) 136
5.14 Monotonic viscosity- Poiseuille flow. Effect of dispersion (a = 1, NTO = 0.01,.
Time = 0.4) 137
5.15 Non-monotonic viscosity at Poiseuille flow conditions. Effect of dispersion (a = 1,
jum= 10, cm = 0.5,NTO=0.01, Time=0.4) 138
5.16 Non-monotonic viscosity- effect of jum (a= 1, //„= 8, cm = 0.5, NTO = 0 ) 139
5.17 Non-monotonic viscosity- effect of jum (a= 1, jum= 20, cm = 0.5, NTO = 0) 140
5.18 Non-monotonic viscosity- effect of cm (a= 1, jum = 10, cm = 0.05, NTO = 0) 141
5.19 Non-monotonic viscosity- effect of cm (a= 1, jum= 10, cm = 0.9, N^ 0) 142
5.20 Schematic of the Experimental Apparatus Setup 143
5.21 Viscosity-concentration profile for a 2-propanol/water mixture 144
5.22 Viscosity-concentration profile for a glycerol/water mixture 145
5.23 Typical image of water displacing 1-propanol 146
5.24 Typical image of water displacing glycerol 147
XI
-------
5.25 Close-up of water displacing 2-propanol 148
5.26 Close-up of water displacing glycerol mixture 149
6.1 The dominant eigenvalue s^Qi) for a non-monotonic viscosity profile with M> 1
(curve a) and M< 1 (curve b) 156
6.2 The dispersion relation a(k~, e) for different values of the dimensionless tune l/4o2
as predicted from equation (6.15) for the conditions of case (a) of Figure 6.1 157
6.3 The dispersion relation a(fc, e) for different values of the dimensionless time l/4o2
predicted from equation (6.15) for the conditions of case (b) of Figure 6.1 158
xn
-------
LIST OF TABLES
3.1 A compilation of mass transfer correlations used in NAPL remediation 58
Xlll
-------
-------
Chapter 1
INTRODUCTION
A commonly used method for in-situ remediation is purnp-and-treat, in which the groundwater
is pumped out of the site, treated on the surface to remove the contaminant and subsequently
reinjected (Mackay and Cherry, 1989). Because of the relatively low solubility of the contaminants,
compounded with mass transfer limitations at the pore scale and from low permeability areas (such
as clay lenses), the overall process efficiency is much lower than desirable. For its enhancement,
various improvements are constantly being developed. The majority involve the injection of remedial
fluids or chemicals which may enhance the solubility of the contaminant (Pennell et al., 1993,
Augustijn et al., 1994) or help to mobilize trapped NAPL (Brandes and Farley, 1993). Injection
of surfactants (West and Harwell, 1992, Pennell et al., 1993), co-solvents, such as alcohols (Final
et al., 1990, Wood et al., 1990), and air (sparging) (Ji et al., 1993, Pankow et al., 1993) in the
groundwater, or injection of steam in the vadose zone (Itamura and Udell, 1993, Ho et al., 1994),
are some of the methods actively under consideration.
A key characteristic of in-situ remediation is the requirement of a high recovery efficiency over
relatively small areas (Mercer and Cohen, 1990, Mayer and Miller, 1993). This is qualitatively
different from enhanced oil recovery, to which it otherwise has many similarities, where recovery
cost is the basic design consideration. Three processes, generic to in-situ groundwater remediation,
require increased scrutiny:
• The emplacement in-situ of the contaminant NAPL;
• The mass transfer of the contaminant from the NAPL to the groundwater or from the ground-
water to a sparging fluid, and of the remedial agents to the NAPL; and
• The development of miscibility and possible microscale flow instabilities.
As with other processes involving porous media, in-situ remediation must be investigated at the
different scales that comprise porous media (Yortsos, 1998). These include
• The pore surface scale (of the order of nm), where adsorption, partitioning between phases,
thermodynamic equilibrium, surface chemistry and wettability are important;
• The pore scale (of the order of tens of /j,m), also called the microscopic scale, where pores can
be defined, fluid-fluid interfaces are confined and capillary phenomena are important;
• The pore-network scale (of the order of mm), where the collective behavior of a process over
an ensemble of pores is investigated;
-------
• The laboratory scale (core scale) (of the order of cm), also called the macroscopic scale, where
laboratory experiments are conducted; and
• The field scale (of the order of m and higher), also called the megascopic scale, which is of
relevance to large-scale spreading.
Phenomena at the different scales have a different manifestation, and input from a lower scale is
necessary for their description at the next higher scale.
With the exception of some recent works, the design of in-situ remediation is currently based on
a description at the macroscopic (and megascopic) scales using computer simulators. Phenomena
at the pore and pore-network scales are typically lumped in terms of averaged quantities, using
empirical or ad hoc expressions. At best, these are adaptations of macroscopic models developed
in enhanced oil recovery (e.g. see Abriola and Pinder, 1985a,b, Brown et al., 1994), and although
a significant improvement over earlied empirical models (e.g. Wilson et al., 1992), they cannot
address fundamental remediation issues that originate at the pore and pore network scales. The
following examples illustrate the inadequacy of the present state of the art.
1. Many aspects of remediation involve immiscible displacements, e.g. the downwards seepage of
NAPL into the subsurface, or the upwards injection of air in air sparging. Displacement patterns,
critically affected by pore-scale heterogeneities, capillary, viscous and gravity forces, dictate the
dynamic or equilibrium configurations of the NA Pi-ground water interface or the air-groundwater
interface (in air sparging), which in turn decisively determine mass transfer. Conventional macro-
scopic simulation based on averages, is inherently unable to provide a detailed description. For
mass-transfer purposes, the NAPL is typically assumed in the form of isolated spherical globules
in an structureless porous medium (Powers et al., 1994b), an assumption that contrasts commonly
accepted displacement patterns, which are often disordered or fractal (Feder, 1988, Lenormand,
1990, Sahimi and Yortsos, 1990). In gas sparging, patterns vary from a fluidized-like regime to a
disordered, channeled flow, depending on the rate of injection and the disorder of the medium (Ji et
al., 1993). Current macroscopic models ignore or lump together these microscale phenomena, with
the consequence that capillary or viscous fingering and channeling at the local scale are grossly
mis-represented. Recent approaches (Celia et al., 1993, Soil and Celia, 1993, Soil et al., 1993, ),
paralleling earlier efforts in oil recovery (Heiba et al., 1982, Parlar and Yortsos, 1988, 1989, Yortsos,
1990), have used Ordinary Percolation theory to relate quantities, such as relative permeabilities
and capillary pressures, to the microstructure. These need to be further advanced by considering
dynamic invasion models that include gravity, viscous forces and correlated heterogeneity.
2. The overall mass transfer coefficient has been estimated by widely different correlations, re-
flecting the different assumptions made by different investigators. Vigorously debated, but presently
inconclusive, are effects of flow velocity, interfacial area and its variation with time (due to dissolu-
tion, mobilization, etc.) (Imhoff et al., 1993). Particularly interesting are the reported differences in
the exponent relating the mass-transfer Sherwood number to the Reynolds or the Peclet numbers,
for which convincing theories have not been provided. This exponent should reflect the relevant
flow regime and the interface geometry (Cussler, 1984). In the most advanced model todate, Powers
et al. (1994) approximated mass transfer by steady-state diffusion from an isolated spherical NAPL
globule (which is akin to a dilute-limit approximation in a Hele-Shaw cell, Li and Yortsos, 1994).
However, such approximations do not adequately represent the true geometrical configurations of
JVj4 Pi-groundwater interfaces in a real system and the concommitant effects on fluid flow and mass
transfer. Similarly poor is the estimation of mass transfer rates from the air-groundwater interface
in air sparging.
-------
3. The development of miscibility between an injected solvent and in-situ NAPLs has also been
addressed only by macroscopic models. The most advanced version (Brown et al., 1994) assumes
equilibrium phase behavior for the average compositions and saturations, similar to analogous mod-
els for surfactant flooding (Pope and Nelson, 1978) or COi flooding (Stalkup, 1983) in oil recovery.
Implied to such equilibrium and first-contact miscibility at the macroscale is the assumption of fast
mass transfer (by diffusion and convection) to lead to spatially uniform profiles. This may not
be generally true in the presence of multiple fluid-fluid interfaces at the the pore scale, however,
where mass transfer will be hindered and the interplay of transport, phase equilibrium and flow is
much more complex than represented in macroscopic models. Local mixing can be accompanied
by changes in fluid viscosity and density, and important viscous or boyancy-driven instabilities at
the pore scale can be induced (Yang and Yortsos, 1997). These are not recognized in macroscopic
simulators. In fact, the current models actually employ a formalism of immiscible displacements,
where miscibility is indirectly accounted for through the dependence of the residual saturation, 5or,
on the capillary number, Ca (in turn indirectly related to concentration). The Sor(Ca) function
used corresponds to a sequence of equilibrium states reached in displacements at a constant Ca,
which may be questionable when the development of miscibility is dynamic. Finally, the REVovei
which averages are defined (Bear, 1972), becomes ill-defined when there exist gradients in average
concentrations or saturations, or when phenomena at the pore-scale are dynamic.
In recent studies, we have addressed some of these issues, although in somewhat different con-
texts. In a series of papers (Li and Yortsos, 1994, 1995a,b, Satik et al., 1995) we have investigated
the pore scale mass transfer towards a growing gas phase from a surrounding liquid in a porous
medium, in the context of solution gas-drive. We have shown by experiments in etched-glass mi-
cromodels and in Hele-Shaw cells, and by simulation using pore network models that the interface
pattern varies from Invasion Percolation (IP) to Viscous Fingering ( VF) (Feder, 1988) and drasti-
cally affects the mass transfer, hence the rate of growth (or dissolution) of a phase. The disorder of
a real porous medium (even if macroscopically homogeneous) was shown to give rise to qualitatively
different phenomena than in effective porous media (such as a Hele-Shaw cell, Li and Yortsos, 1994,
1995a). To our knowledge, these are the first studies, where pore level simulation has been used to
understand mass transfer in multiphase flows in porous media. In another recent study (Yang and
Yortsos, 1997), we have described the miscible displacement in a single capillary or in the gap of
a Hele-Shaw cell using Stokes (rather than Darcy) equations. We have shown that, if the injected
fluid is of a lower viscosity, viscous fingers involving a sharp interface between injected and initial
fluids within a single pore will develop, provided that the transverse Peclet number, which measures
convection over diffusion, is sufficiently large. Similar effects are likely to occur in the context of
remediation as well.
To better understand and to help design improved remediation methods, we present in this
report a pore-scale study of some of the generic aspects of remediation cited above, namely the
emplacement in-situ of the NAPL, the mass transfer characteristics of remediation processes, and
the development of microscale flow instabilities during remedial fluid injection. Experimental, the-
oretical and computational studies are conducted, with emphasis placed at the pore network scale,
although appropriate efforts are also devoted to the single pore scale, whenever appropriate. The
information obtained will be incorporated by up-scaling techniques into macroscopic simulators.
The specific tasks are organized as follows:
Chapter 2 describes the use of pore-network simulation for the characterization of immiscible
displacement patterns expected during NAPL emplacement. We briefly present a pore-netwrok
approach to study the various immiscible displacement patterns that develop between a percolating
NAPL and the groundwater, or between an immiscible remedial fluid and the groundwater. The
pore network simulation includes capillary, gravity and viscous forces and also accounts for pore
-------
scale heterogeneity.
Chapters 3 and 4 describe our study of mass transfer from or to a trapped, stationary NAPL,
during a remediation process. Here, the main objective is the determination of the mass transfer
coefficient as a function of the flow rate, the interfacial area and other physical and geometrical
parameters. Dissolution experiments in etched-glass micromodels are conducted. Then, based on
the particular configuration of the trapped NAPL, a pore-network simulation is carried out of mass
transfer by diffusion and convection in the disordered pore-space occupied by the injected fluid.
Mass tranfer is computed by a pore network simulator, as described below, with local mass-transfer
coefficients determined from single-pore expressions. At sufficiently slow dissolution rates, a quasi-
steady state is likely. Then, an overall mass transfer coefficient is computed and its sensitivity
to the Peclet number and various other properties of the system is determined. Chapter 3 deals
with the etched-glass micromodel experiments, while Chapter 4 presents a general theory of mass
transfer from stationary sources in a porous medium.
The third topic to be studied is the various viscous instabilities during miscible displacements.
Chapters 5 and 6 describe aspects of instability at the single-pore scale (for example in the gap of
a Hele-Shaw cell, which is equivalent to a single capillary) and also at the macroscopic scale for the
case of non-monotonic viscosity behavior, respectively. Some experiments in a Hele-Shaw cell are
also described in Chapter 5. Because viscosities and densities are composition-dependent, spatial
variations in concentration may lead to viscous instabilities or otherwise affect the flow pattern.
The studies in a single pore are similar to Yang and Yortsos (1997) for monotonic displacements.
The results to be obtained can be incorporated into a macroscopic description by applying an
up-scaling technique (Chapter 4). Included among them are homogenization (Bensoussan et al.,
1978), volume-averaging (Quintard and Whitaker, 1994), ensemble averaging (Gelhar and Axness,
1983, Koch and Brady, 1987), and gradient percolation (Gouyet et al., 1988). Depending on whether
scales can be separated or not, these approaches either lead to the formulation of a local boundary
value problem, the solution of which provides the desired coefficients (case of homogenization and
volume averaging), or may lead to an altogether different macroscopic model. When gradients in
the average properties are not large, the application proceeds in a standard fashion, as described in
Quintard and Whitaker (1994). However, problems will arise at larger rates, where concentration
gradients are large, in which case the basic premise of separation of scales breaks down. Asymptotic
descriptions and scaling theories for displacement patterns and saturation profiles are described in
Chapter 2. Corresponding results for the mass transfer coefficients are provided in Chapter 4.
The simulations on mass transfer are compared against experimental data in etched-glass mi-
cromodels. The comparison between theory and experiment will allow a check of the relevance of
the pore network description and to seek corrections, particularly at the single-pore level, if there
exists a mismatch.
REFERENCES
1. Abriola, L.M. and Finder, G.F., A multiphase approach to the modeling of porous media
contamination by organic compounds: 1. Equation Development, Water Resour. Res., Vol.
21, 11-18 (1985a).
2. Abriola, L.M. and Finder, G.F., A multiphase approach to the modeling of porous media
contamination by organic compounds: 2. Numerical simulation, Water Resour. Res., Vol.
21, 19-26 (1985b).
-------
3. Augustijn, D.C.M., Jessup, R.E., Rao, P.S.C. and Wood, A.L., Remediation of contaminated
soils by solvent flushing, J. Envir. Engin., Vol. 120, 42-57 (1994).
4. Bear, J., "Dynamics of fluids in porous media", Elsevier, New York, NY (1972).
5. Bensoussan, A., Lions, J.L. and Papanicolaou, G., "Asymptotic analysis for periodic struc-
tures", North-Holland, Amsterdam (1978).
6. Brandes, D. and Farley, K.J., Importance of phase-behavior on the removal of residual
DNAPLs from porous media by alcohol flooding, Water Envir. Res., Vol. 65, 869-878 (1993).
7. Brown, C.L., Pope, G.A., Abriola, L.M. and Sepehrnoori, K., Simulation of surfactant-
enhanced aquifer remediation, Water Resour. Res., Vol. 30, 2959-2977, (1994).
8. Celia, M. A., Rajaram, H. and Ferrand, L.A., A multi-scale computational model for multi-
phase flow in porous media, Advances in Water Resources, Vol. 16, 81-92 (1993).
9. Cussler, E.L., "Diffusion: Mass transfer in fluid systems", Cambridge University Press (1984).
10. Feder, J., "Fractals", Plenum, New York (1988).
11. Gelhar, L.W. and Axness, C.L., Three-dimensional stochastic analysis of macrodispersion in
aquifers, Water Resour. Res., Vol. 19, 161-180 (1983).
12. Gouyet, J.F., Rosso, M. and Sapoval, B., Fractal structure of diffusion and invasion fronts in
three-dimensional lattices through the gradient percolation approach, Phys. Rev. B, Vol. 37,
1832 (1988).
13. Heiba, A. A., Sahimi, M., Scriven, L.E. and Davis, H.T., Percolation theory of two-phase
relative permeabilities, Paper SPE 11015, presented at the 57th SPE Annual Meeting, New
Orleans, LA, Sept., 26-29, 1982.
14. Ho, O.K., Liu, S. and Udell, K.S., Propagation of evaporation and condensation fronts during
multicomponent soil vapor extraction, J. Contam. Hydr., Vol. 16, 381-401 (1994).
15. Imhoff, P.T., Jaffe, P.R. and Pinder, G.F., An experimental study of complete dissolution of
a nonaqueous phase liquid in saturated porous media, Water Resour. Res., Vol. 30, 307-320
(1993).
16. Itamura, M.T. and Udell, K.S., Experimental clean-up of a dense non-aqueous phase liquid
in the unsaturated zone of a porous medium using steam injection, Vol. 265, Multiphase
Transport in Porous Media, ASME, 57-62 (1993).
17. Ji, W., Dahmani, A., Ahlfeld, D.P., Lin, J.D. and Hill III, E., Laboratory study of air sparging:
Air flow visualization, GWMR, 115-126 (1993).
18. Koch, D.L. and Brady, J.F., A non-local description of advection-diifusion with application
to dispersion in porous media, J. Fluid Mech. Vol. 180, 387-403 (1987).
19. Lenormand, R., Liquids in porous media, J. Phys.: Condens. Matter, Vol. 2, SA79-SA88
(1990).
20. Li, X. and Yortsos, Y.C., Bubble growth and stability in an effective porous medium, Phys.
Fluids A, Vol. 6, 1663-1676 (1994).
-------
21. Li, X. and Yortsos, Y.C., Visualization and simulation of bubble growth in pore networks,
AIChE J., Vol. 41, 214-223 (1995a).
22. Li, X. and Yortsos, Y.C., Bubble growth in porous media, Chem. Eng. Sci., in press (1995b).
23. Mackay, D.M. and Cherry, J.A., Groundwater contamination: Pump-and-treat remediation,
Environ. Sci. Technol., Vol. 23, 630-636 (1989).
24. Mayer, A.S. and Miller, C.T., An experimental investigation of pore-scale distributions of
nonaqueous phase liquids at residual saturation, Transport in Porous Media, Vol. 10, 57-80
(1993).
25. Mercer, J.W. and Cohen, R.M., A review of immiscible fluids in the subsurface: Properties,
models, characterization and remediation, J. Contam. Hydr., Vol. 6, 107-163 (1990).
26. Pankow, J.F., Johnson, R.L. and Cherry, J.A., Air sparging in gate wells in cutoff walls and
trenches for control of plumes of volatile organic compounds (VOCs), Ground Water, Vol.
31, 654-663 (1993).
27. Parlar, M. and Yortsos, Y.C., Percolation theory of vapor adsorption-desorption processes in
porous media, J. Colloid Interface Sci., Vol. 124, 162-176 (1988).
28. Parlar, M. and Yortsos, Y.C., Nucleation and pore geometry effects on capillary desorption
in porous media, J. Colloid Interface Sci., Vol. 132, 425-443 (1989)_.
29. Pennell, K.D., Abriola, L.M. and Weber, W.J., Jr., Surfactant-enhanced solubilization of
residual dodecane in soil columns. 1. Experimental investigation, Environ. Sci. Technol.,
Vol. 27, 2332-2340 (1993).
30. Final, R., Rao, P.S.C., Lee, L.S., Cline, P.V., and Yalkowski, S.H., Cosolvency of partially
miscible organic solvents on the solubility of hydrophobic organic chemicals, Environ. Sci.
Technol., Vol. 24, 639-647 (1990).
31. Pope, G.A. and Nelson, R.C., A chemical flooding compositional simulator, Soc. Pet. Eng.
J., Vol. 18, 339-354 (1978).
32. Powers, S.E., Abriola, L.M. and Weber, W.J., Jr., An experimental investigation of nonaque-
ous phase liquid dissolution in saturated subsurface systems: Steady state mass transfer rates,
Water Resour. Res., Vol. 28, 2691-2705 (1992).
33. Powers, S.E., Abriola, L.M. and Weber, W.J., Jr., An experimental investigation of nonaque-
ous phase liquid dissolution in saturated subsurface systems: Transient mass transfer rates,
Water Resour. Res., Vol. 30, 321-332, Feb. (1994).
34. Powers, S.E., Abriola, L.M., Dunkin, J.S. and Weber, W.J., Jr., Phenomenological models
for transient JVAPi-water mass-transfer processes, J. Contam. Hydr., Vol. 16, 1-33 (1994).
35. Quintard, M. and Whitaker, S., Convection, dispersion, and interfacial transport of contam-
inants: Homogeneous porous media, Advances in Water Resources, Vol. 17, 221-239 (1994).
36. Sahimi, M. and Yortsos, Y.C., Application of fractal geometry to porous media: A review,
paper SPE 20476 presented at the 65th SPE Annual Fall Meeting, Dallas, TX (Oct. 6-9,
1990).
-------
37. Satik, C., Li, X. and Yortsos, Y.C., Scaling of bubble growth in porous media, Phys. Rev. E,
in press (1995).
38. Soil, W.E., Celia, M.A. and Wilson J.L., Micromodel studies of three-fluid porous media
systems: Pore-scale processes relating to capillary pressure-saturation relationships, Water
Resour. Res., Vol. 29, 2963-2974 (1993).
39. Soil, W.E. and Celia, M.A., A modified percolation approach to simulating three-fluid capil-
lary pressure-saturation relationships, Advances in Water Resources, Vol. 16, 107-126 (1993).
40. Stalkup, F.I., Jr., "Miscible Displacement", SPE Monograph, Vol. 8, SPE, New York (1983).
41. West, C.C. and Harwell, J.H., Surfactants and subsurface remediation, Environ. Sci. Tech-
nol., Vol. 26, 2324-2330 (1992).
42. Wilson, D.J., Kayano, S., Mutch, R.D., Jr. and Clarke, A.N., Groundwater cleanup by in-situ
sparging 1. Mathematical modeling, Separation Science and Technology, Vol. 27, 1023-1041
(1992).
43. Wood, A.L., Bouchard, D.C., Brusseau, M.L. and Rao, P.S.C., Cosolvent effect on sorption
and mobility of organic contaminants in soils, Chemosphere, Vol. 21, 575-587 (1990).
44. Yang, Z. and Yortsos, Y.C., Asymptotic solutions of miscible displacement in geometries of
large aspect ratio, Phys Fluids, Vol. 9, 286-298 (1997).
45. Yortsos, Y.C., "Discrete Approach: Percolation Theory, Instabilities: Miscible and Immiscible
Flows, Heterogeneity Description Using Fractal Concepts, Reaction and Transport in Porous
Media", von Karman Institute Lecture Series on Modeling and Applications of Transport
Phenomena in Porous Media, Brussels, Belgium (Feb. 5-Feb. 9, 1990).
46. Yortsos, Y.C., "Fluid Flow and Transport Processes in Porous Media", Springer Verlag (in
progress) (1998).
-------
Chapter 2
DISPLACEMENT PATTERNS
USING PORE-NETWORK MODELS
Baomin Xu, Yuyong Zhang and Yanis C. Yortsos
INTRODUCTION
The determination of immiscible displacement patterns that develop during seepage of NAPLs
in the groundwater or during air sparging of contaminated groundwater is important. Displacement
patterns are critically affected by pore-scale heterogeneities, and by capillary, viscous and gravity
forces, which dictate dynamic or equilibrium configurations of the TV^lPi-groundwater interface
or the air-groundwater interface (in air sparging). In turn, interfacial area and geometry are key
variables to mass transfer in remediation processes. In this chapter we consider a pore-network
study of some generic aspects involved in the emplacement in-situ of the NAPL and displacement
patterns in air sparging.
To simulate displacement patterns in a pore network, we can use either a generally complex pore
network simulator, or a simpler Invasion Percolation in a Gradient (IPG) model, to be described
below, in which gravity, capillarity and viscous forces dictate the displacement. Although simpler,
IPG is not less realistic, and rigorously applies to many problems as will be shown. In either
case, pore-network simulation or simpler statistical physics models conditions can lead to pore-
level configurations for the interface between invading and initial fluid. For example, the spreading
of a contaminant phase over a lens of low permeability, etc., can be readily modeled. An interesting
configuration is that corresponding to an Invasion Percolation IP cluster, where the interface has
the fractal structure of the percolation cluster. Such interfaces are highly disordered, they are
ubiquitous in displacement processes, even in relatively homogeneous systems, as will be shown
below, and cannot be replicated by continuum models. In this chapter we provide a pore-network
level study of the particular patterns expected to develop in immiscible displacements under various
conditions, and will examine the sensitivity to parameters, such as the capillary and Bond numbers,
expressing relative magnitudes of viscous, gravity and capillary forces, and to the disorder of the
medium. The interface configuration obtained will be then used as input for solving the asscociated
mass transfer problem in the next tasks.
This chapter is organized as follows: We first present a new approach for understanding drainage
(namely the displacement of a wetting by a non-wetting phase in a porous medium) in the presence
of viscous forces, based on Invasion Percolation in a Gradient (IPG). This method is detailed
in Xu et al. (1998) and allows the classification of the various patterns expected. Then, we
-------
describe conventional pore network simulations and provide examples of patterns obtained based
on the work of Xu (1995). The effect of gravity on the spreading of DNAPL over an impermeable
boundary is described in a following section, where the displacement characteristics are analyzed
for various parameter values. We close with a description os air sparging patterns obtained with the
use of IPG. Before we proceed, we note that the application of pore-network models to describe
processes in real porous media is a rapidly growing research area. Xu et al. (1997) describe in detail
the specific approach they have followed for simulating laboratory experiments in heterogeneous
carbonate samples.
DRAINAGE WITH VISCOUS FORCES
The displacement (drainage) of a wetting fluid (subscript w) in a porous medium by the injection
of a non-wetting fluid (subscript mo), immiscible to the former, has been analysed in great detail
in past studies. In the absence of viscous or gravity forces, slow drainage is controlled solely by
the capillary pressure, Pc — Pnw - Pw (the difference in pressure between the two fluids), which is
spatially uniform. At the pore network level, this problem can be modeled by Invasion Percolation,
in which the front separating the two fluids advances by penetrating the pore throat at the front
with the largest size (smallest capillary resistance). The properties of IP and its close connection
to Ordinary Percolation (OP) have been extensively studied (Wilkinson and Willemsen, 1983).
In the presence of gravity (Gouyet et al. 1988), or of a gradient in the average pore-size (hence,
of a gradient in the permeability) (Chaouche et al., 1994a), and in the absence of viscous forces, slow
drainage has been modeled with Invasion Percolation in a Gradient, which is a modified version
of Gradient Percolation (GP). Here the capillary pressure varies linearly (or almost linearly)
in the direction of displacement, x. Because of their direct relationship (see also below), this
gradient also results in a gradient in the percolation probability p, usually expressed in terms of
the Bond number, B (where B ~ -dp/dx). For example, for invasion under a hydrostatic gradient,
B = P3*r™, where Ap is the density difference, gx is the gravity component in the direction of
displacement, rm is a typical throat size and 7 is the interfacial tension between the two fluids. For
invasion in a permeability gradient, B = -£^-, where k is the permeability. In IP (5 = 0), the
entire displacement pattern is a percolation cluster (see Xu (1995) for further discussion of these
patterns). However, in IPG, one needs to distinguish between two different cases, depending on
whether percolation is in a stabilizing or a destabilizing gradient.
In the first case (B > 0), for example in the downwards displacement at capillary control of a
heavier fluid by the injection of a lighter fluid, or in drainage in afield of decreasing permeability, the
percolation probability decreases in the direction of displacement. The region where the invasion
has the characteristics of a percolation cluster is only of a finite extent, a, which was shown by
Gouyet et al. (1988) to scale as
a ~ B~^ (2.1)
where v is the correlation length exponent of percolation (Stauffer and Aharony, 1992). Here a
denotes the width of the front (in 2-D) or of the front-tail (in 3-D), where the displacing pattern has
the structure of the percolation cluster and fractal concepts apply. Equivalently, a measures the
maximum extent of the correlation length, which in gradient percolation problems becomes finite
due to the applied gradient. Various properties of GP and IPG have been studied in considerable
detail. During invasion in a destabilizing gradient (B < 0), the percolation probability increases
in the direction of displacement. Then, the displacement proceeds in the form of capillary fingers,
-------
the scaling of the average thickness of which with (the absolute value of) the Bond number also
satisfies (2.1) (see Meakin et al., 1992). Typical examples include the downwards displacement at
capillary control of a lighter fluid by the injection of a heavier fluid, as in DNAPL movement, the
upwards displacement of a heavier fluid by the injection of a lighter fluid, as in air sparging, or
drainage in a field of increasing permeability,
In the presence of both viscous and capillary forces, the displacement is characterized by three
dimensionless numbers: the capillary number, Ca = qpnwhi where q is the injection velocity and
/inw the viscosity of the displacing phase; the viscosity ratio, M — Hw/Hnw, where /j,w is the viscosity
of the displaced phase; and the dimensionless system size L (expressed in units of the average pore
length /). As Co. or L increase, the viscous pressure drop in the two fluids becomes comparable to
capillarity, and one expects that some form of gradient percolation would also describe this process.
The effect of viscous forces on displacements in porous media is of obvious importance to pro-
cess scale-up and large-scale simulation. As length scales increase, viscous effects are increasingly
dominant over capillarity. An understanding of this competition at the pore-network scale is nec-
essary to provide insight on the validity of the conventional continuum description using relative
permeabilities, and to delineate the particular patterns obtained. Here, we consider fully devel-
oped drainage in uncorrelated random media in the presence of viscous forces and proceed by
postulating an analogy with IPG. Depending on the relative magnitude of M, we anticipate the
existence of two different regimes, described by IPG in a stabilizing or a destabilizing gradient,
respectively. These two regimes dictate the development of the saturation profiles in the respective
displacements, as will be shown below.
In the absence of viscous forces, the capillary pressure is spatially uniform and the displacement
proceeds by following the usual rules of IP, namely by successively invading the perimeter pore
with the largest size. In the presence of a viscous pressure drop, a gradient in the capillary pressure
(negative or positive) is generally expected to develop. In view of the relations
and
P-22
* r. —
/oo
a(r)dr
(2.2)
(2.3)
(or, more correctly, p = /r°°. a(r)dr, where rmj-n is the minimum throat size invaded), this, in
turn, implies a gradient in the percolation probability. In the above, a(r] is the probability density
function of the pore size distribution. Problems involving a constant gradient in p are amenable to
GP and IPG, thus we expect that a similar description would also be applicable in the present case
involving viscous forces. Figure 2.1 shows typical patterns of viscous displacements for M = 0.1,
obtained from pore network simulation in the absence of trapping (details of the simulation can be
found in Xu, 1995, and will also be described below). When Ca is low (Figure 2.la), viscous forces
are negligible and the pattern has the fractal structure of the IP cluster. As viscous forces increase
at larger Ca (Figures 2.1b-2.1d), however, the front takes the appearance of a rough (self-affine),
rather than self-similar, curve and has an extent that decreases with increasing Ca. These trends
are consistent with a gradient percolation description.
Following GP notions, we will distinguish two different cases, one in which the percolation
probability p, (hence Pc) decreases in the direction of displacement, and which has features similar
to GP in a stabilizing gradient, and another in which p (hence Pc) increases in the direction of
displacement, with features similar to GP in a destabilizing gradient. Because in our problem, the
10
-------
capillary pressure Pc is controlled by the viscous pressure drop (rather than gravity), the viscosity
ratio M is expected to be an important parameter in delineating these two regimes.
Stabilized Displacement
In the case of a stabilized displacement, we follow Gouyet et al. (1988) and define the location
Xc(t], as the place where the transverse average of the percolation probability is equal to the
percolation threshold, pc, namely
p(Xc) = pc (2.4)
In 2-D lattices, Xc represents the mean front position. The regions on either side of Xc of an extent
a (namely between Xc — a and Xc+cr) have the fractal properties of the percolation cluster (Gouyet
et al., 1988). In 3-D lattices, however, due to the higher connectivity, Xc does not represent the
mean front position (here the front extends far upstream), but rather denotes a mean leading edge.
Nonetheless, a percolation pattern is also expected around Xc. We will focus on the front-tail region
(X > Xc), the extent of which we will also denote by a. As in the corresponding IPG problem,
Xc(t] varies linearly with time in either 2-D or 3-D geometries, with a velocity u, to be determined.
In both cases, the fractal regions are followed by an upstream region, where both invading and
invaded phases are compact and the conventional continuum description is valid. A schematic is
shown in Figure 2.2.
After an in-depth analysis, Xu et al. (1998) (see also Yortsos et al., 1997) showed that the
frontal region of width a satisfies the following scaling with the front capillary number Cap and
the standard deviation of the pore-size distribution E,
(7
2s; (2.5)
where exponent £ is the conductivity exponent of percolation and D denotes the fractal dimension
of the percolation cluster. Equation (2.5) expresses the asymptotic scaling of the front (or front-tail)
width with Cap at small Cap for the case of a stabilized displacement.
Using accepted estimates, the exponent in (2.5) is equal to 0.382 or 0.25 in 2-D or 3-D, respec-
tively. The rather small values suggest a weak sensitivity of a on Cap and S. A plot of a as a
function of the capillary number and the mobility ratio, M, is shown in Figure 2.3 (see also Xu
et al., 1998). It should be kept in mind that the rate dependence in Cap has entered through
the front velocity, v. This can be of some significance in investigating the effect of M, which also
influences v. It is also interesting to note that as the degree of heterogeneity, E, increases the front
width increases. In fact, Equation (2.5) suggests that for such displacements it is more appropriate
to replace, in the definition of the capillary number, 7 by the product 7E. The scaling (2.5) can
be also obtained by using a version of GP, termed Viscous Gradient Percolation (VGP), to be
introduced below for modeling the saturation profile. Although necessary for the saturation profile,
VGP is not required for the derivation or the validity of the scaling (2.5).
11
------- An error occurred while trying to OCR this image.
-------
decreases almost linearly with distance in the region upstream of the front, and the front extent
increases as B decreases. (The different lattice sizes used in the normalized plots of these figures
should be noted). The scaling of the saturation profile is expected to have the general features
of the scaling function (2.6) of the standard GP, with some correction to account for the slightly
different VGP profile.
We summarize this section by emphasizing the different description of percolation processes
involving viscous forces in the two regions, near and away from the front, at least for relatively
low Co,. The different scalings, CaF1+c+"(D~1) and Ca~l (for the continuum case) obtained indicate
that near the front, the continuum description for the profile should be replaced with the more
appropriate VGP Equation (2.6). Either theory suggests an advancing front. In the VGP model,
the profile is a function of X - Xc, Cap and M, the time-dependence entering through Xc, which
varies with time. A travelling state with constant velocity is also contained in the continuum
description. However, the latter predicts a hypodiffusive behavior, namely a profile with a divergent
derivative at the front (a sharp "knee"), in contrast to the tail involved in VGP. Thus, appropriate
caution must be exercised in using the continuum approach in this region.
Capillary-Viscous Fingering
When the nw fluid has a much smaller viscosity (namely M » 1), most of the pressure drop
occurs within the displaced phase. In such cases p increases in the direction of displacement. IPG
problems with spatialy increasing percolation probability involve a negative Bond number and
describe invasion in a destabilizing gradient (Meakin et al., 1992) (see also below). In particular,
capillary invasion in a destabilizing gravity field corresponds to the release of a lighter fluid at the
bottom of a porous column filled with a heavier fluid, for example in air sparging. It was found
that the displacement occurs in the form of distinct capillary (but not DLA-type, namely only
viscous-dominated) fingers. For a sufficiently long column, one single finger emerged. Figure 2.5
reprinted from Chaouche et al. (1994a) shows the structure of such a finger for capillary invasion
in a field of increasing permeability. Scaling arguments similar to the self-consistency arguments of
GP can be used here to show that the finger consists of a string of beads of average width a with
the following scaling behavior
(J
\B\
(2.10)
To apply these findings to the viscous problems of interest here, an expression for B is needed. Xu
et al. (1998) have identified a suitable Bond number for this problem, namely
CaM
2S
(2.11)
Using this definition, a direct comparison with (2.10) gives the following result for the finger width
f T4.
(2.12)
\ '-"-> /
The power law of (2.12) has the familiar GP exponent, with values 0.571 or 0.469 in the respective
geometries (compared to 0.382 or 0.25 of the previous case). Comparison with (2.5) shows that the
scaling exponent almost doubles as the mobility ratio increases from the one limiting regime to the
other, implying a higher sensitivity on the capillary number. Equation (2.12) also shows that the
finger width decreases with an increase in the capillary number Ca and the mobility ratio M, and
13
-------
that it eventually reduces to a single thin finger of the size of a single pore (and where the above
scaling fails and a DLA regime emerges). This behavior is as expected.
We also note that equation (2.12) can be approximated rather well (at least in 3-D) with the
expression
a ~ Ca~°-5 (2.13)
This scaling is consistent with that of the fastest (most dangerous) growing finger predicted by the
linear stability analysis of Chuoke et al. (1959), which suggests an exponent equal to 1/2. However,
the two should not be confused. The present analysis is based on finger widths of the order of the
pore-scale, while the linear stability analysis of Chuoke et al. (1959) is based on a continuum,
large-scale description.
The analysis of the above sections shows that the displacement behavior is different depending
on the relative magnitude of M. This is a direct consequence of the pressure drop in the frontal
region in the two cases of low or high M, respectively. In either case, the pressure drop is associated
with the higher flow resistance. In the first case, it is due to the invading phase, which near the front
occupies a percolation-like cluster. In the second case, it is due to the displaced phase, which near
the front occupies a compact region. The different behavior in these two limits is the origin of the
difference in the scaling exponents in the two limits. These results are helpful in the understanding
of patterns and also in the delineation of conditions where fractal behavior is expected.
In summary, in this section, the effect of viscous forces on drainage displacements in porous
media was studied. We recognized that the process, at least near the front, shares common aspects
with IPG. When M is sufficiently small, the displacement can be modeled by a form of Gradient
Percolation in a stabilizing gradient. We developed the scaling of the front width and the saturation
profile, in terms of the capillary number. These results generalize the theory of Gouyet et al.
(1988). As the stabilized regime is described by the Buckley-Leverett equation, the two share the
same constraints for their validity. In the opposite case, the displacement is described by Gradient
Percolation in a destabilizing gradient and leads to fingering. The particular regime involves a
competition between capillary and viscous forces and was identified for the first time in the context
of viscous displacements. The theory shows that the conventional continuum approach should be
used with caution near the front.
PORE NETWORK SIMULATION
To complement the previous section, we describe here the elements of pore network simulation
for drainage processes. As pointed out, pore network simulators are used to simulate processes
in pore-networks or glass micromodels. A key assumption is the porous medium representation
in terms of a collection of pores, connected to each other by pore throats in a network-like fash-
ion. Conventional networks include regular lattices, such as square or cubic, or Voronoi lattices
(Blunt and King, 1991). Either lattice can be readily constructed, although regular lattices are
computationally more manageable. Pore networks are ideal for replicating experiments in etched-
glass micromodels (see Haghighi et al., 1994, and Li and Yortsos, 1995, for a recent application).
Whether for miscible or immiscible flow, their basic aspects include the following:
* Pores provide volumetric storage, throats control the conductance to flow, heat and mass.
Geometrical characteristics can be statistically distributed.
• Capillary equilibrium at the pore scale is controlled by pore throats during drainage (a non-
wetting (nw) fluid displacing a wetting (w) fluid), and by sites during imbibition (which is the
14
-------
inverse of drainage). The capillary pressure across an interface is given by the Laplace-Young
equation, Pnw — Pw — 2j"H, where 7 is the interfacial tension, and at equilibrium the mean
interface curvature, H, is controlled by the pore geometry.
• Pore bodies can be occupied by one or both of the two fluid phases. Pores occupied by one
phase can become trapped by the other phase, if topologically possible. When a pore is
trapped, liquid phase movement is not allowed, although mass diffusion continues.
• Local pore-scale coefficients are calculated from single-pore studies. Thermodynamic equilib-
rium applies within any given pore, although adjacent pores can be at non-equilibrium.
In the simulations, the governing equations are discretized over the pores (nodes) of the network.
For example, the overall mass balance at node i is Y^j PijQij = 0, where pij, Qij refer to density
and volume flow rate between adjacent sites i and j. The sum is over all neighboring sites j,
and we may use Poiseuille's law, Q^ = ^-(Aff, - pijQij), where the overall fluid conductance,
'
Gij = %[.'}, can be a distributed variable. Here, p, is viscosity, Aij denotes the cross-sectional area
of the pore throat joining sites i and j, r is the throat radius, I is the bond length and g is the
corresponding component of gravity. Subscript ij indicates quantities pertaining to the throat ij.
Fluid properties, such as viscosity, density and interfacial tension can be allowed to vary in different
pores, due to compositional dependence.
In the past, simulators have been developed to model immiscible displacement (Lenormand,
1990, Blunt and King, 1991, Chaouche et al., 1994b, Haghighi et al., 1994). These simulators
provide a pore-by-pore account of the displacement of interfaces during drainage or imbibition, and
can be used to study effects of the viscosity ratio, M, the capillary number, Ca = ^n, the gravity
Bond number, Bg — ^^, and various heterogeneities (including spatial correlations at the pore
scale). Here, q is the flow velocity, k is the medium permeability and Ap the density difference.
As pointed out above, in certain limits, the displacement patterns can be well characterized by
simpler statistical physics models: When Ca 1, or to a
compact (piston-like) displacement, if M < 1. Simulation of the statistical models IP, IPG and
DLA is computationally simple and it is often used to represent displacement patterns. This was
demonstrated in the previous section using VGP and will be shown again in the next section. For
more realistic simulations, the full pore network simulator can be used to explore displacement
patterns under general conditions.
In the following we will illustrate applications of pore network simulations for the description
of displacement patterns in heterogeneous porous media. We present simulation results, using the
full pore network simulator of Xu (1995), carried out in 2-D square lattices of size 100 X 21 and
100 X 60, and for various parameter values. Because of the connection between permeability and
average pore size, the simulation of permeability heterogeneity can be accomplished by a variation
of the pore size. The latter included a sudden increase or decrease in the average pore size in the
direction of displacement by a factor of 2 or 4. Pore.sizes were randomly and uniformly distributed
around their mean value. We shall present results with M — 1. More details can be found in Xu
(1995).
First, we show results for drainage in a homogeneous medium (Figure 2.6). The displacement has
the typical features of a unit mobility ratio displacement. We note that these features are similar to
those obtained in the previous section using the simpler VGP model. A certain amount of capillary
15
-------
trapping occurs at this capillary number. It is important to note the significant fluctuations in the
saturation profile due to the finite size of the pore network, however. These effects must be kept in
mind, when assessing the subsequent simulations in heterogeneous systems. Numerical simulations
for a pore size ramp increase of five lattice spaces are shown in Figures 2.7 and 2.8 for two different
values of the capillary number. A response similar to the one predicted by the continuum model
and also observed in experiments (Chaouche et al., 1994b) is evident, namely the reduction of
the non-wetting phase saturation as the region of permeability increase is approached, followed by
a sharp increase. The downward kink diminishes as the capillary number increases (Figure 2.8),
although not as smoothly as in the experiments, mostly we suspect as a result of finite size effects.
Also, because of the overlap with the boundary end effect at the outlet of the lattice, the saturation
following the jump does not rise to the high level observed in either experiments or continuum
simulations. A better agreement was found in the simulations involving a larger lattice, however.
In all cases, the response diminished substantially as the capillary number increased to levels higher
than 10~3.
A comparison between experiments and pore network simulations is given in Chaouche et al.
(I994b). Both experiments and pore-network models confirmed the validity of the macroscopic
predictions, at least as far as the basic features of the saturation response is concerned. An analysis
of the results based on percolation arguments provided additional support to the adequacy of
continuum models in describing the response of the saturation to the heterogeneity. The sensitivity
of the saturation to heterogeneity at low displacement rates suggests the possibility of heterogeneity
identification from saturation maps.
PATTERNS OF DNAPL SPREADING UNDER THE INFLUENCE OF GRAVITY
The previous sections illustrated the relevance and utility of pore network models, both full-scale
simulation and simpler statistical physics models, in describing immiscible displacement patterns.
In this section, we we will use a similar approach to examine the spreading of DNAPLs.
Spreading of DNAPLs in the subsurface, when they encounter a low-permeability region, is
subject to capillary, gravity and viscous forces. The previous pore-network approach is equally well
applied to describe these displacement patterns. While a full-scale simulation is possible, however,
we shall use here the simpler model of IPG. As discussed above, in a stabilizing IPG, a fractal front
with a width satisfying the scaling (2.5) is followed by a compact front. In a destabilizing IPG,
on the other hand, the invasion process is dominated by the growth of a single finger, assuming
sufficiently large gradients. The corresponding scaling of the finger width is also given by the same
equation. It was shown in the previous sections, that IPG can describe both the effect of gravity
and the effect of viscous forces. In the typical application of IPG, the gradient is specified in the
direction of displacement. For example, this is the case in viscous gradients presented above. In
the case of a DNAPL spreading at the bottom of an impermeable barrier, however, the gradient
due to gravity will be perpendicular to the direction of displacement. This gives rise to a new and
interesting problem of IPG, in which two gradients exist: one in the direction of displacement (for
example, due to viscous forces), and another in the direction perpendicular to the displacement
(due to gravity). To our knowledge, this problem has not been addressed in the literature, so far.
This section explores aspects of the solution of this problem.
Consider IPG in the presence of gradients both transverse and parallel to the direction of
displacement. Simulations were carried out in 2-D (square) and 3-D (cubic) lattices of sizes 800 x
400 and 128x128x128, respectively. The algorithm is an IP algorithm, in which the site with the
lowest threshold among all sites in the advancing front is the next site to be penetrated. Each site
is assigned a threshold according to
16
-------
r = rand + Bx x X + BY x
(2.14)
where rand is a random number in the interval (0,1), and Bx and By are the Bond numbers in the
parallel, X, and transversely, directions, respectively. Coordinates X and Y express dimensionless
distance in terms of pore lengths and the transverse coordinate Y increases upwards. Stabilizing
gradients correspond to a positive Bond number, destabilizing gradients to a negative Bond number.
To simulate DNAPL spreading at the bottom of the boundary, a destabilizing (negative) Bond
number was selected. (The opposite case would correspond to a smaller density fluid, e.g. a gas,
penetrating the medium and collecting at the top boundary). Recall that for the case of gravity,
the Bond number expresses the ratio of buoyancy to viscous forces, while for the case of viscous
forces it is related to the capillary number Ca. In the presence of a gradient in the transverse
direction, we considered injection from the entire face (edge injection) X = 1, all sites of which are
occupied. For the case of only a parallel destabilizing gradient, injection from the middle point of
the face X = 1 was considered. In the latter case, we also used periodic boundary conditions across
the lateral boundaries. A number of different realizations were carried out, most results shown
corresponding to averages over 50 runs.
Results
We first present simulations of a destabilizing IPG, where By — 0. Figures 2.9 and 2.10
correspond to 2-D and 3-D simulations, respectively, for various values of the parallel Bond number,
BX , the physical picture corresponding to spreading due to viscous and capillary (but not gravity)
forces. As expected, the single finger dominates the process, the width of which increases, as the
Bond number decreases. As before, the relation between the width and the (absolute value of) the
Bond number must be the power law of equation (2.10). Results are shown in Figure 2.11. We
obtain the exponent values of 0.6303 for 2-D and 0.5335 for 3-D. These should be compared to
Meakin et al. (1992), where the values of 0.60 and 0.58, respectively, are found. The corresponding
theoretical values from percolation theory are 0.57 and 0.47, respectively. This discrepancy is
probably due to the relatively small number of realizations used and the limited length scale of the
lattice covered by the simulations, as also explained in Meakin et al. (1992).
Figures 2.12 and 2.13 shows displacements using only a transverse gradient (Bx = 0). The
physical picture corresponding to this problem is DNAPL spreading subject to capillary and
gravity (but not viscous) forces. The pattern spreads on the top of the impermeable boundary, as
expected physically, and it has the characteristic percolation features. Of interest is the width, ay,
of the invaded zone along the vertical direction. A theoretical analysis indicates that the problem
belongs to the same class as IPG, thus we expect the same power-law scaling as in (2.10), namely
ay ~ \By\-^ (2.15)
Results plotted in Figure 2.14 confirm the validity of (2.15) in this case, and indicate exponents
equal to 0.6339 and 0.5326 for 2-D and 3-D, respectively. These values compare very well with the
previous case of destabilizing IPG. Thus, it can be concluded that IPG with a single transverse
gradient has the same features as IPG in a destabilizing parallel gradient. This means that in the
absence of viscous forces, the spreading DNAPL has percolation characteristics and a thickness
which decreases as By = Ap^yfc increases, namely as the permeability k increases or the interfacial
tension decreases.
17
-------
Consider, next, the case where both gradients are combined. This would be the case of DNAPL
spreading in the presence of capillary, gravity and viscous forces. A series of simulations were
conducted for stabilizing and destabilizing parallel gradients.
The case of a stabilizing gradient corresponds to viscous-stabilized displacement, which was
described in a previous section. Here, we expect that BX would be related to the velocity through
a relationship of the type Bx ~ (iH°, where a = lD_ (e.g. compare with equation (2.9),
and which takes values equal to 0.66 and 0.53, in the respective geometries. Thus, larger Ca
corresponds to larger BX- Figures 2.15 and 2.16 show 2-D invasion patterns corresponding to a
variety of values of BX and By. The patterns are characterized by a leading front with a fractal
structure, which is expected to scale as (2.5), followed by a compact front. A simple analysis
suggests that this problem is actually IPG in a stabilizing gradient but with a Bond number which
is the geometric sum of the two Bond numbers, namely
l + BY (2.16)
Thus, the size of the front width decreases with an increase in B. For sufficiently large B, the front
is quite sharp and it is well described as a straight line of slope Bx/By (note that this slope is
negative). Figures 2.15 and 2.16 demonstrate indeed the validity of this result. Similar features
were found for invasion in a 3-D lattice. In this case of stabilized displacement, therefore (namely
when M or Ca are sufficiently small, see Yortsos et al., 1997), the slope of the front scales as
slope
(2.17)
namely the front is more horizontal (has more of a tongue-like shape) when the velocity is small
or buyoancy is strong, as expected. However, note the exponent a in the scaling of the front with
the velocity or the heterogeneity, which reflects percolation characteristics. This result should be
contrasted to the classical, effective medium result for the slope of a propagating sharp interface in
a porous medium, which is
slope
qp,
(2.18)
corresponding effectively to o = 1. We conclude that IPG leads to a non-trivial new result on the
pattern structure.
The case of a destabilizing parallel gradient corresponds to the problem of unfavorable mobility
ratio (sufficiently large M). Now, the Bond number is of the type shown in (2.11), namely BX ~
—Qjjf. Patterns obtained from simulations for various values of BX and By are shown in Figures
2.17 and 2.18. We note that in 3-D the width, az, along the lateral (Z) direction of the invading
cluster depends on the parallel gradient. However, the thikness ay along the transverse direction
depends on both gradients. Figure 2.19 shows the dependence of the width along the lateral Z
direction for the 3-D case. As before, the width is a power-law, this time of the parallel gradient,
BX, only, with an exponent estimated at 0.59. This value compares reasonably well with the
previous value 0.5335, above, and 0.58 of Meakin et al. (1992). This result indicates that the
transverse gradient has little influence on CT#, which behaves essentially as in the presence of only
one destabilizing parallel gradient.
To understand the extent of the invasion zone along the transverse direction Y", we consider a
simple analysis based on IPG arguments. This suggests the following scaling
ay = B
x
:/(*)
(2.19)
18
-------
where we denoted z = j£ and the scaling function f(z) is determined with the following reasoning.
In the limit of small transverse gradients, the problem should reduce to that portrayed in Figures
2.9-2.11, hence we should expect the result f(z) ~ 1 as z.->- oo. On the other hand, in the limit of
small parallel gradients, the result pertaining to Figures 2.12-2.14 applies, hence we should expect
the different behavior, f(z) ~ z^+i as z -> 0. Combining the two we have
/(z)~2r"+i as z ->• 0 and f(z) ~ 1 as 2;->• oo (2.20)
To test the validity of these predictions, we plotted the results according to the scaling function
above in Figure 2.20. For the 2-D simulations a total of 220 different combinations is shown (where
the range of BX is (10~5 - 4 X 10"1), and the range of By is (6x 10~6-4xlO~1)), while for the
3-D simulations a total of 44 combinations of BX and By is plotted (where the range of BX is
(10~5-10~1), and the range of By is (5x lO^-lO"1)). For each combination of parameters, the
ensemble average of 50 realizations was obtained. Figure 2.20 shows that the data collapse on a
single scaling curve with the scaling properties described in (2.20), thus supporting well the scaling
result.
We summarize this section by concluding that the spreading of DNAPL on the top of an
impermeable medium may well have strong percolation (namely fractal-like) characteristics, when
the rate of penetration is slow, so that gravity and capillarity dominate over viscous forces (and
also in the case of unfavorable mobility ratio). In such cases, the percolation characteristics must
be carefully considered in the design of remediation processes.
USE OF IPG FOR DESCRIBING AIR SPARGING PATTERNS
We close this chapter with a very brief presentation of the use of IPG for the description of
air sparging patterns. Experiments of Dumore (1970) in oil recovery and of Ji et al. (1993) in
air sparging have showed that the upward migration of gas depends on two conditions related to
the capillary pressure. At conditions of high capillarity, the air-occupied region is essentially a
conical region, while at conditions of low capillarity, only one gas channel develops. These two
conditions were termed dispersive and non-dispersive, respectively. As discussed previously, air
sparging could possibly be described as a destabilizing IPG. However, Figures 2.9 and 2.10 above
show that destabilizing IPG leads to a single finger, regardless of the value of the Bond number.
This constrasts to the experiments for the dispersive case. To model this problem using IPG, we
proceeded with a modification of the traditional algorithm.
Assume that at a given stage, the invading phase penetrates not only the site with the smallest
capillary threshold, but also several sites. The justification is that in reaching "steady-state",
pressure fluctuations in both invading and defending phases develop, so that more than one site
can be penetrated. One algorithm to implement this idea involves penetrating all neighboring sites,
the threshold of which is close to the maximum threshold at the front, at a given stage, namely
\T-Tr,
< P (2.21)
Here the thresholds T are defined as r = rand +BgY, where Bg is a gravity Bond number and Y
indicates the vertical coordinate, while parameter P measures the extent of the fluctuations and
controls the transition from non-dispersive (small P) to dispersive (large P). Using (2.21) we can
aPply IPG except that now multiple sites can be penetrated. The simulations were carried out in
a 200x 200 lattice. Typical results are shown in Figure 2.21 for a given value of the Bond number
and different values of the fluctuation extent, P. We note that as P increases, multiple sites are
penetrated resulting in a conical-shape pattern (the artifact at the top of the pattern should be
19
-------
disregarded), in qualitative agreement with published experiments. Further work is in progress to
improve the modeling of such patterns.
ACKNOWLEDGEMENTS
This research was also partly supported by DOE Contract No. DE-FG22-93BC14899.
20
-------
REFERENCES
10.
11.
12.
13.
14.
15.
16.
Blunt, M. and King, P.R, Relative permeabilities from two- and three-dimensional pore-scale
network modeling, Transp. Porous Media, Vol. 6, 407-433 (1991).
Chaouche, M., Rakotomalala, N., Salin, D., Xu, B., and Yortsos, Y.C., Invasion percolation
in a hydrostatic or a permeability gradient: Experiments and simulations, Phys. Rev. E,
Vol. 49, 4133-4139 (1994a).
Chaouche, M., Rakotomalala, N., Salin, D., Xu, B., and Yortsos, Y.C., Capillary effects
in drainage in heterogeneous media: Continuum modeling, experiments and pore network
simulations, Chem. Eng. Sci., Vol. 49, 2447-2466 (1994b).
Chuoke, R.L., van Meurs, P., and van der Poel, L.B., The instability of slow, immiscible,
viscous liquid-liquid displacements in permeable media, Trans. AIME, Vol. 216, 188-194
(1959).
Dumore, J.M., Development of gas-saturation during solution-gas drive in an oil layer below
a gas cap, SPEJ, Vol. 10, 211-218 (1970).
Gouyet, J.-F., Sapoval, B. and Rosso, M., Fractal structure of diffusion and invasion fronts in
three-dimensional lattices through the gradient percolation, Phys. Rev. B, Vol. 37, 1832-1842
(1988).
Haghighi, M., Xu, B. and Yortsos, Y.C., Visualization and simulation of immiscible displace-
ment in fractured systems using micromodels: I. Drainage, J. Coll. Interface Sci., Vol. 166,
168-179 (1994).
Hulin, J.-P., Clement, E., Baudet, C., Gouyet, J.-F., and Rosso, M., Quantitative analysis of
an invading fluid invasion front under gravity, Phys. Rev. Lett., Vol. 61, 333-336 (1988).
Ji, W., Dahmani, A., Ahlfeld, D.P., Lin, J.D., and Hill, E., Laboratory study of air sparging:
Air flow visualization, GWMR, 115-126 (1993).
Lenormand, R., Liquids in porous media, J. Phys.: Condens. Matter, Vol. 2:SA, 79-88 (1990).
Li, X., and Yortsos, Y.C., Visualization and simulation of bubble growth in pore networks,
AICHEJ, Vol. 41, 214-222 (1995).
Meakin, P., Birovljev, A., Frette, V., Feder, J., and J0ssang, T., Gradient stabilized and
destabilized invasion percolation, Phys. Rev. A, Vol. 46, 227-239 (1992).
Stauffer, D., and Aharony, A., "Introduction to Percolation Theory", Francis-Taylor (1992).
Wilkinson, D., and Willemsen, J.F., "Invasion percolation: A new form of percolation theory",
J. Phys. A, Vol. 16, 3365-3376 (1983).
Xu, B., Ph.D. dissertation, University of Southern California, (1995).
Xu, B., Yortsos, Y.C. and Salin, D., Invasion percolation with viscous forces, Phys. Rev. E,
Vol. 57, 739-751 (1998).
21
-------
17. Xu, B., Kamath, J., Yortsos, Y.C. and Lee, S.H., "Use of Pore Network Models to Simulate
Laboratory Corefloods in a Heterogeneous Carbonate Sample", paper SPE 38879 presented
at the 1997 SPE Annual Conference, San Antonio, TX (Oct. 5-8, 1997).
18. Yortsos, Y.C., Xu, B., and Salin, D., Phase diagram of fully-developed drainage, Phys. Rev.
Lett, Vol. 79, 4581-4584 (1997).
22
-------
(b)
(c)
(d)
Figure 2.1. Displacement patterns from numerical simulation of drainage in a lOOx 100 lattice, with
M = 0.1 and different values of the capillary number: (a) Ca = 3.2 X 10~8 (IP) , (b) Ca = 3.2 X 10~6,
(c) Ca = 3.2 X 10~5 and (d) Ca = 3.2 X 10~4.
23
-------
Figure 2.2. Schematic of the front region and the notation used.
24
-------
io
Figure 2.3. Numerical solution of (2.5) vs. the modified capillary number for three different values
of M.
25
-------
0.1 0.3 0.3 0.4 05 0.6 0.7 0.8 0.8 1
X
0.1 0.2 0.3 0.4 0.5 0.8 0.7 0.8 0.8 1
(a)
(b)
Figure 2.4. Percolation probability (dashed line) and saturation (solid line) profiles vs. normalized
length, for 3-D VGP: (a) B = 10~4 in a lattice of 200 x 200 x 65, (b) B = 10~6 in a lattice of
200 X 200 x 500.
26
-------
Figure 2.5. Single finger in GP in a destabilizing gradient.
27
-------
•i J* •
4- V i-*i?* SJrf^
\^:1"VSF;- .*.^2
•jj^.i !!_,'
^^•v^^^^^r>Jfc! 1
Figure 2.6. Pore network simulations (saturation profile and displacement patterns) for uniform
permeability and Ca = 1.5 x 10~4.
28
-------
10 20 30 40
Figure 2.7. Pore network simulations (saturation profile and displacement patterns) for permeabil-
ity increase (1:2) and for Ca — 1.5 X 10~5.
29
-------
I lnli'
IU-1
Figure 2.8. Pore network simulations (saturation profile and displacement patterns) for permeabil-
ity increase (1:2) and for Ca = 1.5 x 10~4.
30
-------
400
350
300
250
200
150
100
SO
400
350 •
300 •
250-
200
ISO-
100-
50
100 200 300 400 500 600 700 800
100 200 300 400 500 600 700 800
(1)
(2)
200
100 200 300 400 500 600 700 800
(3)
(4)
Figure 2.9. Displacement pattern with parallel gradient in a 2-D 800X 400 lattice: (1) BX
(2) Bx = -0.01, (3) Bx = -0.001, (4) Bx = -0.0001.
= -0.1,
31
-------
120.
100
SO
>• SO-
40.
20
100
(i)
(2)
(3)
(4)
Figure 2.10. Displacement pattern with parallel gradient in a 3-D 128x 128 x 128 lattice: (1)
Bx = -0.03, (2) Bx = -0.003, (3) Bx = -0.001, (4) Bx = -0.0003.
32
-------
-1
slope = • 0.6303
-8 -7-6-5-4 -3 -2 -1
(a)
(b)
-2
Figure 2.11. The dependence of finger width on the Bond number for (a) 2-D, and (b) 3-D lattices.
33
-------
200
IN
180
140
120
20
SO 100 ISO 200 250 300 350 400
X
(I)
(3)
200
ISO
180
140
120
>-IOO
SO
60
40
20
SO
100
150
200
X
250
350
2M250 300 350 400
(2)
(4)
Figure 2.12. Displacement pattern with transverse gradient in a 2-D 800 x 400 lattice: (1) By
-0.1, (2) By = -0.01, (3) By = -0.001, (4) BY = -0.0001.
34
-------
(1)
(2)
120
100
80,
> 80.
(3)
(4)
Figure 2.13. Displacement pattern with gradient in a 3-D 128 x 128 x 128 lattice: (1) By
(2) By = -0.01, (3) By = -0.001, (4) Bx = -0.0005.
= -0.1,
35
-------
4.5
3.5
3
2.5
1.5
-7
-5
(a)
(b)
-2
Figure 2.14. The dependence of finger width ay on By for (a) 2-D, and (b) 3-D lattices.
36
-------
200
> 100
SO 100 1SO 200 250 300 350
50
200
>100
(I)
(4)
200
SO 100 1M 200 250 300 350
350
(2)
(5)
200
100 150 200 250 300 350
200
>100
SO
SO tOO 150 200 250 300 350
(3)
(6)
Figure 2.15. Displacement pattern with both transverse and parallel gradient in a 2-D 800x 400
lattice: (1) Bx = 0.01, BY = -0.001, (2) Bx = 0.01, By = -0.005, (3) Bx = 0.01, BY = -0.01,
(4) Bx = 0.005, By = -0.01, (5) Bx = 0.001, By = -0.01, (6) Bx = 0.001, By = -0.001.
37
-------
200
ISO
>too
SO 100 150 200 250 300 350
200
>100
50
350
(1)
(4)
200
>100
50
200
180
180
140
120
>100
80
60
40
20
(2)
(5)
200
100 1SO 200 250 300 350
SO
(3)
Figure 2.16. Displacement pattern with both transverse and parallel gradient in a 2-D 400 x 200
lattice: (I) Bx = 0.1, By = -0.01, (2) Bx = 0.1, By = -0.05, (3) Bx = 0.1, By = -0.1, (4)
Bx = 0.05, By = -0.1, (5) Bx = 0.01, By = -0.1.
38
-------
200
180
160
140
120
>100
80
60
4O
20
r* -
so
100
150
200
X
(1)
250
300
350
200
180
160
140
120
>100
80
60
40
20
(4)
200
180
160
140
120
>100
80
60 •
40
20
350
400
(2)
300
350
400
(5)
200
180
160
140
120
•100
80
60
40
20
100
150
200
X
(3)
250
300
350
400
350
400
(6)
Figure 2.17. Displacement pattern with both transverse and parallel gradient in a 2-D 400x 200
lattice: (1) Bx = -0.01, By = -0.01, (2) Bx = -0.001, By = -0.01, (3) Bx = -0.0005,
By = -0.01, (4) Bx = -0.01, By = -0.001, (5) Bx = -0.001, BY = -0.001, (6) Bx = -0.0005,
BY = -0.001.
39
-------
129.
100.
*0
29.
1M
...*""": L..--? ""i f"'4. : ''J-. :
>••*": : ..-I"" : :''•••! I'"'•
•"*" • : >» 1 • t :'•-.:
,.- ; ,.'-''*' • ' .-••!""'$• : !••.
••'''"'"
120
100
M
>• tt.
40
20.
120
i ...•>•"'{•. i ':••. :
120.
ICO.
M.
«0.
U.
20.
...
I'*
,» *
."
' 4
J.I-.V
• • •" 1 : : '"••. •
i ... > "'i ! ..-•»•"'•-.. : ''i'-.
:••'"• [-•-'!"""] 1. -J '''••!
.«•' ^ • ^.'''"'1 ; : '"i', '
f "'• \ „.•!••"""'!"•"".[ T •
*.'•""* • ? .*• "J"**"I . * '*'-
.•**"* • ,-'T**' i ! *"••; ; *'
>•••' !"" i ...•;-"'1'""r--j i'"--
** 1 * **•*';"*' '• ' '' ' ' . •
* .. •* . • *""•*.*" • •'* ,
.
.'"H.
:''•• ] «00 .
* . ' 10
•-..
'•-.: >- tt.
•'*••• L u •
;' • .\_ 20 .
iv»'v" '
_,.••
^'
JSi
fee
•t'".
'.'••• :'
.'•*.'.
V'- •»
'"^j
;'
&'
ife
••••"
^a
...•••"' : i "'-i.
• ..»••'•;. '• '''*,, .
..-;• ' " j i""" .;. : *!..
i ..-i' ''I i *"!-.. j "
..••*T"' i i ***!•• i '**••
.i.--!"'k.j '"'I---.!.
i ••••• J t '" •• ;
.•!•**'* * " ***•,• ; '"''i.
' j .-•]•"•[ j'"' !--.!'
:•':."' .'.';-|2H|MJf ;
••-.
"•
100
(5)
Figure 2.18. Displacement pattern with both transverse and parallel gradient in a 3-D 128x 128x
128 lattice: (1) Bx = -0.01, By = -0.1, (2) Bx = -0.01, BY = -0.005, (3) Bx = -0.01, By =
-0.001, (4) Bx = -0.0005, By = -0.1, (5) 5* = -0.0005, By = -0.005, (6) Bx = -0.0005,
By = -0.001.
40
-------
log(
3.5
2.5
1.5
1 -
0.5
-8
o
*
+
X
*
•
Bo=0.1
Bo=0.05
Bo=0.025
Bo=0.005
Bo=0.0025
Bo=0.0005
slope = - 0.5942
-6
-5
-4
-3
-2
Figure 2.19. The dependence of finger width 0% on BX for a 3-D lattice and various values of By.
41
-------
-1
-2
-5
slope = 0.63
-5
0 5
log(.t)
10
15
0.5
-0.5
-1
-1.5
-2
-2.5
O „ O
<§>
6>°
o
o
slope = - 0.53
_1 I I I 1 1 1 L.
-6 -4 -2 0 2 4 6 8 10 12
log(.t)
Figure 2.20. The scaling function f(z) from simulations in: (a) 2-D lattice, (b) 3-D lattice.
42
-------
(3)
Figure 2.21. Modified destabilizing IPG patterns to model air sparging for B = -0.1 and P = 0,
1 and 2, respectively.
43
-------
Chapter 3
VISUALIZATION AND
SIMULATION OF NAPL
SOLUBILIZATION IN PORE
NETWORKS
Chunsan Jia, Katherine Shing and Yanis C. Yortsos
INTRODUCTION
A common source of soil and water aquifer contamination is the accidental release or spill of
organic liquids, commonly known as Non-Aqueous Phase Liquids (NAPL), which percolate to the
subsurface due to gravity. This well-known problem is quite extensive and poses potential health
hazards, particularly if the NAPLs reach and contaminate aquifers. A variety of remediation
techniques have been developed to mitigate the problem.
As pointed out in the Introduction, the most commonly used method is pump-and-treat, in
which the groundwater is pumped out of the site, treated on the surface to remove the contaminant
and subsequently reinjected (Mackay and Cherry, 1989). In-situ remediation involves the injection
of remedial fluids or chemicals which may enhance the solubility of the contaminant (Pennell et
al., 1993, Augustin et al., 1994) or help to mobilize trapped NAPL (Brandes and Farley, 1993).
Injection of surfactants (West and Harwell, 1992, Pennell et al., 1993), co-solvents, such as alcohols
(Final et al., 1990, Wood et al., 1990), and air (sparging) (Ji et al., 1993, Pankow et al., 1993)
in the groundwater, or injection of steam in the vadose zone (Itamura and Udell, 1993, Ho et al.,
1994), are some of the methods actively under consideration.
Of interest to this work is the mass transfer of contaminants from residual NAPLs. As with
other processes, this must be investigated at the different scales that comprise porous media. These
include the pore scale or microscopic scale (of the order of tens of /j,m), where capillary phenomena
are important, the pore-network scale (of the order of mm), where the collective behavior of a
process over an ensemble of pores is investigated, the laboratory (core) scale (of the order of cm),
and the field or megascopic scale (of the order of m and higher), of relevance to large-scale spreading.
Phenomena at the different scales have a different manifestation, and input from a lower scale is
necessary for their description at the next higher scale.
44
-------
With the exception of some recent works, the design of in-situ remediation is based on a descrip-
tion at the macroscopic (or megascopic) scales. Phenomena at the pore and pore-network scales are
typically lumped in terms of averaged quantities, using empirical or ad hoc expressions. At best,
these are adaptations of macroscopic models developed in enhanced oil recovery (e.g. see Abriola
and Finder, 1985a,b, Brown et al., 1994), and cannot address remediation issues originating at the
pore and pore network scales. Recent approaches (Celia et al., 1993, Soil and Celia, 1993, Soil et al.,
1993), paralleling earlier efforts in oil recovery (e.g. Heiba et al., 1982) have used percolation theory
to relate quantities, such as relative permeabilities and capillary pressures, to the microstructure.
However, the problem of mass transfer, which is fundamental to NAPL remediation is currently
based on macroscopic simulation.
The overall mass transfer coefficient has been estimated by different correlations (Table 1),
reflecting the different assumptions made by different investigators. Thus, inconclusive at present
are the effects of flow velocity, interfacial area and its variation with time (due to dissolution,
mobilization, etc.) (Imhoff et al., 1993). Particularly interesting are the reported differences in
the exponent relating the mass-transfer Sherwood number to the local Reynolds number, for which
convincing theories have not been provided. In principle, this exponent should reflect the relevant
flow regime and the interface configuration (Cussler, 1984). In the most advanced model todate,
Powers et al. (1994) approximated mass transfer by steady-state diffusion from an isolated spherical
NAPL globule in a structureless porous medium. However, this approximation does not adequately
represent the true geometrical configuration of NAPL-groundwater interfaces in a real system, which
are disordered (and possibly fractal) (Feder, 1988, Lenormand, 1990, Sahimi and Yortsos, 1990)
and the concommitant effects on fluid flow and mass transfer.
Issues of mass transfer between immiscible phases in disordered porous media at the pore-
network scale were addressed in recent studies, although in a somewhat different context. In a
series of papers, Li and Yortsos (1994, 1995a,b), Satik et al. (1995) and Satik and Yortsos (1996)
used pore-network models to study mass transfer towards a growing gas phase from a surrounding
liquid in a porous medium, in the context of solution gas-drive and/or boiling. Experiments in
etched-glass micromodels and in Hele-Shaw cells and simulations with pore network models showed
that the interface pattern, varying between Invasion Percolation (IP) and Viscous Fingering (VF)
(Feder, 1988), drastically affects the mass transfer, hence the rate of growth (or dissolution) of a
phase. The disorder of a real porous medium (even if macroscopically homogeneous) was shown to
give rise to qualitatively different phenomena than in effective porous media (such as a Hele-Shaw
cell, compare for example Li and Yortsos, 1994, 1995a). We expect that similar effects will occur
in remediation processes as well.
To better understand and to help design improved remediation methods, we present in this paper
a pore scale study of the mass transfer of the remediation of residual phases. From a mass transfer
perspective, the problem is interesting from two points of view: That mass transfer occurs in a
disordered pore space from a disordered stationary source, and that it involves the interplay between
small (pore-size) and large (pore-network size) length scales. The basic problem is described as
follows. An immiscible to water NAPL has percolated through the action of gravity or otherwise to
the subsurface and displaced part of the water in place. As a result of the subsequent displacement
of the NAPL by a displacing fluid, e.g. water, the NAPL phase has been trapped due to capillarity
in various pores in the form of trapped and disconnected ganglia (as shown in the schematic of
Figure 3.1). To remove this trapped liquid, a solubilizing agent, fully micsible to the original fluid
in place, but not to NAPL, to which is partly miscible, is injected. We are interested in evaluating
the mass transfer rates from the trapped NAPL to the solubilization agent and its dependence
on parameters such as the flow rate, in order to be able to upscale the process and to develop
appropriate mass transfer expressions. The particular study in this paper is restricted to mass
45
-------
transfer from residual NAPLs, the interfaces of which remain stationary during the experiment,
thus solubilization by the mobilization of trapped blobs is not considered.
The paper is organized as follows: First, we report on visualization experiments, conducted
in etched glass micromodels. The solubilization of various residual NAPL phases, trapped in the
micromodel following displacement with a driving fluid, is studied. We monitor the variation of
the effluent concentration with time and, after it reaches a steady-state, with the injected flow
rate. Subsequently, a pore network simulator is developed to describe mass transfer. The simulator
is calibrated by comparison against tracer displacement experiments. Then, the experiments are
simulated with the use of the pore network model. It is found that appropriately modified local
mass transfer coefficients need to be introduced in the pore-network model in order to match the
experimental results. Sensitivity studies are then conducted to analyze the dependence of the
response to the microstructure and particularly to the assumed local mass transfer coefficients.
Implications on the upscaling of the results are also discussed.
EXPERIMENTS
To obtain insight on mass transfer processes during NAPL remediation, visualization experi-
ments in an etched-glass micromodel were conducted. Such micromodels are prepared by etching
a 2-D square network pattern of pores and throats on a glass plate (Figure 3.2), which was subse-
quently fused to another glass plate to create a 2-D pore network (see Chatzis, 1982, Buckley, 1990,
also Li and Yortsos, 1995a). One injection port and one production port at opposites sides of the
micromodel were drilled to control the injection and production of fluids. Various patterns, with
randomly distributed pore widths were etched, the typical micromodel size being 29x15 pores. The
pores have a rectangular cross-section with a channel depth of about 0.0312 cm in the typical case.
In all patterns, the lattice coordination number was equal to 4 (square lattices). The pore volume
in the various micromodels ranged between 0.5 cm3 to 2.0 cm3 (see also below).
The schematic of the experimental setup is shown in Figure 3.3. The experimental procedure
was as follows:
(1) First, the micromodel was saturated with distilled water.
(2) Subsequently, using a syringe pump (Harvard Apparatus model 850), water was displaced
at sufficiently high rates (order of 0.053 ml/min) by a NAPL containing a visualization dye,
until all water was completely removed from the micromodel. Often, this required tilting the
micromodel for gravity to assist in the displacement process. Various liquids were used to rep-
resent NAPL, including TCE, PCE and benzonitrile. However, complete sets of experiments
were conducted only for benzonitrile. These are the only reported here.
(3) The next step involved displacing the NAPL by a solvent, usually water or a water-alcohol
mixture, at high rates (of the order of 2.85 ml/min) (see also below). Following a transient
stage, of approximately 1 hour, during which much of the NAPL was displaced, the system
reached a stage in which the NAPL was trapped in well-defined ganglia. During the time of
the subsequent mass transfer experiments in this state, there was a negligible reduction of
the NAPL-water interfacial configuration, suggesting a steady-state as far as mass transfer
is concerned. All processes and patterns were videotaped. Comparison of the patterns at
different times ensured that there was no significant change in the pattern, hence in the flow
distribution.
46
-------
(4) Samples at the production port were collected and analyzed immediately (less than a minute
elapsed) after collection using a Varian 3300 gas chromatograph with a thermal conductiv-
ity detector (TCD) for measuring concentrations. The accuracy of the measurements was
validated by independent comparison of saturated samples with published solubility data.
In this approach, a saturated solution (with solubility obtained from the literature) was an-
alyzed using the GC. The area under the concentration peak was measured and served to
scale the subsequently obtained concentrations. After ascertaining that a steady-state was
indeed reached, the effluent concentration was recorded. Then, the injection rate was re-
duced, another steady-state was established and the effluent samples were analyzed. The
process was repeated for a range of injection rates covering more than three decades (between
0.0007 ml/min and 2.8 ml/min). The transient times required to reach the steady-state varied
between 20 and 30 minutes for the high rates to several hours for the lower rates. All the
experiments were carried out at room temperature.
Miscible displacement experiments of a passive tracer, in the absence of a NAPL, were also con-
ducted to test the validity of the pore-network simulator which will be subsequently used to analyse
the results. In these experiments, a 50% glycerol solution was displaced by the same mixture con-
taining a coloring dye. During fluid injection, the pressure difference across the two ports was
periodically monitored, using water-filled U tubes.
Shown in Figure 3.4 is a typical picture of the configuration of the trapped NAPL at the
onset of a mass transfer experiment. Note the disconnected state of the NAPL, the existence of
ganglia of various sizes, and the varying local interface configuration of the different ganglia. It
should be pointed out that because of the rather high velocity used to isolate and trap the NAPL,
the geometric configuration is not that of a percolation cluster, which one would expect at lower
trapping velocities, where viscous forces are small. However, this is of no significant consequence for
this paper, where the analysis and simulation will be based on the actual, rather than a theoretical,
pattern. All mass transfer experiments were carried out in order of decreasing flow rate in order
to keep the NAPL interfaces fixed, thus allowing for a meaningful sensitivity analysis of the effect
of flow rate on mass transfer. The inverse order, from low to high rates, would have resulted in
the mobilization of some trapped ganglia at increasing velocities (compare with Figure 3.1) with
a resulting change of the interfacial geometry. Such mobilization was indeed observed at high
capillary numbers. This process will be investigated in a separate study.
Some additional remarks are also in order. To keep a relatively fixed interface between NAPL
and solvent, which is necessary for the study of the effect of velocity on mass transfer, requires that
the interface recedes only very slightly during the experiments. Given that the solubility of most
NAPLs in water is generally low, this can be readily accomplished by using water as a solubilizing
agent. On the other hand, low solubility presented resolution problems in the GC analysis. In
addition, during the collection and analysis of the samples, particularly at low rates, which involve
longer transients, the losses of NAPL due to its volatility must be minimal. We sought, therefore, to
work with NAPLs whose solubility is low, but not beyond the resolution of our GC, and with small
volatility. Benzonitrile with a solubility of 0.2 g/100 cm3, which is comparable to TCE solubility,
but with a much higher boiling point (190.7 °F, compared to 86.7 °F for TCE), was one such liquid
used. All results shown here pertain to this specific NAPL.
A typical effluent curve from the 29 X 15 micromodel is shown in Figure 3.5 corresponding to the
solubilization of benzonitrile by water. This curve shows the variation of the effluent concentration,
normalized with respect to its solubility at the given temperature, as a function of the injection
flow rate. The curve has a characteristic erfc-like shape, approaching asymptotically zero at large
rates and 1 at smaller rates, where the actual concentration becomes equal to its solubility. The
47
-------
approach to the latter solubility value was tested independently by measuring the solubility of the
NAPL and also by comparison with published data. The experimental error bars shown in Figure
3.5 are associated with errors in the GC measurement and sample collection, due to the rather low
injection rates used. The analysis of the effluent curve is presented in more detail in a following
section.
MATHEMATICAL FORMULATION
To simulate the experiments and to provide a working model for the process, we proceeded with
a mathematical formulation. The problem considered is mass transfer in the porespace of a porous
medium, such as the micromodel of Figure 3.4, containing a trapped liquid phase, which becomes
solubilized in the flowing solvent. The lateral boundaries of the micromodel are closed to flow, with
the exception of two producing sites at opposites ends. The location of the interfaces was directly
inserted in the simulations from observations of the experiments. (An invasion percolation, or other
displacement, algorithm could also be used to determine this configuration, although it will not
necessarily correspond to the particular experiment, for the reasons described.)
Mass transfer of the solubilized NAPL occurs from the direction of the trapped liquid towards
the effluent end, described in the pore-space by the usual convection-diffusion equation
(3.1)
- + u-VC = DV2C
where C is concentration, u is velocity, and D is the molecular diffusivity. At creeping flow condi-
tions and in the absence of gravity, the injected solvent velocity u obeys in the pore-space Stokes'
law
VPi = /j,V2u (3.2)
where PI is liquid pressure and ft is viscosity. The trapped NAPL is taken to be stationary, and
so is the interface between NAPL-solvent, for the reasons explained above. A receding interface
algorithm could also be used, but the conditions of slow dissolution examined here do not warrant
it. At the interface, the local concentration of the NAPL in the solvent is that of thermodynamic
equilibrium
= CS
(3.3)
These equations are subject to the appropriate flux or concentration conditions at the boundaries,
and to the condition that the concentration on the interface of each cluster is spatially constant.
A detailed simulation of the full problem is possible with the use of lattice gas simulation
(Rothrnann and Zaleski, 1997). Here, we eledrto work, instead, with the simpler, although perhaps
less accurate, approach of pore network simulation. In the pore network simulations, mass and
momentum balances will be discretized, with the latter being approximated using Poiseuille's law.
This is a standard approximation in pore-network models.
Pore Network Simulation
Pore network simulators are used to simulate processes in pore-networks or glass micromodels.
The key approximation is the porous medium representation in terms of a collection of pores,
connected to each other by pore throats in a network-like fashion (Dullien, 1992). Conventional
48
-------
networks include regular lattices, such as square or cubic, or Voronoi lattices (Blunt and King,
1991). Either lattice can be readily constructed, although regular lattices are computationally
more manageable. Pore networks are ideal for replicating experiments in etched-glass micromodels
(see Li and Yortsos, 1995a, for a recent application). Whether for miscible or immiscible flow, their
basic aspects include the following:
• Pores provide volumetric storage, throats control the conductance to flow, heat and mass.
Geometrical characteristics can be statistically distributed.
• Capillary equilibrium at the pore scale is controlled by pore throats during drainage (a non-
wetting (nw) fluid displacing a wetting (w) fluid), and by sites during imbibition (which is the
inverse of drainage). The capillary pressure across an interface is given by the Laplace-Young
equation, Pnw — Pw = 27%, where 7 is the interfacial tension, and at equilibrium the mean
interface curvature, H, is controlled by the pore geometry.
• Pore bodies can be occupied by one or both of the two fluid phases. Pores occupied by one
phase can become trapped by the other phase, if topologically possible. When a pore is
trapped, liquid phase movement is not allowed, although mass diffusion continues.
• Local pore-scale coefficients are calculated from single-pore studies. Thermodynamic equilib-
rium applies within any given pore, although adjacent pores can be at non-equilibrium.
In the past, pore-network simulators have been developed to model immiscible displacement (e.g.
see Lenormand, 1990, Blunt and King, 1991, Haghighi et al., 1994). These simulators provide
a pore-by-pore account of the displacement of interfaces during drainage or imbibition. In mul-
ticomponent displacements, mass transfer of various species to and from interfaces must also be
computed. In recent papers (Li and Yortsos, 1995a,b, Satik et al., 1995, Satik and Yortsos, 1995)
we have described how such a simulator must be constructed for a problem involving phase change.
The governing equations are discretized over the pores (nodes) of the network. For an incom-
pressible fluid, the overall mass balance at node i is
£Qij = 0 (3.4)
3
where Qij refers to volume flow rate between adjacent sites i and j. The sum is over all neighboring
sites j, and we use Poiseuille's law to relate flow rates to pressure drops
Qij —
(3.5)
where the overall fluid conductance, (?»_,•, is a distributed variable. For the square cross-section
channels of the glass micromodel, we have in particular,
Gij = iiL^ (3-6)
where w,j is the width of the pore throat joining sites i and j, h is the pore depth (taken to be
spatially uniform due to etching) , / is the bond length (uniform for this problem), /j, the fluid
viscosity, and kij is a geometric correction factor given by the following correlation of data from
Purday (1949)
= -0.033124 + 0.80114 -
f,- + 30.3497r?- - 3.0355ri7- + 67.6754
y
(3.7)
49
-------
where r = Wij/h.
The mass balance for solute in site i not neighboring the NAPL-solvent interface reads
Vs
AC,-
At
(3.8)
3 " 3
where Vs is the site volume, D is diffusivity and A;J — Wijh is the cross-section of the pore throat.
Here, concentrations were assigned to sites only, and the mass transfer between adjacent sites
was assumed to occur by both diffusion and convection. If a site neighbors a NAPL, however, a
modification to account for local mass transfer is necessary, to be described below. Published values
were used to estimate fluid properties, such as viscosity and diffusivity. Contaminant diffusivity can
be approximated by infinite-dilution coefficients, in the case of low concentrations, using methods
such as Tyn and Calus (1975), or by more complex models (Taylor and Krishna, 1993) in the
case of concentrated solutions. Solubility parameters for various systems were estimated from
currently available data (Mercer et al., 1990, Edwards et al., 1994, Glinski et al., 1994). The set
of concentration and pressure equations was simultaneously solved by a numerical method using
SOR.
Local mass transfer coefficients
As previously noted, the NAPL-water interface could reside in throats or in bodies (sites).
Discretization (3.8) for mass transfer applies only to cases in which interfaces are not within the site,
or within pore throats adjacent to the site over which the mass balance is written. For sites which
are nearest neighbors of the NAPL-solvent interface, however, the mass transfer expression needs
to be more carefully examined. Indeed, near a stagnant interface, the local mass transfer cannot be
captured simply by diffusion and convection, where the flow field obeys Poiseuille's (or effectively,
Darcy's) law, even if a more refined grid was to be used (for example, if we were to also assign
concentrations to throats, which has also been considered here). Consider, for example, typical
configurations of the NAPL-solvent interface, shown in Figure 3.6, obtained from photographs of
the micromodel. To more properly account for the mass transfer rate from the interface under
these conditions requires accounting for Stokes flow dynamics, near the interface, which, however,
are not inherent to our pore network.
To proceed, we will keep the hydrodynamic field as dictated by Darcy's law (3.4) and (3.5), but
modify the terms for mass transfer from the interface to an adjoining site by introducing a local
mass transfer coefficient k and expressing the corresponding flux to site i as
Nsi = kAsi(Cs - d} (3.9)
Thus, the local hydrodynamic effect to mass transfer is kept through the local (pore-scale) mass
transfer coefficient, while the global effect is to be computed by solving the flow field, subject to
Darcy's law in the entire pore network.
Various mass transfer correlations (analytical and numerical) have been developed over simpli-
fied geometries, e.g. for flow over a sphere, a flat plate, etc. For example, for flow over a flat plate
and typical conditions, we have (Cussler, 1984)
Sh = 0.664JRe1/25c1/3 (3.10)
where Sh = ^ is the Sherwood number, Re — ^ is the Reynolds number and Sc = jjj is the
Schmidt number. For the more complex geometries of the type shown in Figure 3.6, the literature is
50
-------
more limited. A close analogue of this geometry is mass transfer over a closed cavity. This problem
has been addressed in the context of pit corrosion in electrochemistry (for example, see Alkire and
Deligianni, 1988, Alkire et al., 1990). Experimental and numerical work has led to correlations of
the mass transfer coefficient of the form
Sh = 0.3 |'^n
M J
0.83
Pe1
,0.33
(3.11)
where we introduced the Peclet number, Pe = 3jf-, and the cavity aspect ratio W/H. Occhialini and
Higdon (1992) proceeded with another calculation of the same quantity using a spectral technique.
Their numerical results were fitted using different correlations depending on the cavity aspect ratio.
For example, we have found that
Sh = 0.1647Pe°-33
f W
for _ =
Sh = 0.1778Pe°-4
W
for - = 2
(3.12)
We have proceeded with adopting (3.12) to represent the local mass transfer rates to sites adjacent
to the interface for various geometries of this type. Wherever appropriate as dictated from a
visual observation of the pattern, expression (3.10) is also used. Of course, in real porous media,
correlations of similar form have to be developed to account for the more complex local geometries
expected.
SIMULATION OF THE EXPERIMENTAL RESULTS
Using the above pore-network simulator, we proceeded to simulate the mass transfer exper-
iments performed at various rates. Because of the 2-D nature of the micromodel, the throats
were approximated as equivalent rectangular channels. Even though the etching depth of the glass
micromodel is unknown, it can be estimated by measuring the pore volume or by the system per-
meability (see below). The geometry of the network and the width of the pores in the simulation
are the same with the computer pattern used to construct the glass micromodel, although the true
sizes of the micromodel pores are likely to be somewhat different due to various defects during the
fabrication process. Before attempting to simulate and match the mass transfer experiments, we
considered a partial validation of the model using a simple miscible displacement experiment.
(a) Validation
The validity of the pore-network simulator was tested by comparing simulated and measured
pressure drops and a simple miscible displacement experiment. First, we matched the measured
pressure drop at various rates and at conditions of single-phase flow in the absence of NAPL. This
comparison is a test of the ability of the simulator to match the overall hydraulic conductivity.
In the simulations, we used as input the measured pore widths and a uniform depth of 0.0312
cm. Measured and predicted results are shown in Figure 3.7a. As expected from Darcy's law, the
curves are linear at these flow rates. Agreement between theory and experiment is good, except
for low flow rates, where the experimental error is relatively large. A second flow rate experiment
was also conducted, in which we compared theoretical and experimental pressure drops for single-
phase flow around the trapped NAPL (for example, for the geometry corresponding to Figure 3.4).
This is a rather more stringent test, due to the more tortuous pore-space involved. The results,
shown in Figure 3.7b, also indicate a good agreement. We conclude from these two independent
51
-------
mesaurements that the flow conductances are well represented (at least on the average) with the
hydraulic expression (3.6) and (3.7).
Next, we compared experimental and theoretical predictions for the evolution of the front in
a miscible displacement process at unit mobility ratio, in the absence of NAPL, which is a test
of the ability of the pore-network simulator to simulate mass transfer. A 50% glycerol solution
was displaced by the same mixture containing a coloring dye. This mixture has a relatively higher
viscosity than water and was also used for pressure drop calculations. The following values were
assigned to the two physical parameters, fj, = 1 cp and D = 1.0 x 10~5 cm?/s. Figures 3.8 and 3.9
show a comparison of the position of the miscible front in the glass micromodel and as predicted
by the pore-network simulator, respectively, at two different times corresponding to 0.091 and 0.21
injected pore volumes, respectively. The agreement in the front positions between the two models
is good and shows that the simulator represents well the flow dynamics in the micromodel.
(b) Simulation of mass transfer experiments
The simulator was next used to simulate the mass transfer experiments. In each experiment,
the actual position of the interface of the trapped NAPL was used as input to the numerical
model. Local mass transfer coefficients were used, as appropriate as discussed above. Figure
3.10 shows the comparison between simulated and measured effluent concentration profiles in an
experiment with benzonitrile, as a function of the Peclet number, P€L. The latter was based on
the overall system size and is defined as PSL = 77, where u is the flow velocity at the entrance
port (u = Q/A, where A is the cross-sectional area of the entrance pore). We must note that in
the experiments (and the simulations) we only varied the overall injection rate, which also affects
the local mass transfer coefficient through expression (3.11)). No adjustable parameters were used
in the simulation. The agreement between experimental and numerical results is quite good for
the range of injection rates tried, which spans more than two decades. Figure 3.11 shows another
comparison between experiments and simulations using a different pattern. Here, we used the same
parameters as in the previous simulation, except that the position of the NAPL-liquid interface is
different. The agreement is equally good. Certainly, a more stringent test of the simulator would
involve comparison of local concentration profiles. However, such data were not collected. Overall,
the agreement in matching the overall behavior is encouraging of the ability of the simulator to
match important aspects of the mass transfer process under the conditions of the experiments.
A key finding from the comparison between experiments and simulation was the importance
of properly accounting for the local mass transfer rates. This is illustrated with the simulation of
Figure 3.12, where the local mass transfer from the interface to the site was taken to be by diffusion
only (using film theory with constant film thickness (Cussler, 1984)), namely here we used locally
5/i = 1. Experiments and simulations are in reasonable agreement in the region of low Peclet
numbers, where diffusion is still predominant. However, they fail to match each other at higher
flow velocities, where local convection is progressively more important, the deviation increasing
with the Peclet number. In fact, this inadequacy motivated the use of the local mass transfer
coefficients. Appropriate accounting of the local geometry in the mass transfer coefficients led to
the much more satisfactory agreement shown in Figures 3.10 and 3.11.
(c) Sensitivity study
The pore-network simulator was subsequently used to explore two issues: The derivation of
effective mass transfer coefficients for use in a macroscopic model, and the effect of the local mass
52
-------
transfer on the effluent concentration. In this paper we will focus on the second issue for the
particular experimental setup described here, the first belonging to the more general problem of
scale-up, which will be considered in a subsequent chapter.
To explore the sensitivity of the effluent concentration on the local mass transfer coefficient, we
assigned the following expression
Sh = a + bPec
(3.13)
uniformly for all pores. Parameters a, b and c were independently varied, while the Peclet number
in (3.13) is based on the microscale (which can be one to two orders of magnitude smaller than
that based on the macroscale). All simulations were performed for the micromodel and interface
geometries of Figure 3.11. The corresponding effects are shown in Figures 3.13-3.15.
Figure 3.13 shows the sensitivity to the parameter a, with values in the range 0.001 to 2. Because
the term containing a is independent of Fe, its effect reflects the importance of diffusion. We note
that the effect of a is minimal, despite the significant range in its variation (the various curves
essentially coincide in the figure). Some effect exists at relatively low velocities, where diffusion
predominates, as expected from an asymptotic study of the low rates regime, which will be published
elsewhere. At relatively large flow rates, however, the effect is negligible, as convection dominates
the local mass transfer expression (3.13) and all curves approach a limiting curve, which in a log-log
plot eventually becomes a straight line. Figure 3.14 shows the sensitivity to the parameter 6, which
is varied in the range 0.0001 to 2. In contrast to the previous, the effect is significant, and results
in a parallel shift of the effluent curve to the right, as 6 increases. This effect essentially reflects the
rescaling of the Peclet number by the factor &1/0, and corresponds to a shift of the axis in the logPe
coordinate. Finally, Figure 3.15 shows the effect of c, which was varied in the physically realistic
range (0,1). We note that in the log-log plot of the Figure, the effluent concentration approaches
a straight line asymptote, with a negative slope, — s, which decreases as c increases. In fact, the
variation of s is almost linear and obeys the relation
s — I — c
as shown in Figure 3.16. These results can be analytically explained as follows:
At steady state, the effluent concentration Ce is obtained from a mass balance
qCe = - D
(3.14)
(3.15)
where the integral is evaluated over the NAPL interface A, n is the outer normal and q is the
injection rate. Using the local mass transfer coefficient k, the local mass transfer rate under the
integral is given by (3.9). Now, at sufficiently large Peclet numbers, the site concentration is small,
namely d ,e — Ce/Cs.
At large Peclet numbers, we also have
k
DbPec
which inserted in (3.16) leads to
(3.17)
(3.18)
53
-------
as also found numerically.
The asymptotic dependence found can be used to infer an estimate of the value of the effective
parameter c in the experiments. By fitting the late part of the effluent concentration vs. Peclet
number curve in a log-log plot to a straight line we find the values of 0.45 for the experiments in
Figure 3.10 and Figure 3.11. These numbers reflect the distribution of c in the local mass transfer
coefficients used in the experiments (where the exponent c in the local coefficient varies between 0
(pure diffusion) and 0.5 (flow over a flat interface).
CONCLUSIONS
In this chapter, we conducted visualization experiments and numerical simulations in pore net-
works to understand basic aspects of mass transfer during the solubilization of residual NAPLs.
The experiments were carried out in 2-D etched-glass micromodels with randomly distributed pore
sizes. We monitored the concentration of the effluent at steady-state as a function of the Peclet
number to study the dependence of the overall mass transfer rate on injection velocity. A pore
network numerical model, based on the convection-diffusion equation using appropriate modifica-
tions for the local mass transfer coefficients, was developed to simulate mass transfer during the
solubilization of a residual phase. The pore network simulator was found to match well the exper-
imental results, provided that the local mass transfer was properly accounted for. Expressions for
mass transfer over cavities were used to represent the latter. Sensitivity studies were subsequently
conducted to investigate the dependence of the overall solubilization rate on parameters, such as the
Peclet number. A scaling law was developed for the resulting overall mass transfer rate assuming
a uniform local mass transfer coefficient. The results obtained here can be used in a macroscopic
model to predict mass transfer rates, as shown in the next chapter. Issues of scale-up and the effect
of macroscopic heterogeneities, are further discussed there.
REFERENCES
1. Abriola, L.M. and Finder, G.F., A multiphase approach to the modeling of porous media
contamination by organic compounds: 1. Equation Development, Water Resour. Res., Vol.
21, 11-18 (1985a).
2. Abriola, L.M. and Finder, G.F., A multiphase approach to the modeling of porous media
contamination by organic compounds: 2. Numerical simulation, Water Resour. Res., Vol.
21, 19-26 (1985b).
3. Alkire, R., and Deligianni, H., The role of mass transport in electorchemical pattern etching,
J. Electrochem. Soc.: Electrochemical Science and Technology, Vol. 135, No. 5, 1093-1100
(1988).
4. Alkire, R., Deligianni, H. and J.-B. Ju, Effect of fluid flow on convective transport in small
cavities, J. Electrochem. Soc., vol. 137, No.3, 818-824 (1990).
5. Augustijn, D.C.M., Jessup, R.E., Rao, P.S.C. and Wood, A.L., Remediation of contaminated
soils by solvent flushing, J. Envir. Engin., Vol. 120, 42-57 (1994).
6. Blunt, M. and King, P.R., Relative permeabilities from two- and three-dimensional pore-scale
network modeling, Transport in Porous Media, Vol. 6, 407 (1991).
7. Brandes, D. and Farley, K.J., Importance of phase-behavior on the removal of residual
DNAPLs from porous media by alcohol flooding, Water Envir. Res., Vol. 65, 869-878 (1993).
54
-------
8. Brown, C.L., Pope, G.A., Abriola, L.M. and Sepehrnoori, K., Simulation of surfactant-
enhanced aquifer remediation, Water Resour. Res., Vol. 30, 2959-2977, (1994).
9. Buckley, J.S., in "Interface Phenomena in Petroleum Recovery", Morrow, N.R. (ed), Marcel
Dekker, Inc., New York (1990).
10. Celia, M. A., Rajaram, H. and Ferrand, L.A., A multi-scale computational model for multi-
phase flow in porous media, Advances in Water Resources, Vol. 16, 81-92 (1993).
11. Chatzis, I., PRRC Report No. 82-12, New Mexico Petroleum Recovery Research Center,
Socorro, NM (1982).
12. Cussler, EX., "Diffusion: Mass transfer in fluid systems", Cambridge University Press (1984).
13. Dullien, F.A.L., "Fluid transport and pore structure", Academic Press, New York (1992).
14. Edwards, D.A., Liu, Z. and Luthy, R.G., Experimental data and modeling for surfactant
micelles, HOCs, and soil, J. Envir. Engin., Vol. 120, 23-41 (1994).
15. Feder, J., "Fractals", Plenum, New York (1988).
16. Friedlander, S.K., "Mass and heat transfer to single sphere and cylinders at low Reynolds
numbers, AIChEJ, vol. 3, No. 1, 43-48 (1957).
17. Geller, J. T. and Hunt, J.R., Mass transfer from nonaqueous phase organic liquids in water-
saturated porous media, Water Resour. Res., Vol. 29, 833-845 (1993).
18. Glinski, J., Chavepeyer, G., Flatten, J-K., An empirical relation between mutual solubilities
and interface tension for two partially miscible liquids, Physica B, Vol. 193, 154-160 (1994).
19. Guarnaccia, J.F., Imhoff, P.T., Missildine, B.C., Oostrom, M., Celia, M.A., Dane, J.G.,
Jaffe, P.R. and Pinder, G.F., Multiphase chemical transport in porous media, Environmental
Protect Agency report, contract EPA/600/5-92/002, Kerr Environ. Res. Lab., Ada, Okla.
(1992).
20. Haghighi, M., Xu, B. and Yortsos, Y.C., Visualization and simulation of immiscible displace-
ment in fractured systems using micromodels: I. Drainage, J. Colloid Interface Sci., Vol. 166,
168-179 (1994).
21. Heiba, A. A., Sahimi, M., Scriven, L.E. and Davis, H.T., Percolation theory of two-phase
relative permeabilities, Paper SPE 11015, presented at the 57th SPE Annual Meeting, New
Orleans, LA, Sept., 26-29 (1982).
22. Ho, C.K., Liu, S. and Udell, K.S., Propagation of evaporation and condensation fronts during
multicomponent soil vapor extraction, J. Contam. Hydr., Vol. 16, 381-401 (1994).
23. Imhoff, P.T., Jaffe, P.R. and Pinder, G.F., An experimental study of complete dissolution of
a nonaqueous phase liquid in saturated porous media, Water Resour. Res., Vol. 30, 307-320
(1993).
24. Itamura, M.T. and Udell, K.S., Experimental clean-up of a dense non-aqueous phase liquid
in the unsaturated zone of a porous medium using steam injection, Vol. 265, Multiphase
Transport in Porous Media, ASME, 57-62 (1993).
55
-------
25. Ji, W., Dahmani, A., Ahlfeld, D.P., Lin, J.D. and Hill III, E., Laboratory study of air sparging:
Air flow visualization, GWMR, 115-126 (1993).
26. Lenormand, R., Liquids in porous media, J. Phys.: Condens. Matter, Vol. 2, SA79-SA88
(1990).
27. Levich, V.G., Physicochemical Hydrodynamics, Prentice-Hall, Englewood Cliff's, NJ (1962).
28. Li, X. and Yortsos, Y.C., Bubble growth and stability in an effective porous medium, Phys.
Fluids A, Vol. 6, 1663-1676 (1994).
29. Li, X. and Yortsos, Y.C., Visualization and simulation of bubble growth in pore networks,
AIChE J., Vol. 41, 214-223 (1995a).
30. Li, X. and Yortsos, Y.C., Bubble growth in porous media, Chem. Eng. Sci., Vol. 50, 1247-
1271 (1995b).
31. Mackay, D.M. and Cherry, J.A., Groundwater contamination: Pump-and-treat remediation,
Environ. Sci. Technol., Vol. 23, 630-636 (1989).
32. Mayer, A.S. and Miller, C.T., An experimental investigation of pore-scale distributions of
nonaqueous phase liquids at residual saturation, Transport in Porous Media, Vol. 10, 57-80
(1993).
33. Mercer, J.W. and Cohen, R.M., A review of immiscible fluids in the subsurface: Properties,
models, characterization and remediation, J. Contam. Hydr., Vol. 6, 107-163 (1990).
34. Miller, C. T., Poirier-McNeill, M.M. and Mayer, A.S., Dissolution of trapped nonaqueous
phase liquids: Mass transfer characteristics, Water Resour. Res., Vol. 26, 2783-2796 (1990).
35. Occhialini, J.M., and Higdon, J.J.L., Convective mass transport from rectangular cavities in
viscous flow, J. Electrochem. Soc., vol. 139, No. 10, 2845-2855 (1992).
36. Pankow, J.F., Johnson, R.L. and Cherry, J.A., Air sparging in gate wells in cutoff walls and
trenches for control of plumes of volatile organic compounds (VOCs), Ground Water, Vol.
31,654-663 (1993).
37. Parker, J.C., Katyal, A.K., Kaluarachchi, J.J., Lenhard, R.J., Johnson, T.J., Jayaraman, K.,
Unlu, K. and Zhu, J.L., Modeling multiphase organic chemical transport in soils and ground
water, Rep. EPA/600/2-91/042, U.S. Environ. Protect. Agency, Washington, D. C. (1991).
38. Pennell, K.D., Abriola, L.M. and Weber, W.J., Jr., Surfactant-enhanced solubilization of
residual dodecane in soil columns. 1. Experimental investigation, Environ. Sci. Technol.,
Vol. 27, 2332-2340 (1993).
39. Pfeffer, R., and Happel, J., "An analytical study of mass and heat transfer in multiparticle
systems at low Reynolds numbers, AIChEJ, vol. 10, No. 5, 605-611 (1964).
40. Pinal, R., Rao, P.S.C., Lee, L.S., Cline, P.V., and Yalkowski, S.H., Cosolvency of partially
miscible organic solvents on the solubility of hydrophobic organic chemicals, Environ. Sci.
Technol., Vol. 24, 639-647 (1990).
56
-------
41. Powers, S.E., Abriola, L.M. and Weber, W.J., Jr., An experimental investigation of nonaque-
ous phase liquid dissolution in saturated subsurface systems: Transient mass transfer rates,
Water Resour. Res., Vol. 30, 321-332, Feb. (1994).
42. Sahimi, M. and Yortsos, Y.C., Application of fractal geometry to porous media: A review,
paper SPE 20476 presented at the 65th SPE Annual Fall Meeting, Dallas, TX (Oct. 6-9,
1990).
43. Satik, C. and Yortsos, Y.C., A pore-network study of bubble growth in porous media driven
by heat transfer, J. Heat Transfer, Vol. 118, 455-462 (1996).
44. Satik, C., Li, X. and Yortsos, Y.C., Scaling of bubble growth in porous media, Phys. Rev. E,
Vol. 51, No. 4, 3286-3295 (1995).
45. Soil, W.E., Celia, M.A. and Wilson J.L., Micromodel studies of three-fluid porous media
systems: Pore-scale processes relating to capillary pressure-saturation relationships, Water
Resour. Res., Vol. 29, 2963-2974 (1993).
46. Soil, W.E. and Celia, M.A., A modified percolation approach to simulating three-fluid capil-
lary pressure-saturation relationships, Advances in Water Resources, Vol. 16, 107-126 (1993).
47. Stalkup, F.I., Jr., "Miscible Displacement", SPE Monograph, Vol. 8, SPE, New York (1983).
48. Taylor, R. and Krishna, R., "Multicomponent Mass Transfer", John Wiley and Sons (1993).
49. Tyn, M.T. and Calus, W. F., Diffusion Coefficients in dilute binary-liquid systems, J. Chem.
Eng. Data, Vol. 20, 106-109 (1975).
50. Wakao, N., and Kaguei, S., Heat and mass transfer in packed beds, Gordon and Breach
Science, New York (1982).
51. West, C.C. and Harwell, J.H., Surfactants and subsurface remediation, Environ. Sci. Tech-
nol., Vol. 26, 2324-2330 (1992).
52. Wood, A.L., Bouchard, D.C., Brusseau, M.L. and Rao, P.S.C., Cosolvent effect on sorption
and mobility of organic contaminants in soils, Chemosphere, Vol. 21, 575-587 (1990).
57
-------
Table 1: A compilation of mass transfer correlations used in NAPL remediation.
Reference
Equation •
Valid Conditions
Friedlander, 1957"
Levich, 1962°
Pfeffer and Happel, 19646
Wakao and Kaguei, 1982°
Miller et al., 1990d
Parker et al., 1991°
Powers et al., 1992d
Guarnaccia et al., 1992d
Imhoffet al., 1993d
Geller et al., 1993d
Powers et al., 1994d
Sh = 0.89Re°-33Sc°-33
Sh = ^=Re°-5Sc°-5
Sh = bRe°-33Sc°-33
b=b(n)
S7i = 2 + l.lJ?e°-6Sc0-33
Sh = 425JRe°-750n°-6°
Sh = 1240Re°-75en°-60
Sh = 57.7[(^-0n)Re]0-6ld500-64UiOA1
Sh' =
Sh = 3400n°-87JRe°-71(^)-0-31
Sh =
Shi =
a = 4.13 ±1.01
A = 0.598 ± 0.073
Pz = 0.673 ±0.156
/?3 = 0.369 ±0.119
Pe > 1000
Re« 1
allPe
3 < Re < 3000
0.016 < 6n < 0.07
0.0015 < Re < 0.1
0.02 < en < 0.03
0.1 < Re < 0.2
On ^constant
0.012 < (-On)Re<0.2
0 < Re < 0.14
0.002
-------
Figure 3.1. Schematic of the configuration of trapped NAPL in a porous medium.
59
-------
•»«•••• •
• • ••••••••• ••••••••••••••••• « f •••
Figure 3.2. Typical pattern used for constructing a 2-D pore network.
60
-------
Camera
Syringe Pump
Figure 3.3. Schematic of the experimental apparatus.
61
-------
Figure 3.4. Photograph showing the steady-state NAPL distribution in the glass micromodel for
typical experimental conditions.
62
-------
1.2
-o.a
o>
o
o
o
§0.6
"33
CD
0.4
0.2 -
0
2 3
log(Pe)
Figure 3.5. A typical effluent curve for the 29x 15 micromodel corresponding to the solubilization
of benzonitrile by water.
63
-------
Figure 3.6. Typical NAPL interface configurations in the porespace of the micromodel. Note the
cavity-like characteristics.
64
-------
0.1 0.2 0.3 0.4 0.5
Flow rate(cc/min)
14
12
^
Q.
O
T3
1
(0
(O
ol
8
1234
Flow rate cc/min
Figure 3.7. Overall pressure drop vs. injection rate in single-phase flow in the micromodel for two
different cases: (a) In the absence of trapped NAPL, and (b) in the presence of trapped NAPL (for
example as shown in Figure 4). The solid line is the theoretical prediction, experimental data are
bracketed by error bars.
65
-------
Figure 3.8. Comparison between experimental (top) and numerical (bottom) concentration profiles
during miscible displacement of 50% glycerol solution for 0.091 pore volumes injected.
66
-------
Figure 3.9. Comparison between experimental (top) and numerical (bottom) concentration profiles
during miscible displacement of 50% glycerol solution for 0.21 pore volumes injected.
67
-------
o
•f
•s
CD
o
o
0.6
CD
"§0.4
en
0.2
0
234567
log(Pe)
Figure 3.10. Comparison between experiments and simulations for mass transfer from trapped
benzonitrile for pattern 1. The concentration at the effluent is normalized with its theoretical
solubility value. The solid line is the theoretical prediction, experimental data are bracketed by
error bars.
68
-------
11
0)
20.8
o
o
§0.6
"§0.4
_N
"5
I 0.2
0
3 4
log(Pe)
Figure 3.11. Comparison between experiments and simulations for mass transfer from trapped
benzonitrile for pattern 2. The concentration at the effluent is normalized with its theoretical
solubility value. The solid line is the theoretical prediction, experimental data are bracketed by
error bars.
69
-------
J 1
I
§0.8
I
"S
§0.6
-------
g
I
o>
o
§
•*-*
0>
^
0)
•a
CD
,N
la
10
-1
10'
10C
10°
Pe
Figure 3.13. Sensitivity analysis of the effect of parameter a in the local mass transfer coefficient
equation (3.13) on the effluent concentration curve. The effect of a is minimal at large values of
the Peclet number.
71
-------
2
"c
-------
.g
'"£
«*•*
0) _i
g 10
o
u
"c
(U
lio"2
0)
"co
10
-3
4
Pe
10°
Figure 3.15. Sensitivity analysis of the effect of parameter c in the local mass transfer coefficient
equation (3.13) on the effluent concentration curve. Note the asymptotic approach to a straight
line at large values of the Peclet number.
73
-------
Figure 3.16. Numerical simulation results for the dependence of the slope s of the effluent con-
centration curve on the parameter c of equation (3.13). The curve essentially coincides with the
theoretical result (3.14).
74
-------
Chapter 4
CONVECTIVE MASS TRANSFER
FROM STATIONARY SOURCES IN
POROUS MEDIA
C.Jia, K. Shing and Y.C. Yortsos
INTRODUCTION
Many applications involve the mass transfer of a component from a stationary source to a flowing
liquid in a porous medium. A particular example of interest to environmental remediation is the
dissolution of trapped organic phases, commonly known as NAPLs (Non Aqueous Phase Liquids),
by the injection of a solvent. The NAPL may be a single or a distributed source, for example in
the form of an organic liquid phase trapped in a low-permeability region or as disconnected ganglia
trapped by capillarity. The removal of this residual phase is typically attempted by the injection
of another fluid, which solubilizes the NAPL (for example, see Powers et al., 1994, Augustijn et al.,
1994). In the case of partial solubility, the mass transfer of the solute from the NAPL is the rate-
controlling process and needs to be determined. Then, the basic problem becomes mass transfer
by convection and diffusion in a porous medium with distributed sources. Similar problems also
arise in oil recovery from subsurface rocks, if the oil phase is trapped and components partition in
a liquid or gas phase flowing around it.
The geometric configuration of the trapped ganglia depends on a variety of factors, including
the method by which the NAPL first percolated to the subsurface, and the heterogeneity of the
soil. Figure 4.1 shows a typical configuration from the experiments reported in Chapter 3 on
the dissolution of a trapped NAPL in a glass micromodel. If a percolation process, in which
capillary forces solely control the fluid displacement is in effect, the geometry of the NAPL is
that of an invasion percolation cluster, which for a random and spatially uncorrelated medium
has a fractal structure with well-studied properties (Wilkinson and Willemsen, 1983, Feder, 1988).
More generally, the medium heterogeneity affects significantly the configuration of the source. The
resulting disorder in the source geometry and transport introduce important elements, which are
absent from classical mass transfer problems.
A common engineering approach to modeling mass transfer in this context is by assuming that
the source is extensively distributed, so that the rate of transfer from the stationary source to the
flowing liquid is similar to a first-order reaction. For unidirectional flow, for example along the
75
-------
z-direction in a cylindrical core, the following volume-average mass balance of solute in the flowing
phase is typically postulated
(4.1)
where is the porosity of the porous medium, C is a volume-averaged concentration of the solute
in the liquid, cs is the equilibrium NAPL solubility concentration, t is time, U is the Darcy velocity,
De/f is an appropriate dispersion coefficient and K is an effective mass transfer kinetic coefficient
with dimensions of inverse time. The latter is, in principle, related to the interfacial area between
NAPL and solvent per unit total volume and the local mass transfer kinetics. Equation (4.1) is
routinely used to estimate mass transfer coefficients from steady-state experiments in laboratory
cores. A partial compilation of the various results obtained with this approach is given in Table 3.1.
The wide-spread variation in the dependence of the mass transfer coefficient on various transport
parameters, and the flow rate in particular, should be noted. Given that the exponent in the flow
rate dependence reflects the underlying flow geometry, the variation shown needs to be further
explored.
Although widely used, equation (4.1) was not rigorously derived until recently, when Quintard
and Whitaker (1994) applied a volume-averaging approach in periodic media to obtain an effective
equation of a similar type. In this approach, the effective reaction-rate coefficient is obtained from
the solution in a unit cell of a closure problem involving a variant of Stokes' equation for fluid
momentum and a convection-diffusion equation for mass transfer. Quintard and Whitaker (1994)
computed the variation of this coefficient with the Peclet number for the simplified unit of Figure
4.2. In the same general context are also the studies of Koch and Brady (1987) and Mauri (1991).
In many applications, features of the microstructure at the pore-network level, such as pore-
size distribution and connectivity, may dominate the process and need to be accounted for. This
is particularly true for sources emplaced by a percolation process, as is often the case in NAPL
contamination, where the configuration can be a self-similar fractal spanning many length scales
(Feder, 1988). In such cases, a detailed analysis of the effects of the pore microstructure and the
source configuration is needed in order to obtain a better understanding of the various empirical
mass transfer correlations.
A detailed simulation of the problem in lattices of some degree of complexity is possible, for
instance with the use of lattice gas simulation (Rothman and Zaleski, 1997). Here, we will take
a simpler, although perhaps less accurate, approach, that of pore network simulation, where mass
and momentum balances are discretized across neighboring pores. This discretization allows one
to explore large-scale effects, related to the medium or the source disorder, over reasonably large
lattices. Pore-network models have been successfully used in the simulation of immiscible displace-
ments in porous media, for example in drainage and wetting (for a recent work see Bakke and Oren,
1997). They have also been successfully applied in porous media pocesses where diffusive transport
of mass or heat is important, for example in the evolution of a gas phase from a supersaturated or
a superheated liquid (Li and Yortsos, 1995, Satik and Yortsos, 1996), or in drying (Nowicki et al.,
1992, Prat, 1993). In these applications it is diffusion of mass (or heat) that forms the predominant
transport mechanism. For the case of mass transfer from a stationary source to a flowing fluid,
however, convection by the flowing fluid is controlling and must be also considered.
In Chapter 3 we reported on experimental investigations of convective mass transfer from a
trapped NAPL in the geometries shown in Figure 4.1. The effluent concentration was measured
as a function of the Peclet number. The process was simulated successfully with a pore-network
model. A systematic study of the effect of the various source geometries and of the particular
mass transfer characteristics was not attempted, however. In addition to establishing the viability
76
-------
of pore-network simulation in such problems, an important finding from the study of Chapter 3
was that the accurate matching of the experimental results required the use of appropriate mass
transfer coefficients near the source interface. For example, it was found that if the local mass
transfer from the trapped ganglia to the adjacent pores was approximated using diffusion only,
the overall transport rates were underpredicted as the Peclet number increased. Good agreement
was reached, on the other hand, when local mass transfer coefficients appropriate to the local
geometries and flow conditions were used. Recognizing that a diffusive boundary condition at the
source is associated with an "effective" medium description in the absence of microstructure, its
inadequacy to provide a good match of the experiments reflects a breakdown of the corresponding
"effective porous medium" description near the source interface (namely in places of porous media
"discontinuities").
Consider, for example, flow over a "flat plate" embedded in a porous medium (the flat plate
representing in this analogy a "flat" NAPL interface), which is a problem encountered in various
applications in porous media (e.g. see Vafai and Thiyagaraja, 1987). Although macroscopically
flat, however, the variability of the interface at the pore-scale cannot be neglected. There, both
Darcy's law for momentum transport and the conventional effective mass transfer equation lose
validity. To resolve the difficulty with the momentum transport, two different approaches have
been proposed: (i) Replacing Darcy's law by Brinkman's equation, in which a macroscopic viscous
shear is retained throughout the porous medium (done for example by Vafai and Thiyagaraja,
1987); and (ii) Retaining Darcy's law, but introducing an appropriate boundary condition at the
discontinuity. The latter is obtained from an analysis of the local problem (e.g. see Tio and
Sadhal, 1994, and references therein). The two approaches are similar in that Brinkman's equation
essentially reduces to Darcy's law everywhere except within a narrow interface boundary layer,
which scales with the square root of permeability (and it is of the order of a typical pore size, Katz
and Thompson, 1986).
Analogous problems should exist with the transport of heat or mass near a discontinuity, where
the effective equation formalism cannot capture the detailed transport process, particularly at high
flow rates. However, the approach typically taken in the literature is to modify only the equation
for momentum transport, while retaining an effective equation for the transport of heat (or mass)
(e.g. see Vafai and Thiyagaraja, 1987). This approach is commonly used in the study of natural
convection in porous media (e.g. see Gouyeau et al., 1996, for a recent application). Given the
predominant contribution of the heat conduction in the solid, this approximation may lead to
acceptable results in heat transfer problems. However, it is unclear whether it also applies for mass
transfer, particularly when convection is dominant. Here, the local mass transfer near the source
interface must be properly accounted, particularly at larger values of the Peclet number.
The objective of this chapter is to study mass transfer in porous media from stationary sources
of various geometries in which the effects of the porous medium and the source geometry are
explored using a pore-network representation, and where to account for the Stokes flow dynamics
near the interface, we will use a local modification of the mass transfer coefficient. Our motivation
is to improve the understanding of the dependence of overall mass transfer rates on the various
parameters shown in Table 3.1, and particularly on the flow rate. We must remark, in advance,
that some of these correlations may reflect the combined effects of dissolution and mobilization
of residual phases. In this chapter, only dissolution from stationary sources will be considered.
Deriving an effective equation similar to (4.1) will not be attempted. However, equation (4.1) will
be used to define effective Damkohler and Sherwood numbers for the case of distributed sources,
and to investigate their dependence on various parameters.
The chapter is organized as follows: First, we present the mathematical formulation, which
includes a modification of the local interface condition. Its effect is explored by analyzing two
77
-------
different cases: Plow over an isolated source of macroscale size, such as flow over a flat plate, a
sphere and over distributed sources. The first geometry allows for an analytical solution, while for
the other geometries, asymptotics at small and at large Peclet numbers are derived. Subsequently,
we consider the more general problem of mass transfer over distributed sources of microscale size.
Asymptotic and numerical solutions are derived. Pore-network simulations are conducted to test
the theoretical predictions. The flow over self-similar surfaces is particularly analyzed. Using
overall mass balances, the effective Sherwood number corresponding to equation (4.1) is computed
for various geometries. The results are then compared to the various experimental correlations.
MATHEMATICAL FORMULATION
Consider mass transfer from a stationary source to a fluid flowing in the pore-space of a porous
medium. As pointed out previously, the source geometry may be an effective plat plate, or an
effective sphere, a self-similar percolation cluster or a distributed, disconnected interface, such as
shown in Figure 4.1. Mass transfer in the pore-space is described by the usual convection-diffusion
equation, the steady-state of which reads
u-Vc = DmV2c
(4.2)
where Dm is molecular diffusivity and u is the local velocity. At the creeping flow conditions
assumed, Stokes' equations describe the momentum balance. The solute is assumed at a sufficiently
small concentration to not affect properties of the mixture, such as density or viscosity. At the
source interface, the solute concentration is equal to its thermodynamic solubility value
c = c.
(4.3)
These equations are subject to the appropriate flux or concentration conditions at the far-field
boundaries. Typically, we assume uniform flow with zero net flux at the inlet, a Dankwerts boundary
condition at the outlet and no-flux boundary conditions at the lateral boundaries. For flow over a
flat plate, a sphere or a self-similar surface, the medium is assumed infinitely large.
As previously noted, our objective is to study the solution of the above problem using a pore-
network analogue or its corresponding effective equation analogue, where appropriate. In general,
the first approach must be taken for the case of sources of a size comparable to the microscale /,
while an effective medium description is possible for macroscale sources. For the latter, and at
conditions of low volumetric content of the source (for example, in the case of flow over a flat plate,
etc.) the effective equation corresponding to (4.2) is a volume-averaged mass balance
U-VC = V(De//-VC) (4.4)
Here De// is a velocity-dependent dispersion tensor, and U is the Darcy velocity satisfying
1?
U=--VP
(4.5)
where k is the permeability of the porous medium and fj, represents fluid viscosity. Alternatively,
the full problem can be discretized by assuming a pore-network description.
Pore network simulation is a useful tool for processes in pore-networks or in experiments in-
volving glass micromodels. The key approximation is the porous medium representation in terms
of pores, connected to each other by pore throats in a network-like fashion (Dullien, 1992). In the
network, sites represent pores and provide volumetric storage, while bonds represent throats and
78
-------
control the conductance to flow, heat and mass. The various geometrical and transport charac-
teristics can be spatially distributed. The governing equations are discretized over the sites of the
network. For example, velocity and pressure fields are determined using the overall mass balance
£Qy=0 (4.6)
3
where Qij refers to volume flow rate between adjacent sites i and j, and the sum is over all sites j
neighboring site i, and the use of Poiseuille's law
Qij=Gij&Pij (4.7)
where the overall conductance, Gij, can be distributed. For mass balance purposes, concentrations
are assigned to sites only, and the mass transfer between adjacent sites (except those neighboring
the interface, see below) is assumed to occur by conventional modes of diffusion and convection.
Then, the mass balance for solute in site i reads
Vt
AC,-
At
DmAi
AC;.
(4.8)
3 L ' J j L J
where Vs denotes the site volume and AJJ is the cross-section of the pore throat. If site i neighbors
a source, however, a modification for local mass transfer is necessary, as noted above.
Local mass transfer coefficients
In a boundary layer near the interface (of the order of Vk) , the above description breaks down
and needs to be replaced. Indeed, local mass transfer cannot be described simply by diffusion and
convection using Poiseuille's law for fluid flow (or Darcy's law, in the respective effective continuum
case), but requires coupling Stokes flow to convection and diffusion. This is not inherent to the pore
network approach (or to an effective medium) used here. One alternative is to keep Poiseuille's law
(or Darcy's law in the effective problem) for the hydrodynamics outside this boundary layer and
modify the mass flux, by using a mass transfer coefficient h[,
Nsi = hi(cs-Ci) (4.9)
where, Ns{ is the mass flux from the interface to an adjoining site i. The corresponding modification
of the boundary condition for the effective medium problem is then
hi(c, - C) = -De//-VC-n (4.10)
where n is the unit normal vector to the interface pointing outwards. We note that in the effective
medium problem, equation (4.10) would be evaluated at the edge of the boundary layer, which for
all practical purposes we may take to coincide with the source interface itself. (This assumption
works well provided that the source volumetric density is not very high).
The coefficient hi reflects the local mass transfer appropriate to the particular geometry. Ex-
amples of various local configurations expected in a NAPL-solvent system are shown in Figure 4.3
reprinted from the experiments of Chapter 3. In general, hi is related to the local flow velocity at
the edge of the boundary layer, given here by the outer (or effective medium) solution, using mass
transfer correlations appropriate to the local geometry. For simple interface geometries, various
mass transfer expressions (for example, flow over a flat plate, etc.) have been developed (Cussler,
79
-------
1984). For the more complex geometries of the type shown in Figure 4.3, however, the literature is
rather limited. A close analogue to some of these is mass transfer over a cavity, which is a problem
in electrochemical pit corrosion studied in some detail by Alkire and Deligianni (1988) and Alkire
et al. (1990) (see also Occhialini and Higdon, 1992). These authors developed correlations of the
following type
Sh; = 0.3 (
0.83
pe0.33
(4.11)
where Sh; = ^- denotes the local Sherwood number, Pe, = ^ denotes the local Peclet number,
and w/h is the"cavity aspect ratio. This analogy was successfully used in the pore-network simu-
lations in Chapter 3 to match the etched-glass micromodel experiments. More generally, however,
we will take in this chapter the general dependence (see also Browning and Fogler, 1996)
Shj = A + bPef
(4.12)
This relationship contains the expected asymptotic behavior of the local mass transfer rates in the
two limits, namely the approach to a constant at small Pe/ and to a power-law at large Pej. In (4.12)
the dependence of mass transfer on the local flow velocity is through exponent c, parameters A and
6 being functions of the geometrical and physical properties only. In general, all these parameters
can be spatially distributed, reflecting the pore geometry disorder (Figures 4.1 and 4.3), hence Sh;
is in general position-dependent. Using a characteristic macroscale velocity U0, we can rewrite
(4.12) as
Shi = A + b
u
U~o
Pef
(4.13)
where an obvious notation was implied for Pei0. The new boundary conditions (4.9) and (4.10),
respectively, will be used to simulate the mass transfer problem.
EFFECTIVE POROUS MEDIA: MACROSCALE SOURCES
To obtain an insight on the effect of the new interface condition on the overall problem, we
will first analyze the problem in a homogeneous effective medium with sources of a macroscopic
size. Pore-network simulations will be described in a later section for mass transfer with distributed
sources of a microscale size. Using for characteristic length the macroscale of the source, L, for char-
acteristic velocity the free-stream velocity, U0, and for characteristic concentration the equilibrium
value, cs, the steady-state mass transfer equation for an effective medium reads in dimensionless
notation
V(DD,e/rVC)
(4.14)
where subscript D denotes dimensionless quantities. The velocity Up satisfies potential flow (for
the case of a homogeneous medium), Do,eff is a dimensionless velocity-dependent dispersion tensor
and we introduced the macroscale Peclet number Pei, = ^fo, where £>e//,o is a measure of the
transverse-dispersion coefficient at the velocity U0. We make the following remarks:
1. For an effective continuum description to be valid, the ratio of length scales, j, must be
sufficiently large.
2. The dependence of De//,0 on the velocity can be approximated for the present purposes by
80
-------
DB//t0 =
\Pei]
(4.15)
where A is an O(l) positive constant. Captured in (4.15) are the two asymptotic limits of molecular
diffusion, at low velocities, and of mechanical dispersion, at high velocities. In the latter, we used
the well-known result that for an uncorrelated and random porous medium, the dispersivity is
proportional to the microscale / (e.g. see Koch and Brady, 1987).
3. As a result of (4.15), the macroscale Peclet number can be expressed as
UnL
Pe;
(4.16)
Thus, for macroscale problems (L > /), Pe^ is generally large, unless Pe; is sufficiently small (and
of the order of ^).
Finally, in terms of this dimensionless notation, the boundary condition corresponding to (4.10)
becomes
• -CD) = -]
where we introduced the macroscale Sherwood number,
L
I
DmSh,
(4.17)
(4.18)
- Deff,0 \IJ De/f,0
Because of its dependence on the local Peclet number Pe/, Sh^ is in general spatially variable.
This introduces a significant complexity in the analysis of the problem. However, certain source
geometries, and specifically a flat plate, allow for an analytical solution. In this section, we will
discuss in order flow over a flat plate and flow over a sphere and develop analytical and asymptotic
results. In these three geometries, the source is embedded in an infinite medium. An asymptotic
analysis of the more general problem, in which the source is distributed throughout the medium,
will also be provided subsequently.
Flow over a flat plate
Consider flow over a flat plate embedded in an infinitely large porous medium. Because the
velocity profile is uniform (Ur> = 1) the effective dispersion coefficient and the local Sherwood
number are spatially constant. Of interest here is mass transfer at relatively large Pez,, thus we
will keep convection in the streamwise direction, £, and dispersion in the transverse direction, y.
According to (4.16), large Pe^ occurs under two conditions: When Pe/ is O(l) or larger, or when
Pe/ < 1 and simultaneously £ < -^ < } (and where the effective medium condition L > / was
tacitly implied). Then, the equation reads
dCD
d2CD
dxD ~ PeL c
which is to be solved subject to the interface condition
ShL(l -
dCD
at
(4.19)
(4.20)
and the boundary conditions
81
-------
— 0 at XD — 0 and yD = co
The solution of this problem is readily obtained using a Laplace transform. We find
(4.21)
CD = -exp(ShiyD + « ED) erfc I CI
+ erfc
(4.22)
where we introduced the parameter a = Sh^Pe^1 . From this we can evaluate the effective mass
transfer coefficient over a flat plate of length L,
Shav =
avg
= ShLf(expa2erfca - 1 + -
"
(4.23)
ef/i0 o oyD
Equation (22) shows that the normalized overall mass transfer rate (namely the ratio
is a function of the parameter a only. This function is plotted in Figure 4.4. One distinguishes two
limits:
(i) At small a, we have the asymptotic behavior
= ShL
fl -
\_
0(a2)
(4.24)
1
Using (4.16) and (4.18), we find that this is possible under the conditions Pe/ > 1 and y yPef
1, a necessary constraint for which is c < 1. The latter is more difficult to be satisfied as c
approaches 1. In this limit, the effective mass transfer coefficient has to leading-order in a the
same scaling with the Peclet number as the local mass transfer coefficient, and reflects only local
geometry effects. The overall mass transfer is dominated by the local mass transfer at the source
interface and follows the velocity scaling of the local coefficient, namely havg — h\ to leading order.
Indeed, here, the interface concentration is to leading-order proportional to a, thus it vanishes as
a-^-0.
(ii) At large a, on the other hand, we obtain from (4.23)
QI ^ -p 1/2 , x-j / —IN (A oc\
•V/H"
This result arises when Pe/ is generally finite. Now, mass transfer is dominated by the external
problem. Indeed, in this limit, the interface boundary condition reduces to CD — 1, and the
problem admits the classical erfc solution. Here, the overall mass transfer is dominated by the
macroscale geometry and is insensitive to the local mass transfer.
These two limits correspond to two different mass transfer regimes in the macroscale problem,
namely where local or global characteristics dominate, respectively. We note that the local mass
transfer characteristics dominate the overall mass transfer rates only in the case of sufficiently large
Pe/. Then, the coefficient havg has the following scaling
havg ~ VI , Pe/ > 1 (4,26)
In all other cases, the global problem masks the mass transfer characteristics at the microscale, and
the conventional effective continuum description with boundary condition (4.3) applies. Now, the
overall coefficient has the different scaling
^Deff>0U0 , Pe/finite
(4.27)
82
-------
which, if approximated by a power-law in U0 would lead to an exponent ranging between 0.5 and 1.
These conclusions also reflect the dependence of the dispersion coefficient on velocity. For example,
in the absence of the latter (A = 0) , the first regime would require the more stringent condition
Pep1/2 < 1, namely c < 1/2.
Numerical results for Shavg for flow over a flat plate are shown in Figure 4.5 for various depen-
dences of Shi, on the local Peclet number. It must be pointed out that in this Figure both Shavg
and Pe£ are normalized with Dm instead of Defft0. Essentially, the vetical axis scales with havg,
while the horizontal with U0, hence we should expect scalings of the form (4.26 and (4.27) to arise.
Indeed, it is shown that there exists a segment common to all three curves, with a scaling close to
(4.27), while at larger values of Peclet number, equation (4.26) is obeyed, as expected.
Flow over a sphere
Consider, next, flow over a sphere embedded in an effective porous medium. Working in r and
9 coordinates, the potential flow solution is (see for example, Lamb, 1932)
1 _
r3
TD
and
(4.28)
(4.29)
where we have normalized lengths with the sphere radius R. The interface boundary condition
— 1) now reads
(1 - CD)ShL(smO) = -VD,effVCD-nr (4.30)
where we have explicitly expressed the dependence of Shi, on the interface velocity UD,B = — fsin#.
An analytical solution to this problem is generally difficult. Instead, we will focus on the asymptotic
behavior of the problem in the limits considered above.
At relatively large Pe^, (namely for any non-negligible value of Pe;), a boundary layer of or-
1/2
der PeL' develops. Neglecting longitudinal dispersion, as well as the variation of the dispersion
coefficient with velocity inside the boundary layer, the mass transfer equation reduces to
37/y
dCD
dY
(4.31)
in terms of the streamwise coordinate 77 = cos#, —1 < 77 < 1, and the boundary layer coordinate
Y =?&£ (r — I). In these variables, the interface condition at Y = 0 becomes
(4.32)
Even under these conditions, the exact solution of this problem is difficult since parameter a =ShiPe£ '
is position-dependent. However, the leading-order dependence can be extracted in a relatively
straightforward fashion, as follows.
(i) If a is uniformly small, the problem will be dominated by the local mass transfer, as in the
case of the flat plate. Then, the leading scaling is expected to be
83
-------
i /
2. J-i
(4.33)
which is the equivalent of the flat plate expression (4.24). Using equation (4.13) for the evaluation
of Shjt,, we obtain after some manipulations
Sh
DmL
atw«-=r=-r(A + &'Pe?0) (4.34)
whereft' = 6(|)C/-l(vT^
and F(z) is the Gamma function. We note that
the local power-law dependence on the velocity (exponent c) has survived, although now with a
modified prefactor. The previous conditions for a flat plate (but now in terms of Pe/0) also hold for
the applicability of this result.
(ii) In the opposite case of large a, the overall mass transfer will be dominated by the external
problem. Now, the boundary condition becomes CD = 1 and the problem can be solved with
V 4 A _(_ §.« _ i
-
V
classical methods, namely using the similarity variable £ = -rr> where g(ij) =
find
z\
CD = erfc£
thus, the average mass transfer coefficient is
Shavo = -Per
'•avg
1/2 f1 9CD
'/
J-i
dY
Y=0drj =
\X2~7r
We
(4.35)
(4.36)
This power law has the 1/2 scaling, which is identical to the equivalent expression (4.25) for the
flat plate problem.
Numerical results for Shavg for flow over a sphere are shown in Figure 4.6 for various dependences
of Shi on the local Peclet number. As in the simulations of Figure 4.5, both Shavg and Pej, are
normalized with Dm instead of Defft0. We should expect, therefore, a common segment with a
slope of about 1/2, in the log-log plot, followed by a slope equal to c, at large velocities. These
results are well obeyed in the respective limits, as expected.
For completeness, we also present the asymptotic behavior in the other limit of small Peclet
numbers, where the problem is diffusion-dominated. Because this geometry (as well as that of the
self-similar surface to follow) corresponds to a finite source embedded in an infinite medium, the
Sherwood number should approach a constant. Now, the concentration profile becomes simply
CD =
Shi.,0
where
= limpei_).o
(ShL,0 + l}rD
, and the dimensionless mass transfer coefficient is
., 2ShL,0
(4.37)
(4.38)
This expression reduces to the classical diffusion-only result for a sphere, Shavg — 2, when we invoke
the continuum limit L ~^> I, in which case Sh£,i0 » 1 (compare with (4.18)).
Asymptotic Behavior With Distributed Sources
84
-------
In the previous section we discussed flow over a finite source embedded in an infinite medium.
In the more general problem, the source is distributed throughout the volume. We should expect a
similar behavior at large Pe;, but a different behavior at small Pe/. We note, again, that although
distributed, the source volume fraction cannot be taken too large if the present method is to be
valid. Indeed, the assumption was made in this chapter that we can capture the effect of Stokes
dynamics in boundary layers around the source. For this to be valid, layers from adjacent sources
should not overlap, which requires a sufficiently small volume fraction. Again, we will examine the
mass transfer behavior in the two limits of small and large Peclet numbers, respectively.
Consider, first, the asymptotic solution at small Peclet numbers, where diffusion predominates,
and take the expansion
C = C0
(4.39)
where n is a positive exponent to be determined and we have used the extent of the domain L,
as the characteristic length scale. Here, and in the following, we will suppress subscript D for
simplicity in notation. Note that from physical considerations it follows that C\ < 0. Substitution
of the expansion into equation (4.14) gives
U-[VC0
(4.40)
where, for simplicity, we omitted the dependence of the dispersion coefficient on velocity at this
level. At the interface boundary, the condition reads as
(ShL,0 + PeLShLil
Co -
8(C0
dn
(4.41)
where we expanded Sh^ in powers of Pe^, while at the upstream boundary x — 0, we have
C0
e
Pon-J
PeL Ox L
dx
= 0
(4.42)
where x is the streamwise coordinate. Finally, at the downstream and the lateral boundaries
dn
1 = 0
(4.43)
The exponent n is determined by applying in (4.40)-(4.43) a dominant balance. It is easily shown
that n = 1, hence the leading-order term satisfies Co = 1, everywhere, which is expected at
conditions of diffusion control. The next term in the expansion satisfies the steady-state diffusion
equation
(4.44)
V2Ci - 0
with upstream boundary condition @J£- = 1, with interface condition at the source boundary
Sh^oCi = ^- and with the no-flux boundary conditions ^- = 0 on all other boundaries. The
solution of this problem depends on the source configuration.
Substitution of the above results in the definition for Shavg leads to the following for the mass
transfer coefficient in this limit
- f Un
A.e JAe
(4.45)
85
-------
where a\ = *§*• is the ratio of the cross-sectional (entrance or exit) area to that of the area of the
sources, Un is the normal component of the velocity at the exit and we made use of the relation
- / d-jj-dA = PeL f Un
JAi on JAe
CedA
(4.46)
We conclude that in the small Fez, limit, Shavg is to leading-order proportional to Pe^, reflecting
the finite size of the domain used. This is to be contrasted to the previous problem of transport
over a source of finite size, embedded in an infinite medium, where Sha,,5 approaches a constant.
To obtain the next-order term in the expansion requires solving the boundary value problem for
Ci, which must be obtained numerically, in general. Because C\ < 0, however, we can infer a priori
that the second-order correction would be negative.
In the other limit of large Pe/, the concentration field practically vanishes everywhere in the
domain, except in narrow boundary layers near the source interfaces, where the local mass transfer
dominates. The problem becomes analogous to flow over a finite source in an infinite medium,
discussed in the previous sections, hence in accordance with the previous analysis, we should expect
the scaling
Sh
avg
V
or
-Pef,
'lo
(4.47)
(4.48)
for the average mass transfer coefficient. These predictions will be tested numerically below.
PORE NETWORKS: MICROSCALE SOURCES
The pore-network simulator was subsequently used to evaluate mass transfer coefficients and
their dependence on the flow parameters and the geometrical configuration of the source. A de-
scription of the simulator was given in Chapter 3. Here, both Shavg and Pe^ were defined based
on a macroscale L (equal to the size of the domain) and on the molecular diifusivity (rather than
dispersion), essentially corresponding to A = 0 in (4.15), (4.16) and (4.18). Dispersion is inherent
in the solution of the problem from the variation of the velocity in the pores of the pore network.
We also note that the pore-network approach can also be viewed as a finite-difference analog of
the effective equations (although it stands on its own right as a model of real physical system).
As a result, some of the asymptotic predictions of the continuum problem may also apply here as
well. The following source geometries were investigated: a Koch surface, a percolation cluster and
a periodically distributed source.
Flow over a self-similar surface
The study of self-similar objects is important in the context of porous media, since the source
geometry may very well have self-similar characteristics as a result of the method of emplacement
(for example, if a percolation process has controlled the spatial distribution of the source). Self-
similar objects are characterized by a cascade of length scales, ranging from an upper cut-off, Lu, to
a lower cut-off, L;. As a result, length-scale dependent quantities, for example the Peclet number
86
-------
or the Reynolds number need to be carefully defined. In accordance with the previous, we define
the average mass transfer coefficient as
oiio.'!//} —
Lavg
-
A
dCD
dA
(4.49)
where Ai refers to the total surface area of the source and lengths have been dimensionalized with
the upper cut-off Lu (which here coincides with the linear extent L). In this section, we consider the
solution of the mass transfer problem over the geometry of a Koch curve (Figure 4.7) (Mandelbrot,
1983, Feder, 1988). A corresponding problem involving a percolation cluster, which is quite relevant
to remediation problems, and where LI coincides with I, will be discussed later in conjuction with
the effective equation (4.1). As in the previous geometries, we will discuss the asymptotic behavior
of the problem in the two limits of small and large Fez,, respectively.
At small Pei, (namely small Pe;), diffusion predominates and at all scales the problem is steady-
state diffusion from a source of constant concentration, although of a self-similar geometry, into
a quiescent medium. The particular problem of diffusion from a self-similar surface was studied
by Gates and Witten (1987), who showed that the overall mass transfer rate is equivalent to that
for diffusion from an equivalent homogeneous source of size equal to the effective radius of the
self-similar source. Their conclusions were confirmed in a related study of mass transfer-driven
bubble growth in porous media by Satik et al. (1995). We should expect, therefore, that at small
Peclet numbers, the overall mass flow rate from the Koch curve will be equivalent to that from a
flat plate, that from a percolation cluster to be equivalent to diffusion from an equivalent sphere,
etc. Using definition (4.49), it follows that in this limit we should have the scaling
(4.50)
(4.51)
if the source is embedded in an unbounded medium, or
Sh
avg
Fez. f- (l + 0(PeL))
if the bounding medium is finite (for example, compare with (4.45)). Here D is the surface fractal
dimension (and E is the Euclidean dimension embedding the object). For the Koch curve of Figure
4.7, we have E — 2 and D — In5/ln3=1.4645, while for a 3-D percolation cluster (Stauffer and
Aharony, 1992) D « 2.5. Note that because the exponent in (4.50) is negative, the so defined ShOV5
decreases as the surface roughness (namely L/L{) increases. Equation (4.50) states that in the limit
of small Pei,, the average Sherwood number is a power law of the ratio of the characteristic length
scales. This means that under otherwise identical conditions, Shavg is size dependent.
As Pei, increases, the situation is more complex. Figure 4.8 shows a plot of the velocity vectors
corresponding to flow over the Koch curve of Figure 4.7. The decrease of the velocity magnitude
as cavities of lower accessibility are probed is apparent. The velocity inside the cavities is smaller,
as the cavity size decreases, thus even though mass transfer may be controlled by convection at
the external surface, it becomes progressively diffusion-controlled inside smaller size cavities. At
sufficiently large Pet, however, we expect that the system would eventually approach a state in
which the local mass transfer coefficient controls the mass transfer rate. This is illustrated in Figure
4.9 which shows a plot of the local mass fluxes along a coordinate tangential to the Koch curve for
various values of the Peclet number. For most of the cases shown, the local flux is quite low in the
inaccessible areas, but increases rapidly in the part of the surface which is exposed. At sufficientlly
large Pei,, however, all fluxes increase and may eventually equalize, reflecting the increased role of
87
-------
convection inside the surface. The corresponding dimensionless concentration profile is plotted in
Figure 4.10: along the streamwise coordinate x on the top of the surface of Figure 4.7, in Figure
4.10a, and along the coordinate tangential to the perimeter in Figure 4.10b. Profiles at different
Peclet numbers are shown. The concentration at the top reflects the existence of the various surface
cavities, but it is otherwise relatively smooth. On the other hand, the profile along the cavities
of the self-similar surface fluctuates significantly, particularly at larger Pe/,. The concentration is
smaller near the surface exposed to flow, but increases substantially inside the cavities, reflecting the
fact that mass transfer by convection is substantially slower inside inaccessible cavities. However,
for a sufficiently large Peclet number, given a fixed ratio of cutoff scales, the interface concentration
becomes small almost everywhere, reflecting the fact that mass transfer is convection-controlled
even at the small scale. We should note that these simulations correspond to the case 6=0,
namely the local Sherwood number is independent of the local velocity.
The transition from a diffusion-dominated to a convection-dominated regime occurs when the
minimum local Peclet number inside the surface, namely the minimum in ^-, becomes sufficiently
large, where u\ denotes the tangential velocity at the surface. For a fractal surface, we expect this
minimum to scale as a power of the ratio of the length scales ^-, and more specifically
mm-
(4.52)
where a is a negative exponent. Then, the transition between the two regions should be controlled
by a scaling function /(z), in terms of the variable z = U0 (j^f -^ = Pe;0 (zr)"- The properties
of this function are discussed below.
At sufficiently large Pe^, the local mass transfer dominates, and using equation (4.49) we have
Sh
avg
(4.53)
For the same reasons cited above we should expect that the geometric factor above also scales as
a power of the ratio of the two scales, namely
(4.54)
^ '
-J-/ (%.
AiJAi\U0
where J3(c) is another exponent that depends on c (/3(0) = 0). Substitution in (4.53) gives the
result
-J (4.55)
The transition from diffusion- to convection- control is modulated by the scaling function. For
the case of a bounded embedding medium, equation (4.51), small and large Pe^ behaviors may be
combined using the following scaling ansatz
Sh^ = PeL - fc(z)
\J
where the scaling function fc(z) has the behavior
(4.56)
/c(0) ~ const and fc(z) ~ zc l as z —>• oo (4-57)
This scaling works if a = —D and /3(c) = etc. Equation (4.57) reflects the fact that given a ratio
of upper to lower cutoff scales, a sufficiently large Peclet number exists, such that the problem
88
-------
becomes convection-controlled. Inversely, for a fixed macroscale Peclet number, self-similar objects
of sufficiently many generations eventually become diffusion-controlled. Although it appears to
collapse the two behaviors correctly, however, a more robust scaling theory must be developed.
This requires further investigation.
Figure 4.11 shows a plot of Shavg as a function of Pe^, for b = 0 and for three different Koch
curves corresponding to different generations. The curves have the general behavior expected from
the previous analysis. Indeed, the asymptotic behavior shown is consistent with (4.51) and (4.55).
The low Pe£ scaling with the scale ratio was tested and found to obey (4.51), as anticipated. Figure
4.12 shows corresponding results for c = 1/3, which also support the previous scalings. We are stil
in the process of investigating the validity of the scaling function, however.
Numerical results for flow around the percolation cluster of Figure 4.13 are shown in Figure 4.14
for three different values of the local exponent c. The results confirm the validity of the previous
asymptotic expressions in the two limits of small and large Pe^,. The scaling for mass transfer on
the available surface area similar to (4.51) was also verified for the percolation cluster problem.
The numerical results for flow over distributed sources at the microscale, will be analyzed in
relation to the effective equation (4.1). For this it is necessary to establish a connection between
Shavg and the corresponding mass transfer coefficients. To proceed we note that Shavg was calcu-
lated as defined in (4.49). The overall mass flow rate, J, is then, J = DeffL 'CsSha,,g. We also note
that Shavg is related to the dimensionless effluent concentration CD,S through
Shavg = ^CD>ePeL (4.58)
j\i
a result obtained from a mass balance. Note that the ratio Ae^efj is essentially independent of L.
Equation (4.58) will be used to evaluate the Damkohler number and the effective mass transfer
coefficient corresponding to equation (4.1).
Effective Thiele Modulus and Mass Transfer Coefficient
In dimensionless notation with length scaled with the extent of the domain L and concentration
with the solubility value cs, the steady-state boundary value problem at constant fluid injection
rate U0 reads as
with boundary conditions
and
CD -
1 rif^1
D
dxD
=0 at XD = 0
= 0 at xD =
(4.60)
^ '
(4.61)
Here we introduced the macroscopic Peclet number Pei, — ^^- and we also defined the Damkohler
number Da= -jj-. Parameter K can also be expressed in terms of an effective mass transfer
coefficient, he/f, such that K = hef/avi, where ay; is the source surface area per total volume.
This leads to an effective Sherwood number, She// = —e£f _ . Note that here we used the source
e//
-------
size for characteristic length scale. For example, Quintard and Whitaker (1994) have used an
analogous definition in their study of mass transfer. This coefficient is not equal to the Sherwood
number, Shot,5 defined previously, although the two can be related as shown below.
We will use the results of the simulation of the full problem, as well as the asymptotic dependence
of Shat/5 discussed previously, to relate Da and She// to Pe£, and the other parameters of the
problem. This step is important for comparison against experimental correlations. To do so we
will use equation (4.58), substitute numerical or asymptotic results for Shavfl from the full problem
and an expression for the effluent concentration, CD,S = CD(XD — 1), by solving (4.59)-(4.61). The
solution of this problem can be obtained analytically. The effluent concentration is readily found
to be
CD>e = 1 -
2mexp
m2)sinh
2mcosh
(4.62)
where m =
-------
9. Feder, J., "Fractals", Plenum, New York (1988).
10. Friedlander, S.K., Mass and heat transfer to single sphere and cylinders at low Reynolds
numbers, AIChEJ, vol. 3, No. 1, 43-48 (1957).
11. Geller, J. T. and Hunt, J.R., Mass transfer from nonaqueous phase organic liquids in water-
saturated porous media, Water Resour. Res., Vol. 29, 833-845 (1993).
12. Gouyeau, B., Songbe, J.-P. and Gobin, D., Numerical study of double-diffusive natural con-
vection in a porous cavity using the Darcy-Brinkman formulation, Int. J. Heat Mass Transfer,
Vol. 39, No. 7, 1363-1378 (1996).
13. Guarnaccia, J.F., Imhoff, P.T., Missildine, B.C., Oostrom, M., Celia, M.A., Dane, J.G.,
Jaffe, P.R. and Pinder, G.F., Multiphase chemical transport in porous media, Environmental
Protect Agency report, contract EPA/600/5-92/002, Kerr Environ. Res. Lab., Ada, Okla.
(1992).
14. Imhoff, P.T., Jaffe, P.R. and Pinder, G.F., An experimental study of complete dissolution of
a nonaqueous phase liquid in saturated porous media, Water Resour. Res., Vol. 30, 307-320
(1993).
15. Katz, A.J. and Thompson, A.H., Quantitative prediction of permeability in porous rock,
Phys. Rev. B, Vol. 34, 8179-8181 (1986).
16. Koch, D.L. and Brady, J.F., A non-local description of advection-diffusion with application
to dispersion in porous media, J. Fluid Mech. Vol. 180, 387-403 (1987).
17. Lamb, H., "Hydrodynamics", Dover, New York (1945).
18. Levich, V.G., "Physicochemical Hydrodynamics," Prentice-Hall, Englewood Cliffs, NJ (1962).
19. Li, X. and Yortsos, Y.C., Visualization and simulation of bubble growth in pore networks,
AIChE J., Vol. 41, 214-223 (1995).
20. Mandelbrot, B., "The fractal geometry of nature", Freeman, San Francisco (1983).
21. Mauri, R., Dispersion, convection and reaction in porous media, Phys. Fluids A, Vol. 3, No.
5, 743-756 (1991).
22. Miller, C. T., Poirier-McNeill, M.M. and Mayer, A.S., Dissolution of trapped nonaqueous
phase liquids: Mass transfer characteristics, Water Resour. Res., Vol. 26, 2783-2796 (1990).
23. Nowicki, S.C., Davis, H.T. and Scriven, L.E., Microscopic determination of transport param-
eters in drying porous media, Drying Technology, Vol. 10, 925-946 (1992).
24. Occhialini, J.M. and Higdon, J.J.L., Convective mass transport from rectangular cavities in
viscous flow, J. Electrochem. Soc., vol. 139, No. 10, 2845-2855 (1992).
25. Parker, J.C., Katyal, A.K., Kaluarachchi, J.J., Lenhard, R.J., Johnson, T.J., Jayaraman, K.,
Unlu, K. and Zhu, J.L., Modeling multiphase-organic chemical transport in soils and ground
water, Rep. EPA/600/2-91/042, U.S. Environ. Protect. Agency, Washington, D. C. (1991).
26. Pfeffer, R. and Happel, J., An analytical study of mass and heat transfer in multiparticle
systems at low Reynolds numbers, AIChEJ, vol. 10, No. 5, 605-611 (1964).
92
-------
27. Powers, S.E., Abriola, L.M. and Weber, W.J., Jr., An experimental investigation of nonaque-
ous phase liquid dissolution in saturated subsurface systems: Steady state mass transfer rates,
Water Resour. Res., Vol. 28, 2691-2705 (1992).
28. Powers, S.E., Abriola, L.M. and Weber, W.J., Jr., An experimental investigation of nonaque-
ous phase liquid dissolution in saturated subsurface systems: Transient mass transfer rates,
Water Resour. Res., Vol. 30, 321-332, Feb. (1994).
29. Prat, M., Percolation model of drying under isothermal conditions in porous media, Int. J.
Multiphase Flow, Vol. 19, 691-704 (1993).
30. Quintard, M. and Whitaker, S., Convection, dispersion, and interfacial transport of contam-
inants: Homogeneous porous media, Advances in Water Resources, Vol. 17, 221-239 (1994).
31. Rothman, D. and Zaleski, S., "Lattice gas automata", Cambridge University Press (1997).
32. Satik, C., Li, X. and Yortsos, Y.C., Scaling of single-bubble growth in a porous medium,
Phys. Rev. E, Vol. 51, 3286-3295 (1995).
33. Satik, C. and Yortsos, Y.C., A pore-network study of bubble growth in porous media driven
by heat transfer, J. Heat Transfer, Vol. 118, 455-462 (1996).
34. Stauffer, D., and Aharony, A., "Introduction to Percolation Theory", Taylor and Francis,
London (1992).
35. Tio, K.K. and Sadhal, S.S., Boundary conditions for Stokes flows near a porous membrane,
Applied Sci. Res., Vol. 52, I (1994).
36. Vafai, K. and Thiyagaraja, R., Analysis of flow and heat transfer at the interface region of a
porous medium, Int. J. Heat Mass Transfer, Vol. 30, No. 7, 1391-1405 (1987).
37. Wakao, N. and Kaguei, S., Heat and mass transfer in packed beds, Gordon and Breach
Science, New York (1982).
38. Wilkinson, D., and Willemsen, J.F., Invasion percolation: A new form of percolation theory,
J. Phys. A, Vol. 16, 3365-3376 (1983).
93
-------
Figure 4.1. Photograph showing the steady-state NAPL distribution in the glass micromodel for
typical experimental conditions (from Jia et al., 1997)
94
-------
f$8LCT-phase
«.^^*r^;
Figure 4.2. Schematic of the unit cell used by Quintard and Whitaker (1994).
95
-------
Figure 4.3. Typical NAPL interface configurations in the porespace of the micromodel. Note the
cavity-like characteristics (from Jia et al., 1997).
96
-------
10
10
Figure 4.4. The normalized Sherwood number as a function of parameter a.
97
-------
10'
10'
Figure 4.5. Numerical results of Shavg vs Pei, for mass transfer over a flat plate and for three local
exponents (c = 0, c = 1/3 and c = 1/2, from bottom to top, respectively).
98
-------
io
10'
10'
10"
10'
10"
PeL
10°
Figure 4.6. Numerical results of Shavg vs Pe£, for mass transfer over a sphere and for three local
exponents (c = 0, c = 1/3 and c = 1/2, from bottom to top, respectively).
99
-------
.""•!
Figure 4. i. A Koch curve of 3 generations used for the simulation of mass transfer over a self-similar
surface.
100
-------
35
30
g
'o
° 20
•^
S
jg
3
| 15
CD
S-
CD
Q.
5 -
S
f t
/ /1
/ / t
ill
i I I
J_
10
20 30
x—flow direction
40
50
60
Figure 4.8. Velocity profile for flow over the Koch curve of Figure 4.7.
101
-------
X10~
20
60 80
X—Flushing direction
100
120
140
Figure 4.9. The local mass flux profile over the perimeter of the Koch curve of Figure 4.7 for various
values of the Peclet number.
102
-------
0.9
0.8
0.7
| 0.6
I
8o.5
-a
I
|0.4
o
Z
0.3
0.2
0.1
Pet=8.3886e+,06
10
20
30
40 50
X-flushing direction
60
70
80
90
60 80
X—Flushing direction
140
Figure 4.10. The local concentration profile over (a) the streamwise coordinate, (b) the perimeter
of the Koch curve of Figure 4.7, for various values of the Peclet number.
103
-------
10"
10'
10"
10'
10"
ID'2 1C'1 10° 101
10"= 10a
PeL
10'
10°
10
10
Figure 4.11. Numerical results of Shavg vs Pejr, for mass transfer over the Koch curve of Figure 4.7
for three different generations and for c = 0.
104
-------
10
10
10'
Figure 4.12. Numerical results of Shavg vs Fez, for mass transfer over the Koch curve of Figure 4.7
for three different generations and for c = 1/3.
105
-------
Figure 4.13. Invasion percolation cluster used for the simulation of flow over a source with perco-
lation characteristics.
106
-------
10
10°
Figure 4.14. Numerical results of Shavg vs Pe£, for mass transfer over the percolation cluster of
Figure 4.13 and for three local exponents (c — 0, c = 1/3 and c — 1/2, from bottom to top,
respectively).
107
-------
10'
10'
Figure 4.15. Numerical results of She// vs P&L for mass transfer over the percolation cluster of
Figure 4.13 and for three local exponents (c = 0, c = 1/3 and c — 1/2, from bottom to top,
respectively).
108
-------
20
10
20
10
0 0
Figure 4.16. Uniformly-distributed sources to simulate flow over a distributed source.
109
-------
10'
10'
10'
10'
10
10"
10'
10'
PeL
10°
10'
10°
10°
Figure 4.17. Numerical results of She// vs Pe£, for mass transfer over a uniformly-distributed source
and for three local exponents (c = 0, c = 1/3 and c = 1/2, from bottom to top, respectively).
110
-------
Chapter 5
PORE-SCALE INSTABILITIES
DURING SOLVENT INJECTION
Maryam Shariati, Kathy Shing and Yarns C. Yortsos
INTRODUCTION
Solvent injection is a currently studied remediation method. It involves the injection of a mixed
solvent (frequently a mixture of alcohols) at a site contaminated with organic chemicals (Wood et
al., 1990). Typically, in situ solvent injection consists of the following steps: (1) The solvent mixture
is injected upstream of the contaminated zone; (2) the solvent with the dissolved contaminants is
extracted downstream and treated above ground to recover the solvent and (3) the recovered solvent
may be re-injected (Wood et al., 1990). Compared to pump and treat method, the extraction rate
of contaminants is enhanced because the addition of remediation fluids leads to an increase in
solubility of the nonpolar organic contaminants (Final et al., 1990, Wood et al., 1990).
Addition of alcohols alters the physical properties such as viscosity and density of the fluid phase.
Other physicochemical effects can occur when the remediation fluid comes into contact with NAPL.
For example, as NAPL dissolves in the injected fluid, some alcohol(s) from the remediation fluid
may also preferentially dissolves into the NAPL causing swelling. Such changes in the properties
of the contacting phases can have profound effects on the effectiveness of remediation treatment.
Development of solvent injection processes and the proper selection of remedial fluids to clean up
particular organic pollutants are complex tasks. Their success depends to a large extent on our
ability to understand (1) the development of miscibility at various scales, (2) the prediction of
migration of liquid, vapor, and dissolved organics and (3) the mass transfer and thermodynamic
partitioning of these compounds.
It is a well-accepted fact that addition of co-solvents such as alcohols has a significant effect on
the aqueous solubility of organic contaminants (Pinal et al., 1990). The solvent phase, in this case,
might be quite complex. It may consist of various completely miscible organic solvents, as well as
partially miscible organic solvents. The modification of solubility of organic contaminants in the
presence of these complex mixtures is frequently referred to as cosolvency. Recent studies by Wood
et al. (1990) indicate that with increasing cosolvent content the contaminant elutes at higher
concentrations, leading to an improvement in the contaminant removal efficiency. Khodadoust
et al. (1994) evaluated the removal of pentrachlorophenol (POP) from contaminated soils using
continuous solvent flushing and batch solvent washing. They used two different solvents: acetone
and ethanol. It was found that acetone was less effective than ethanol in the removal of POP.
Furthermore, the lowest solvent flow rate removed more POP, which is believed to be due to longer
111
-------
hydraulic residence time. Final et al. (1990) found that a partially miscible organic solvent (or
PMOS) present in a mixed solvent can alter the solubility of a hydrophobic organic chemical (or
HOC). In an aqueous environment, PMOSs with strong polar functional groups are most likely to
exhibit significant cosolvency effect. Nonpolar PMOSs such as TCE, octanol, toluene and other
hydrocarbons on the other hand, are not expected to have appreciable cosolvency. The cosolvency
of CMOSs(completely miscible organic solvents) increases with decreasing solvent polarity, whereas
the opposite is true for PMOSs.
While the thermodynamic partitioning of contaminants between the solid and fluid phases and
the desorption kinetics play key roles in determining the efficiency of extraction (or solubility) based
remediation processes, for in situ solvent injection processes, issues such as the mechanisms involved
in miscibility, the development of miscibility and miscroscale flow instabilities are equally important
because these impact both the efficiency of delivery of injected solvents to the contaminated sites as
well as the surface recovery of the solvents. These issues, however have not received the same level
of research interest. Because of the relative lack of understanding of these issues, in this chapter we
will examine microscale flow instabilities in mixtures with non-monotonic viscosity- concentration
profiles as they are typical of water/alcohol mixtures.
When a solvent system is injected for in situ remediation of NAPL trapped in a porous medium,
we must consider the issue of how the injected fluid penetrates the saturated porous medium, either
prior to or after contacting the NAPL contaminated regions. The issues of miscible displacement,
how miscibility develops and the possibility of microscale instabilities (fingering at the pore scale)
naturally arise. These issues are important because the extent of mixing of the injected solvents
with water in the saturated medium affects the concentration of the remediation fluid when it comes
in contact with the contaminant (often adsorbed on porous medium surfaces or collected in the
pores), thus influencing both the thermodynamic partitioning as well as the rate of contaminant
transfer into the remediation fluid. The further dilution of the remediation fluid as it flows towards
extraction points for collection and treatment prior to recycle also involves miscible displacement
and affects the total remediation cost. The development of miscibility and microscale instability
depend on a number of factors, including the fluid species present, the geometry of the medium as
well as the flow conditions. The system could be either locally fully miscible or it could consist of
immiscible and miscible regions depending on the aqueous and NAPL composition and properties
of different phases. In case of complete miscibility, a single phase exists, however in case of partial
miscibility, two distinct phases may exist. A brief review of previous work on miscible displacement
and flow instability is given in the following.
LITERATURE REVIEW
The literature in the study of miscible displacements in porous media can be classified into
two categories: Studies dealing with porous media, or Hele-Shaw cells, and studies in single-pore
geometries, such as the gap of a Hele-Shaw cell or a capillary. These are reviewed separately below.
Comprehensive reviews of viscous fingering can be found in Saffman (1986) and Homsy (1987).
One of the pioneers in the field of miscible displacement was Hill (1952). He studied viscous
fingering in miscible displacement. Saffman and Taylor (1958) determined experimentally that when
a fluid displaces a more viscous fluid in a Hele-Shaw cell, the interface between them is unstable
leading to viscous fingering. The properties of a single penetrating finger were determined. Slobod
and Thomas (1963) applied an X-ray technique to visualize displacements in a two-dimensional
porous slab. They found that the viscous fingering pattern is determined by the mobility ratio
and the velocity of the displacing fluid. Wooding (1969) studied gravity-driven fingering in a
112
-------
Hele- Shaw cell, a phenomenon analogous to viscous fingering. His experiments at low velocities
showed fingering patterns similar to Slobod and Thomas (1963). At a high inclination angle,
tip-splitting phenomena were observed, suggestive of the 'nonlinear interaction between fingers.
Pitts (1980) gave another account of the Saffman-Taylor problem. Kopf-Sill and Homsy (1988)
conducted experiments in Hele-Shaw cells to study nonlinear viscous fingering. They observed
tip splitting, shielding, and spreading. Bacri et al. (1992) studied viscous fingering flow in a
rectilinear, homogeneous, and isotropic porous medium in the presence of gravitational forces.
Their experimental work indicated the existence of a crossover between a diffusive and a convective
regime.
To interpret the behavior of miscible displacement and viscous fingering, a number of stability
analyses have been carried out. Saffman and Taylor (1958) analyzed the shape of a single steady
translating finger in immiscible displacement in the limit of zero surface tension. Chouke et al.
(1959) developed a basic stability theory and showed that if the displacing fluid is less viscous than
the displaced fluid, the unfavorable mobility profile will lead to the well-known fingering instability.
This causes the displacing fluid to channel through the displaced zone, thereby reducing the effi-
ciency of the displacement process. Hickernell and Yortsos (1986) studied displacement processes in
the absence of dispersion and showed that mobility profiles with any segment of decreasing mobility
are unstable in the linear regime. Tan and Homsy (1986) analyzed the linear stability of miscible
displacements in the presence of dispersion. They computed initial growth rates and wavelengths
of the fingers and reported good agreement between their analysis and previous experimental work
(Slobod and Thomas, 1963). Tanveer (1987) performed an analytic theory of the Saffman-Taylor
finger in the limit of zero surface tension and showed why a finger with half-spacing is the finger
selected. Chikhliwala and Yortsos (1988) presented linear and nonlinear stability analysis results
for immiscible displacement. Yortsos and Zeybek (1988) analyzed the effect of velocity-dependent
anisotropic dispersion on the stability of miscible displacements. They identified a short wave in-
stability for short times and steep concentration profiles, the properties of which depended on the
displacement velocity, the mobility ratio, and the velocity-dependent dispersion.
A number of simulation studies have been conducted to predict the behavior of miscible dis-
placement in porous media. Tryggvason and Aref (1983) simulated the interactions between fingers
in immiscible displacement in a Hele-Shaw cell. Tan and Homsy (1988) studied the nonlinear behav-
ior of viscous fingering in miscible displacements. They concluded that the spreading and shielding
effects are caused by a span wise secondary instability, and are aided by the transverse dispersion.
Zimmerman and Homsy (1991) examined the effect of anisotropic dispersion on nonlinear viscous
fingering in miscible displacements. Their simulations for a Hele-Shaw cell produced shielding,
spreading, tip splitting, and pairing of viscous fingers. Tchelepi et al. (1993) analyzed stable and
unstable displacements experimentally and performed simulations of these displacements. Roger-
son and Meiburg (1993) simulated the nonlinear evolution of the interface between two miscible
fluids of different densities and viscosities. Finally, Manickam and Homsy (1995) studied fingering
instabilities in vertical miscible displacement flows in porous media driven both by viscosity and
density contrasts.
The study of displacements at the pore-scale is restricted to flows in the gap of a Hele-Shaw cell
or in a cylindrical tube. Reinelt and Saffman (1985) studied theoretically immiscible displacements
in the gap of a Hele-Shaw cell in the presence of capillary forces. Paterson (1985) carried out
miscible experiments in Hele-Shaw cells. Joseph (1990) performed experiments in which a drop
of one fluid moved through another fluid that was miscible in all proportions. Petitjeans and
Maxworthy (1996) carried out experiments in a capillary tube using miscible fluids of different
viscosities, while Chen and Meiburg (1996) simulated these experiments. Rakotomalala et al.
(1996) simulated miscible displacement in the gap of a Hele-Shaw cell using a two-dimensional
113
-------
BGK lattice gas method. They solved the full Navier-Stokes equation coupled with the convective-
diffusive equation. Their findings are similar to those of Yang and Yortsos (1997). The latter
authors performed an asymptotic analysis of miscible displacement in the limits of large aspect
ratio under conditions of Stokes flow and identified the characteristics of displacement in the gap
of a Hele-Shaw cell or a single capillary. Finally, Lajeunesse et al. (1997) conducted experimental
work on the downward miscible displacement of a fluid by a lighter and less viscous one in the
gap of a Hele-Shaw cell. Their findings were interpreted successfully using the theory of Yang and
Yortsos (1997).
In the analysis of miscible displacements, a key difficulty is that the fluid viscosity is a func-
tion of concentration, which varies as miscibility develops. As a result, the flow field cannot be
decoupled from concentration. In the context of soil remediation, a further complication may arise
because several mixtures involved in in-situ solvent injection (1-propanol-water, 2-propanol-water,
and methanol-water, for example), have non- monotonic viscosity concentration profiles. There-
fore, as miscibility develops, locally varying viscosity results in locally varying mobility ratios which
may change from favorable to unfavorable, thus greatly complicating the issue of flow instability.
This brings to the forefront, the issue of non-monotonic viscosity profiles. So far, the effect of
nonmonotonicity of viscosity profiles on the dynamics of transport has been largely unexplored.
It is postulated that non- monotonic profiles create the potential for propagation of nonlinear fin-
gers (Manickam and Homsy, 1993). These authors provided the only published work so far on
non-monotonic viscosity profile in miscible displacements using Darcy's law. However, for miscible
displacements at the pore scale, Darcy's law does not apply. At this scale, no experimental or
theoretical work on fluids with non-monotonic viscosity profiles has been conducted at the pore-
scale. In addition, viscous fingering experiments involving non-monotonic profiles have also not
been reported.
In view of their importance on microscale instabilities, and the relative scarcity in fundamental
studies of miscible displacements, we consider in this chapter miscible displacements in systems
with non-monotonic viscosity-concentration profiles, which are of interest in alcohol-based solvent
injection remediation processes.
First, we consider the application of the methodology of Yang and Yortsos (1997) to develop
solutions of Stokes' flow in the gap of Hele-Shaw cell for non-monotonic viscosity profiles. This
configuration is a limiting case of the more general problem of miscible displacement in a pore
network, which will be explored in a subsequent study. As Yang and Yortsos (1997) have shown,
this problem is also qualitatively similar to displacement in a capillary. By taking advantage of the
large aspect ratio, we develop a simple asymptotic analysis of miscible displacement. We show that
the leading-order term in an expansion with respect to the inverse aspect ratio satisfies a single
integro-differential equation for the concentration profile in the gap or in a capillary, the passive
case limit of which yields the Taylor-Aris problem. The equation obtained is numerically solved to
give concentration profiles at various mobility ratios, a, non-monotonicity features and for different
values of the Peclet number and the viscosity ratio. The results are used to delineate the conditions
for which transverse averaging and the conventional CDE formalism for the associated Hele-Shaw
miscible displacement is valid.
Subsequently, we report on some visualization experiments in Hele-Shaw cells involving non-
monotonic viscosity profiles, for example using fluid pairs of water, 2-propanol and glycerol. The
preliminary experiments show aspects consistent with the stability analysis of Manickam and Homsy
(1993) as well as its simplification reported in the next chapter. Microscale (single-pore) instabilities
are also indicated. These can be addressed using the approach of Yang and Yortsos (1997) as
outlined below.
114
-------
MISCIBLE DISPLACEMENT IN THE GAP OF A HELE-SHAW CELL FOR
NON-MONOTONIC VISCOSITY
Consider the miscible displacement of fluid 1 by fluid 2 between two parallel plates or in a
cylindrical capillary of the geometries shown in Figure 5.1. The two fluids are Newtonian and
the injection rate q is sufficiently low for the displacement to be in the Stokes regime. We follow
the conventional approach for the description of miscible flow, in which, in the absence of volume
of mixing, the velocity is divergence-free and the difFusivity is constant (Manickam and Homsy,
1993). The alternative formalism of Joseph (1992), in which an extended divergence-free velocity is
defined, is also applicable. However, we ignore Korteweg stresses (Joseph, 1992). For simplicity, we
restrict this study to displacements in the absence of gravity. The assumption of constant difFusivity
does not change qualitatively the results to be obtained (particularly if it is recalled that we are
interested in the high velocity limit, where convection predominates).
Mathematical Formulation
In displacement between parallel plates, we are interested in the concentration profile in the
gap (X - Y) between the two plates (Figure 5.la), thus we take the displacement independent of
the third coordinate Z. At large aspect ratio, this is the geometry of a Hele-Shaw cell, in which.
the 2-D planar features of miscible displacements in the X - Z plane have been extensively studied
(e.g. see Homsy, 1987, for a review). Of interest to this paper are profiles in the gap, for which
little is currently known. Under the previous assumptions, we have
v-u = o
0= -VP + V-T
(5.1)
(5.2)
(5.3)
where C denotes the concentration of the injected fluid, T is time, U is the two-dimensional velocity
vector with components U and V, P is pressure, and the deviatoric stress tensor, T, is
(5.4)
(5.5)
The viscosity is assumed to be a function of concentration, /i(C). Using the following dimensionless
notation
dV dU
X
x = —;
C
Y
y--^'i
* = 77-; * = —;
u
u=-—;
Q
_ H_
6~ L'
V
„=-;
P* =
where subscripts 1 and 2 denote initial and injected fluid respectively, and A is a normalized mobility,
the system of equations can be rearranged as follows
115
-------
dc
d(uc) i d(wc] ^
dx ' dy
du dw
1
dx dy
( 2 (Pc
\
— o
> d /. idu\ 9 d
V A~ a- +€ 7T
ax \ da/ ay
_£^.4
da;"1
dp
^r t o ^ - l \—l ~ " l j A ~
dx
d_
dy
^ ~
dy.
L^ +
2\j i i u w
i \ -i
fc -^— I A -^~
dx
= 0
(5.7)
(5.8)
(5.9)
(5.10)
In the above, we introduced the rescaled transverse velocity w = e 1v and defined the modified
inverse Peclet number, NTD — ^2 = (e^6)"1) where Pe = a§-.
Consider next the asymptotic limit e —)- 0, NTD finite, arising at large values of the aspect ratio
and large Pe. In the context of flow in porous media, this is known as the Vertical Flow Equilibrium
approximation, (e.g. see Lake 1989, Yortsos, 1995). We proceed with a formal regular asymptotic
expansion in powers of e. At zeroth order, we obtain
dc d(uc) d(wc]
dt^~ dx "*" dy TD
d2c\
du
Hfx
dp d (,-ldu\_
—5 r ~^~ I A -5— I —
dx dy \ dy)
9P _„
(5.11)
(5.12)
(5.13)
(5.14)
where all variables denote zeroth-order in the corresponding expansion. Equation (5.14) shows
that to first-order, pressure changes only in the direction of the displacement, as in the lubrication
approximation (Ockendon, 1995). This does not imply absence of transverse convective mixing,
which, as is apparent in (5.11)-(5.12), is retained to leading order. Direct integration of (5.13) and
(5.14), and use of the no-slip condition at the walls and the incompressibility constraint /Oa udy =1,
leads to a compact result for the streamwise velocity component
u =
SoGdy
where we introduced the function
rv /•! ry yi
G(y\ A) = / Xdy I Xydy - I Xydy I Xdy
JO Jy JO Jy
(5.15)
(5.16)
the notation G(y, A) implying the functional dependence on the viscosity mixing rule. Equations
(5.15)-(5.16) show that in this approximation, the longitudinal velocity is explicitly related to
the concentration field. Then, the problem reduces to the solution of a single integro-differential
equation (5.11), where u is given by (5.15)-(5.16) and w follows from the continuity equation
w
ry
=-i
d_u_
dx
dy
(5.17)
116
-------
Equation (5.17) shows that transverse convective mixing exists in all places where u is x— depen-
dent, namely where fronts are not parallel to the streamwise direction. In the passive tracer limit,
A = 1, we recover the Poiseuille problem, with the velocity profile
u - Qy(l - y)
Numerical Results
(5.18)
The integro-diiferential equation was numerically solved using a Total Variation Diminishing
method with flux limiter, which is a dynamic weighting between lower- and higher- order difference
schemes (see Blunt and Rubin, 1992). In this approach, time derivatives are approximated by a
first-order implicit scheme and the convective term by
m
- tt
-
-l/2i-l/2
where
dx
CT+i/2 =
? + -
-------
= TT — sin_i(-
r~~ — ??™
^773. *ITTI
a =
(5.28)
(5.29)
(5-30)
The three parameters a, /Lim, cm characterize the shape of the viscosity-concentration profile. Pa-
rameter a is the ratio of the endpoint viscosities, while fj,m is the maximum value of the viscosity,
which occurs at the given concentration cm. Several sets of simulations were performed to examine
the effects of these parameters, as well as the effect of dispersion.
Results
Miscible displacements with monotonic and non-monotonic viscosity profiles, to be referred to
subsequently as MN, and NMN, respectively, at various endpoint mobility ratios were simulated.
It is known from previous work on MN miscible displacements in porous media, that end-point
viscosity ratios corresponding to a > 1 are unfavorable, while cases with a < I are favorable.
Unfavorable viscosity contrast leads to instablities in porous media, as reviewed previously.
In their analysis of miscible displacements in the gap of a Hele-Shaw cell using MN profiles,
Yang and Yortsos (1997) reported that depending on the endpoint viscosity ratio, two different
types of displacements exist. The displacement pattern in the gap of the Hele-SHaw cell viewed
from the side of the cell has the characteristics of a finger developing at the center of symmetry
when a > 3/2 (Figure 5.2a) whereas the displacement has a parabolic shape with a protruding tip
(a bullet-like shape) when a < 3/2 (Figure 5.2b) (note that in Figure 5.2, M is used in place of a).
For the case of non-monotonic viscosity concentration behavior considered here, the picture is
complicated by the fact that the viscosity increases first to a maximum value and then monoton-
ically decreases. As a result, both (locally) favorable (bullet shape) and unfavorable (finger like)
displacement patterns may be expected to exist in the same system. Manickam and Homsy (1993)
reported that the NMN displacement in the Darcy flow regime has an unstable region adjacent to
a stable one. In the unstable region, the viscosity increases with concentration in the flow direction
to its maximum, and it is adjacent to a stable region where the viscosity decreases in the flow di-
rection. They also reported that the viscous fingering instabilities become evident in the unstable
region. These findings will be compared to experiments below. However, in the case of flow in
the gap (or, equivalently, in single-pore level displacements) examined here, there are no available
results.
We first simulated the miscible displacement of a fluid by a less viscous one for Stokes flow for
both MN and NMN cases. The concentration profile for the monotonic case with a — 5 in the
absence of diffusion is shown in Figure 5.3. (It must be noted that the spreading apparent in the
concentration profile in these and other similar figures is due to numerical dispersion, even though
the physical dispersion is set to zero, NTD — 0). Corresponding results for the non-monotonic
viscosity profile are shown in Figure 5.4. The displacement patterns in both have the general
characteristics of a finger developing at the center of symmetry, which is typical of unfavorable
displacements. However, in the NMN case there is a slight protrusion at the tip, which is a
characteristic of a favorable displacement (a < 1.5, compare with Figure 5.2b). This is consistent
with the fact that for the NMN case considered here, the viscosity passes through a maximum
so that there is a region near the leading edge of the front where the viscosity is decreasing in
118
-------
the flow direction resulting in locally favorable displacement. For the MN case in Figure 5.3, the
viscosity always increases in the flow direction, so the displacement pattern has a finger like shape
typical of unfavorable displacement (as in Figure 5.2a). The MN front also propagates slightly a
longer distance than in the NMN case. The reverse would be expected if the NMN case involved
a viscosity minimum. In the NMN case, the front spreads further in the y-direction, probably
a reflection of the partial "bullet-like" shape of favorable displacement patterns associated with
regions of decreasing viscosity. The effects of dispersion are shown in Figures 5.5 and 5.6. As
expected, dispersion generally has the effect of broadening the mixing zone and to overshadow
the differences in local concentration resulting from the differences in viscosity. Consequently, the
protrusion in the NMN tip displacement vanishes as dispersion is enhanced.
Figures 5.7 and 5.8 show concentration profiles for a typical "favorable displacement", in which
a fluid displaces another of lower viscosity a — 0.2. Figure 5.7 is the MN case while Figure 5.8
shows corresponding results for the NMN case. The displacements have an overall similar pattern
but there are some differences in the tip area. The NMN case has a more compact nose which is
a characteristic of the unfavorable displacement surfacing as a result of non-monotonicity in the
viscosity concentration profiles. The effect of dispersion is shown in Figures 5.9 and 5.10. Again,
the broadening of the mixing zone as result of dispersion leads to a reduction in the differences in
the tip area between the MN and NMN cases. There is an unusual curvature in the tails of the
NMN displacement patterns, which is not present in the monotonic displacement, presumably due
to the increased complexity of the flow due to a greater viscosity gradient in the NMN case. More
pronounced effects of non-monotonicity are shown in the figures to follow.
Figure 5.11 shows the displacement pattern for the case a = 1, with a MN viscosity- concentra-
tion profile. This profile matches well the classical Poiseuille flow results, as expected in this case
of small dispersion. Figure 5.12 shows for comparison purposes the NMN case, also with a — 1 but
with a maximum (dimensionless) viscosity of 10, resulting from non-monotonicity. Some unusual
features are observed here. As a result of the favorable local mobility ratio in regions where the
viscosity is decreasing in the flow direction, the NMN displacement pattern has some of the features
of the "bullet-like" shape characteristic of favorable displacement, compared to the MN case. Fur-
thermore, the NMN patterns are nonsymmetric indicating the development of a possible instability.
Since both systems have end-point viscosity ratio of unity, this unusual behavior in the NMN case
could only be explained in terms of the details of the non-monotonic viscosity- concentration profile
where the viscosity first increases to a maximum and then decreases to one again. In the region of
increasing viscosity, flow instabilities can arise due to the locally unfavorable viscosity ratios. The
instabilities are then propagated through the mixing zone (see for example Figure 5.13 where the
displacement patterns at 5 successive time intervals are shown). The growth of the unstable region
at the tip is very clear. The effect of dispersion is shown in Figures 5.14 and 5.15, for MN and NMN
cases, respectively. As before, dispersion broadens the mixing zone, and the extent of asymmetry
in the NM case is diminished. However, there is a clear difference between the two profiles, even in
the presence of dispersion, in contrast to the previous problems (Figures 5.5 and 5.6, and 5.9 and
5.10, for example). The NMN profile exhibits an expansion of concentration contours at the tip,
and a compression at the sides of the tip, compared to the MN case.
Figures 5.16-5.17 show the effect of varying the maximum viscosity on the NMN displacement
patterns. The end-point mobility ratio is kept at unity, but parameter m, which is the maximum
viscosity relative to the pure component viscosity, was varied (from 8 to 20). Each set of concen-
tration contours corresponds to a time step. The displacement patterns are essentially the same,
although the instability is more pronounced at the largest maximum viscosity (Figure 5.17). Fur-
ther investigations should be carried out to determine if a threshold limit for this behavior exists.
Figures 5.18-5.19 depict concentration profiles for various values of cm, which is the concentration
119
-------
at which maximum viscosity occurs. The displacement patterns seem to indicate that instabilities
are more pronounced as cm increases. In their study of Darcy flows, Manickam and Homsy (1993)
pointed out that the larger the value of cm, the more efficient it is in stunning the downward
propagation of instabilities leading to viscous fingering. They also reported that for non-monotonic
viscosity concentration profiles, cm determines the length of the stable zone. These effects are not
present in our analysis, which is a single-pore study and involves Stokes flows, instead. The analysis
of these features is still underway.
PRELIMINARY EXPERIMENTS
To assess the effect of non-monotonicity on miscible displacement features, preliminary exper-
iments were also carried out in Hele-Shaw cells. We must note that these experiments were not
designed to test the previous theory, although such tests would be desirable and will be conducted
in the future. Rather, the experiments probed the predictions of Manickam and Homsy (1993) as
well as those presented in the next chapter.
The experiments were carried out in a Hele-Shaw cell. The cell was constructed of two sheets
of 0.5-inch thick glass with 11X4" dimensions. The two glass plates were separated with a 0.09-cm
thick rubber gasket. The gasket sealed the cell and also maintained a constant thickness between
the two glass plates. The plates were clamped together using ten C-clamps. Vacuum grease was
applied to the edges of the cell, to prevent air leakage. A schematic diagram of the experimental
apparatus is shown in Figure 5.20. The Hele-Shaw cell was first filled with the resident fluid. The
lines were de-aired by pumping the fluid at a slow rate. To distinguish the advancing front, water
was dyed with Bromocresol dye coloring it dark blue. A Harvard Apparatus 850 model syringe
pump was used to inject the invading fluid into the cell. The flow rate of the injected fluid was
then set. The experiments were conducted at two different flow rates: (1) 7.48 ml/sec and (2)
1.91 ml/sec. A Javelin digital camera, which was connected directly to a television and VCR set,
was used to capture the images. After the breakthrough of the injected fluid was complete, the
syringe pump was turned off. The experiments were performed using four different fluids paired
with water. Their physical properties are listed in Table 1.
Table 1: List of Miscible Fluids Used
Component
1-propanol
2-propanol
Methanol
Glycerol
Molecular Weight
60.01
60.01
32
92.10
Density (g/cm3)
0.802
0.912
0.785
1.158
Viscosity (cp)
2.223
2.428
0.585
1759.6
Although a variety of experiments were conducted, only two visualization experiments will
be reported, water displacing 2-propanol, and water displacing a water-glycerol mixture. The
viscosity-concentration profiles for the two problems are shown in Figures 5.21-5.22. The Figures
show that 2-propanol-water mixtures have nonmonotonic viscosity concentration profiles, while the
glycerol/water pair has a monotonic viscosity-concentration profile. Non-monotonicity also exists
in 1-propanol-water and methanol-water mixtures, as well (not shown here)..
The displacement of 2-propanol by water is shown in Figure 5.23. This miscible displacement
has a mobility ratio of 2.223, which is a nominally unfavorable displacement because the invading
fluid has a lower viscosity leading to the possibility of flow instability and viscous fingering. The
120
-------
displacement shows characteristics associated with non-monotonocity. Thus, the leading part of
the mixing zone, where the concentration of water is low,.appears to be rather stable. By contrast,
a fingered zone follows upstream, where the concentration of water is relatively high, and the
viscosity-concentration curve of Figure 5.21 indicates a region of increasing viscosity. The existence
of these two zones is consistent with the viscosity profiles and the stability analysis of Manickam
and Homsy (1993) (see also next chapter).
In order to compare the nonmonotonic viscosity to its monotonic counterpart, a 28-weight
percent glycerol/water mixture with a viscosity of 2.274 cP was prepared. Figure 5.24 shows the
displacement of this mixture by water at an injection flow rate of 7.48 ml/sec. The mobility ratio
for this case is 2.274. Comparison with Figure 5.23 indicates a more unstable displacement in
the monotonic case. This is particularly true at the leading edge of the displacement. This is
presumably due to the fact that in the non-monotonic case, there is a region in the mixing zone
where the displacement is favorable because viscosity is decreasing in the flow direction resulting
in a more compact front. Finally, a close-up of the two displacements is shown in Figures 5.25
and 5.26, respectively. Filaments are clearly identified in the close-ups of both displacements.
However, for the case of water-glycerol displacements, which corresponds to a monotonic profile,
these filaments also exist at the leading edge of the front, in contrast to the non-monotonic case.
In either case, the filaments indicate structures of a size of the order of the gap of the Hele-Shaw
cell, thus reflecting single-pore level phenomena. We are in the process of analysing these by using
the approach outlined at the beginning of this chapter, namely the use of Stokes equations in
constricted geometries. An overall analysis of the large-scale instabilities is also presented in the
next chapter.
CONCLUSIONS
In this chapter, it was shown that non-monotonic viscosity concentration profiles have an im-
portant effect on miscible displacements. As a general observation, the displacement profiles have
the characteristics of both favorable and unfavorable displacements. This can be explained by the
fact that the non-monotonic viscosity concentration profiles possess two regions, one in which the
viscosity increases with concentration (unfavorable local mobility ratios), and the other in which
the reverse is true. Instabilities were shown to exist both at the large-scale (Darcy scale) and at
the small (single-pore) scale, provided that a sufficiently strong non-monotonicity was imposed.
The instabilities exist even when the two pure fluids have equal viscosities. This rather interesting
behavior needs to be analyzed further at the pore-scale. A further discussion of the large-scale
instabilities is presented in the chapter to follow. In either case, the instability that ensues may
lead to a reduction of the contact efficiency, hence of the removal efficiency of the process. The
degree, however, depends on the maximum mobility contrast that develops.
REFERENCES
1. Bacri, J. C., Rakotomalala, N., Salin, D. and Woumeni, R., Miscible viscous fingering: ex-
periments versus continuum approach, Phys. Fluids A, Vol. 4, 1611-1619 (1992).
2. Blunt, M., and Rubin, B., Implicit flux limiting schemes for petroleum reservoir simulation,
J. Comp. Phys., Vol. 102, 194-206 (1992).
121
-------
3. Chen, C.Y., and Meiburg, E., Miscible displacements in capillary tubes. Part 2: numerical
simulations, J. Fluid Mech., Vol. 326, 57-90 (1996).
4. Chikhliwala, E.D., and Yortsos, Y.C., Theoretical investigation on finger growth by linear
and weakly nonlinear stability analysis, SPERE, 1268-1278 (1988).
5. Chouke, R.L., van Meurs, P., and van der Poel, C., The instability of slow, immiscible, viscous
liquid-liquid displacements in permeable media, Trans. AIME, Vol. 216, 188-194 (1959).
6. Hickernell, F.J., and Yortsos, Y.C., Linear stability of miscible displacement processes in
porous media in the absence of dispersion, Stud. Appl. Math., Vol. 74, 93-115 (1986).
7. Homsy, G.M., Viscous fingering in porous media, Ann. Rev. Fluid Mech., Vol. 19, 271-311
(1987).
8. Joseph, D.D., Fluid dynamics of two miscible liquids with diffusion and gradient stresses,
Eur. J. Mech. B/Fluids, Vol. 9, 82-88 (1990).
9. Joseph, D.D., and Renardy, Y.Y., "Fundamentals of Two Fluid Dynamics, Part II: Lubricated
Transport Drops and Miscible Liquids", Springer-Verlag (1992).
10. Khodadoust, A.P., Wagner, J.A., Mackram, T.S., and Steven, I.S., Solvent washing of POP
contaminated soils with anaerobic treatment of wash fluids, Water Environ. Research, Vol.
66, 692-697 (1994).
11. Koch, D.L., and Brady, J.F., A non-local description of advection-diffusion with application
to dispersion in porous media, Phys. Fluids, Vol. 31, 242-249 (1988).
12. Lajeunesse, E., Martin, J., Rakotomalala, N., and Salin, D., 3D instability of miscible dis-
placements in aHele-Shaw cell, Phys. Rev. Lett., Vol. 79, 5254-5257 (1997).
13. Lake, L.W., "Enhanced Oil Recovery", Prentice Hall, New York (1989).
14. Manickam, O., and Homsy, G.M., Stability of miscible displacements in porous media with
nonmonotonic viscosity profiles, Phys. Fluids A, Vol. 5, 1356-1367 (1993).
Manickam, O., and Homsy, G.M., Fingering instabilities in vertical miscible displacement
flows in porous media, J. Fluid Mech., Vol. 288, 75-102 (1995).
15.
16. Ockendon, H., and Ockendon, J.R., " Viscous Flow", Cambridge University Press (1995).
17. Patterson, L., Fingering with miscible fluids in a Hele-Shaw cell, Phys. Fluids, Vol. 28, 26-30
(1985).
18. Petitjeans, P., and Maxworthy, T., Miscible displacements in capillary tubes. Part 1: exper-
iments, J. Fluid Mech., Vol. 326, 37-56 (1996).
19. Final, R., Rao, P.S.C., Lee, L.S., Cline, P.V., and Yalkowsky, S.H., Cosolvency and sol-
ubility of hydrophobic organic chemicals in mixed solvents: evaluation using partially-and
completely-miscible organic solvents, Environ. Sci. Tech. Vol. 24, 639-647 (1990).
20. Rakotomalala, N., Salin, D., and Watzky, P., Saffman-Taylor finger in parallel Stokes flow, J.
Fluid Mech., Vol. 338, 277-297 (1997).
122
-------
21. Reinelt, D.A., and Saffman, P.G., The penetration of a finger into a viscous fluid in a channel
and tube, SIAM J. Sci. Stat. Comput., Vol. 6, 542-561 (1985).
22. Rogerson, A., and Meiburg, E., Numerical simulation of miscible displacement processes in
porous media flows under gravity, Phys. Fluids A, Vol. 5, 2644-2660 (1993).
23. Saffman, P.G., Viscous fingering in Hele-Shaw cells, J. Fluid Mech., Vol. 173, 73-94 (1986).
24. Saffman, P.G., and Taylor, G.I., The penetration of a fluid into a porous medium or Hele-Shaw
Cell containing a more viscous fluid, Proc. Roy. Soc. A., Vol. 245, 312-329 (1958).
25. Slobod, R.L., and Thomas, R.A., Effect of transverse diffusion on fingering in miscible-phase
displacement, Soc. Pet. Eng. J., Vol. 3, 9-13 (1963).
26. Tan, C.T., and Homsy, G.M., Stability of miscible displacements in porous media: rectilinear
flows, Phys. Fluids, Vol. 29, 3549-3556 (1986).
27. Tan, C.T., and Homsy, G.M., Simulation of nonlinear viscous fingering in miscible displace-
ment, Phys. Fluids, Vol. 31, 1330-1338 (1988).
28. Tanveer, S., Analytic theory for the linear stability of the Saffman-Taylor finger, Phys. Fluids
Vol. 30, 2318-2342 (1987).
29. Tchelepi, H.A., and Orr Jr., P.M., Dispersion, permeability heterogeneity, and viscous fin-
gering: acoustic experimental observation and particle-tracking simulations, Phys. Fluids A
Vol. 5, 1558-1574 (1993).
30. Tryggvason, G., and Aref, H., Numerical experiments on Hele-Shaw flow with a sharp inter-
face, J. Fluid Mech., Vol. 136, 1-30 (1983).
31. Wood, A.L, Bouchard, B.C., Brusseau, M.L., and Rao, P.S.C., Cosolvent effect on sorption
and mobility of organic contaminants in soils, Chemosphere, Vol. 21, 575-587 (1990).
32. Wooding, R.A., Growth of fingers at an unstable diffusing interface in a porous medium or
Hele-Shaw cell, J. Fluid Mech., Vol. 39, 477-495 (1969).
33. Yang, Z., and Yortsos, Y.C., Asymptotic solutions of miscible displacements in geometries of
large aspect ratio, Vol. 9, 286-298 (1997).
34. Yee, H.-C., and R. F. Warming, Implicit total variation diminishing (TVD) schemes for
steady-state calculations, J. Comp. Phys., Vol. 57, 327-338 (1985).
35. Yortsos, Y.C., A theoretical analysis of vertical flow equilibrium, Transp. Porous Media Vol
18, 107-129 (1995).
36. Yortsos, Y.C., and Zeybek, M., Dispersion driven instability in miscible displacement in
porous media, Phys Fluids, Vol. 31, 3511-3518 (1988).
37. Zimmerman, W.B., and Homsy, G.M., Nonlinear viscous fingering in porous media with
anisotropic dispersion, Phys. Fluids A, Vol. 3, 1859-1872 (1991).
123
-------
(a) Geometry for Flow Between Parallel Plates
(b) Geometry for Flow in Cylindrical Capillary
Figure 5.1: Typical displacement geometries (from Yang and Yortsos, 1997).
124
-------
1
OB
0.6
O4
O2
1
0.8
0.6
0.4
O2
M.IOO.NTD.O MCanc-ttrtonProlito
1
OB
06
02;
O2 O4 O6 OB
M.100. NTD»0.01
02 a4 O6 as
M.IOO.NTCfc.1
1
OB
06
O4
O2
1
as
ae
a4
02
O2 0.4 O6 OB
02 . NTO*O.Ot
OJJ 0.4 O6 OS
1
OS
O6
O4
O2
02 04 0.6 OB «
Mri.01.NTOa I
O4 O6 OB
MaO.OI.NTOalO
1
0.8
0.6
0.4
02i
O2 O4 0.6 OB
O2 0.4 0.6 OS
Figure 5.2-b: Concentration profiles for Miscible Displacements with a Monotonic
Viscosity-Concentration Profile Between Parallel Plates at Time Step=0.3 and M=0.01
for Various Values of the Dispersion Coefficient NTD- M is the endpoint viscosity
mobility ratio, which is defined as the ratio of the viscosities of the initial to the injected
fluid. NTD^S Pe)"1, where Pe=qH7D, q is the injection rate, H is the gap width, D is the
diffusion coefficient [Yang and Yortsos, 1997].
125
-------
Figure 5.3: Monotonic viscosit- unfavorable mobility ratio (a: = 5, NTD=0, Time=0.5).
126
-------
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Figure 5.4: Non-monotonic viscosity- unfavorable mobility ratio (a = 5, /J,m = 10, cm = 0.5,
NTD=0, Time=0.5).
127
-------
Figure 5.5: Monotonic viscosity. Effect of dispersion (a = 5, NTD=0.01, Time=0.3).
128
-------
Figure 5.6: Non-monotonic viscosity. Effect of dispersion (a = 5, /j,m = 10, cm — 0.5, NTD=0.01,
Time=0.3).
129
-------
0.1 0.2 0.3 0.4
0.5 0.6 0.7 0.8
Figure 5.7: Monotonic viscosity- favorable mobility ratio (a = 0.2, NTD=0, Time=0.5).
130
-------
Figure 5.8: Non-monotonic viscosity- favorable mobility ratio (a = 0.2, fim = 10, cm = 0.5,
NTD=0.01, Time=0.5).
131
-------
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Figure 5.9: Monotonic viscosity- favorable mobility ratio. Effect of dispersion (a = 0.2, NTD=0.01,
Time=0.5).
132
-------
Figure 5.10: Non-monotonic viscosity- favorable mobility ratio. Effect of dispersion (a = 0.2,
l*m - 10, cm = 0.5, NTD=0.01, Time=0.5).
133
-------
0.1 0.2 0.3 0.4 0.5 0.6 0.7
Figure 5.11: Monotonic viscosity- Poiseuille flow (a = 1, NTD=0, Time=0.4).
134
-------
Figure 5.12: Non-monotonic viscosity at Poiseuille flow conditions (a = 1, /j,m = 10, cm — 0.5,
NTD=0, Time=0.4).
135
-------
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Figure 5.13: Effect oft on non-monotonic viscosity for a = 1 (fj,m = 7, cm = 0.5, NTD=0).
136
-------
0.1 0.2 0.3 0.4 0.5 0.6 0.7
Figure 5.14: Monotonic viscosity- Poiseuille flow. Effect of dispersion (a = 1, NTD=0.01, Time=0.4).
137
-------
Figure 5.15: Non-monotonic viscosity at Poiseuille flow conditions. Effect of dispersion (a
Hm = 10, cm = 0.5, NTD=0.01, Time=0.4).
138
-------
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Figure 5.16: Non-monotonic viscosity- effect of/zm (a = 1, /zm = 8, cm = 0.5, NTD=0).
139
-------
0 0.1 0.2
0.3 0.4 0.5 0.6 0.7 0.8
Figure 5.17: Non-monotonic viscosity- effect of ^m (a = 1, /J,m = 20, cm = 0.5, NTD=0).
140
-------
Figure 5.18: Non-monotonic viscosity- effect of cm (a = 1, \im = 10, cm = 0.05, NTD=0).
141
-------
<> 0.1 0.2 0.3
0.4 0.5 0.6 0.7 0.8
Figure 5.19: Non-monotonic viscosity- effect of cm (a = 1, fj,m = 10, cm = 0.9, NTD=0).
142
-------
Image Processing
Digital Camera
00
IT
Heie-Shaw Cell
Syringe Pump
Figure 5.20: Schematic of the Experimental Apparatus Setup.
143
-------
Figure 21: 2-Propanol/Water Mixture
A
•T
Q C _
O.«J
0 _
o
ST
O O c: -
*•* /t.o
•*
£ 2-
8 -
o 1 5 _
w ''°
">
1 -
I
n *=; -
u.o
n -
U T
0 20 40 60 80 100
Weight per cent
Figure 5.21: Viscosity-concentration profile for a 2-propanol/water mixture.
144
-------
Figure 22: Glycero I/Water Mixture
10.0
I UU
qn -
ou
fto
ou
•— - 70
Q. 'u
*-* 60 -
vJVJ
& 50 -
w ou
R Ar\
8 4U
> ?n
*^ ou
on
£,\J
10
1 VJ
••>•>• > > • • i
» & ^ A i 1
_*-— - ^^
» — > — * — ^^
/
/
/
/
/
/
/
/
- 1 1
0 20 40 60 80 100
Weight Percent
Figure 5.22: Viscosity-concentration profile for a glycerol/water mixture.
145
-------
Figure 5.23: Typical image of water displacing 1-propanol.
146
-------
Figure 5.24: Typical image of water displacing glycerol.
147
-------
Figure 5.25: Close-up of water displacing 2-propanol.
148
-------
Figure 5.26: Close-up of water displacing glycerol mixture.
149
-------
Chapter 6
A NOTE ON THE EFFECT OF
DISPERSION ON THE STABILITY
OF NON-MONOTONIC MOBILITY
PROFILES IN POROUS MEDIA
Yanis C. Yortsos (with the collaboration of Didier Loggia and Dominique Salin)
INTRODUCTION
In the previous chapter we discussed the relevance of non-monotonic behavior at the single-
pore scale. The final chapter of this report is a short discussion of the effect of non-monotonicity
in promoting viscous instabilities at the larger scale.
In a recent paper by Manickam and Homsy (1993), to be referred to subsequently as MH, the
stability of non-monotonic mobility profiles in miscible displacements in porous media was analysed.
This problem involves the miscible displacement of one fluid (at dimensionless concentration C = 0)
by another fluid (at C = 1) at constant rate q. Non-monotonicity in mobility may result either
from a non-monotonic viscosity function /J-(C), when gravity effects are not considered, or from
a monotonic, but non-linearly dependent on density, viscosity function, when gravity is included
(Manickam and Homsy, 1995). In this chapter we focus on the first case, the extension to the gravity
problem being discussed elsewhere (Loggia et al., 1995). The base state profile, C, propagating in
the direction of displacement x, satisfies the equation
(6.1)
where f = a; - qt is the moving coordinate, b = 1^/(D\\i) is the length of the transition region, and
D\\ denotes the longitudinal dispersion coefficient. In stability problems of this kind, it is important
to differentiate between the two cases M > I and M < 1, respectively, where M = ^jj is the ratio
in the mobilities at the two end points.
In the absence of dispersion, the linear stability of general mobility profiles (which include non-
monotonic and gravity cases) was analysed by Hickernell and Yortsos (1986), to be referred to
subsequently as HY. Using normal modes, they derived the following eigenvalue problem
d
(6.2)
150
-------
where a — 2?r/A and w are the wavenumber and the rate of growth of the disturbance, respectively.
They showed that (6.2) admits a discrete spectrum of infinitely many eigenvalues u with the
following properties:
(i) When M > 1, the dominant (more unstable) eigenvalue has the long wavelength (LW) expansion
M-l\ „, ^
WTi)a + 0(a] <6-3)
the leading term of which can be recognized as the Saffman-Taylor (1958) term, and the short
wavelength (SW) limit
as a —> oo
(6.4)
All other eigenvalues are subdominant.
(ii) When M < 1, the dominant eigenvalue has the different LW expansion
•fO (a3) (6.5)
the SW limit remaining the same as in (6.4). Here, u;2 is a non-negative constant (positive in the
non-monotonic mobility case and zero, otherwise) with dimensions of length. For example,
(6.6)
•K
where we assumed ^ > 0 for f > & and ^ = 0 at £ = £0. A mode with the Saffman-Taylor
behavior (6.3) also exists, but it is subdominant. For this case, therefore, the ratio jjf=± is not the
relevant instability index.
In MH the linear stability analysis includes dispersion following a quasi-steady-state approxima-
tion, as in previous works (Gardner and Ypma, 1982, Tan and Homsy, 1986, Yortsos and Zeybek,
1988, and Bacri et al., 1992). Analytical results were derived (when A > 0, see below), for the LW
expansion of an arbitrary base state. At leading-order, the dominant eigenvalue was found to have
the same LW expansion as (6.3), except that the role of the index $=± is here played by
^ (6.7)
Thus, in the MH analysis, it is the difference in the derivative of the mobility at the end-points,
rather than the mobility itself, that dictates LW instability. MH also analyzed numerically the sta-
bility of more general, diffuse profiles, parametrized by time, and reached the following conclusions:
(i) When M > 1, the displacement is LW unstable, as long as A > 0. However, if A < 0,
a critical time, to be subsequently denoted by tC)1, was found, below which the displacement is
unconditionally stable. tCjl is of the order of £>||/g2, at which time the base state has spread over
a distance of the order of a pore (assuming that D\\ ~ aq, where a is a characteristic measure of
the microstructure). This result, predicting stability, even though the end-point mobility ratio is
unfavorable, is unexpected.
(ii) When M < 1, another critical time, to be denoted by tCj2, was found such that the displacement
is unconditionally stable for t < £c>2. For appropriate viscosity functions, this critical time can be
several orders of magnitude larger than iCjl, thus corresponding to a base-state spread over a
sufficiently large number of pores for a continuum description to be acceptable. The existence
of such a critical time, above which a displacement with M < 1 may become unstable, is a key
151
-------
contribution of MH. The authors explained the onset of this instability in terms of the spreading
of the base-state due to diffusion, hence of the development of unfavorable mobility regions, using
an extension to small wavelengths of their LW instability cirterion (6.7) based on A.
In this note, we follow a different approach to this problem. First, we provide a heuristic method
which explains the main findings of MH, particularly those of (ii) above, but in which the main
parameter is M - 1 rather than A. The validity of this method in the large time limit, e = -^- < 1
is subsequently discussed. Here, we note the relationship e - 5^, where tD is the dimensionless
time of MH.
THEORY
Consider the linear stability of
— +
dt q
where y is the transverse direction, and
"-dy*
with
V-q =
(6.8)
(6.9)
(6.10)
Assume, next, that the effect of longitudinal dispersion is only in the development of the base state
(for example, as in (6.1), and which will be used to scale the profile's dependence on £) and not in
the stability analysis, in which we will keep only the effect of transverse dispersion. This separation
of the effect of dispersion is consistent with the qualitative interpretation of the main effect of
dispersion given in MH. By taking normal modes c' = Sexp(wt + iay), p' = riexp(w£ + iay), and
proceeding as in HY, we obtain the following eigenvalue problem
d
dtt
u
(6.11)
Define, now
a' =
(6.12)
and make the further substitution \P = S. Upon integration, the resulting problem becomes
identical to (6.2), but with a' substituted in place of w. Consider now the rescaled variables n = ba
an{j s _ i£l) and C = C(z) = ^erfcz, where z = £/&, in terms of which the eigenvalue problem
reads
dz L dz \ V s dz
(6.13)
n
Given, therefore, the C(z) dependence (e.g. the erfc solution) and a fj,(C) functional relation (here,
non-monotonic), the solution of (6.13) provides the dimensionless growth rate, s = SHy(
terms of SHY (n), the properties of which were studied in detail in HY.
We recall that the dominant eigenvalue has the asymptotic behavior (6.3)-(6.5)), namely
152
-------
(i)s~
we have
n when M > 1, and s ~ s2n2 when M < 1, at small ra, where in parallel with (6.6)
-
7T2
(6.14)
(ii) s -4- smax = maxs^ at large n.
For a non-monotonic viscosity profile, SHY (n) is positive (and so are s2 and smax) and monotonically
increasing (see Figure 6.1). Back substitution in (6.12) gives the following result for the rate of
growth
( . f.\ r+j 7~) 2 i " (hrv\ (ft 1 ^
b
This relation is the main result of this note. It shows that in this approximation the growth rate
consists of the additive contributions of a viscous instability (second term in the RHS of (6.15))
and of a transverse dispersion stability (first term in the RHS of (6.15)). Longitudinal dispersion
enters indirectly through b, which scales with the square root of time and affects the magnitude
of the destabilizing contribution. From the properties of SHY, it is apparent that for a fixed
wavenumber a, the destabilizing contribution increases with time (increasing 6), until a maximum
is reached, following which the instability decays asymptotically. The following results follow from
(6.13)-(6.15):
(i) When M > 1, the problem is LW unstable (but SW stable), regardless of the particular mobility
details. This prediction differs from the findings of MH, which rely on the sign of A (see also below).
(ii) When M•< 1, the displacement is stable at sufficiently small b (or t}. As t increases, however,
b also increases, as a result of longitudinal dispersion, and a critical time tc may be reached when
the displacement becomes unstable. To show this more concretely, consider the LW expansion of
(6.15), which reads (e.g. compare with (6.5))
u --
s2qb) (6.16)
For the non-monotonic case, we have s2 > 0, thus, at a sufficiently large time, t = tc, (sufficiently
large 6), the quantity in parenthesis becomes positive and instability will set in. This result is
in qualitative agreement with the effect of the longitudinal dispersion in promoting instability,
discovered by MH. Comparison with (6.14) shows that tc increases as the degree of non-monotonicity
decreases. However, in our approach, it is the aggregate behavior of the mobility profile, for example
as reflected in the integral of (6.14), rather than the end-point mobility derivatives, which controls
the onset of this instability. One should also note that because of the infinite family of eigenvalues
of (6.13), the rate of growth contains a similar multiplicity (evidence of multiplicity was also found
numerically in MH). Because in the non-monotonic case SHY > 0, this multiplicity is restricted to
eigenvalues s > 0, namely to u + DI_C? > 0. This is also consistent with the numerical findings of
MH.
Typical dispersion curves obtained from (6.15) are shown in Figures 6.2 and 6.3 for two char-
acteristic viscosity profiles taken from the representation used in MH and for the case of isotropic
dispersion. In the Figures we plot the dimensionless rate of growth, a = -^, vs. the dimensioless
wavenumber, k — €ab, as defined in MH. The case M > 1 (Figure 6.2) shows instability for all
profiles (all times t). This is in contrast to Figure 6.7a of MH. On the other hand, the case M < 1
(Figure 6.3) clearly shows the effect of longitudinal dispersion (time) in making conditionally un-
stable an initially stable profile (although the rather small values of the rate of growth should be
noted). This is consistent with the main findings of MH.
153
-------
Although correctly predicting the main quantitative features, equation (6.15) cannot be ex-
pected to be in strict quantitative agreement with MH, certainly not in the case bq/Dn = O(l).
For example, the mismatch with the step profile predictions of MH is due to the fact that' the latter
actually correspond to the long-wave limit, but with the product D\\a/q fixed (large diffusion). The
same problem arises in the analysis of the Orr-Sommerfeld equation using step profiles (e.g. see
Drazin, 1981) and has also been noted in the stability of immiscible displacements in porous media
(see Huang and Yortsos, 1984). A closer agreement should exist when time is large compared to
D||/92, namely in the limit c < 1. It should be pointed out that in the present context, this limit
is actually the most physically relevant, as the base state has time to spread to a sufficient extent
over the rnicrostructure for the continuum equations used in the present description to be used
with confidence. (The same cannot be said when time is of the order of D\\/q2, as is the case with
£ci, where the profile varies sharply over distances of the order of the microscale. This may be the
origin of the intuitively unexpected prediction of MH when M > I and A < 0, see above).
Consider, now, the full eigenvalue problem including dispersion. In our nomenclature, we obtain
the fourth-order system
and
-E s + ^ +A ^
A2
= E
(6.17)
(6.18)
where we defined the mobility A = l//i. In the relevant limit e -* 0, a numerical analysis of (6.17)
and (6.18) showed that s —> s#y(^)> although non-uniformly in n. At small n, of main interest to
the present study, this limit is singular and a boundary layer, <5(e) (where 5(e) —)• 0, as e —>• 0) must
be inserted, namely
-sHy(n)g (
n
6(6}
(6.19)
iFrom our numerical results, the cross-over function g(u) approaches 1 as u —>• oo, although its
particular behavior at u -> 0 depends on the mobility profile details at the end points. In this limit,
therefore, the quantitatively correct expression for LJ should read
ba
(6.20)
which should replace (6.15) for accurate quantitative predictions. The main effect of this correction
is to shift the critical wavenumber at which instability first sets in to a larger value, in agreement
with the numerical results of MH. The specific features of 6(e) and of the small u behavior of g(u)
require further study.
REFERENCES
1. Bacri, J.-C., Salin, D., and Yortsos, Y.C., Analyse lineaire de la stabilite de 1' ecoulement de
fluides miscibles en milieux poreux, C.R. Acad. Sci. Paris, Vol. 314, 139-144 (1992).
2. Drazin, P.G., Discontinuous velocity profiles for the Orr-Sommerfeld equation, J. Fluid Mech.,
Vol. 10, 571-586(1961).
154
-------
3. Gardner, J.W., and Ypma, J.G.J., SPE paper 10686, presented at the 1982 SPE Enhanced
Oil Recovery Symposium, Tulsa, OK (1982).
4. Hickernell, F.J., and Yortsos, Y.C., Linear stability of miscible displacement processes in
porous media in the absence of dispersion, Stud. Appl. Math., Vol. 74, 93-115 (1986).
5. Loggia, D., Rakotomalala, N., Salin, D., and Yortsos, Y.C., Evidence of new instability
thresholds in miscible displacements in porous media, Europhys. Lett., Vol. 32, 633-638
(1995).
6. Manickam, O., and Homsy, G.M., Stability of miscible displacements in porous media with
nonmonotonic viscosity profiles, Phys. Fluids A, Vol. 5, 1356-1367 (1993).
7. Manickam, O., and Homsy, G.M., Fingering instabilities in vertical miscible displacement
flows in porous media, J. Fluid. Mech., Vol. 288, 75-102 (1995).
8. Saffman, P.G., and Taylor, G.I., The penetration of a fluid into a porous medium or Hele-
Shaw cell containing a more vicous liquid, Proc. Roy. Soc. London., Vol. A 24, 312-329
(1958).
9. Tan, C.T., and Homsy, G.M., Stability of miscible displacements in porous media: Rectilinear
flow, Phys. Fluids, Vol. 29, 3549-3556 (1986).
10. Yortsos, Y.C., and Huang, A.B., Linear stability analysis of immiscible displacement: Part
1- Simple basic flow profiles, SPE Reservoir Engineering, 378-390 (1986).
11. Yortsos, Y.C., and Zeybek, M., Dispersion-driven instability in miscible displacement in
porous media, Phys. Fluids, Vol. 31, 3511-3518 (1988).
155
-------
100
Figure 6.1. The dominant eigenvalue s#y(n) for a non-monotonic viscosity profile with M > 1
(curve a) and M < 1 (curve b) .
156
-------
0.1
0.05
a
-0.05
-0.1
.0001
t . . .
0.2
0.4
0.6
0.8
k
Figure 6.2. The dispersion relation a(k; c) for different values of the dimensionless time
predicted from equation (6.15) for the conditions of case (a) of Figure 1.
as
157
-------
2x10-*
a
-1x10-*
1000
6000
0 0.01 0.02 0.03
k
Figure 6.3. The dispersion relation cr(fc;e) for different values of the dimensionless time
predicted from equation (6.15) for the conditions of case (b) of Figure 1.
as
158
-------
Chapter 7
CONCLUSIONS
In this report we presented studies on NAPL immiscible displacement patterns using pore-network
simulation, NAPL dissolution experiments using micromodels, the theory of mass transfer in porous
media using distributed sources, the development of single-pore flow instabilities in non-monotonic
viscosity fluids, expected in solvent flushing, and the stability of displacements in porous media
using the same fluids.
The effect of viscous forces on drainage displacements in porous media was studied. We rec-
ognized that the process, at least near the front, shares common aspects with IPG. When M is
sufficiently small, the displacement can be modeled by a form of Gradient Percolation in a stabiliz-
ing gradient. We developed the scaling of the front width and the saturation profile, in terms of the
capillary number. As the stabilized regime is described by the Buckley-Leverett equation, the two
share the same constraints for their validity. In the opposite case, the displacement is described
by Gradient Percolation in a destabilizing gradient and leads to fingering. The particular regime
involves a competition between capillary and viscous forces and was identified for the first time in
the context of viscous displacements. The theory shows that the conventional continuum approach
should be used with caution near the front.
We used pore network models, both full-scale simulation and simpler statistical physics models,
for describing immiscible displacement patterns. In particular we used a simple version of IPG with
both transverse and parallel gradients, to simulate the spreading of DNAPLs. It was shown that
the spreading of DNAPL on the top of an impermeable medium may well have strong percolation
(namely fractal-like) characteristics, when the rate of penetration is slow, so that gravity and
capillarity dominate over viscous forces (and also in the case of unfavorable mobility ratio). In
such cases, the percolation characteristics must be carefully considered in the design of remediation
processes.
Subsequently, we reported on visualization experiments and numerical simulations in pore net-
works conducted to understand basic aspects of mass transfer during the solubilization of residual
NAPLs. The experiments were carried out in 2-D etched-glass micromodels with randomly dis-
tributed pore sizes. We monitored the concentration of the effluent at steady-state as a function of
the Peclet number to study the dependence of the overall mass transfer rate on injection velocity. A
pore network numerical model, based on the convection-diffusion equation using appropriate mod-
ifications for the local mass transfer coefficients, was developed to simulate mass transfer during
the solubilization of a residual phase. The pore network simulator was -found to match well the
experimental results, provided that the local mass transfer was properly accounted for. Expressions
for mass transfer over cavities were used to represent the latter. Sensitivity studies were subse-
quently conducted to investigate the dependence of the overall solubilization rate on parameters,
159
-------
such as the Peclet number. A scaling law was developed for the resulting overall mass transfer rate
assuming a uniform local mass transfer coefficient.
In the subsequent chapter, we developed a theoretical analysis and a corresponding numerical
simulation of convective mass transfer from stationary distributed sources in a porous medium. The
key assumption made in the analysis is the use of a traditional, effective porous medium description
in the porespace outside of the sources, coupled with a flux condition (rather than a constant
concentration condition) at the source interface, to reflect the effect of the local microstructure on
the mass transfer rate, as well as the fact that the effective description breaks down in such places.
We developed theoretical predictions for the dependence of the overall mass transfer coefficient on
the flow velocity and the source geometry. In general, we derived the asymptotics of the problem
at both small and large Pe/> We expect that the asymptotics at low P&L should be particularly
useful to macroscopic simulations. It was shown that the large P&L behavior reflects the local
mass transfer exponent c, thus offering an interpretation on the various experimental correlations.
Particular emphasis was placed on mass transfer from a self-similar surface, such as a Koch curve or
a percolation cluster. A scaling theory was proposed for these surfaces but not proved. Additional
work is needed in this direction. Finally, we established a connection between the macroscopic
equation typically used for analyzing data, and the microstructure of the problem, and used pore-
network simulation to compute the effective parameters. This methodology should be useful for
macroscopic simulation.
Finally, it was shown that non-monotonic viscosity concentration profiles have an important
effect on miscible displacements. As a general observation, the displacement profiles have the
characteristics of both favorable and unfavorable displacements. This can be explained by the
fact that the non-monotonic viscosity concentration profiles possess two regions, one in which the
viscosity increases with concentration (unfavorable local mobility ratios), and the other in which
the reverse is true. Instabilities were shown to exist both at the large-scale (Darcy scale) and at
the small (single-pore) scale, provided that a sufficiently strong non-monotonicity was imposed.
The instabilities exist even when the two pure fluids have equal viscosities. This rather interesting
behavior needs to be analyzed further at the pore-scale. A further discussion of the large-scale
instabilities was also presented, emphasizing the importance of diffusion in the onset, growth and
decay of the instabilities.
160
------- |