PB88-131006
EPA/600/2-87/101
November 1987
PHYSICS OP IMMISCIBLE PLOW IN POROUS MEDIA
by
J. C. Parker, R. J. Lenhard and T. Kuppusamy
Virginia Polytechnic Institute and State University
Blacksburg, Virginia 24061
Project No. CR-812073
Project Officer
Thomas Short
R. S. Kerr Environmental Research Laboratory
Ada, Oklahoma 74820
R. S. KERR ENVIRONMENTAL RESEARCH LABORATORY
OFFICE OF RESEARCH AND DEVELOPMENT
U.S. ENVIRONMENTAL PROTECTION AGENCY
ADA, OKLAHOMA 74820
-------
FOREWORD
EPA is charged by Congress to protect the Nation's land, air and water
systems. Under a mandate of national environmental laws focused on air and
water quality, solid waste management and the control of toxic substances,
pesticides, noise and radiation, the Agency strives to formulate and imple-
ment actions which lead to a compatible balance between human activities and
the ability of natural systems to support and nurture life.
The Robert S. Kerr Environmental Research Laboratory is the Agency's
center of expertise for investigation of the soil and subsurface environment.
Personnel at the Laboratory are responsible for management of research pro-
grams to: (a) determine the fate, transport and transformation rates of
pollutants in the soil, the unsaturated and the saturated zones of the
subsurface environment; (b) define the processes to be used in character-
izing the soil and subsurface environment as a receptor of pollutants; (c)
develop techniques for predicting the effect of pollutants on ground water,
soil, and Indigenous organisms; and (d) define and demonstrate the applica-
bility and limitations of using natural processes, indigenous to the soil
and subsurface environment, for the protection of this resource.
This report presents the conceptual formulation, numerical implementa-
tion, and experimental validation of a model for the movement of organic
chemicals in the subsurface environment which are introduced as nonaqueous
phase liquids through spills or leaks.
Clinton W. Hall
Director
Robert S. Kerr Environmental
Research Laboratory
i i i
-------
ABSTRACT
This project addresses the conceptual formulation, numerical implementation
and experimental validation of a model for the movement of organic chemicals
which are introduced into soils as nonaqueous phase liquids via surface spills
or leakage from subsurface containment facilities. Numerical procedures were
investigated and implemented for solution of the governing equations for flow
in three phase systems under the assumption of constant gas phase pressure,
coupled with the equations for three phase component transport with
equilibrium phase partitioning. A two-dimensional finite element code was
developed which was used to evaluate a proposed constitutive model for
properties governing multiphase flow. The code was also applied to
investigate various hypothetical problems involving leakage from underground
storage facilities.
Continuum-based mathematical models for fluid flow in porous media require
knowledge of relationships between fluid pressures (P), saturations (S) and
permeabilities (k) for the fluids and porous media of concern to enable
prediction of convective velocities and saturations. We developed a concise
parametric description for S-P relations based on a scaling procedure
formulated to separate porous medium-dependent and fluid-dependent effects.
Relative permeability-saturation relations are predicted from the scaled S-P
functional, which is assumed to reflect the pore size distribution of the
medium. The general form of the model provides a unified theoretical
framework for the description of three phase k-S-P relations subject to
arbitrary saturation paths, including effects of hysteresis and nonwetting
fluid entrapment. Model calibration is predicated on the availability of two
phase saturation-pressure measurements alone. Data requirements may be
further reduced by employing interfacial tension data to determine scaling
coefficients.
Experimental studies were carried out to validate the constitutive model for
cases involving monotonically draining water and total liquid saturation paths.
Procedures were developed for the direct measurement of fluid pressures and
saturations in three phase systems. Measurements of S-P relations for static
two and three phase systems provided model validation and demonstrated the
feasibility of using interfacial tension data to simplify model calibration. Model
validation under transient flow conditions was investigated by comparison of
numerically simulated and experimentally measured results. Two phase
investigations were conducted involving measurement of water displacement
during constant pressure NAPL injection. Three phase studies entailed
measurements of fluid saturations and pressures in soil columns during a
NAPL infiltration and redistribution event. These validation studies indicated
the k-S-P model provides a reasonable description of the system behavior
under the imposed experimental conditions and may be satisfactorily calibrated
using, relatively simple procedures.
iv
-------
CONTENTS
Abstract lv
Acknowledgements vi
1. Introduction jj
2. Model for nonhysteretic multiphase flow 4
Formulation for nonhysteretic multiphase flow 5
Characterization of S-P functions 5
Characterization of fluid permeability functions 14
Summary and conclusions 16
3. Measurement and estimation of static two phase S-P relations 18
Theoretical 19
Methods and materials 21
Results and discussion 23
Summary and conclusions 30
4. Measurement and estimation of static three phase S-P relations.. 31
Methods and materials 32
Results and discussion 36
Summary and conclusions 41
5. Finite element analysis for multiphase flow 42
Model formulation 42
Hypothetical two dimensional problems 48
Summary and conclusions 52
6. Model validation and calibration for two phase flow 53
Calibration methods 53
Results and discussion 55
Conclusions 59
7. Model validation and calibration for three phase flow 60
Methods and materials 60
Results and discussion 67
Summary and conclusions 80
8. Model for hysteretic k-S-P relations 81
Basic notation 83
Description of scaled S-P relations 86
Analysis of fluid entrapment 90
Two phase permeabili ty model development 97
Three phase permeability model development 10°
Model calibration procedures 109
Predicted Je-S-h relations for hypothetical scenario 115
Summary and conclusions 121
9. Model for multiphase component transport 123
Model development 123
Numerical approach 129
Hypothetical problem 131
Summary and conclusions 132
References • 133
Publications and presentations associated with project 139
-------
ACKNOWLEDGEMENTS
We would like to express sincere thanks to Carl Enfield, who served as project
officer during the first year of this study and who continued to provide
advise on many occasions afterwards, and to Tom Short who served as project
officer for the last two years of the study. Special thanks are also due to
Jamea McNabb whose assistance in leveling administrative obstacles was keenly
welcome on numerous occasions. Successful laboratory investigations were
made possible by the very capable and dedicated efforts of Marshall McCord,
who helped develop laboratory protocols and carry out the laboratory studies,
and Ron Alls who designed and built the specialized equipment which was
needed. Jacob Dane of Auburn University provided expertise needed to
develop procedures for simultaneous measurement of three phase liquid
saturations and kindly provided us with access to the Auburn dual gamma
attenuation facility while a unit was under construction at VPI. His input into
Chapter 7 of this report was substantial. Assistance with the development of
numerical solutions to the multiphase flow and transport problem was provided
by Jopan Sheng and Bon-Hsiang Lien who served as graduate research
assistants. Fruitful discussions over the years of this project with colleagues
from a variety of fields and from many different institutions have subtlely but
surely shaped our thoughts in many ways. Though impractical to name all
such individuals, we note their contributions in the aggregate.
Vi
-------
CHAPTER 1
INTRODUCTION
Contamination of groundwater by organic chemicals has become a serious
threat to subsurface water resources. Among the most widespread and
hazardous contaminants are organic liquids of low water solubility introduced
into the subsurface environment via surface spills, leaks from underground
storage facilities, or seepage from improperly designed or managed landfills or
land disposal operations. Numerous contamination incidents have been
documented resulting from petroleum hydrocarbons and more recently to an
increasing extent from non-petroleum sources (e.g., see review of case
histories by Abriola, 1984). Such contamination problems generally involve
complex mixtures of multiple organic constituents which may move in aqueous
and nonaqueous liquid phases and in the gas phase and which may undergo
chemically and biologically mediated degradation with time. Modeling of these
systems requires consideration of multiphase fluid flow and multicomponent
transport and reaction within each phase. The development of efficient
numerical algorithms for such problems presents many difficulties. An even
more fundamental impediment, however, has been the substantial lack of
information concerning constitutive relationships governing multiphase flow
and transport.
Under assistance agreement CR-812073 with R. S. Kerr Laboratory, we have
focused attention on the characterization of properties governing multiphase
flow with the following specific goals:
*
• develop a parametric model for hysteretic permeability-saturation-
pressure relationships in three phase porous media systems,
• develop experimental methods and specialized apparatii to directly
measure fluid pressures and saturations in three phase systems for
purposes of model validation and refinement,
• develop experimental and computational procedures for routine model
calibration from two phase system measurements,
• implement a two-dimensional finite .element model for liquid flow and
single component transport in a three phase system at constant gas
pressure and experimentally validate flow model for selected saturation
paths,
Substantial progress was made under CR-812073 and the above objectives have
been met. Results of studies pertaining to this project have been published,
or are currently in'press or in review for publication, in the open scientific
literature and have been reported at a number of meetings (see section
"Publications and presentations associated with project"). This report
presents a compendium of these various publications, edited to achieve
continuity and uniformity with a minimum of redundancy. In the remainder of
this section, we will briefly review the most important results of the study.
To model the transport of organic contaminants introduced into the subsurface
as nonaqueous phase liquids (NAPLa), it is first necessary to accurately
-------
decribe the convective movement of coexisting fluid phases, i.e., water, air and
NAPL. Continuum-based mathematical models for fluid flow in porous media
require knowledge of relationships between fluid pressures (P), saturations (S)
and permeabilities (k) for the fluids and porous media of concern to enable
prediction of convective velocities and saturations. Unfortunately, such
relationships are quite complex in three phase (air-NAPL-water) systems and
very difficult to measure directly. Therefore, much of our effort has
addressed the development of methods to describe k-S-P relations with
maximum accuracy, while keeping procedures for model calibration to a
minimum to facilitate practical application to real problems given reasonable
resources. Our approach to this problem was to develop a concise parametric
description for S-P relations based on a scaling procedure formulated so as to
separate porous medium-dependent and fluid-dependent effects. Relative
permeability-saturation relations are predicted from the scaled S-P function
which is regarded as an index of the pore size distribution of the medium
(Chapt. 2). The formulation permits calibration to be performed by measuring
only two phase system S-P relations. Furthermore, we have found that
scaling coefficients may be estimated with good accuracy from fluid interfacial
tension data, thus enabling model calibration from only air-water S-P data and
interfacial tension measurements. Validation of such calibration procedures is
demonstrated for static two phase systems in Chapter 3 and is extended to
transient and three phase systems in subsequent chapters (4, 6, 7).
A critical assumption in the k-S-P model is that two phase S-P data can be
used to predict three phase S-P relations. Although such assumptions have
been commonly invoked in the field of petroleum engineering for many years,
no direct validation has been reported due to a lack of experimental
procedures for measuring three phase S-P relations. We therefore developed
a new experimental device expressly»for the purpose of measuring three phase
S-P relations and have utilized this apparatus to test the three phase S-P
model. These tests provided direct validation of the S-P submodel for cases
of monotonically .draining water and total liquid saturation paths (Chapt. 4).
In lieu of evaluating the k-S submodel by direct measurement of three phase
k-S-P relations, which is a very tedious and arduous task, we have taken the
approach of validating the model by comparing results of transient flow
experiments with numerical simulations based on the parametric model.
Numerical procedures were investigated and implemented for the solution of
the governing equations for flow in two and three phase systems under the
assumption of constant gas phase pressure in two dimensional spatial domains
using finite element methods which have been verified numerically and
employed to investigate a number of hypothetical problems (Chapt. 5).
Experimental validation of the multiphase flow model was achieved first for
transient two phase experiments involving measurement of water displacement
by TCE and a benzene-derivative hydrocarbon (Chapt. 6).
In order to validate the k-S-P model for transient flow conditions in three
phase systems, procedures had to be developed for the measurement of fluid
saturations and pressures in three phase porous media systems. To
accomplish this, a dual energy gamma attenuation device was designed, built
and installed in a laboratory dedicated to this function. The apparatus
enables simultaneous determinations of air, water and NAPL saturations to be
made with spatial resolution of less than 1 mm on soil columns or flumes. Soil
columns were designed and constructed for use in the dual energy gamma
-------
system and were equipped with specially fabricated sensors capable of
measuring NAPL and water phase pressures concurrently. Instruments for
measuring NAPL pressures were adopted from the three phase S-P apparatus
design utilizing a special technique for producing stable hydrophobia sensors
by silanization of ceramic tensiometers. Experiments were performed involving
simultaneous measurement of liquid pressures and saturations during transient
flow in three phase air-NAPL-water systems which provided successful
validation of the three phase k-S-P model and has demonstrated the feasibility
of employing model parameters calibrated using only two phase air-water S-P
relations and interfacial tension data (Chapt. 7).
The studies mentioned to this point have all been constrained to cases in
which hystersis and nonwetting fluid entrapment is ignored or eliminated
experimentally by considering only situations involving monotonic water and
total liquid drainage. In field situations, effects of hysteresis and especially
fluid entrapment may be expected to have significant effects on flow and mass
transport. The volume that a NAPL plume will eventually occupy before
becoming effectively immobilized due to saturation path reversals will exert a
major influence on the long term behavior of aqueous and gaseous phase
transport from a contamination event, and this behavior cannot be modeled
without explicitly considering entrapment effects on the k-S-P relations.
Previous analyses of hysteresis in three phase systems have been incapable of
dealing with arbitrary saturation paths and have generally considered
hysteresis in k-S relations only.
To deal with the problem of hysteresis and fluid entrapment, we have
developed a unified theoretical framework for the description of three phase
k-S-P relations subject to arbitrary saturation paths (Chapt. 8). Thw theory
extends our previous scaled parametric model, which it includes as a special
case. Model calibration is still predicated on the availability of two phase data
alone and requires only slightly more information than the nonhysteretic
model. We have experimentally investigated calibration procedures for the
hysteretic S-P submodel which has also provided model validation for two
phase main drainage and primary imbibition paths.
Numerical studies were initiated in the final stages of the project dealing with
tho analysis of coupled three phase flow and single component transport
controlled by local phase equilibrium. The mathematical formulation for the
problem was derived and implemented numerically for a two dimensional spatial
domain (Chapt. 9).
-------
CHAPTER 2
MODEL FOR NONHYSTERBTIC MULTIPHASE FLOW
To model the movement of nonaqueous phase liquids (NAPLs) in the vadose
zone, consideration must be given to possible simultaneous convective
movement of gas, water and NAPL. A major difficulty in these modeling
efforts has been a substantial lack of information concerning the constitutive
relationships governing multiphase contaminant movement — in particular the
functional relationships between fluid pressures, saturations, and
permeabilities of coexisting phases. Because of the great difficulty of directly
measuring permeability-saturation-pressure relationships of three fluid-phase
systems, most researchers have relied upon various methods- of predicting
three-phase relationships from two-phase measurements (Leverett and Lewis,
1941; Snell, 1962; Stone, 1973; Aziz and Settari, 1979). Even with these
methods the experimental burden remains substantial. Further potential for
reductions in experimental effort are provided by procedures for estimating
relative permeabilities from comparatively easily measured saturation-pressure
data (Corey et al., 1956; Nakornthap and Evans, 1986); however, these methods
have been constrained by the use of saturation-pressure relationships having
rather narrow applicability.
Our purpose in this chapter is to present a general but concise parametric
formulation for relative permeability-saturation-pressure relationships in three
fluid phase systems which can be calibrated from a minimum of experimental
data involving measurements of two phase saturation- pressure relations. In
order to obtain a parsimonious parametric description of the constitutive
relations and to simplify generalization of experimental results to systems with
different fluid compositions, a scaling method is introduced which conceptually
permits fluid-dependent and porous medium-dependent factors to be explicitly
distiguished.
FORMULATION OF THE MULTIPHASE FLOW PROBLEM
Modeling of multiphase contaminant transport in general requires the solution
of n flow equations and p component transport equations where n is the
number of mobile phases and p is the number of mobile components. Our
focus here is on convective phase movement without regard to heterogeneous
phase compositions. Consideration will be given to convective-dispersive
transport of individual phase components in Chapter 9. Here we simply note
that the only effect .on the flow equations will be to add terms for rates of
interphase mass transfer.
For an incompressible porous medium, the water phase flow equation may be
written as
] (2.1)
-------
where * is the porous medium porosity, Sw is the water saturation, Pw is the
water pressure, k is the intrinsic permeability tensor, which we consider to be
unique to the porous medium; krw is the relative permeability to water; i\w is
the water viscosity; pw is the water density; and g is the scalar magnitude of
gravitational acceleration; and z is elevation. For organic liquid flow, we have
similarly
(2-2)
where subscript o signifies variables for NAPL. The gas phase flow equation
is written as
(2.3)
where subscript a denotes the gas phase. Owing to their dimensional
convenience, it is prefered to utilize heads rather than pressures so that
subsequently we employ, in lieu of fluid pressures, water-equivalent heads
defined by
(2.4)
(2.5)
(2.6)
To solve (2.1) - (2.3) for phase pressures and saturations, the porous medium
properties + and k and the fluid properties i*,, ifc, i±, pwt /*,, and 4, (and
their pressure dependence where significant) must be known, as well as the
functional relating JrPW» kro kro Sm So and Sa to h*, ho and hff Our
attention in this chapter will focus on formulation of the latter functional
relations with the proviso that hysteresis is ignored.
CHARACTERIZATION OP SATURATION-PRESSURE FUNCTIONS
Parameterization for Two Fluid-Phase Systems
At the microscopic scale it is apparent that pressure differences between
adjacent fluid phases will, in the absence of significant fluid-solid surface
interactions, reflect the local curvature in the fluid interfaces. For example,
consider a two fluid-phase system of water and oil. If an experiment were
-------
conducted with oil monotonically displacing water, we may expect each
incremental decrease in water saturation to cause an increase in oil pressure
relative to water as the curvatures of oil-water interfaces increase and water
is displaced. If the porous medium is rigid, it is reasonable to suppose that
the relationship between interface geometry (hence interfacial pressure drop)
and wetting fluid saturation depends uniquely on the pore geometry for a
monotonic process. Accordingly, it is assumed that in a two-phase system the
interfacial pressure difference is a function of wetting fluid saturation.
Assuming that wettability follows the order: water > oil (NAPL) > air, then this
assumption may be stated as
h0- hw= h^) (2.7)
where how will be referred to as the oil-water (or NAPL-water) capillary head
and superscripts are used to denote the fluid phases present in the system
(we follow the convention of listing subscripts and superscripts in order of
increasing wettability). In two phase air-oil and air-water systems, the
corresponding saturation-capillary head relations may be defined analogously
by
ba- h0= Aao( (2.8)
ha - bu = ha^^) (2.9)
If, as we have assumed, the porous medium is rigid and fluid-solid surface
interactions are minimal, then at a specified fluid saturation in an arbitrary
two-phase system following a monotonic saturation path, the fluid-fluid
interface geometry will depend upon the pore size distribution of the medium
and the interfacial tension of the fluid pair. In order to simplify description
of the system, we assume that a unique scaled saturation-capillary head
function S*(h«) may be obtained for a porous medium subjected to monotonic
displacement of an arbitrary wetting fluid by an arbitrary nonwetting fluid
starting from wetting fluid saturation. Physical rationale for scaling
saturation-capillary pressure functions in geometrically similar porous media
have been reviewed by Corey (1986). Here we simply pose the problem in an
empirical fashion and seek to scale two-phase saturation-capillary head
functions via the following transformations
(2-10)
(2.11)
6
-------
where the /Jjj (i,j=a,o,tf, i*j) are fluid pair-dependent scaling factors, and
are effective wetting fluid saturations in the two phase system defined by
1 - S,.
(2.13)
in which SQ, is an apparent minimum or "irreducable" wetting fluid saturation.
For a rigid porous medium with negligible fluid-solid surface interactions, we
will regard SH to be independent of fluid properties or saturation history.
Due to nonwetting phase entrapment associated with nonmonotonic displacement
histories, we may anticipate hysteresis in S*(h*). In the present analysis we
constrain our attention to monotonic wetting phase drainage from saturation.
Many empirical mathematical functions have been presented in the literature to
relate capillary pressure to wetting fluid saturation for monotonic displacement
in two phase systems (e.g. Brooks and Corey, 1964; Campbell, 1974; Su and
Brooks, 1975; van Genuchten, 1980). Here we employ an adaptation of van
Genuchten's (1980) expression to describe the scaled saturation-capillary head
relation as
For h* * 0
S* = [1 + (a/*)*]"" (2.14a)
For h* * 0
S* = 1 (2.14b)
where a and n are curve shape parameters and m = 1-l/n. Prom (2.10)-(2.14)
the general two phase saturation-capillary head relationships are given by
For hij * 0;
Sjj = S» + (1 - S,)[l+ («ijAij)n] ~" (2.15a)
For AIJ * 0:
SJj = 1 (2.15b)
for i,j = a,o,w with i*j and where for brevity we define a{jfaft{j.
-------
Experimental Analysis of Scaling Theory
To demonstrate the application of (2.15), we investigate saturation- pressure
relationships measured on two phase air-water, air-NAPL, and NAPL-water fluid
systems for a compacted well-graded sand with p-cymene as the NAPL and for
a clayey soil material (dominantly kaolinitic mineralogy) with o-xylene as the
NAPL. All experiments, performed in triplicate, involved incremental
displacement of the wetting fluid by the nonwetting phase subject to
monotonically increasing capillary pressures in a manner described in detail in
Chapter 3 which also discusses details of data analysis and interpretation.
A comparison of the scaled and unsealed saturation-capillary head data for the
two soils are presented in Figures 2.1-2.2. Instead of seeking to
nondimensionalize S*(h*), we take the air-water system as a reference and set
ftw • I (hence a • a&w) so S*(h«) is defined by (2.10). Parameters for the
scaled retention function were obtained by nonlinear least-squares fitting of
data for all three fluid pairs and replicates in each porous medium to (2.15)
(Table 2.1). Smooth curves shown on figures of both scaled and unsealed data
correspond to model predictions from these parameters. The results indicate
that the scaling procedure is remarkably successful in separating fluid-
dependent and porous medium-dependent effects on saturation-pressure
relationships, thus permitting a parsimonious and accurate description of two
phase system behavior for arbitrary fluid pairs.
Table 1. Parameters of the multiphase saturation-pressure function fit to
the main drainage branches of the sandy and clayey porous media.
Parameters1 Sand Clay
« aw
«ao
«ow
n
s»
0.052
0.099
0.110
1.84
0.00
0.032
0.108
0.060
1.86
0.36
Units of ocjj in ca"1; n and S^ diaensionless.
8
-------
ICO
Q
N
I
E
U
U
I
75 -
25 -
Q.
<
U
AIR-WATER
— • AIR-OIL
.* OIL-WATER
0. O 0. 2 O. 4 0. 6 O. 8
SATURATION
1.0
16O
E
U
120 1
UJ
I
< eo -
Q.
<
U
a
u
u
U)
4O
A AIR-WATER DATA
• OIL-WATER DATA
• AIR-OIL DATA
0.0 0.2 0.4 0.6 0.8 1.0
EFFECTIVE SATURATION
Figure 2.1. (a) Unsealed capillary head versus wetting fluid saturation for
two phase air-water, air-NAPL and NAPL-water systems in a sandy porous
medium with p-cymene as NAPL. Bach data point is a mean of triplicate
experiments. (b) Corresponding scaled capillary head versus wetting fluid
saturation data. Solid lines are model predictions from fitted function.
-------
300
a
IM
i
E
U 200
u
I
K 100 •
_l
t-i
Q.
U
o
0.0
AAIR-WATER
—— • AIR-OIL
— • OIL-WATER
•I I
0. 2 0. 4 0. 6 0. 8
SATURATION
l.O
40O
r\
E
U
300 •
UJ
X
200 •
a.
<
u
u
(/I
100 -
A AIR-WATER DATA
• AIR-OIL DATA
• OIL-WATER DATA
O. 0 0. 2 0. 4 0. 6 0. 8 1.0
EFFECTIVE SATURATION
Figure 2.2. (a) Unsealed capillary head versus wetting fluid saturation for
two phase air-water, air-NAPL and NAPL-water systems in a clayey porous
medium with o-xylene as NAPL. Bach data point is a mean of triplicate
experiments. (b) Corresponding scaled capillary head versus wetting fluid
saturation data. Solid lines are model predictions from fitted function.
10
-------
Extension to Three Fluid-Phase Systems
The analysis of three phase fluid saturation-pressure relationships is markedly
more difficult than that of two phase systems due to the much greater
potential complexity in the manner in which three fluids may interact in a
porous medium. These problems are compounded by the difficulty of directly
measuring fluid pressures and saturations simultaneously in
three phase systems. In practice, three phase behavior is commonly estimated
from two phase measurements. The extension of two phase saturation-capillary
pressure relationships to three phases was first suggested by Leverett (1941).
Regarding interfacial curvatures for a monotonic displacement process to be
fixed by the pore-size distribution of the medium, Leverett assumed total
liquid saturation, Sf= S^SQ, in a multiphase system to be a function of the
interfacial curvature at the gas-liquid interface, independent of the number or
proportions of liquids contained in the porous medium.
In three phase systems of air, water and NAPL in which water is the dominant
wetting fluid, the total liquid saturation would thus be a function of air-NAPL
capillary head. It is also commonly assumed on similar rationale that water
saturation in a three-phase system depends on the NAPL-water capillary head
when water is the wetting phase and the NAPL has wettability intermediate
between water and air (Leverett and Lewis, 1941; Corey et al., 1956; Aziz and
Settari, 1979; Eckberg and Sunada, 1984). The foregoing assumptions may be
written as
(2.16)
(2.17)
where
and superscript aow denotes the three-phase system. Combining (2.16)-(2.17)
with the two phase parametric formulation (2.15) yields for water saturation
in a three phase system
For hgif *• 0:
4*°"= [l+(«
-------
For h * 0:
ao
(2.19b)
where effective saturations in the three phase systeo are defined by
1 - S,, (2.20)
(2.21)
Implicit in (2.16)-(2.19) is the assumption that no air-water interfaces occur in
the three phase system until NAPL saturation diminishes to a level at which
the organic liquid exists only as discontinuous blobs or pendular rings at
particle contacts. Thus, air-water capillary pressures are presumed to have
no physical relevance as long as a continuous organic liquid phase exists.
Equations (2.16)-(2.17) indicate that for all conditions in which NAPL occurs,
the Sw and St functions given by (2.18) and (2.19) are independent of
whether a two or three phase system is described. At locations in the porous
medium where no NAPL occurs, i.e. in an air-water system, (2.18) is replaced
by
For hau ± 0:
s£*= [l+(aawAajy)n]~" (2.22a)
For hau * 0:
1 (2.22b)
Since Sepl-S^Sf, (2.18)-(2.22) fully define the saturation-pressure relations
in any two phase or three phase system during monotonic drainage.
Analysis of Saturation Derivatives
In order to solve the How equations given by (2.1)-(2.3), it is necessary to
eliminate the saturation time derivatives. This may be done for the NAPL flow
equation, noting that So=So(/i,wJioJia) by applying the expansion
12
-------
'A0 »hw »ha
= c°° IT * °"~n~ * °"1T (2'23)
or, in general, for the saturation time derivative of fluid i (i=a, o,w)
•Si *h\
= £ i Ci i — i__ _ .„ /p OA\
_. ^JX^J-i ^^ Cf, C/, W \ f*. *»Y /
where the fluid pair capacities are defined by
'Si
Cij* *— r (2.25)
«/
Owing to the form of the saturation-pressure relationships given by
(2.18)-(2.22), expressions for the fluid capacities Cy may be obtained in closed
form as follows:
System with organic liquid present
(2.26a)
(2.26b)
Coo = GW+ Caa (2.26c)
GIW = Co» = -Ciw (2.26d)
Cao = Cog = -Caa (2.26e)
C»o= Cw= 0 (2.26f)
phase air-water system
" (2.27a)
C«a= -Cw= C,^ (2.27b)
Cw = -CW (2-27c)
other terms - 0
For notational convenience (2.26) and (2.27) are written in terms of effective
fluid saturations rather than pressures. Using (2.18)-(2.20) the saturation
derivatives may be expressed directly in terms of fluid pressures, however,
13
-------
computationally it will be more efficient to first compute effective fluid
saturations and then evaluate saturation derivatives as well as fluid
permeabilities discussed in the next section which are also more conveniently
expressed in terms of fluid saturations.
It may be noted that for the analysis of compressible fluid flow as for the gas
phase, a term for the fluid density time derivative arises also on the left hand
aide which can be evaluated very simply given an appropriate constitutive
relationship for fluid density as a function of fluid pressure (e.g. ideal gas
law).
CHARACTERIZATION OP FLUID PERMEABILITY FUNCTIONS
Because of the considerable difficulty in measuring relative permeabilities in
three phase systems, much attention has focused on the development of
procedures to estimate three phase relative permeabilities from two phase
system measurements. In early work on three phase system behavior,
Leverett and Lewis (1941) observed that permeability to water in water-wet
gas-oil-water systems was essentially dependent only on water saturation,
while permeability to gas was a function only of gas saturation. Accordingly,
these functions may be readily determined from appropriate two phase system
measurements. Relative permeabilities to oil in three-phase systems were
observed to depend on both gas and water saturation. While this fact
precludes the possibility of directly measuring relative permeabilitity functions
for fluids of intermediate-wettability from a single two phase system test, a
number of methods have been devised to enable their estimation from suitable
combinations of measurements from two phase systems with different fluid
pairs (e.g., air-oil and oil-water).
Two methods for predicting three-phase oil relative permeabilities from two
phase measurements, which are widely uaed in petroleum reservoir modeling,
were proposed by Stone (1970, 1973) (see also Aziz and Settari, 1979; Payers
and Mathews, 1984). The methods involve the use of various empirical
formulae to predict three phase oil permeabilities from measured water and oil
permeabilities in two phase oil-water systems and air permeabilitites in two
phase air-oil systems. An alternative method for evaluating three phase
relative permeabilities from two-phase data proposed by Corey et al. (1956) is
based on the estimation of mean flow channel' diameters from
saturation-capillary pressure relationships using a modification of Burdine's
(1953) two phase relative permeability model. Comparisons of theoretically
predicted and measured oil permeabilities of an oil-brine-gas fluid system in
sandstone cores yielded good agreement between experimental and predicted
permeabilities. However, as Corey et al. point out, the applicability of their
method to a broad range of porous media would be limited by the lack of
generality of the assumed saturation-capillary pressure function.
Saturation-capillary pressure functions with greater generality than that used
by Corey et al. (1956) have been employed with flow channel distribution
models by Brooks and Corey (1964), Nakornthap and Evans (1986) and van
Genuchten (1980) to predict relative permeabilities of two phase systems with
14
-------
considerable success (e.g. van Genuchten and Nielsen, 1985). Here, we propose
extending the method of van Genuchten to the prediction of three phase
system relative permeabilitites.
Assuming water relative permeability to be a function only of water saturation
regardless of the occurrence of other fluid phases, the analysis of van
Genuchten (1980) for two phase air-water systems may be employed directly.
Using Mualem's (1976) model to evaluate effective hydraulic radii from the
saturation-capillary pressure function yields the expression
*rv = S*1/a f r(0,^) / (F(0,l) ]J (2.28)
where the term in brackets represents the ratio of the mean hydraulic radius
of unsaturated soil at effective saturation S^ to the mean hydraulic radius of
saturated soil and we now omit superscripts on the effective water saturation
since (2.28) is assumed to hold for any two or three phase water-wet system.
Taking l/h* versus S* as an index of the pore size distribution allows T(a,b)
to be written as
Integration of (2.29) with the expression for h»(S*j obtained by inverting
(2.14) gives
F(a,b) = « {[l-a/«]« - [l-bV«]»} (2.30)
which yields, upon substitution into (2.28), the desired expression for water
relative permeability written in terms of effective water saturation
(2.31)
As noted previously, this relation is presumed to be valid both in three phase
systems as well as two phase systems (e.g., air-water and NAPL- water).
To evaluate air relative permeabilities, we assume dependence only on air
saturation. Extension of Mualem's model to the gas-phase case and taking into
consideration the so-called gas slippage or Klinkenberg effect (Corey, 1986),
which may be important in fine-grained porous media, yields
Sal/' ( T(Sttl) I r(0,l) ]* (2.32)
where C = 1 + */ha is a gas slippage correction with * a porous
medium-dependent parameter which is expected to diminish as grain-size
increases. Upon substitution of (2.30) for F, (2.32) yields
15
-------
*ra = C 5a* (1 - Stl/m)*n . (2.33)
Thus, ignoring gas slippage effects, kra is predicted to be a function only of
Sa since Sfl-S9 Again, we omit any superscripts to denote phase
compositions since this is assumed to have no effect on air permeability (note
that St=Sw when So=0, i.e. in an air-water two phase system).
For the fluid phase of intermediate wettability to air and water in a three
phase system (i.e. NAPL), the appropriate formulation of Mualem's model is
kro ~ (St~su)1 f T(Si»St) / r(0,l) ]a (2.34)
Froa (2.30) for T we obtain
1/1 (['-*•"•]•-
*ro - (^r^) l - *, - l - (2.35)
indicating that in general kro is a function of both water and oil phase
saturation. In addition to representing the three phase case, (2.35)
degenerates to the appropriate relations for any two phase case except for
gas slippage correction terms. In particular, it 'may be noted that if Sa=0 and
NAPL ia the nonwetting phase in a water-wet system, then kro=kro(S0) has the
same form as kra(Sa) given by (2.33) with 0=1; if 8^=0 and organic liquid ia
the wetting phase with air then k^S^) has the same form as kfW(Sw) given
by (2.31).
SUMMARY AND CONCLUSIONS
In this chapter, we have presented a parametric model for relative
permeability-saturation-pressure relationships in porous media which contain
up to three coexisting fluid phases following monotonic saturation paths. The
model provides simple closed-form expressions for fluid saturation-pressure
functions and their derivatives and for relative permeability-saturation
functions.
By scaling saturation-capillary head relations, two fluid-phase system behavior
for an arbitrary fluid pair ia described by a scaled function involving three
Porous medium-specific parameters (a, n, Sm) and a fluid pair-dependent
scaling factor (ft). Experimental results for two porous media with different
two-phase fluid combinations indicate the scaling approach provides a good
representation of two phase saturation-pressure data. Extension of the scaled
two phase relations to three phase systems based on the assumption that
wettability follows the sequence water > NAPL > air provides a model for three
Phase saturation-pressure relations involving only five parameters (
-------
By interpreting the scaled saturation-capillary head function as an index of
the pore size distribution of the medium, expressions for the relative
permeabilites of air, water, and NAPL in a three phase system are derived
which also degenerate to appropriate functions for any subset two phase case.
No additional parameters are introduced over those involved in describing the
saturation-capillary head with the exception of a gas slippage correction factor
for gas permeabilites. In fact, several parameters in the saturation-capillary
pressure function do not arise in the permeability-saturation functions (i.e.,
the ojj ). Direct verification of the relative permeability model will require
comparison of model-predicted permeabilities with experimentally observed
values on two phase and three phase systems. Alternatively, model validation
may be sought by comparison of observed transient multiphase flow behavior
with numerical simulations based on the proposed model as will be presented
in Chapters 6 and 7.
The proposed parametric model for constitutive relationships governing
multiphase convection is expected to provide an adequate representation of
fluid-porous media systems subject to their meeting certain criteria implied in
deriving the modeL One difficulty to contend with is the assumption of a
rigid porous medium with no significant fluid-solid phase interactions. The
data shown for a clayey soil, albeit of relatively low swelling potential,
indicates that the proposed scaling procedure may not deteriorate seriously in
fine-grained materials. However, even rather small changes in the pore size
distribution, which in certain instances may not be clearly evident in the
saturation-pressure relations (e.g., development of small fissures due to clay
shrinkage in the presence of low dielectic organic fluids), may markedly affect
the permeability functions. In some instances such problems may be a
substantial concern and research is warranted to investigate methods of
dealing with their effects quantitatively.
Another proviso on the model which should be emphasized is that we have
assumed saturation paths correspond to monotonically decreasing wetting
phase saturations in all instances, and for three phase flow that total liquid
saturation also monotonically decreases. When these conditions are not met,
hysteresis in the constitutive properties is expected to be of importance,
particularly when imbibing and draining saturation-capillary head curves do
not close at zero capillary head. For example, following an initial injection of
oil into a water-saturated core, reflooding with water will commonly fail to
achieve full water saturation due to entrapment of some oil. Extension of the
model to consider such effects will be considered in Chapter 8.
17
-------
CHAPTER 3
MEASUREMENT AND ESTIMATION OP STATIC
TWO PHASE SATURATION-PRESSURE RELATIONS
Very little experimental research has been conducted to elucidate fluid-porous
medium properties governing multiphase contaminant flow. Experimental
determinations of permeability-saturation-pressure relations of two phase
porous media systems have been well documented and are relatively
•™«htforward to perform (Brooks and Corey, 1964; Corey, 1986; Scheidegger,
iy74). Measurement of three phase permeability-saturation-pressure relations
is much more difficult and is hampered by a lack of reported procedures for
the simultaneous measurement of organic liquid and water phase pressures in
Porous media. In practice, investigators normally utilize functional
relationships measured in two phase systems {i.e., air-organic and organic-
water) to estimate fluid behavior in three phase systems.
The extension of two phase saturation-capillary pressure relations to three
Phases was first suggested by Leverett (1941) who postulated that in
gas-oil-brine systems, total liquid saturation will be a function of interfacial
curvature at gas-liquid interfaces which, for a monotonic displacement process
m the absence of specific fluid-solid phase interactions, is fixed by the pore
size distribution of the medium. It is also commonly assumed that water
saturation is a function of the interfacial curvature at the organic-water
interface.
When modeling air-organic-water flow phenomena, investigators typically use
separate parametric expressions to relate water and total liquid saturations to
differences in pressures between the appropriate contiguous fluid phases. In
Chapter 2, we described a technique to scale two phase saturation-pressure
relations so that a single parametric expression can be used. When applied to
the van Genuchten (1980) function, a scaled equation was obtained which
describes saturation-pressure relations of two phase systems or, by utilizing
assumptions of extending two phase measurements to three fluid phases, can
oo extended to describe three phase systems. The scaling method is general
and can be applied to numerous retention functions that exist in the
literature. In this chapter we will evaluate the proposed scaling technique by
applying it to systems with various organic liquids. The scaling procedure
will be applied to the Brooks-Corey (1964) and van Genuchten (1980) retention
functions and a theoretical procedure for the prediction of fluid-dependent
scaling factors from ratios of fluid interfacial tensions will described and
tested using directly measured and handbook tabulated interfacial tension
data. Details of experimental methods and precautions which ahould be
followed to obtain physically consistent parameters to describe
saturation-pressure relationships for three phase systems in rigid porous
media will be discussed.
18
-------
THEORETICAL
We wish to consider the measurement of fluid pressures and saturations in two
phase porous media systems for monotonic displacement of a wetting fluid j by
nonwettmg fluid i such that a unique saturation-capillary pressure function,
Sj(fy), is obtained where Sj is the wetting fluid saturation and fli=fl-ft is
referred to as the capillary pressure with ^ and Pj denoting the respective
phase pressures. In lieu of pressures, it is often more convenient to employ
heads and we accordingly define capillary heads on an equivalent water height
basis by hjj^j/p^g, where p^ is the density of water and g is the scalar
magnitude of gravitational acceleration. In terminology more common in
hydrology, the negative of the air-water capillary head would be denoted as
the pressure head in a two phase air-water system.
Let us consider a hypothetical fluid-porous media system in which interfacial
pressure differences are due dominantly to capillary effects. Then at the
microscopic scale, the capillary pressure for two contiguous fluid phases may
be related to the curvature of the fluid-fluid interface by Laplace's equation
of capillarity
where 2a)
= S(R) (
S0(2rao/J>ao) = S(R) (3.2c)
k
where subscripts a, o and w refer to air, organic liquid and water phases,
respectively. Equation (3.2) indicates that two phase saturation-capillary
pressure functions may be scaled by interfacial tensions as suggested by
Leverett et al. (1942), Miller and Miller (1956), Dumore and Schols (1974) and
others.
In the scaling procedure proposed in Chapter 2, saturation-capillary head
relations of two phase air-water, air-organic and organic-water systems are
adjusted so as to obtain a unique scaled function for a given porous medium
19
-------
after applying a linear transformation to capillary heads such that
(2.10)
(2.11)
(2.12)
where S*(h*) is a scaled saturation-capillary head function, and ftij are fluid
pair-dependent scaling factors. At a given wetting phase saturation, (2.10) -
(2.12) indicate that
= ftao^ao (3.3)
and from (3.2), upon converting from pressure to head units, we observe that
(3.4)
Since S*(h*) is arbitrary in absolute magnitude, the scaling coefficient of a
reference system may be arbitrarily specified. We choose the native air-water
system (i.e., uncontaminated by anthropomorphic organic components) as a
convenient reference and let ftaw • 1. Employing this convention and utilizing
(3:3) and (3.4) enables the remaining scaling factors to be estimated as
(3.5a)
ftao = *au/*ao- (3.5b)
Since the presence of relatively small amounts of organic chemicals dissolved
in the aqueous phase or volatilized into the gaseous phase may substantially
alter air-water interfacial tensions in contaminated systems, it is useful to
define an additional scaling coefficient pertinent to two phase air-water
systems in which organic concentrations in aqueous and vapor phases
correspond to unit activity (i.e., phase saturation with contaminant) as
0aw* = ^aw/'aw' (3.5c)
where
-------
and the scaling factor definitions given by (3.5), it is seen that
^ao + *au = »aw* (3.7)
which via (3.5c) indicates that
(3.8)
thus enabling estimation of 0aw' from air-oil and oil-water interfacial tension
data.
MATERIALS AND EXPERIMENTAL METHODS
Saturation-capillary pressure relations were measured, in triplicate, for
air-water, air-organic and organic-water two phase systems in a sandy porous
medium containing 0.92, 0.05 and 0.03 kg kg~l sand, silt and clay,
respectively. Four different organic liquids were used. They were benzene
o-xylene [CfeHifCIfcte], p-cymene [ClfcCfeClKCHate], and benzyl alcohol
Both p-cymene and o-xylene are essentially insoluble in water
while benzene and benzyl alcohol have water solubilities of 0.7 and 40 g kg"1 ,
respectively {Dean, 1979).
Retention cells were constructed from stainless steel, teflon tubing and brass
fittings to minimize retention cell-organic liquid interactions. The porous
medium of 112.6 raL volume (51.8 mm length x 52.6 mm diameter) was bounded
at the upper surface by a coarse sintered brass disc and at the lower surface
by a fine porous ceramic disc, the cells were designed so that there would
be no leakage or fluid for pressures exceeding 1 MPa. On the top of the cell
an inlet line supplied nonwetting fluid under regulated positive pressures and
an outflow buret connected to the bottom maintained the wetting fluid
pressure near atmospheric pressure and permitted volumetric determination of
the wetting fluid outflow.
An oven-dry equivalent mass of air-dry sandy porous medium was packed into
the retention cells to obtain a prescribed bulk density of 1.55 Mg m~3 using
procedures to ensure uniform packing of the cells and hence similar pore size
distributions for all replicates. A vacuum was applied to the packed cell
through the top inlet to saturate the porous medium with wetting fluid which
entered the cell through the lower surface. The porous ceramic disc was
vacuum-saturated with wetting fluid be/ore placing soil into the cells.
A pressure difference was imposed between the fluid phases comprising a
particular two phase system by applying a positive pressure to the nonwetting
fluid and maintaining the wetting phase pressure at slightly above atmospheric
pressure at the lower soil surface. Pressure differentials were applied in a
manner to produce monotonic wetting fluid drainage paths to measure the main
drainage branches of the saturation-capillary head relations. Fluid saturations
were allowed to equilibrate with the applied pressure differentials as indicated
21
-------
by outflow rate measurements before incrementing the pressure again. Fluid
saturations were determined by gravimetric measurements of the cells, after
equilibrium was obtained, for the air-water and air-organic systems, and by
recording the water outflow volume for the organic-water systems after
correcting for evaporation from 'blank' buret measurements. Six pressure
increments were applied to each replicate to yield 54 saturation-capillary head
data pairs for each triad of two phase air-water, air-organic and
organic-water measurements.
The scaling procedure outlined in Chapter 2 was applied to the Brooks-Corey
(1964) and van Genuchten (1980) retention functions to obtain multiphase
retention functions. Model parameters were determined by nonlinear
least-squares regression analyses of the measured air-water, air-organic and
organic-water saturation-capillary head data for each of the four organic
liquids subject to various constraints which will be outlined later. The scaled
van Genuchten (1980) function as described in Chapter 2 is
f -,-M-l/n
5* = [ 1 + (« /*)n] A* > 0 (2.14a)
5* = 1 » * 0 (2.14b)
where a and n are empirical parameters. Applying the same scaling procedure
to the Brooks-Corey (1964) function yields
A* > h
-------
RESULTS AND DISCUSSION
Parameters in the Brooks-Corey and van Genuchten models were determined by
fitting the respective functions to the experimental saturation-capillary head
data by nonlinear regression analyses via three methods: (i) fit a0jj, n and
Sjn (van Genuchten model) or hd/ftj, \ and Sn, (Brooks-Corey model)
separately to each two phase system (18 data points each) with no constraints
on parameter values, (ii) fit model parameters to the pooled set of air-water,
air-organic and organic-water data for each organic fluid (54 data points)
subject to the constraint that n or X and Sm are constant for porous medium
and that paw-l» and (Hi) same procedure as Method ii but with S,n=0.
Parameter estimates, 95% confidence intervals of parameter estimates and
regression R2 values are given for each method in Tables 3.1 and 3.2 for the
van Genuchten and Brooks-Corey models, respectively. Confidence intervals
for aftij and h^/ftij values determined by the various methods did not differ
significantly at the 0.05 probability level. Confidence intervals for n and X
values in general also overlap and appear to be independent of the specific
fluid pair as anticipated from the scaling theory. Exceptions to the
overlapping of confidence intervals for n occur in the benzene-water and
benzyl alcohol-water systems (Table 3.1) and in the benzyl-alcohol systems for
X (Table 3.2). In each of these cases, the estimated values of §„, and n or X
for Method i are observed to be much larger than those for Method ii.
Inspection of the parameter covariances indicated high correlations (r= 0.69 -
0.98) between S,n and n or X which will result in a marked widening of the
hyperellipsoids defining the parameter confidence regions (Beck and Arnold,
1977). Thus we conclude that Sm has poor identiflability when n or X is also
unknown, as evidenced by the rather wide confidence intervals for these
parameters when all are estimated.
*
As a result of the poor identif lability of Sm, constraining Sm to be a constant
for the porous medium had little effect on the accuracy of the model (see ff
values, Tables 3.1 and 3.2). Imposing such a constraint also precludes
possible inconsistencies which may arise in the description of three phase
systems from two phase data when Sm for the air-organic system exceeds that
for the organic-water system. Several examples of this undesirable trend in
irreducable saturations fitted from individual two phase data (Method i) are
evident in Tables 3.1 and 3.2. We note that it is feasible to fit variable Sm
values to the pooled data subject to the constraint that values for the
air-organic system are not less than those for the organic-water system while
a and n or hd and X are held constant. However, if this is done one finds an
average, improvement in the regression Ra of less than 0.003 at the cost of
increasing the number of coefficients by two.
It is observed that estimation uncertainties for constant SQ, vlaues (Method ii)
are such that their, 95% confidence intervals include zero. If SJQ = 0 is
assumed (Method iii), the decreases in R' over the constant but nonzero Sm
cases are less than 0.007 for the van Genuchten model and 0.001 for the
Brooks-Corey model (Tables 3.1 and 3.2). Thus, for sandy media such as
employed here, taking SH, - 0 may be a justifiable simplification in the
parametric model.
23
-------
Table 3.1. Parameters of the van Genuchten retention function for four organic liquid systens with /?aw«l.
SYSTEM "PBW
air-water (Method j) .058 * .019*
benzene
two-phase (Method j)
air-organic -
organic-water — -
pooled Multiphase
S^O (Method H) .053 * .008
S.=0 (Method Hi) .054 * .008
benzyl alcohol
two-phase (Method i)
air— organic -
organic-water -
•ultiphase (pooled)
S.*0 (Method H) .050 * .008
&V=0 (Method Hi) .050 * .010
two-phase (Method jf)
air-organic -
organic— water -
pooled Multiphase
5^*0 (Method H) .051 * .006
Sm-0 (Method Hi) .051 * .008
o-xyleoe
two-phase (Method j)
air-organic -
organic-water -
pooled Multiphase
SB'0 (Method Jj) .045 * .006
SB=O (Method Hi) .043 ' .006
«*ac
.112 *
-
.115 «
.117 *
.061 *
-
.064 *
.067 *
.088 *
-
.100 *
.103 *
.084 *
-
.093 *
.094 '
i
.019
.026
.026
.009
.016
.020
.011
.018
.020
.009
.018
.020
*&»
-
.079 *
.097 *
.099 *
-
.229 *
.279 *
.298 *
-
.109 *
.113 *
.115 *
-
.076 *
.084 *
.085 *
f
.004
.022
.022
.026
.072
.084
.009
.021
.022
.016
.016
.018
n
1.746 *
1.835 *
3.476 *
1.962 *
1.802 *
2.332 *
4.267 *
2.317 *
1.858 *
2.335 *
2.339 *
2.025 *
1.837 *
2.602 *
3.207 *
2.358 *
2.003 *
.343
.045
.051
.338
.101
.701
1.63
.527
.145
.575
.313
.308
.092
.485
.448
.372
.127
s.
10-W , .
10-1° * .
.244 * .
.075 * .
0.0
.120 * .
.261 * .
.151 * .
0.0
.107 * .
.163 * .
.082 * .
0.0
.104 * .
.186 * .
.177 * .
0.0
161
251
034
123
178
073
107
133
062
103
099
041
080
R2
0.933
0.968
0.992
0.952
0.951
0.962
0.897
0.916
0.909
0.962
0.990
0.955
0.953
0.985
0.994
0.959
0.956
1 units of <*
-------
Table 3.2. Parameters of the Brooks-Corey retention function for four organic liquid systems with law'1*
N)
SYSTEM l»d//W
air-water (Method j) 9.224 * 2.14'
benzene
two-phase (Method i)
air— organic —
organic-water -
pooled Multiphase
SM*0 (Method jj) 9.296 * 1.63
SB=O (Method jjj) 9.296 * 1.61
benzyl alcohol
two-phase (Method i)
air-organic -
organic-water -
pooled Multiphase
S*»0 (Method ii) 9.550 * 1.52
S.=0 (Method Hi) 9.545 * 1.51
p-cy»ene
two-phase (Method j)
air-organic -
organic-water -
pooled Multiphase
Sm*0 (Method jj) 9.534 * 1.28
SB=0 (Method Hi) 9.531 * 1.28
o-xylene
two-phase (Method 2)
air-organic -
organic— water —
pooled Multiphase
Sm*0 (Method jj) 9.288 * 1.69
Sa=0 (Method Hi) 9.265 * 1.63
hd/Pao
4.219 *
4.
4.
5.
7.
-
184 *
184 *
017 *
-
726 *
7.719 *
5.
420 *
5.288 *
5.286 *
bd/low
.734
1.45
1.43
2.908 *
4.355 *
4.355 *
.963
1.55
1.53
1.50
1
1
1
.80
.79
.48
1.05
1
4.763 * .
-
4.643 *
4.631 *
1
1
.05
819
.25
.21
3.114 *
1.520 «
1.519 *
4.320 *
4.425 *
4.424 *
—
4.565 *
4.719 *
4.700 *
.386
.403
.403
.579
.876
.873
.888
1.30
1.24
A
.532 *
.521 *
.372 *
.536 *
.536 *
.401 *
1.825 *
.552 *
.552 *
.576 *
.534 *
.551 *
.551 *
.574 *
.507 «
.539 *
.534 «
.103
.094
.427
.077
.075
.105
.978
.069
.068
.513
.069
.049
.049
.517
.101
.060
.054
s.
1Q-16 « 1Q-32
1Q-15 « .032
10~5 ,
10-16 ,
0
10-16
.255
10-16
0
f*
10-7
10-16
10-"
0
f
10~5
10-16
10-15
fc .336
t 10-32
.0
* 10-31
* .109
» 10-32
.0
* .206
ft 10-32
< 10-32
.0
* .228
* 10-32
* .019
0.0
R2
0.912
0.932
0.841
0.901
0.901
0.848
0.885
0.894
0.894
0.924
0.959
0.931
0.931
0.924
0.908
0.915
0.916
* units, of hcj/0ij in on; other variables dimensionless
* 95X confidence interval from regression analyses and by method of Rao et al. (1977) for the pooled data
-------
Saturation-capillary head data are shown for the air-water system (Figure 3.1)
and for the air-ortfanic and organic-water systems for benzene and benzyl
alcohol (Figures 3.2-3.5) along with the curves corresponding to the
Brooks-Corey and van Genuchten functions fit to the pooled air-water,
air-organic, and organic-water data with Sm = 0. It is noted that the
pooled-fit models sacrifice very little accuracy in describing the individual two
phase systems. The scaled van Genuchten and Brooks-Corey models both
describe the data reasonably well, although the former was in all caaea
slightly superior (see Ra values, Tables 3.1 and 3.2). The curves are in most
instances encompassed by the experimental scatter arising from the difficulty
of packing triplicate cores to obtain exactly the same pore size distribution.
It is interesting to note that there was more scatter in the air-water
measurements than in the other fluid pair measurements. This may be a
consequence of the greater difficulty in achieving a uniform initial degree of
saturation for the air-water system due to the much higher interfacial tension
between air and water than any of the other fluid pairs. Among the organic
liquid systems, the greatest relative error occurred for the benzyl
alcohol-water system (Figure 3.5) which was observed to drain at very low
capillary heads due to the very low interfacial tension for the system. A
consequence of the low interfacial tension is that a large volume of fluid may
drain with a small change in the capillary head leading to experimental
difficulties and reflected in the lower Ra values.
160
a
CM
I
E
0
LJ
I
cr
<
_i
_i
i—i
a.
<
u
120 -
SO -
40 -
WATER-AIR SYSTEM
—^— van Ganuclltan
—— BrooK«-Coray
0. O O. 2 0.4 0.6 0.8 1.0
SATURATION
Figure 3.1. Brooks-Corey and van Genuchten retention functions fit to two
phase air-water data.
26
-------
1QQ
/•s
o
(M
I
E
U
LU
X
a
J
I—I
a.
<
u
75 -
100
a
(M
I
E
U
a
<
u
i
d
HI
Q_
<
U
75 -
SO -
BENZENE-AIR SYSTEM
—— r«gra«sion
_— pradictad
Q. 0 0. 2 O. 4 0. 6 Q. 8
SATURATION
1.0
BENZENE-WATER SYSTEM
—— ragra««lon
—— predicted
>8rooh«-Cor«y
O. Q O. 2 O. 4
0. 8
O. 8
1. 0
SATURATION
Figure 3.2. (a) Measured benzene-air data and curves for Brooks-Corey and
van Genuchten multiphase retention functions using parameters
fitted to the pooled data set and those determined from air-water
data with scaling factors predicted from interfacial tension data.
(b) same for benzene-water system.
27
-------
BENZYL ALCOHOL-AIR
—— ragrammion
pradie-tad
van Ganuchtan
Sc-ooK«—Cer-a
O. 2 Q. 4 0.6 0.8 l.O
0.0
O. Q
SATURATION
BENZYL ALCOHOL-WATER
ragi-o«»lon
— — pradletad
0.2 0.4 0.6 0.8 1.0
SATURATION
Figure 3.3. Measured (a) benzyl alcohol-air and (b) benzyl alcohol-water data
and multiphase retention functions as in Figure 3.2.
28
-------
Table 3.3. Measured interfacial tensions by drop method for the different
fluid systems; values in mN m~l as means * one standard deviation
at 22.5 * 1.5*0.
Liquid
benzene
o-xylene
p-cymene
benzyl alcohol
water
Air-Liquid
35.0 t 1.8
30.6 * 1.1
34.3 * 1.1
36.3 * 1.4
66.4 * 1.3
Liquid-Water
32.0 * 0.2
28.6 * 0.5
31.8 * 2.6
7.8 * 1.1
-
Table 3.4. Scaling factors obtained from regression analyses of pooled
two-phase saturation-capillary head data and from measured
interfacial tensions.
Fluid System
benzene
0OW
o-xylene
lao
p-cymene
benzyl alcohol
Paw
Pij from pooled regression1 0jj from measured
van Genuchten Brooks-Corey interfacial tensions2
2.18 * 0.34
1.98 * 0.29
2.16 * 0.32
1.96 * 0.29
2.01 * 0.26
2.25 * 0.30
1.35 * 0.26
5.96 * 1.09
2.22
2.13
2.00
1.99
1.80
2.16
1.24
6.29
* 0.66
* 0.66
* 0.39
* 0.39
* 0.26
* 0.31
* 0.21
* 1.34
1.90
2.07
2.17
2.32
1.94
2.10
1.83
8.51
t 0.27
* 0.11
* 0.23
* 0.16
t 0.18
* 0.46
* 0.21
* 3.09
1 Mean * 95% confidence limits from estimation uncertainty in regression.
2 Mean * 95% confidence limits from first-order analysis of experimental error
by method of Rao et al. (1977).
29
-------
Scaling factors calculated via (3.5) using measured interfacial tension data are
compared in Table 3.4 with values determined by nonlinear regression
analyses. Scaling factors estimated from measured interfacial tensions compare
closely with those determined from the pooled regression analyses of the van
Genuchten and Brooks-Corey functions except for the benzyl alcohol system.
In the latter case, the discrepancies reflect wide confidence limits on
parameter values due to difficulties in taking precise measurements of
interfacial tensions and saturation-pressure relations for this fluid. In all
cases, confidence regions of measured and predicted ftj overlap well within
95% probability limits. Predicted van Genuchten and Brooks-Corey multiphase
retention functions using parameters determined from air-water S-P data and
interfacial tensions are compared with the measured saturation-capillary head
data for the benzene and benzyl alcohol systems, which typify best and worst
case agreement between measured and estimated parameters, in Figures 3.2
and 3.3. The results indicate that air-organic and organic-water S-P curves
may be predicted with reasonable precision using scaling factors estimated
from interfacial tension data. As a result, the material characterization
problem is greatly simplified, since one need only directly measure the
saturation-capillary head function of a single two phase system to predict that
of any other two phase system. By the extension of Leverett (1941), S-P
relations of three phase systems may be similarly predicted.
SUMMARY AND CONCLUSIONS
The format proposed in Chapter 2 to scale saturation-capillary pressure
relations of two phase air-water, air-organic and organic-water porous media
systems was evaluated by applying the procedure to relations measured for
four organic liquid systems in a sandy porous medium. Multiphase versions of
the Brooks-Corey (1964) and van Genuchten (1980) retention functions were fit
to the experimental data using nonlinear regression analyses. The results
validate the proposed scaling procedure in which it is assumed that the
paramters a, n and Sg, or hj, X and SB, are porous medium constants while
Pao and Pow depend on the organic fluid properties (when 0aw=1 » assumed).
Relatively good fits of both functions to the experimental data were obtained,
although the van Genuchten function yielded slightly higher correlation
coefficients in all cases.
A procedure was outlined to obtain parameters of retention functions which
are used to predict saturation-pressure relations in three phase systems that
preserves a unique pore size distribution for rigid porous media. In addition,
it was shown that fluid-dependent scaling factors used to scale saturation-
capillary pressure relations of different two phase systems in the sandy
porous medium investigated are well approximated from ratios of interfacial
tensions provided the latter are measured using fluids subject to any
impurities occurring' in the actual porous medium system. Using scaling
factors from interfacial tension data permits prediction of saturation-capillary
pressure relations for any two phase fluid system in a porous medium from
direct measurements of a single two phase system and appropriate interfacial
tension data. Via the assumptions of Leverett (1941), three phase system
behavior may likewise be predicted.
30
-------
CHAPTER 4
MEASUREMENT AND ESTIMATION OP THREE PHASE
SATURATION-PRESSURE RELATIONS
The extension of two-phase S-P relations to predict three phase behavior was
first suggested by Leverett (1941). Regarding interfacial curvatures for a
monotonic displacement process to be fixed by the pore-size distribution of
the medium, Leverett (1941) assumed effective total liquid saturation, Sj, to be
a function of interfacial curvature or capillary pressure at the gas-liquid
interface, independent of the number or proportions of liquids contained in
the porous medium. In three phase air-NAPL-water systems in which water is
the dominant wetting fluid, the total liquid saturation would thus be a
function of air-NAPL capillary pressure, Poo-P^-Pot where capillary pressure
is defined as the difference in pressures between contiguous nonwetting and
wetting fluids [i.e., S£=f(Pao)]. It is also commonly assumed (Aziz and Settari,
1979), where fluid wettabilities follow the order water>NAPL>air, that effective
water saturation in an air-NAPL-water system is a function of the NAPL-water
capillary pressure, POvrpNAPL>air until oil saturation diminishes to a level at
which the oil phase exists only as discontinuous blobs or pendular rings at
particle contacts. Thus, air-water capillary pressures are presumed to have
no physical relevance as long as a continuous oil phase exists.
Although (2.16) and (2.17) represent the most common method of depicting
three phase S-P relations (e.g., Sheffield, 1969; Kazemi et al., 1978; Coats, 1980;
Aziz and Settari, 1979), other formulations have been suggested (e.g., Abriola,
1984; Peery and Herron, 1969; Shutler, 1969). Very little experimental research
has been conducted to directly test the assumptions implied by (2.16) and
(2.17). Eckberg and Sunada (1984) concluded from an air-oil-water flow
experiment that water saturation was a function only of oil-water capillary
pressure in agreement with (2.16). Total liquid saturation was inferred not to
be exclusively a function of the air-oil capillary pressure in contradiction to
31
-------
(2.17); however, oil phase pressures were not directly measured and possible
effects of hysteresis were ignored. If water and total liquid saturation path
histories in two and three phase experiments do not correspond, effects of
hysteresis and fluid entrapment will substantially complicate evaluation of
(2.16) and (2.17) as will be discussed in detail in Chapter 8.
Our purpose in the present chapter is to present a technique that can be
employed to directly measure three phase air-NAPL-water S-h relations in
unconsolidated porous media. Details of the experimental apparatus and
procedures used to measure three phase S-h relations are discussed. Results
of three phase measurements will be compared to two phase data to test
assumptions of extending two phase S-h relations to three phase systems as
postulated by Leverett (1941) and others. The scaling format described in
Chapter 2 will be evaluated to test whether a single scaled function can
adequately describe S-h relations in two and three phase systems.
METHODOLOGY AND MATERIALS
Fluid and porous medium properties
Characteristics of the three soil materials employed in this study are given in
Table 4.1. The NAPL was Soltrol 170 which is a mixture of branched alkanes
with carbon numbers ranging from CIO to CIS having very low solubility in
water and a liquid density of 0.785 Mg m~3. Air-water, air-oil and oil-water
interfacial tensions, as measured by the drop method (Adamson, 1967), were
67.0, 24.2 and 44.2 mN m'^x-espectively.
Table 4.1. Particle size distributions and "packed bulk densities of porous
media used in experiments.
Mass fractions in size classes (MI) Packed
Soil —————————————~~~~—•—~~~ density
no. 2-1 1-.5 .5-.25 .25-.! .1-.05 .05-.002 <<.002 Mg aT3
1 3.3 40.3' 37.0 14.0 2.6 0.8 1.9 1.65
2 5.4 10.1 72.6 10.4 0.7 0.4 0.3 1.60
3 0.6 38.0 10.0 14.9 9.9 21.1 5.5 1.55
32
-------
Three phase S-h measurement apparatus
Before describing the three phase retention cell which is capable of
characterizing three phase Sa(ha,h0,hw)t ^(ha^ho.h^) and ^(ha/ho.hwr) or two
phase Sa^h^hw), Sw**(hah^, S^h^b^ and So^h^ho) relations of
unconaolidated porous media, a brief review of theory utilized in two phase
S-h measurements is apropos. In two phase S-h measurements, a porous
material which acts as barrier to the nonwetting fluid phase is utilized to
obtain a continuous wetting fluid phase extending from the porous medium to
an external reservoir where the wetting fluid pressure can be measured or
regulated at a particular level. Measurement or control of the nonwetting
phase is achieved via direct contact between nonwetting phase in the porous
medium and an external fluid supply. For measurements of three phase S-h
relations for systems in which fluid wettabilities for solid surfaces are time
invariant and follow the order water>NAPL>air, continuity of two liquid phases
and a gas phase needs to be established with meaaurment/control devices.
To establish continuous water and oil (Soltrol 170) phases between the porous
medium and external reservoirs in the three phase measurement apparatus, we
employ untreated and treated porous ceramics, respectively. Employing
untreated porous ceramic barriers enables maintenance of a water saturated
condition in the barrier providing continuity with the water phase so long as
the difference in pressures between oil and water phases does not exceed the
oil entry capillary pressure of the ceramic which is a function of the maximum
pore size of the ceramic. To obtain a continuous oil phase extending from the
porous medium to a fluid pressure measuring device located outside the
measurement cell, we silanized porous ceramic rings with chlorotrimethysilane
to create hydrophobic ceramic surfaces. Ceramic sections were immersed in
chlorotrimethysilane for approximately 2 hours followed by thorough rinsing
with toluene and then methanol. Reaction of chlorotrimethylsilane with silica
ions on the ceramic surface forms durable Si-O-Si bonds, and the trimethyl
group, which is attached to one of the silica ions, repells water molecules in
preference to oil and air (Ogden, 1985). Caution needs to be exercised when
working with chlorotrimethysilane because it will react violently with water.
The three phase measurement cell was constructed of alternating plexiglas
sleeves containing treated and untreated porous ceramic rings on their inner
surfaces. A schematic of the retention cell is shown in Figure 4.1. Porous
ceramic rings with an internal diameter of approximately 4 cm were affixed in
machined plexiglas sections with an epoxy resin. Two hydrophobic (treated)
and two hydrophilic (untreated) ceramic plexiglas sections were bolted
together in an alternating pattern with silicons rubber in the joints. The
total height of the retention cell was 8.5 cm. The NAPL we employed did not
interact with the plexiglas so as to prohibit use of the latter. For organic
liquids that might react with plexiglas, stainless steel or some other inert
machinable material could be employed.
Brass fittings and tubing connected the plexiglas sections of the three phase
retention cell to respective pressure transducers, burets and vacuum/pressure
regulators (Figure 4.1), This configuration, which is a modification of a two
phase technique developed by Su and Brooks (1980), permitted NAPL and
aqueous phase pressures to be independently controlled through the
33
-------
vacuum/pressure regulators and permitted monitoring the flow of liquid phases
into or out of the porous medium packed in the retention cell. Liquid phase
pressures were measured with pressure transducers after the 3-way stopcocks
below the burets are positioned to seal the b'quid phases in the porous
medium from the burets which are regulated at negative pressures. The three
phase retention cell is placed in a chamber with containers of free NAPL and
water to minimize evaporation of liquids from the porous medium. Narrow
pathways connect the interior of the chamber to the ambient atmosphere to
ensure the gas phase remains at atmospheric pressure.
IS MOM.
© MMIUMtOIMCI
0 VACUUM IOIMCI
(Q- J. WAT ITOrCOCK
O mittUM TDAMIOUCM
© MfllUMflUCUUU MOUUTM
II mL WATCH
t
tNMC FHAM
MUNTttN CfU
SCHEMATIC OF THREE-PHASE FLUID
PRESSURE-SATURATION APPARATUS
FITTINGS
To TRANSDUCER.
BURET AND
hESSURE/VACUUM
REGULATOR *•:
POROUS MEDIA
DETAIL OF CELL SEGMENT
FITTING
TO TRANSDUCER,
BURET ANO
PRESSURE/VACUUM
REGULATOR »-;
CERAMIC RING
» TO TRANSDUCER,
BORETANO
PRESSURE/VACUUM
REGULATOR
PLEXIGLASS CAVITY
figure 4.1. Three phase saturation-pressure apparatus: system configuration
(top) and enlargement of cell design (bottom).
34
-------
Two and three phase S-h measurement procedures
Procedures used to measure three phase air-oil-water S-h relations involved
draining water from an initially water saturated porous medium to an arbitrary
water content, adding oil so the porous medium was resaturated with liquid,
i.e., Sw + S0 = 1, and then draining either water or oil. This technique
yielded two phase air-water S-h measurements corresponding to draining
water from the initially water saturated condition [Swaw(haw)] and a single
two phase oil-water S-h data point lSwow(how)1 when the porous medium was
saturated with oil and water. Water saturation as a function of oil-water
capillary head [Swfhow)) and total liquid saturation as a function of air-oil
capillary head [St^ao)) in the three phase system were obtained as water or
oil subsequently drained. Water or oil were independently drained by
adjusting the vacuum/pressure regulators.
Water initially was allowed to imbibe into the packed air-dry porous media
slowly through the untreated ceramic section located at the bottom of the
retention cell. We adopted this procedure to minimize gas entrapment on
wetting fluid imbibition saturation paths. Air could easily be displaced from
voids in the porous media as the water front- moved upwards. The apparently
water saturated packed retention cell was left undisturbed for approximately
18 hours to permit entrapped air to diffuse out of the sample.
Assuming the air phase to be everywhere at atmospheric pressure,
measurements of two phase air-water system saturations were determined by
recording the outflow volume with a buret as water drained and water
pressure was measured by recording the water transducer output voltage
after the 3-way stopcock was positioned to close off the buret from the cell
arresting water drainage. Liquid phase pressures stabilized relatively quickly
after the stopcock was turned. To ensure a monotonic drainage saturation
path, the time period between initiating liquid drainage and turning the
stopcock needs to be long enough so that redistribution of fluids in the
porous medium after turning the stopcock is minimized (i.e., static equilibrium
is approached before positioning the stopcock). This sequence was repeated
numerous times along a saturation path.
The "porous media were drained of water to arbitrary water saturations at
which point Soltrol 170 was slowly added to the air-water-porous media system
through the lower treated ceramic section to minimize air entrapment. Organic
liquid was added until the porous media were apparently saturated with liquid
(i.e., water and organic) and a period of approximately 12 hours was allowed
to elapse before measurements were taken to facilitate diffusional outgassing.
The fluid-porous media system at this juncture is only a two phase oil-water
system in which the water content corresponds to the water saturation when
drainage of water in the two phase air-water system concluded.
Measurements of S^hyh^h^, S^h^hoh^ and Sdha.ho.h*) relations were
accomplished by monotonically and incrementally draining oil or water from the
total liquid saturated condition. Oil and water saturations were determined by
recording the volumetric outflow of the liquids between each measurement
point. Oil and water phase pressures were recorded after the stopcocks were
turned and after volt readings of the oil and water pressure tranducers
35
-------
stabilized. Experimental procedures used in three- phase measurements were
similar to procedures used for the air-water measurements previously outlined.
This methodology produced a monotonic drainage total liquid saturation path.
Two phase air-oil measurements were also conducted with the three phase
retention cell. Procedures used in two phase air-oil measurements were the
same as those for two phase air-water and three phase measurements. Two
phase air-oil S-h measurements were necessary to evaluate the assumption of
extending Soac^ao) relations to predict total liquid behavior in three phase
air-oil-water systems as proposed by Leverett (1941). A total of 108
measurements from 7 packings of the retention cell (5 packings for two phase
air-water and oil-water and three phase measurements, 2 packings for two
phase air-oil measurements) were conducted for Soil 1 of which 48 were two
phase air-water, oil-water or air-oil measurements and 60 were three phase
measurements. A total of 95 measurements from 5 packings of the retention
cell were conducted for Soil 2 and 37 measurements from 3 packings were
completed on Soil 3. Of the 95 and 37 S-h measurements with Soils 2 and 3,
46 and 14 were three phase measurements, respectively.
RESULTS AND DISCUSSION
The validity of the assumptions embodied in (2.16) and (2.17) for predicting
three phase behavior from two phase data may be tested by analysis of the
experimental results from the three phase measurement cell. Assumption (2.16)
states that the relationship between water saturation and oil-water capillary
head in a two phase oil-water system can be utilized to predict water
saturation in a three system. We may evaluate this assumption by comparing
measurements of two phase Swow(how) relations to three phase Sw(hoW)
measurements. The 5 measured two phase oil-water S-h data points and 30
three phase Sw(JiOMr) points from 5 packings of Soil 1 are compared in Figure
4.2. Scatter in the experimental data is about the same for both two and
three phase measurements and correspondence is close. Also shown is a
retention curve obtained by fitting the measured two phase oil-water S-h data
to the van Genuchten (1980) function. Measured three phase water saturations
as a function of oil-water capillary head closely parallel the two phase
oil-water retention curve. This indicates that the assumption of extending
swow(how) relations to predict water saturations in three phase porous media
given by (2.16) is valid under the experimental conditions imposed.
The second assumption in extending two phase S-h data to three phase
systems is that liquid saturations are solely a function of air-oil capillary
head as stated by (2.17). Figure 4.3 shows two phase oil saturations and
three phase total liquid saturations as a function of air-oil capillary head
measured for 7 packings of Soil 1. Also shown is a retention curve obtained
by fitting the van Genuchten (1980) function to the two phase air-oil S-h
data. Three phase total liquid saturation measurements closely bound the
retention curve fitted to two phase data. Experimental scatter in three phase
measurements is not significantly different than scatter in two phase
measurements and three phase data are approximately equally distributed on
both sides of the two phase air-oil retention curve. The experimental data in
36
-------
Figure 4.3 corroborate the assumption advanced by Leverett (1941) that two
phase S0ao(hao) data can be employed to predict three phase St(hao) relations
in rigid porous media for monotonically draining total liquid saturation paths.
To evaluate the applicability of the scaling procedure introduced in Chapter 2
to the combined set of two and three phase experimental S-h data of Figures
4.2 and 4.3, scaling coefficients were fit to the pooled data with the van
Genuchten equation to describe S*(h*). The fitted S*(h*) function and the
scaled two and three phase data are shown in Figure 4.4. The scaling
procedure appears to produce a unique S*(h*) function which can be employed
to describe S-h relations in two phase air-water, air-oil and oil-water systems
and in three phase air-oil-water systems.
In Chapter 3 we demonstated that scaling coefficients may be estimated from
ratios of interfacial tensions. This suggests that three phase air-oil-water
S-h relations may be predicted from measurements of two phase air-water S-h
relations and interfacial tension data, which would significantly reduce the
experimental effort required for model calibration. The van Genuchten
retention function fit to two phase air-water data only is shown in Figure 4.5
along with scaled two and three phase S-h data for Soil 1 using fcj s
calculated from measured interfacial tensions. Comparing Figure 4.5 to Figure
4.4 indicates that using air-water S-h relations to characterize parameters of
the van Genuchten (1980) function and ratios of measured interfacial tensions
to estimate ftj's results in a less accurate description of the system as
evidenced by greater scatter in Figure 4.5. On the other hand, considering
the great reduction in data requirements for the second calibration method,
thin loss in precision may well be an acceptable tradeoff.
Two and three phase S-h data for Soil 2 are presented in Figure 4.6 in scaled
form along with the fitted S*(h*> function. Model parameters were determined
by fitting to the entire S-h data set in the same manner employed in Fl*ufe
4.4 for Soil 1. The greater scatter in the experimental data m Figure 4.6 IB
attributed to difficulty in reproducably packing this porous medium due to its
narrow grain size distribution (see Table 4.1). Small deviations in the packing
procedure can result in significant differences of measured S-h relations
between replicates. Despite this greater scatter, the mean behavior is well
described by the scaled function indicating that the scaling assumptions and
the assumptions pertaining to extending two phase S-h relations to predict
three phase behavior are both closely obeyed for this data.
Figure 4.7 shows S-h data for Soil 3 scaled by best-fitting model parameters
to the pooled data set. Scatter in the measured data is low and agreement
with the fitted S*(h*) function is good, again corroborating the model. It may
be observed that the range in fluid saturations spanned in the experiments on
Soil 3 were rather narrow. At a scaled capillary head of 150 cm in Figure 4.7,
over half of the void volume still retains liquid. This reflects, of course, the
finer particle size in Soil 3 than in Soils 1 and 2, but illustrates a limitation
in the experimental retention apparatus to relatively low pressures. At room
temperatures, gas captation will generally occur when the pressure ^posed
on water is less than -300 to -400 cm H,0. These air bubbles become lodged
in the tubing and affect liquid levels in the bureta.
37
-------
u
LJ 24 -
I
16 •
CL
<
U
tc
UJ
0
TWO PHASE < • )
OIL-WATER AND
THREE PHASE C O >
WATER SATURATION
MEASUREMENTS
0. O 0. 2 0. 4 0. 6 O. 8
WATER SATURATION
l.O
Figure 4.2. Comparison between water saturation in two phase oil-water
systems (closed squares) and water saturation in three phase
air-oil-water systems (open squares) as functions of oil-water
capillary heads for Soil 1.
Ul
I 16
Q_
U
_l
f*
O
I
(T
8 -
TWO PHASE
AIR-OIL. (•>
AND THREE PHASE
TOTAL LIQUID <• >
MEASUREMENTS
O. O O. Z O. 4 O. 6 O. 8
LIQUID SATURATION
1. O
Figure 4.3. Comparison between oil saturation in two phase air-oil systems
(closed diamonds) and total liquid saturation in three phase
air-oil-water systems (open diamonds) as functions of air-oil
capillary head for Soil 1.
38
-------
O. O O. 2 O. 4 Q. B O. 8
SCALED SATURATION
1. o
Figure 4.4. Scaled wetting fluid saturations of two phase air-water (stars),
air-oil (open diamonds) and oil-water (open squares) systems, and
water (closed squares) and total liquid (closed diamonds)
saturations of three phase air-oil-water systems as functions of
scaled capillary heads for Soil 1 with best fit model parameters.
0. O 0.2 0.4 0.8 0.8 1.0
SCALED SATURATION
Figure 4.5 Scaled wetting fluid saturations of two phase systems, and water
and total liquid saturations of three phase systems as functions
of scaled capillary heads for Soil 1 with model parameters
estimated from two phase S-b data and interfacial tensions. Same
legend as Figure 4i4.
39
-------
o. o o. a o. 4 o. a o. a
SCALED SATURATION
i.o
Figure 4.6. Scaled wetting fluid saturations of two phase systems, and water
and total liquid saturations of three phase systems as functions
of scaled capillary head in Soil 2 with best fit model parameters.
Same legend as Figure 4.4.
jao
E
O
a
UJ 120
(£
<
_J
a.
<
u
o
ta
_i
5
so -
a. a Q.Z a.4 a. s a. a
SCALED SATURATION
i.o
Figure 4.7. Scaled wetting fluid saturations of two phase systems, and water
and total liquid saturations of three phase systems as functionu
of scaled capillary head in Soil 3 with best fit model parameters.
Same legend as Figure 4.4.
40
-------
SUMMARY AND CONCLUSIONS
An experimental apparatus was developed to measure three phase air-oil-water
S-h relations in unconsolidated porous media. The apparatus is also capable
of measuring S-h relations of two phase air-water and air-oil systems.
Measurements of three phase S-h relations were conducted and compared to
two phase S-h measurements for three porous media to directly test
assumptions involved in the prediction of three phase S-h relations from two
phase measurements. Good agreement between total liquid saturations in three
phase air-oil-water systems and oil saturations of two phase air-oil systems as
functions of air-oil capillary head was observed. Close agreement was also
observed between water saturation versus oil-water capillary head relations
measured in two and three phase systems.
Agreement between two and three phase S-h relationships for monotonic
saturation paths imply that researchers modeling flow of continuous oil and
water phases in the vadose zone may employ more readily available two phase
S-h measurements to predict flow phenomena in three phase air-oil-water
systems if fluid saturations are assumed to be unique functions of capillary
head (i.e., no hysteresis). We have not considered here the behavior of
systems subjected to nonmonotonic water and total liquid saturation paths
which will be discussed in Chapter 8.
In addition to the proviso that results have only been presented for
monotonically draining water and total liquid saturation paths, the procedures
for predicting three phase S-h relations are expected to be valid only for
strongly water-wet and rigid porous media. Soils and aquifer materials are
generally expected to meet the criterion of being strongly water-wet, and at
least in relatively coarse grained materials the rigid medium constraint may be
sufficiently closely met as well. In finer textured porous media, the rigid
medium assumption is most prone to cause deviations from theory, and more
extensive investigations for materials with a wide compositional range is
warranted in the future.
The scaling procedure that we have advanced combined with a suitable
parametric representation of the scaled retention function enables development
of a multiphase retention model which can be used to describe two and three
phase S-h relations. The results indicate that calibration of the model by
fitting to two phase air-water, air-oil and oil-water S-h data yield an accurate
description of two and three phase drainage S-h relations. If two phase
air-oil and oil-water S-h relations are not available, model parameters can be
estimated from two phase air-water S-h data and relatively simple
measurements of air-water, air-oil and oil-water interfacial tensions with only
moderate deterioration in model precision.
41
-------
CHAPTER 5
FINITE ELEMENT MODEL FOR MULTIPHASE FLOW
Much of our current knowledge of the physics of multiphase flow and of
pertinent numerical solution methods has come from studies in petroleum
reservoir engineering (e.g., Aziz et al., 1979; Lewis et al., 1984). However,
substantial differences in specific conditions (e.g., displacing fluid pressures,
geometry) pertinent to petroleum reservoir and water resources problems often
precludes cross-application of numerical codes between these fields. Solutions
for simultaneous flow of air and water in the vadose zone have been
presented by Morel-Seytoux (1973) and Touma and Vauchlin (1986) among
others. Numerical analyses of the movement of NAPL arising from organic
waste disposal or leaking underground storage facilities have been presented
recently by several authors. Faust (1985) analyzed the simultaneous flow of
water and NAPL in a three fluid-phase system under the assumption of a
static gas phase at atmospheric pressure using a finite difference method in a
vertical two-dimensional domain. A similar flow problem, coupled with
component transport for a slightly soluble and volatile NAPL, was analyzed by
Abriola and Finder (1985a,b) also using a finite difference method. Finite
element methods have been applied to two-phase flow problems by Huyakorn
and Pinder (1978) and Osborne and Sykes (1986).
In the present paper, we present a solution for liquid flow in a three
fluid-phase porous medium system in two space dimensions with gas at
constant pressure using a variational finite element method. The k-S-P
relations will be described by the^model which was introduced in Chapter 2
a"d experimentally validated for static S-P data in Chapters 3 and 4.
Hypothetical field-scale problems will be presented in the present chapter and
in following chapters model validation for transient flow conditions will be
investigated.
MODEL FORMULATION
Governing equations for multiphase /low
In this paper we restrict our attention to the analysis of multiphase flow with
no interphase transfer. Extensions to consider concomitant flow and
component transport with interphase transfer will be presented later (Chapter
9). For an incompressible porous medium with constant fluid densities and
viscosities, the flow equation for liquid phase i can be written as
»o .
(5.1)
*t
+ is the porosity, Si is the saturation of phase i (i = w,o where w
denotes water phase and o denotes NAPL), kri is the relative permeability of
42
-------
phase i, TJJ is the fluid viscosity, p{ is . the fluid density, Pj is the phase
pressure, k is the instrinsic saturated permeability tensor, g is the scalar
gravitational acceleration, z is elevation, and t is time.
It is convenient to write (5.1) in the form
* 7T = v ' (Ki Hi) (5.2)
where H{ is the total head in equivalent height of water for fluid i defined by
H = h +
in which pw is the density of water, hj = P{/gpw is the fluid pressure head in
equivalent water height) and the fluid conductivity tensor is defined here by
Ksi =
We consider here the occurrence of three immiscible phases, namely air, water
arid NAPL, in a two-dimensional domain under the assumption that the air
phase remains at atmospheric pressure and the principal axes of the
permeability tensors are aligned with the Cartesian coordinates. Noting that
Sj = Si(hyhj) for ij = o,w and applying the chain rule to the left-hand side of
(5.2) yields the coupled differential equations for movement of two mobile
fluids as
CWO 7j- + GWW — = — (Kxw — ) + — (Kyy, — ) (5.3)
coo g * Cow a. .•- (Kxo a,, + •- (Kyo a. , „.„
where Kxj = kriKxai and Kyj = kriKySj(i = o,w) are the components of Kj in the
x and y directions with* saturated values Kxsi and Kyai and Cij (ij = o,w) are
fluid capacities defined by
*s
To describe the relationships for fluid capacities and conductivities as
functions of water and NAPL pressures the parametric model descibed in
Chapter 2 is employed.
43
-------
Finite element formulation
A variety of methods may be used to derive finite element formulations for the
multiphase flow problem posed including variational and weighted-residual
methods (e.g., Galerkin's). Here we utilize the variational method and derive
the quadratic functionals corresponding to (5.3) and (5.4) as
Kv
dxdy
(5.5a)
-------
where Nj are the bilinear shape functions and Hwj and Hoj are the nodal
heads. The derivative of the shape function with respect to the spatial
coordinates is
[B] = ['Ni/ax *Ni/«y]T (5.8)
Substituting (5.7) and (5.8) into (5.5) and taking variations, yields the element
equations in matrix form as
[K]w{Hn}w + [Ktw]w{Hn}w + [Ktw]o{Hn}o = Ww (5.9a)
[K]0{Hn}0 + [Kto]W{Hn}w + [Ktolo{Hn}o = Wo (5.9b)
where the overdot represents differentiation with respect to time. The
conductivity matrix [K] and capacitance matrix [Ktl for the water phase are
given by
"" [B]dxdy (5.10a)
[Ktwlw = JJyC^tNjTtNldxdy (5.10b)
[Ktolw = fJyCwo[N]T[N]dxdy (5.10c)
and those for the NAP1 phase by
[K]0 = Jx}y[B]T X° [B]dxdy (S.lla)
[Ktwlo = Jx}y(WN]T[N]dxdy (5. lib)
[Ktolo = (Jc0o[N]T[N]dxdy (5.lie)
, '"' y
arid load vectors for the water and NAPL phases are ffiven, respectively, by
(5.12)
{Q}0 =
45
-------
By coupling the water and NAPL phase equations we can write
[K]{Hn} + [Kt]{Hn} =
(5.13)
where
[K] =
'[K]w
.[0] [K]<
[Kt] =
{Hn} =
HOt Hoa HO3 Ho»}
••••—,
^01 Hoa HO3 HO4}
Qoi
Evaluation of the integrals in (5.10) and (5.11) is performed using four-point
Gauss quadrature.
Time integration and nonlinear iteration techniques
Time derivatives are analyzed by the finite difference method to provide a one
step time integration scheme. Applying the general '8' algorithm (Belytschko
and Liu, 1983) to (6.13) yields
- e)[Kt]t+At[Kt]t-'[K]t
OUt
(5.14)
where 8 is a time-weighting factor, 0 * 8 * 1, At is the time increment, and
subscripts denote the time at which variables are evaluated.
For transient flow problems, [KtJt * [KtH+At and (Q)t is not necessarily equal
Jt° CQJt+At' To simplify the calculation 8 = 1 is used in this study leading to
the fully implicit backward scheme given by
46
-------
t + Wt-nt (5-15)
In (5.15), [Kth+At and [Klt+At are functions of {Hn}tfAt« Thus (5.15) is
nonlinear and requires an iteration technique to handle the nonlinearity
numerically. Direct iteration (Picard method) or tangent methods
(Newton-Raphson) are commonly employed. Direct iteration methods are simple
in numerical operation, but for severe nonlinearities may be less efficient than
tangent methods due to slower convergence (Huyakorn and Finder, 1983). In
this study, a direct iteration method is used with modifications to improve
convergence (modified Picard).
The known terms in (5.14) are {Hn)t (Q)t+At and *t. The unknown terms
[Kth+At anc' tKH+At are required in order to solve for (Hn}t+Af Tne iteration
starts from the initial estimates [Ktlt+At* and [Klt+at1 based on {Hn)t where
superscript 1 indicates the first iteration. After {HnJt+At1 *a solved from
(5.15), [KtJt+At3 and [K]t+*ta are calculated based on a relaxed head vector
{Hn}*= X{Hn}t+Atl + (1 - X) {Hn}t (5.16)
where X is a relaxation factor. [Ktlt+Ata and [KJu^t2 are then employed in
(5.15) to solve for {Hn)t+Ata- Iteration continues with the relaxed head vector
at the jth iteration evaluated as
{Hn}*= A{Hn}tJ-nt + (1 - X) {HnJtJ'1 (5.17)
until meeting the convergence criterion
|HJJ - HiJ~l|t+At< e|Hi|t for i = o,w (5.18)
at all nodes, where t is a small number.
In the two chapters to follow, model validation will be investigated by means
of various one-dimensional transient laboratory studies. In the remainder of
the present chapter, we will investigate two hypothetical numerical problems
involving leakage of immiscible contaminants from underground storage
facilities.
47
-------
HYPOTHETICAL TWO-DIMENSIONAL PROBLEMS
Problem 1
To illustrate a field-scale application of the finite element code, hypothetical
problems involving a leaking buried storage tank are analysed. The storage
tank is 10 m wide and 2 m deep below the ground surface. In the first
scenario, a lighter-than-water benzene derivative, p-cymene (ib/pw = 0.86),
leaks from the tank which is assumed to be continuously filled to the ground
surface with negligible impedance to flow through the tank relative to hat of
the soil. The soil domain is 210 m wide and 16 m deep as shown in Figure
5.1. The base of the soil domain is chosen as the datum in this problem.
Initially, a steady water flow regime is assumed for a water table height of 16
B on the left-hand boundary and 8 m on the right with no flow across
boundaries in the vadose zone and normal to the lower surface. No NAPL is
present in the system initially for either scenario, which is achieved
operationally by making the initial NAPL total head Ife = - 1 m. The assumed
material parameters are given in Table 5.1. A mesh with 140 elements and 169
nodes was used for the analysis (Figure 5.1). Experimentation with finer
meshes and smaller time steps indicated that spatial and temporal
discretization employed was suitable.
The predicted p-cymene plume at 166, 572, and 1590 days is shown in Figure
5.2 with contours denoting the zone in which NAPL saturation is greater than
zero. The general behavior of the plume is as anticipated for an immiscible
organic fluid having a density less than that of water, with lateral spreading
occurring in a high organic content zone 'floating' above the water-saturated
zone. The horizontal velocity of the NAPL plume is observed to reach a
quasi-steady state rate of approximately 0.07 m/day. Closer inspection of the
vertical displacement behavior of the plume reveals some interesting and
unanticipated phenomena. In Figure 5.3 variations in NAPL saturation with
time are shown at two nodes A and B slightly downstream of the tank and
about 5 m beneath the original water table level (Figure 5.1). Due to the
weight of injected organic fluid, the plume displaces the original water table
downward reaching node A at 300 days. NAPL saturation continuously
increases to 0.56 at 750 days and then decreases. That is, initially p-cymene
displaces the water phase at point A, and then the water phase starts to
displace p-cymene. The response at node B is delayed but similar in nature.
This saturation rebound phenomena may be explained as follows. As the water
table is displaced downward, a hydraulic constriction occurs as the aquifer
thickness diminishes and an increased hydraulic head is forced to
develop on the upstream side of the plume. If the rate at which the upstream
water table rises is less than the initial rate of vertical displacement by the
NAPL plume, then the plume displacement will eventually be reversed when the
hyrdrauHc head builds up sufficiently. The occurrence of this rebound
behavior will be highly dependent on boundary and initial conditions and soil
and fluid properties which control the relative rates of vertical and horizontal
NAPL and water phase flow. For example, we find that if the NAPL and water
conductivities in the preceding hypothetical case are reversed (i.e., NAPL less
48
-------
viscous than water) or if a horizontal to vertical anisotropy ratio of 2 is
assumed, then no such rebound behavior is observed. In both of the latter
instances, the diminished resistance to downstream movement .of the NAPL
plume enables the upstream water head increase to keep pace with the the
vertical displacement. We note finally concerning; this phenomenon that when
the three-dimensional nature of the contaminant source is considered, the
damming effect of the NAPL plume on upstream water movement may be
markedly diminished with significant effects on overall system behavior.
Table 5.1. Parameters in the multiphase retention function, saturated
conductivities and fluid densites for two-dimensional leaking
storage tank simulations for Problem 1 with p-cymene and Problem
2 with TCE.
p-cymene
TCE
a
n
SB
Paw
0ao
Pan
Kxao=Kyso
Kxgw=Kysw
Po/Pw
5.0 nTl
1.9
0.0
1.0
2.0
2.2
1.5 in d-'
2.0 m d~l
0.86
6.4 m-1
1.7
0.0
1.0
2.3
3.0
7.7 m d~l
12.0 n hr~>
1.46
-60M-
-HTANK
1
16 M
A
^
^-8M
Figure 5.1. Finite element mesh used for Problems 1 and 2.
49
-------
V\>i%\^_->/l
168 doys ^^572 doys
Figure 5.2. NAPL plume at three times for Problem 1. Contours demark zone
containing nonzero NAPL saturation.
400 800 1200 I60O 2000
TIME (days)
Figure 5.3. NAPL saturation versus time at two locations for Problem 1.
50
-------
25 day*
55 dayi
91 days
134 Jays
lib days
248 day.
649 days
ffp?^:.^;::^:^^^
Figure 5.4. Model predicted TCE plume at aelected timea after initiation of 100
day TCE leakage event from storage tank for Problem
-------
ProbJem 2
In the second problem, the NAPL is taken to be trichloroethylene (TCE) with
soil-fluid properties as given in Table 5.1. The physical domain and initial
conditions are selected to be identical to those of Problem 1. TCE is assumed
to be released from the leaking tank for a period of 100 days, after which no
flow occurs along the atmospheric boundary. The predicted TCE plume at 25,
55, 91, 134, 186, 248, and 649 days after onset of leakage is shown in Figure
5.4 as contours denoting the zone in which organic fluid saturation is greater
than zero. As anticipated for an organic fluid more dense than water
(/»o/Pw=l-46), the plume rather quickly sinks to the bottom of the aquifer
where it spreads horizontally in both the upstream and downstream directions.
Since the lower surface of the aquifer is horizontal, gravitional forces exert
little influence on movement along the bottom and the direction of TCE flow is
controlled by viscous drag exerted by the flowing water phase. The
horizontal velocity of TCE is observed to be about 1/2 that of the water
Phase. It is important to note, however, that with small changes in the slope
°* the lower aquifer boundary, the rate of TCE flow could be substantially
affected ~ even to the extent of reversing the flow direction to move counter
to that of the ground water flow. We point out also that in the wake of the
TCE front after the 100 day injection period, water is predicted to completely
displace TCE. This behavior is a result of the fact that the k-S-P
relationships have been assumed to be reversible and nonhysteretic. In
reality, hysteresis will be observed and entrapment of nonwetting phases (air
and TCE) will occur during water imbibition. Thus, we emphasize that
consideration of hysteretic behavior will likely be of marked importance for
the accurate prediction of immiscible organic movement under conditions of
nonmonotonic saturation paths. Procedures for incorporating hysteresis into
the constitutive model will be discussed in Chapter 8.
SUMMARY AND CONCLUSIONS
A finite element formulation based on the variational approach has been
successfully employed to solve the coupled immiscible flow problem involving a
three fluid-phase system of air, water, and NAPL in soil with air at constant
Pressure. For two hypothetical fluid storage tank problems, the leakage of
organic fluids was analysed and the NAPL plume movement predicted. The
results show the potential of applying the finite element method developed
h«re to the design of monitoring systems and remediation schemes for leaking
underground storage tanks. Extension of the numerical model to consider
interphase mass transfer and convective-dispersive transport of components in
multiple phases will be considered in Chapter 9. Before addressing these
complications, we first turn our attention to validation of the multiphase flow
model under transient conditions.
52
-------
CHAPTER 6
MODEL VALIDATION AND CALIBRATION FOR
TWO PHASE TRANSIENT FLOW
In this chapter, we will investigate three methods of calibrating the k-S-P
model developed in previous chapters involving: (1) direct measurement of S-P
relations for air-water, air-oil and oil-water systems, (2) direct measurement of
air-water S-P relations with scaling coefficients estimated from interfacial
tension data, and (3) numerical inversion of transient oil-water displacement
experiments using a nonlinear optimization procedure. The methods will be
evaluated by analysis of results for two NAPL-porous media systems.
CALIBRATION METHODS
Method 1
In this method two phase S-P relations are determined directly by incremental
equilibrium drainage of wetting fluid controlled by stepwiae adjustment of
capillary pressure in short soil columns (52.6 mm diameter x 51.8 mm length).
Details of the procedure are as described in Chapter 3. Measurements were
carried out for air-water, air-oil and oil-water porous media systems using two
sandy porous media with different NAPL. Studies with Soil 1 utilized
p-cymene as the NAPL and those with Soil 2 involved trichloroethylene (TCE).
Bach two phase stepwise drainage experiment was performed in triplicate on
similarly packed soil cores. To avoid soil dispersion, 0.01 M CaCla was used in
all experiments as the aqueous phase.
Model calibration is achieved in Method 1 by simultaneously fitting parameters
in (2.15) to the combined set of replicated two phase air-water, air-oil and
oil-water saturation-capillary head data using an unconstrained nonlinear
least-squares procedure. Preliminary analyses indicated that for both soils
studied, Sm could not be statistically distinguished from zero. Hence, Sm = 0
was assumed in all analyses thus reducing the unknown parameters to: a, n,
flao, and £,w to be estimated by the regression. The nonlinear optimization
problem is solved using the Marquardt-Levenberg method (Kool et al., 1987).
*
Method 2
In the second method, use is made of the theoretical analysis of scaling
presented in Chapter 3 which indicates scaling factors in a rigid porous may
be estimated from interfacial data. To estimate /Jow and 0ao in this method, we
use interfacial tensions measured by the drop-weight method (Adamson, 1967).
Interfacial tension measurements were performed in quintuplicate and had
coefficients of variation generally within 5X. Measured air-NAPL interfacial
tensions were 34.3 and 27.3 mN m~l for p-cymene and TCE, respectively, and
53
-------
NAPL-water values were 31.8 and 19.7 mN oT1. Air-water measurements
yielded craw = 66.4 mN m"1 for uncontaminated 0.01 M CaCla extracts of the
soils studied.
Values of a and n were obtained by analysis of air-water saturation-capillary
head data for systems uncontaminated by p-cycmene or TCE in a manner
analogous to that employed in Method 1, but using only air-water data. As in
Method 1, preliminary analyses indicated SQ to be indistinguishable from zero
for either soil, so Sjn = 0 waa assumed in the regression analyses.
Method 3
In Method 3, we adopt the same procedure as in Method 2 to estimate scaling
coefficients ftjO and &w» Dut now estimate a and n from transient flow
experiments rather than from static S-P data. The transient experiments
involved soil cores initially at equilibrium close to full water saturation and
subsequently subjected to a step change in NAPL pressure at the upper
surface inducing water displacement. Water saturated ceramic plates at the
lower surfaces of the cores prevented NAPL egress from the system. Core
dimensions were the same as for the static S-P tests. Tests for each soil
were performed on duplicate coes and mean cumulative water outflow versus
time during NAPL displacement was used to infer parameters « and n by
nonlinear regression of observed versus numerically simulated outflow. The
linear finite element solution of the coupled water and oil flow equations of
Chapter 5 was combined with a nonlinear least squares regression code to
analyze the data. The procedure is essentially analogous to the method
reported by Parker et al. (1985) for the analysis of transient gas-driven water
outflow experiments, except that the direct problem formulation in the
air-water case adopts Richard's approximation to avoid explicit consideration of
nonwetting phase flow in the direct problem solution.
Additional input needed to solve the direct problem includes boundary and.
initial conditions pertaining to the experiments and water and oil saturated
conductivities. Boundary and initial conditions for each experiment were
directly measured and thus input into the inversion code. Initial
investigations were conducted to evaluate the feasibility of treating oil and/or
water saturated conductivities as additional variables in the inversion. It was
f°und that sensitivities to these parameters were low, resulting in very large
confidence limits for the estimated values and a general deterioration of the
^ell-posedness of the inverse problem. Accordingly, measured saturated
conductivities were used in the inversions as obtained from falling head water
and oil conductivity tests performed on triplicate soil cores. Means of the
^Plicate measurements were used as known values in the inverse solution.
Mean water and oil saturated conductivities (as defined in Chapter 5) were 23
and 41 cm h'1 for Soil 1 with p-cymene, and 50 and 32 cm h'1 for Soil 2 with
TcE. Coefficients of variation for conductivity measurements were 16% for Soil
1 and 30% for Soil 2, reflecting the greater difficulty in reproducably packing
c°lumns for Soil 2.
54
-------
RESULTS AND DISCUSSION
Results of the various parameter estimation procedures for the two soils are
summarized in Table 1. In addition to the scaling factors &o and £,w, values
of coefficients fow* for organic contaminated air-water systems given by (3.8)
are shown in Table 6.1. The results indicate that the presence of aqueous
and gaseous TCE has a substantial effect on air-water interfacial tension and
thus on air-water saturation-capillary head relations in contrast to the minor
effects for p-cymene. The higher aqueous solubility and volatility of TCE is
expected to induce greater reductions in air-water interfacial tensions than
p-cymene, although theoretical calculations based on methods of Lyman et al.
(1982) predict
-------
Q
U
c
120
•0
4O
AIH-MAICII OAT*
940
a
5 i>o
X
i
-! M
a. o a. a 0.4 a. a a. •
SATURATION
i.o
0.2 a.4 a. a a. a
SATURATION
t.o
100
D
U
-I
•4
0.
u
75-
so •
F-CTNCNC-HATU OATA
Q.O 0.2 0.4 aa 0.0
SATURATION
J
•^
a.
u
40
TCC-MTU OATA
t.O
0.0 0.3 0.4 0.0 o. a
SATURATION
I.O
IOO
I
E
U
w
a
x
cc
a.
u
3-
a
a.
u
0.0 0.2 a.4 ,a.« o.«
SATURATION
OATA
i.a
0.0
1.0
a. 2 0.4 a. a a.•
SATURATION
Figure 6.1. Measured equilibrium two phase S-P data and predicted relations
using; Method 1 (solid lines), Method 2 (dashed lines) and Method
3 (dashed-dotted lines) for: (a) Soil 1 air-water, (b) Soil 1
p-cycmene-water, (c) Soil 1 air-p-cymene, (d) Soil 2 air-water,
(e) Soil 2 TCE-water, and (f) Soil 2 air-TCE systems
56
-------
Model parameters determined by inversion of transient flow experiments via
Method 3 deviate more substantially from those of Methods 1 or 2 (Table 6.1),
although the corresponding differences in predicted saturation-capillary
pressure relations are not very great, especially for Soil 2 (Figure 6.1). In
contrast, observed cumulative water outflow during transient flow experiments
corresponds more closely with numerical simulations using Method 3 parameters
than those for Methods 2 and 3 (Figure 6.2). It is, of course, not surprising
that Methods 1 and 2 agree more closely with equilibrium data and Method 3
with transient data, since the objective functions in the respective
optimization problems are formulated in terms of those data. However,
considering these fundamental differences in parameter estimation methodology,
the relatively small effects of estimation method on predicted transient flow
behavior are quite encouraging.
Deviations between actual and assumed forms of the k-S-P relations could be
responsible, at least in part, for differences in parameter values determined
from the equilibrium and transient data. Comparison of the observed
equilibrium data and Method 1 and 2 predictions suggest that the assumed
form of the S-P relations provides a reasonable representation of the actual
behavior. Since transient outflow results depend on k-S as well as S-P
relations, errors in the theoretical k-S model will affect Method 3 results, but
not those of Methods 1 or 2. Comparisons of Method 3 parameters with those
of Methods 1 and 2 indicate that the greatest deviations occur for the factor
«flow (note that etfao has no effect on the two phase oil-water experiments),
while n values deviate less markedly — especially for Soil 1. Prom equations
(2.31), (2.33) and (2.35), it is noted that relative permeabilities are independent
of «, 0ao and 0OW and depend only on parameter n. Thus, predicted two
phase relative permeability functions corresponding to parameters from
Methods 1-3 differ negligibly for Soil 1 and only slightly for Soil 2 (Figure
6.3).
\
Another possible source of discrepancy between Method 3 parameters and
those for Methods 1 and 2 is error in initial and boundary conditions and in
saturated conductivity values which are assumed to be known exactly for
transient flow inversions. In fact, these values are subject to some
uncertainty due to finite experimental error and, in the case of conductivity
values, to uncertainty due to performing measurements on samples other than
those used for transient displacement experiments. Errors in these input data
for the flow inversion analyses will bias estimated parameters from their true
values. To minimize such problems, care should be exercised in the design of
experiments to be utilized for parameter estimation by inversion of the
transient flow equations. Especially for soils with narrow pore size
distributions, boundary conditions should be chosen to achieve low final
wetting phase saturations. Use of air-water rather than oil-water experiments
would probably be advisable for routine calibration analyses, since air-water
displacement experiments are simpler to conduct and are less prone to error
in the regulation of initial and boundary conditions yet provide the same
information as oil-water experiments if interfacial tensions are employed to
determine scaling coefficients. Also, simulations for air injection can commonly
be performed without explicit consideration of nonwetting phase flow which
reduces the computational effort and avoids errors associated with uncertainty
in saturated oil conductivity.
57
-------
to - 10 - 10 • 10- 10
10
TIME >
Figure 6.2. Measured water outflow versus time for duplicate cores during
constant pressure NAPL injection (data points) and predictions
using Method 1 (solid lines), Method 2 (dashed) and Method 3
(dashed-dotted): (a) Soil 1 w/ p-cymene, and (b) Soil 2 w/ TCE.
_l
^
m
10
10
u
g 10-
0.
SI io-
Ui
a.
10
10
0.0 a 2 0.4 O. 8 0.0 1.0
WETTING FLUID SATURATION
to
0.0 0.2 a.4 a. a a. a i.o
WETTING FLUID SATURATION
Figure 6.3.
Model predicted two phase water and oil relative permeability
functions: (a) for Soil 1 using Method 1 (solid lines), Method 2
(dashed) or Method 3 (dashed-dotted), and (b) Soil 2 in which all
methods yielded indistinguishable curves.
58
-------
Pertinent to the question of properly specifying oil and water conductivities
for solution of the flow inversion problem, we note that saturated
conductivites are commonly assumed to be related to fluid viscosities as
(6.1)
where KaW and KgQ are water and oil saturated conductivities as employed in
Chapter 5 and T»W and TJO are water and oil viscosities. Measured KSO/K8W
ratios for Soils 1 and 2 with p-cymene and TOE, respectively, were 1.78 and
0.64. Using published pure fluid viscocity data (Dean, 1985; Weiss, 1980), we
calculate rjw/7»p to be 0.29 for p-cymene and 1.73 for TCB, which clearly do not
agree well with the measured conductivity ratios. Measured conductivity
ratios are subject to some uncertainty due to variability of replicate cores,
although calculated confidence limits on measured conductivity ratios exclude
the corresponding viscosity ratios at probability levels greater than 95X.
One possible explanation is that (6.1) simply does not adequately describe the
dynamic fluid behavior. This explanation seems improbable, at least for the
sandy porous media used in these studies. A more likely explanation is that
small differences in air entrapment during "saturated" conductivity
measurements lead to errors.
In any event, it was also observed that outflow simulations which employ
saturated oil conductivities predicted from measured water conductivities and
viscosity ratios via (6.7) did not adequately describe the experimental results.
Simulated transient outflow results using Method 1 and 2 retention function
parameters were found to deviate more severely from observed outflow data
when predicted rather than measured oil conductivities were used. Also,
parameter estimates obtained by inversion of the transient How data exhibited
greater deviations from Method I and 2 parameters when predicted rather
than measured oil conductivities were used.
CONCLUSIONS
Of the methods investigated here for calibrating the proposed multiphase
k-S-P model, all have yielded reasonably good descriptions of both static S-P
relations and of transient two phase displacement experiments, although the
equilibrium-based procedures (Methods 1 and 2) naturally describe S-P data
with greatest accurately, while the transient inversion procedure (Method 3)
predicts transient flow behavior more closely. The substantial agreement
between various calibration methods not only facilitates the selection of
calibration procedures to ease experimental effort involved, but provides
validation of the form of the proposed k-S-P model for the conditions involved
in the present study.
59
-------
CHAPTER 7
MODEL VALIDATION AND CALIBRATION FOR
THREE PHASE TRANSIENT PLOW
In this chapter, we discuss experimental procedures for simultaneously
measuring transient fluid saturations and pressures in three phase porous
media systems in the laboratory and compare measured liquid saturations and
pressures during a transient three phase flow experiment with numerical
simulations based on the proposed Jc-S-P model. Water and oil saturations are
determined with a dual energy gamma radiation apparatus and liquid pressures
are measured using hydrophobia and hydrophilic ceramic tensiometers. The
experimental regime imposed monotonically draining water and total liquid
saturation paths to avoid hysteretic effects in order to evaluate the k-S-P
model of Chapter 2 under conditions for which it should be valid.
MATERIALS AND METHODS
Fluid and porous medium properties
The porous medium employed in this study was an unconsolidated sandy
material comprised of 4.8, 94.7, 0.4 and 0.3 X fine gravel, sand, silt and clay
sizes, respectively. The average bulk density of the sand packed in the flow
column as measured with the dual energy gamma radiation apparatus was 1.60
S cm'3 which corresponds to a porosity of 0.40. The NAPL was a 1:9 mixture
DX volume of 1-iodoheptane and Soltrol-170, a hydrocarbon solvent produced
by Phillips Chemical Company consisting of a mixture of branched alkanes with
carbon numbers ranging from CIO to 016 and having very low solubility in
w«*ter. The aqueous phase was a 7X solution of NaBr by weight. Addition of
NaBr and 1-iodoheptane to water and Soltrol, respectively, was employed to
enhance differences in gamma attenuation coefficients (Oak and Bhrlich, 1985).
The aqueous phase and NAPL mixture densities were 1.05 and 0.83 g cm' ,
respectively. Air-water, air-oil and oil-water interfacial tensions were
measured by the drop method (Adameon, 1967) in quintuplicate using the
mixtures after they were in contact with the sandy porous medium for an
extended period (in excess of 90 h) yielding values of 73.8, 26.0 and 43.3 mN
N'1, respectively. The liquids were mixed with the porous medium before
interfacial tension measurements were conducted to obtain values reflective of
in the porous medium.
To calibrate the Jr-S-P model employed in the multiphase flow code, two phase
air-water S-P relations were measured in triplicate according to procedures
outlined in Chapter 3 on 5.26 cm diameter x 5.18 cm length cores of the sandy
Porous medium compacted to the same density measured in the transient flow
experiment. Details pertaining to calibration of the S-P model will be
discussed in. the section on numerical simulation. Saturated conductivities
*ere measured by a falling head method on triplicated cores having the same
60
-------
dimensions used for air-water S-P analyses. Replicate means of water and oil
saturated conductivites were 77.4 and 37,7 cm h"1, respectively, using fluid
pressure heads defined in' water-equivalent heights (Chapter 5).
Flow cell design and liquid pressure measurements
A plexiglas flow cell (Figure 7.1), 100 cm in height and 7.5 cm square in cross
section, was equipped with specially designed porous ceramic tensiometera to
enable independent measurement of NAPL and water phase pressures. Soltrol
does not interact with plexiglas in a manner to prohibit its use. Ceramic
tensiometera connected to pressure transducers (Figure 7.1) were placed at 10
cm intervals on both sides of the cell with water and NAPL phase pressures
being measured on opposite sides at 5 elevations. Drainage of liquid from the
column was accommodated by an outlet at the base of the cell which was large
enough not to restrict liquid outflow during drainage.
Measurements of liquid phase pressures were accomplished by establishing
continuous liquid phases between the porous medium and pressure transducers
located outside the flow cell. Continuous water and oil phases were obtained
by employing untreated and treated ceramic tensiometers, respectively. When
fluid wettabilities follow the order water>oil>air, a continuous water phase will
extend from pores in the soil to those in untreated ceramic plates provided
the oil-water capillary head does not exceed the oil-entry capillary head of
the ceramic. The latter is a function of the largest continuous pore size in
the ceramic and the interfacial tension between oil and water. A continuous
oil phase may be obtained by treating the ceramic to make it hydrophobia.
This was achieved with chlorotrimethylsilane which forms durable Si-O-Si
bonds with the ceramic surface, while attached trimethyl groups repell water
molecules in preference to oil and air (Ogden, 1985). Ceramic sections were
immersed in' chlorotrimethylsilane for approximately 2 hours and then rinsed
thoroughly with toluene followed by methanol. Caution needs to be exercised
when working with chlorotrimethylsilane because it may react violently with
water.
HYDROPHOBIC
CERAMIC INSERTS
OIL
PRESSURE
TRANSDUCERS —
'
X.
[ ^\
^^J
pw
ls^
rs^
^3
^^
^1
^
r\
1 ^s
V\l
.^
L/i
[^
^
^
X"
TN
|r«x
1 1 ^fc
II 1
•^
L
|l. 1
1 x
{t^l
1 ^
{•sl
nx<»»
%
| ^^
..
[^
iw^
&
^J
C3
x
— HYOROPHILIC
CERAMIC INSERTS
— WATER
PRESSURE
TRANSDUCERS
1
u— DRAINAGE OUTLET
Figure 7.1. Experimental three phase flow column design for simultaneous
measurement of oil and water pressures.
61
-------
Simultaneous determination of water and oil saturations
A dual-energy (Am-241 and Cs-137) gamma radiation system, as described by
Hopmans and Dane (1985, 1986), was used to simultaneously determine
volumetric , fluid contents of a three phase gas-NAPL-water porous medium
system. The general attenuation equation (Beer's Law) for a radioactive beam
passing through a porous medium containing three fluids, viz. water (w), NAPL
(n) and gas (g) is:
= I expx - ~ ~ (7.1)
in which I is the exiting radiation intensity expressed as a count rate (counts
per second, cps), !„ is the incident count rate (cps), Mi are mass attenuation
coefficients (cm2 g-1), ^ are maBB densities (g cm'3), ^ are volumetric fluid
contents (cm-* cm'3), and x is the path length through the porous medium
(cm), where , denotes soil (s), water (w), NAPL (n) or gas (g). Assuming the
mass density of soil (Pa), i.e., the bulk density of the porous medium, to be
constant and the contribution of the gas term to be negligible, (7.1) can be
rewritten as
= I
where l' = I
Following Nofziger and Swartzendruber (1974) and differentiating between
energy* sources Am-241 (subscript a) and Cs-137 (subscript c), (7.2) yields
ha(tr- U^e^) (7.3a)
oc exp-fccO^- U^e^ (7.3b)
where Uij = ^j/jjjx in which subscripts i and j indicate the material and
energy source, respectively, and LJJ r LJJ exp(-p8JPsx).
Solving (8) for Bw and 9n yields
(7.4a)
(7.4b)
«rhere k = UwcUna - U^U^. To solve (7.4) U^, Unc, Uwa, U^c, ^ and r-
Twer'e dTtlrmi d 1 ana, sa wc cPnc «n
use were determined in preliminary experiments with a plexiglas cell of similar
62
-------
length and width as the flow cell which was divided 'into 4 vertical
compartments of equal thickness transverse to the gamma radiation beam.
Radiation counts were recorded for the empty cell and after each compartment
was filled with the material of concern (i.e., water, NAPL or soil) to obtain 5
measurements of exiting radiation intensity as a function of path length. This
procedure was replicated four times yielding a total of 20 data pairs for each
material and radiation source from which HjPij were determined by regression
to the linearized attenuation equation (see equation 7.3)
ln(Ij) = ln(I0j) - Pijpijx (7.5)
where x is the path length through the material, which is a multiple of the
thickness of a single compartment, and Ioj is the incident count rate for
energy source j. To determine the mass attenuation coefficients of dry soil
0%j) under conditions in which the mass density of porous medium (fa)
inevitably varies somewhat between compartments, ^\ are determined by
regression to .
(7.6)
where volumetric water content, 9^, and bulk density, p^, in each compartment
were determined gravimetrically; x is the path length through the porous
medium; and j again refers to the energy source. Values for jiwaPwa* PwcPwc.
VnaPna, VncPno Psa» and Psc so obtained were 0.337596, 0.088411, 0.735462.
0.071888 cm'1 and 0.242413 and 0.076085 cmj g-1, respectively.
Gamma attenuation measurements during the transient flow experiment were
taken at the same elevations as pressure measurements. Path lengths
corresponding to the 5 measurement locations along the column needed to
calculate Uj/s were determined by first taking count rates at the selected
measurement positions along the empty flow column to yield values of 1^ and
loc, followed by counts on the flow column filled with water of known /^p^
and A»wa/>wa providing values of IQ and 1^. Path lengths corresponding to the
measurement locations, x, for each energy source could then be calculated as
x =
(7.7)
During the calibration procedure all counting times were 5 min in duration.
To determine path lengths, 4 cycles were completed on the empty column and
4 on the column filled with water, resulting in a total counting time of 20 min
at each position for t,he empty and water filled column, respectively.
Path lengths calculated for Am-241 and Cs-137 using (7.7) differed by a
maximum of 1.17X. Volumetric water and oil contents for the transient flow
experiment calculated using path lengths determined from Am-241 and Cs-137
sources differed by less than IX. Thus, this discrepancy in path lengths was
considered to be acceptable and mean values were employed in practice.
63
-------
All count rates were corrected for resolving time as outlined by Fritton (1969).
Radiation emitted by the Cs-137 source, but detected in the Am-241 window,
was accounted for by the method of Nofziger and Swartzendruber (1974). All
count rates in the above equations refer, therefore, to corrected count rates.
Upon packing and saturating the porous medium with water, 4 cycles of 5 min
counting periods were completed to determine the initial water content and
bulk density at the five locations according to equations (15) and (16) of
Hopmans and Dane (1986). Bulk density values were needed for loa and IQC.
With all parameters known in (7.3), 9W and 6n could be determined during
subsequent transient conditions. Volumetric fluid contents were converted
into fluid saturations from the measured bulk density at each of the 5
measurement positions (i.e., saturation = volumetric fluid content/porosity).
During the transient flow experiment count rates~ were obtained over 1 min
periods.
Imposed experimental initial and boundary conditions
The sandy porous medium was packed on top of a 4 cm thick coarse gravel
layer to a height of 56.5 cm from the bottom of the flow column. The gravel
layer ensured that drainage of liquid would not be impeded and that liquid
flow would be one dimensional throughout the column. The packed column was
gently flushed with CO3 before imbibing water slowly through the bottom of
the cell to minimize gas entrapment. In addition, the water saturated column
was left undisturbed for approximately 12 hours to enable entrapped gas to
diffuse out of the porous medium. Locations of the ceramic tensiometers (i.e.,
measurement locations) were at 5.2, 15.2, 25.2, 35.2 and 45.2 cm below the soil
surface.
Before drainage of liquid from the column commenced, 200 mL of the NAPL
mixture was ponded on top of the water saturated soil. The initial depth of
oil on the soil surface was 3.56 cm. Drainage of liquid from the column at t =
0 was permitted after the outflow elevation was positioned at 25.2 cm below
the soil surface. Measurements of liquid saturations and pressures were
initiated shortly after drainage commenced. The outflow point was reset to 1.5
cm above the lower soil boundary at t = 41.5 min to allow additional drainage.
This stepwise drainage procedure was used to prevent excessively large flow
rates during the time period employed for gamma counts. The imposed initial
and boundary conditions ensured monotonically draining water and total liquid
saturation paths.
Numerical simulation
A one dimensional version of the finite element multiphase flow code described
in Chapter 5, which solves the coupled partial differential equations governing
flow in a three fluid phase porous medium system with assumed constant air
phase pressure, was used to simulate the experiment. The k-S-P model of
Chapter 2 was employed to describe the constitutive relations of the porous
medium. Parameters in the model consist of «, n, Sjg, ftaw, 0ao» anc' Pow> as
defined by (2.10) - (2.14), and additionally the porosity, +, and water and oil
64
-------
saturated conductivites, KaW and KBO, respectively. Values of +, KgW and Kg0
were determined as previously discussed. Due to the coarse grain size of the
porous medium, we assumed S^fO. Since contiguous air and water phases are
absent in the experiment considered here, the factor atftaw is not required for
the simulations. The remaining parameters needed to calibrate the model are
n, «0ao and <*Pow (note parameters a and 0H occur only as a product in the
constitutive equations).
Two methods were utilized to evaluate the latter parameters. In Method 1,
independently measured air-water S-P data were used to fit parameters n and
a ( the air-water system was arbitrarily taken as the scaling reference, i.e.,
Paw - 1). and scaling factors ftao and ftaw were obtained from ratios of
mterfacial tensions as described in Chapter 3. In Method 2, the parameters n,
a/?ao and <*0OW were obtained by fitting water saturation vs. oil-water capillary
head data and total liquid saturation vs. air-oil capillary head data actually
measured during the transient flow experiment to (4.L) and (4.2). Fluid
saturations were determined with the dual gamma apparatus and capillary
heads were obtained from the oil and water pressure measurements; air
pressure was assumed constant at that of the ambient atmosphere. Values of
the parameters employed in the simulations for both calibration methods are
listed in Table 7.1.
Table 7.1. Parameters of the three phase Jr-S-P model employed in numerical
simulations as obtained by two calibration methods.
Parameters
°Paw (cm-l )
<*0ao (cm"1 )
<*ftoff (cm"1 )
n
Sm
*
1(319 ( cm h~ * )
KSO ( cm h~ * )
Method 1
0.054
0.154
0.092
3.248
0.0
0.40
77.4
37.7
Method 2
0.104
0.091
2.729
0.0
0.40
77.4
37.7
65
-------
z = 3.56(Poil/Pwater)
L.5 min by
A*(52.5,t) = 27.3 en
q0(52.5,t) = Q
for tMl.5 min by
A*(52.5,t) = 1.5 on
q0(52.5,t) = 0
W.hr S^L".™±± l?*'^- ** «».«PP«1 boundary «, = 0), . zero
for t > o
A0(0,t) = (po/p^tQtot - Q(t)] for Q(t) <
Q(t) , Qtot
• ( ,o(0
,t) dt
Qtot = 3.56 cm is the initial height of oil above the soil surface at t = 0.
66
-------
RESULTS AND DISCUSSION
Experimental results
Measured water, oil and total liquid saturations using the dual energy gamma
radiation apparatus for the five measurement locations are shown in Figures
7.2-7.6. Also shown are water and oil pressure heads expressed in centimeters
of water height. As drainage of water from the bottom of the flow column was
initiated, water saturation decreased at a very rapid rate near the upper
surface (Figure 7.2) which was mirrored by an increase in oil saturation. All
of the applied oil imbibed into the porous medium within 8 min.
With the outflow level positioned at z = 25.2 cm insufficient volume of water
drained from the column to allow air infiltration at z = 5.2 cm (Figure 7.2),
i.e., measured total liquid saturations remain at unity. A rapid approach to
equilibrium is evidenced by the small slope of saturation with respect to time
following the initial rapid changes in liquid saturation. Recorded cumulative
water outflow volume from the column at t = 36.5 min was only 210 cm3,
indicating only 10 cm3 of air had imbibed in the porous medium.
At t = 41.5 min the water table position was lowered again to its final position
at z = 51.0 cm. As soon as the water table was lowered, there were rapid
decreases in oil and water saturations (e.g., total liquid saturation) followed
by very gradual decreases at later times. Shortly after lowering the water
table, air was present at z = 5.2 cm. Changes in measured water and oil
pressure heads at z = 5.2 cm also reflect lowering of the water table.
Only a small volume of oil penetrated to 15.2 cm below the soil surface (Figure
7,3) when the water table was positioned at z = 25.2 cm. After the water table
was lowered to z s 51.0 cm, there were rapid decreases in water saturation
and pressure head with corresponding increases in oil saturation. Air did not
penetrate to z = 15.2 cm until 85 min after commencement of drainage. At z =
25.2 cm (Figure 7.4), the elevation of the first water table position, oil was not
detected until after the water table was lowered to z = 51.0 cm. Since
Soltrol-170 is less dense than water, oil should not penetrate to z = 25.2 cm
when the water table is positioned at the same elevation. There was an
anomaly in measured liquid saturations at about 50 - 60 min when it appeared
that water saturation decreased without a corresponding increase in oil
saturation. At t = 65 min, however, oil became evident in the porous medium
at this depth. The anomaly may be a consequence of the rapidly changing
liquid saturations relative to the duration of gamma counts. At later times
water and oil saturations decreased gradually as equilibrium conditions were
approached corresponding to the final water table elevation.
At z = 35.2 cm (Figure 7.5), which was 15.8 cm above the final water table
position, the fluid system in the sandy porous medium changed from a single
fluid phase system (water) at early times to a two fluid phase system (water
and oil) at later times. When oil appeared at this depth at t « 100 min, there
was a dramatic change in the oil pressure head. Treated (hydrophobia)
ceramics are wet by organic liquids and air in preference to water, therefore,
a small volume of air may occur between the ceramic surface and the water
67
-------
phase when oil is not present. Pressure readings from hydrophobic
tensiometers in the absence of oil in the porous medium will reflect pressures
of the water phase transferred by the gas film and modified due to air-water
capillary effects in the porous medium and air-oil capillary effects at the
tensiometer surface. When a continuous oil phase is present in the porous
medium, oil will wet the hydrophobic ceramic surface in preference to air and
oil pressure head measurements will reflect the pressure of the oil phase at
that elevation. The conspicuous change in oil pressure head between t = 92
to 100 min in Figure 7.5 corresponded to an oil saturation change from 0.01 to
0.05 cm3 cm~s. However, since at z = 35.2 cm the porous medium was liquid
saturated, i.e., Sj = 1, the oil pressure heads should be positive or equal to
zero when oil has imbibed to this elevation. An explanation for the anomolouH
oil pressure heads when oil is present at z = 35.2 cm may be experimental
error in calibrating the oil pressure transducer at his location. At the fifth
measurement position, z =- 45.2 cm (Figure 7.6), the porous medium remained
essentially water saturated throughout the duration of the experiment and
measured oil pressure heads were zero when only a small volume of oil was
present (So = 0.01).
Note that oil pressure heads at z = 15.2 cm (Figure 7.3) changed very little
(0.5 cm) as the total liquid saturation decreased from 0.81 to 0.55 whereas at z
- 5.2 cm (Figure 7.2) the oil pressure heads changed from -10.5 to -12.4 as
tne total liquid saturation decreased from 0.84 to 0.54 and at z = 25.2 cm
(Figure 7.4) the oil pressure heads changed from -8.3 to -10.2 as the total
liquid saturation decreased from 0.67 to 0.55. This suggests that either the
oil pressure transducer at z = 15.2 cm was not functioning properly at later
times (^ 64 min) or that the porous medium was not packed tight enough
against the chemically treated ceramic insert and when the total liquid
saturation decreased below approximately 0.8 contact between the treated
ceramic surface and the liquids contained in the porous medium was severed
as Uquid drainaged from the relatively large flow channel between the porous
medium and treated ceramic tensiometer. Consequently, we do not expect good
agreement between simulated and measured oil pressure heads at this
elevation.
Model validation
Calit>ration °t tne three phase Jr-S-P model was achieved either by use of two
phase air-water S-P data in conjunction with interfacial tension measurments
(Method 1) or by directly fitting parameters to three phase S-P data obtained
from the transient flow experiment (Method 2). To evaluate the accuracy of
the*»e procedures, three phase S,y(J)OMf) and St(hao) data, obtained from oil and
water saturation and pressure measurements during the transient flow
experiment, are compared in Figure 7.7 with predictions by the two methods.
Differences in S^h^ predictions for the two methods are minor and both
^gree well with the observed data. Scatter in the data is fairly substantial,
reflecting composite effects of measurment errors and heterogeneity in packing
from depth to depth.
68
-------
TOTAL LIQUID
O. O
75 ISO 225
TIME CmirO
3QO
O
150 225
TIME CmirO
3QQ
Figure 7.2.
Results of transient flow experiment at 5.2 cm depth showing (a)
measured total liquid, water and oil saturations and (b) measured
water and oil pressure heads as fuctiona of time.
69
-------
1. 0
TOTAL LIQUID
O. O
170 255
TIME CmirO
340
D
<
UJ
x
LU
o:
D
in
(/)
u
o:
o_
-3O
-20 -
-10 -
85 17O 255
TIME (man)
34O
Figure 7.3. Results of transient flow experiment at 15.2 cm depth showing (a)
measured total liquid, water and oil saturations and (b) measured
water and oil pressure heads as fuctions of time.
70
-------
TOTAL LIQUID
MM4
7O 140 21O
TIME CmirO
2SO
-28
U
I
UJ
tr
D
0)
0)
LJ
(T
OL
-21 -
-14 -
-7 -
O
OIL
O 70 14O 21O 28O
TIME CmirO
Figure 7.4. Results of transient flow experiment at 25.2 cm depth showing (a)
measured total liquid, water and oil saturations and (b) measured
water and oil pressure heads as fuctions of time.
71
-------
TOTAL LIQUID
O. O
85 17O 255
TIME CmirO
34O
Q -30
<
LJ
LJ
o:
w
ui
Li)
o:
£L
B
65
170
255
340
TIME CmirO
Figure 7.5. Results of transient flow experiment at 35.2 cm depth showing (a)
measured total liquid, water and oil saturations and (b) measured
water and oil pressure heads as fuctions of time.
72
-------
1. O
_ O. 8 -
O
£ O. 6 -
D O. 4 -
I-
(/) O. 2 -
O. O
T \
WATER
TOTAL LIQUID
OIL
85 170 255
TIME CmirO
340
-10
LJ
I
LJ
CC
D
CO
CO
u
IT
Q_
85 170 255
TIME Cmin)
340
Figure 7,6. Results of transient flow experiment at 45.2 cm depth showing (a)
measured total liquid, water and oil saturations and (b) measured
water and oil pressure heads as fuctions of time.
73
-------
Discrepancies between measured and predicted St(Jiao) relations are greater
than for Sw(J]Ow) results. The scaling factor estimated from interfacial tension
data (i.e., Method 1) is significantly greater than that inferred directly from
the three phase data (Method 2). Scatter about the Method 2-fitted curve is
more severe than for the Sw(haHr) data. Since deviations between observed
and fitted St(Jiao) and S^(how) for the same observation depths were not
significantly correlated, porous medium heterogeneity seems an unlikely cause
of the large scatter for St(hao) data. Since errors in oil pressures would
effect Sw(hOMr) as well as St(hao) results, this source of error is also
improbable. Inspection of the St(hao) data points which lie above the Method
2-fitted curve indicates that these values largely correspond to measurements
at the 35.2 cm depth where there was only slight change in pressure readings
(i.e, 0.5 cm) as the total liquid saturation decreased from 0.81 to 0.55. If the
latter points are ignored, we note that agreement between Method 1 and
Method 2 St(hao) curves would be substantially improved. Deviations in the
assumption of atmospheric gas phase pressure could also cause errors in the
"observed" capillary head values.
Fluid saturations and liquid pressures were simulated with the multiphase flow
code for the period in which measurements were conducted using Method 1
end Method 2 Jr-S-P model parameters. Lowering the water table from z = 25.2
to 51.0 cm at t = 41.5 min resulted in some numerical difficulties evidenced by
oscillations in computed oil pressures. In order to obtain more stable
solutions, we modified the computational lower boundary conditions so that
instead of a single large drop in the water table position, the water table was
lowered incrementally over a longer time to a slightly higher final elevation.
These modifications dampened the oscillations and are not expected to have
major effects on predicted fluid saturations and pressures.
J*igures 7,8-7.11 show comparisons between predicted and measured water and
oil saturations and pressure heads as a function of time for elevations z = 5.2,
15.2, 25.2, and 35.2 cm, respectively. At z = 45.2 cm, experimental data and
numerically simulated saturations (data not shown) both reflected that at this
elevation the porous medium was water saturated during the duration of the
experiment. Solid lines in these figures represent predicted saturations and
pressure heads obtained when the model was calibrated by Method 1 and
dashed lines when the model was calibrated by Method 2. Predicted water
saturations correspond favorably to the measured water saturations, except at
z - 25.2 cm. Some of this error may be attributable to lowering the water
table more slowly in the numerical simulation resulting in a greater water
Content. However, since the discrepancy is substantial only at this depth, a
tuore likely explanation is that local heterogeneity in soil packing at this depth
has resulted in deviations of average soil properties from those inferred from
tjie S-P measurements used to calibrate the model. There is good agreement
between predicted and measured oil saturations at all elevations. The
numerical simulation predicts three phase air-oil-water systems at locations z =
$.2, 15.2 and 25.2 cm, a two phase oil-water system at z s 35.2 cm and a
predominantly water saturated system at z = 45.2 cm. This sequence matched
that experimentally observed. The average difference between predicted and
Measured liquid saturations for all times was less than 10% of the pore volume
if discrepancies between predicted and measured water saturations at z = 25.2
are excluded.
74
-------
0.2 0.4 O. 6
WATER SATURATION
l. 0
Figure 7.7.
0.0 0.2 0. 4 O. 6 0. S 3.0
TOTAL LIQUID SATURATION
Comparison of observed saturation-pressure data in transient
three phase flow experiment (points) and predictions using
Method 1 parameters (solid lines) and Method 2 (dashed lines):
(a) water saturation vs. oil-water capillary head, and (b) total
liquid saturation vs. air-oil capillary head.
75
-------
75 ISO 225
TIME
3OQ
1.0
o. o
75 ISO 225
TIME
3OO
ISO
TIME GnirO
225
30Q
Figure 7.8. Comparison of measured water saturations and heads (solid
symbols) and oil saturations and heads (open symbols) with
numerically simulated values at 5.2 cm depth using model
parameters from Method 1 (solid lines) and Method 2 (dashed
lines): (a) measured versus predicted water saturations, (b)
measured versus predicted oil saturations* and (c) measured
versus predicted water and oil pressure heads.
76
-------
1. O
o. o
170 255
TIME CmirO
340
-3D
U
I
u
o:
D
If)
(/)
LJ
CE
Q.
-2D -
-1O -
B
OIL
i
85
17O 255
34O
TIME CmirO
Figure 7.9. Comparison of (a) measured versus predicted water and oil
saturations and (b) measured versus predicted water and oil
pressure heads at 15.2 cm depth. Legend as Figure 7.8.
77
-------
1.0
0. 0
7O 140 210
TIME CmirO
28O
-3D
LJ
LJ
(T
m
V)
LJ
•* *
7O 14O 21O
TIME CmirO
28O
Figure 7.10. Comparison of (a) measured versus predicted water and oil
saturations and (b) measured versus predicted water and oil
pressure heads at 25.2 cm depth. Legend as Figure 7.8.
78
-------
o. o
Q5 170 255
TIME CmirO
34O
Q -OU -
UJ
I
-20 -
LU
o:
(/) -10 -
(/)
LU
o:
n o -
B ^««.
•
WATER
****** /
<**"" ^^o* « o oo « o
l^":\^^—'^ ^-OIL
1_ ^ —
0 85 170 255 34.
TIME CmirO
Figure 7.11. Comparison of (a) measured versus predicted water and oil
saturations and (b) measured versus predicted water and oil
pressure heads at 35.2 cm depth. Legend as Figure 7.8.
79
-------
Predicted and measured water pressures generally corresponded closely, while
°il pressures agreed less well. Deviations between predicted and measured
Pressure heads appeared to increase at lower elevations. Method 2
Parameters, which were obtained from measured Sv(how) and St(hao) relations
including oil pressure head data at z = 15.2 cm, predicted water and oil
Preaaure heads as well as or more accurately than Method 1 parameters at all
elevations. This was expected because Method 2 parameters were calibrated
from S-P relations measured during the flow experiment. In practical
circumstances, Method 1 parameters would be much more accessible and,
considering the simplicity of the calibration method, the accuracy of oil and
*ater saturations and pressures predicted with Method 1 parameters is very
encouraging.
SUMMARY AND CONCLUSIONS
A three phase air-oil-water flow experiment was designed and conducted
involving simultaneous measurements of liquid saturations and pressures
d"ring monotonic water and total liquid drainage. Oil and water saturations
w«re measured with a dual gamma radiation apparatus after doping the liquid
Phases to enhance separation of gamma attenuation coefficients. Oil and water
Pressures were measured with hydrophobic (treated) and hydrophilic
Untreated) ceramic tensiometers, respectively.
e water saturation vs. oil-water capillary head data during the three
How experiment agreed well with predictions using the Pr°P°Jfd k'^
. of Chapter 2 calibrated from static air-water S-P relations and
^terfacial tension data (Method 1). Measured total liquid saturation vs. air-oil
c*PUlary head data deviated more severely from the air-water data predictions
* exhibited marked scatter, which was inferred to reflect, in part, error, , m
Perimental measurements or perhaps due to deviations from the -^J^*™
Constant gas phase pressure. Model parameters were also evaluated by
rectly fitting to the observed three phase S-P data (Method Z).
g<* agreement was found between measured liquid saturation
I?* and space during the transient flow experiment and those
£« finite Element multiphase flow code using the P^^
2'ibrated using either Method 1 or Method 2. Measured liq
£re«* more closely with simulations employing Method 2 .
, considering the practical advantages of using Method 1 1 tor mode
the accuracy of the simulations based on these estimates is very
the experin,ental regime considered m
A
of the
no.b.o
to evaluate field scale multiphase flow behavior.
80
-------
CHAPTER 8
MODEL FOR HYSTERETIC
PERMEABILITY-SATURATION-PRESSURE RELATIONS
Hysteresis in k-S-P relations is often ignored in numerical analyses of flow in
three phase air-nonaqueous liquid-water-pouous media systems (e.g. Peery and
Herron, 1969; Casulli and Greenspan, 1982; Abriola and Finder, 1985; Faust,
1985; Kuppusamy et al., 1987), despite the fact that saturation paths are often
markedly nonmonotonic. Disregarding hysteresis in numerical simulations of
flow in two phase systems has been shown to result in significant errors in
predicted fluid distributions (Hoa et al., 1977; Gillham et al., 1976; Kool and
Parker, 1987; Kaluarachchi and Parker, 1987). Whereas two directions of
saturation change are feasible in two phase systems ~ i.e., wetting phase
drainage or imbibition — in three phase systems, 12 path directions are
possible, i.e., IID, IDI, IDD, DII, DID, DDI, GDI, CID, ICD, DCI, IDC, and DIC
where I, D and C designate imbibition, drainage or constant saturation,
respectively, of water, NAPL, and gas (Saraf et al., 1982). As a result,
disregarding hysteresis in simulations of flow in three phase systems subject
to complex boundary conditions is likely to produce greater errors than is
observed for two phase systems. Killough (1976) has demonstrated major
differences in hyateretic and nonhysteretic simulations of hydrocarbon
recovery from three phase petroleum reservoirs, even for simple boundary
conditions associated with continuous oil and gas extraction or waterflooding.
Furthermore, in groundwater contamination problems, the phenomenon of
nonwetting phase entrapment associated with certain saturation paths (Land,
1968; Wilson and Conrad, 1984) may be expected to have extremely important
effects on predictions of transport of organic components which can partially
dissolve into the water phase or vaporize into the gas phase. This will be
particularly so when occluded immobile NAPL acts as a source of contaminants
which are gradually released into mobile aqueous or gaseous phases.
Prediction of the source distribution, and hence of resulting aqueous and
gaseous phase plumes, will not be possible without consideration of the
process of nonwetting fluid entrapment.
Hysteresis in S-P relations of two phase systems has been well documented
and various theoretically and empirically baaed models have been developed to
describe two phase hysteretic S-P functionals from primary drainage and
imbibition path measurments (e.g., Mualem, 1974; Scott et al., 1983). A
procedure for incorporating air entrapment in the description of hysteretic
air-water saturation-pressure relations has been given by Kool and Parker
(1987). Killough (1976) presented a method for describing hysteretic
saturation-capillary pressure relations in three phase systems, but did not
consider effects of fluid entrapment.
Snell (1962) seems to have been the first to consider possible effects of
saturation history on three phase relative permeability-saturation
relationships. For the highly nonmonotonic saturation histories imposed by
Snell, relative permeabilities of all fluids were observed to depend markedly
81
-------
°n saturations of all phases. Oil relative permeabilities were found to be
consistently lower when oil and water were simultaneously draining than for
other saturation histories. Extensive three phase relative permeability
fceasurements by Saraf et al. (1982) for various saturation paths indicated
water relative permeability to be insensitive to saturation history, while gas
and oil relative permeabilities exhibited marked path dependence. Oil relative
Permeabilities were dependent on saturations of all three phases and were
8"ghtly hysteretic with respect to gaa saturation path.
UG to the considerable difficulties in performing and interpreting three phase
Relative permeability measurements, considerable effort has been given to the
evelopment of various predictive models. Based on experimental results of
^everett and Lewis (1941) and others, it is often assumed that three phase
water and gas relative permeabilities can be satisfactorily estimated directly
rom two phase system measurements. Various theoretical and semi-empirical
methods have been devised to estimate oil relative permeabilities in three
?JLase systems from two phase air-oil and oil-water measurements (Stone, 1970,
J3'3; Aziz and Settari, 1979; Payers and Ma thews, 1984). Hysteretic effects are
°t considered explicitly in these models, although they may be accomodated
by using appropriate hysteretic two phase permeability data as
by Killough (1976).
*ernative methods for evaluating three phase relative permeabilities from two
80 C'ata nave been proposed which are based on the estimation of flow
c^str*Dut")n8 from saturation-capillary pressure relations. Corey et al.
derived an expression for oil relative permeabilities in three phase
under conditions in which water and total liquid saturations are both
. The method, based on Burdine's (1953) pore distribution model and
power function for the saturation-pressure function, yielded generally
at ** a£reement with drainage permeabilities for Berea sandstone cores, except
low total liquid saturations. Other formulations for two and three phase
"^»nage relative permeabilites, based upon more general parametric models for
**tu ration-pressure relations, have been presented by Nakornthap and Evans
u"6), Corey (1986), and ourselves (see Chapter 2).
°difications of the Corey-Burdine approach have been proposed by Naar and
(1961) and Land (1968) to estimate three phase imbibition relative
(i.e., water and total liquid saturations both increasing) on the
that hysteresis in permeability relations is due largely to
^etting fluid entrapment, which is estimated via various empirical formulae.
°nwetting fluid entrapment affects wetting fluid conduction by obstructing
of ine p0re space which would otherwise be occupied by wetting fluid
by displacing some of the wetting fluid into larger pores. It also affects
ting fluid conduction by reducing the volume of continuous nonwetting
in the larger pores.
models for nonhysteretic S-P relations combined with predictive
k~S relations nave been developed for the description of two phase
8 ^van G«nuchten, 1985) as well as three phase systems (Corey et al.,
j Chapter 2 of this report). The combination of such models with
Co!^et»cally or empirically based hysteresis schemes can provide general yet
ncise descriptions of k-S-P relations for arbitrary saturation paths,
82
-------
resulting in major reductions in data requirements for model calibration, as
demonstrated by Kool and Parker (1987) for flow in air-water systems
described by Richard's equation.
Our purpose in this chapter is to develop a general method to describe two or
three phase Jc-S-P relations for arbitrary saturation paths, including effects
of nonwetting phase entrapment, based on extensions of the nonhysteretic
three phase Jc-S-P model of Chapter 2. The result will be a parametrically
parsimonious model which can be calibrated from measurements performed on
two phase systems and which can be easily implemented in numerical
multiphase flow codes.
BASIC NOTATION
For notational consistency with Chapter 2, we utilize water height equivalent
air-water, air-oil and oil-water capillary heads, respectively, defined by
(Pa ~
(Pa ~
(2.4)
(2.5)
(2-6)
where Pw, Po, and Pa are water, "oil" (i.e., NAPL) and air phase pressures,
respectively, g is the scalar gravitational acceleration, and pw is the density
of water.
It will also simplify subsequent analyses to introduce effective saturations of
water, oil, air and total liquid in three phase systems given respectively by
Sw -
1 - S
l -
Sa =
+ S -
1 - 5,
I -
(8. Id)
where Sw," SQ, and Sa are actual water, oil and air saturations, respectively,
St (zSw+So) is total, liquid saturation, and Sm is a minimum or irreducable
wetting phase saturation at which capillary heads appear to extrapolate to
infinity. We assume that wettability follows the order: water * oil *> air. It
may also be noted that Si = l-Sa =
In order to predict k-S-P relations for arbitrary saturation paths, it will be
necessary to describe the manner in which fluids are distributed within the
pore space. This will necessitate distinguishing between a continuous oil
83
-------
Phase which is free to move correctively ("free" oili «n^ *•
°r ganglia of oil occluded within thTwIJn «K ' ! d dlslcontln«°u8 islands
conditions ("trapped" oil) Under *?£ pha?4.aBd immobile under ambient
emulsified by the water phase in SSS" COndltlon8' «* °1°°* -nay become
Phases" trapped air occluded within continuous water or oil
o/1 0* (8.2a)
= Saf + Sat (8.2b)
= Satv+ Sat0 C8.2c)
P*
total Uquid ^turationa for two or
fluds fluld 8atu~tions within the
(8>3a)
Sr = ^ + S0 + Sat (8.3b)
.,
S«r®J>re8enl8 a saturation equivalent to the effective saturation of
"ld! occluded within the water phase and ST represents a
^t.10 the effective "Cation of water, free and trapped ofl
within the water and oil phases.
WlU *enerally not «»«" •* *U locations or times in the porous
gro,undwater contamination scenarios, it is necessary to consider
wa , ° phase air-water systems as well as that of three phase
modS 8^,f ?v j;urthermore» calibration of the proposed three phase
ai? ,*U1 UtlllZe data which ""'y be obtained from measurements of two
"•-water, air-oil and oil-water systems. In order to clearly indicate
84
-------
which fluid system is under consideration, we will explicitly denote saturations
in two phase systems by the superscripts aw, ow, or ao to represent the
corresponding fluid pair. Absence of superscripts will be taken to imply a
three phase system. Effective wetting and nonwetting phase saturations are
defined for two phase air-water, oil-water and air-oil systems, respectively, by
Sao =
1 - SB
w0" - S
1 - 8
S0ao _
1 - S
(8.4a)
(8.4c)
(8.4e)
Saaw
1 - Sn
So0**
1 - S,,,
Saao
1 - S
(8.4b)
(8.4d)
(8.4f)
The foregoing assumes SQ, to be independent of the fluid system. This may
not be strictly true. However, since it is three phase water-wet systems
which are of concern, such difficulties may be largely circumvented by
performing air-oil experiments on media initially near the irreducable water
saturation.
We introduce free and trapped phase effective saturations for two phase
systems in a manner analogous to the three phase case by
(8.5a)
(8.5b)
(8.5c)
and apparent wetting phase saturations by
= Sw+ Sa
(8.6a)
Saat°
(8.6b)
(8.6c)
where the nototion follows that for three phase variables.
85
-------
DESCRIPTION OP SCALED S-P RELATIONS
Definition of Scaling Relations
Wetting fluid saturations in two phase systems are normally regarded as
unique functions of capillary head in the absence of nonwetting fluid
entrapment or hysteretic effects. We propose accomodating effects of fluid
entrapment by employing apparent rather than actual or effective saturations
in a manner similar to that used for air-water systems by Payer and Hillel
(1986). Hysteresis will be accomodated by treating apparent
saturation-capillary head relations as history dependent functional.
Furthermore, we assume that two phase S-P functional for arbitrary fluid
pairs in a given rigid porous medium may be scaled to eliminate fluid
dependent effects via
(8.7a)
(8.7b)
(8.7c)
where 0aw, 0ao and ftow are fluid dependent scaling factors as introduced in
Chapter 2 and which may be estimated from fluid interfacial tensions as
demonstrated in Chapters 3, 4, 6 and 7.
It is commonly assumed that in three phase systems Sw will be a function only
of how and St will be a function of /^, for monotonic saturation paths in the
absence of fluid entrapment (Leverett, 1941; Leverett and Lewis, 1941; Corey
et el., 1956; Aziz and Settari, 1979). We propose generalizing these
assumptions by using Sw and ST in lieu of S^ and St to eliminate fluid
entrapment effects and by treating Sw(hoW) and Sr(Jfe0) as independent
hysteretic functional. Furthermore, we assume that both hysteretic
functional may be represented by a single scaled functional defined as for
two phase systems by
(8.8a)
(8.8b)
where 0^ and 0OW are fluid dependent scaling factors assumed identical to
those for the corresponding two phase systems.
Equations (8.8a) and (8.8b) pertain so long as a free oil phase exists, that is
if ST * Syf. If SW = ST and Sot = 0, then a simple two phase air-water
system is indicated and (8.7a) pertains. If Sw = ST and Sbt * 0, then a more
subtle case exists in which we assume air-water interfaces control water
86
-------
saturation, although the scaling factor may differ from the simple two phase
case due to the presence of soluble oil components in the water phase. For
the latter case we designate the scaling relation as
where few' is a scaling factor for the system with trapped oil but no free oil,
the value of which should differ from 0aw in relation to the ratio of air-water
interfacial tensions for organic-free and organic-contaminated water.
It should be noted that the use of apparent saturations leads to scaled
hysteretic functional which are single valued at zero capillary head and at
infinite capillary head. This behavior is in contrast to that of effective
wetting fluid saturation-capillary head relations, which will in general not
yield closure at zero capillary head due to fluid entrapment during imbibition.
In the remainder of this section we will outline procedures for describing the
hysteretic S*(h*) functional. Analyses of fluid entrapment, necessary to
transform between actual fluid saturations and apparent or scaled saturations,
will be taken up in a subsequent section.
Parametric Description of Main Branches
Various empirical forms have been proposed for the description of S-P
relations which may be adopted to parameterize the description of S*(h ).
Here, we employ van Genuchten's (1980) function and describe the main
imbibition and drainage branches, respectively, of the scaled retention
function for h* » 0 by
(8.9)
= (l + (4* #)a}~* (8.10)
where 4«, ^a and n are curve shape parameters and m=l-l/n. For h* ^ 0,
>S*(Ji*) = dS*(h*) - 1. The assumption of constant n for imbibition and
drainage has been shown by Kool and Parker (1987) to be a reasonable
approximation for air-water systems and will be invoked here for simplicity.
Preforatory superscripts i and d are used to denote main imbibition and
drainage branches, respectively. Note that when S^(^ow) or Syfaw(haw) are
under consideration, imbibition and drainage refer to increasing and
decreasing apparent water saturations, respectively, whereas increasing and
decreasing apparent total liquid saturation are signified for Sq>(hao). We will
use the terms imbibition (or wetting) and drainage (or drying) in this sense
throughout this report.
87
-------
Scanning Curve Interpolation
To describe S*(h*) scanning curves, we employ an empirical procedure which
rescales main branch functions to pass through appropriate reversal points.
The procedure is similar to the methods of Scott et al. (1983) and Kool and
Parker (1987), except that in the present method closure of scanning loops is
enforced. Consider an arbitrary imbibition scanning curve beginning at point
(*S$jf*h§d, corresponding to the most recent reversal from drainage to
imbibition, and ending at point (AS*idiAMd)» corresponding to the reversal
from imbibition to drainage on the preceding drainage curve. Rescaling the
main wetting branch of S*(h») to interpolate between the reversal points is
achieved via
(8.11)
where *S*(Ahjd) and *S*(Ah$j) are values of the %*(A*) function given by (8.9)
with h* = A$d and Ah£ii respectively.
For drainage scanning paths, an analogous procedure is used to rescale the
main drainage curve through appropriate reversal points. For an arbitrary
drainage scanning curve beginning at point (AS^d»AMd)» corresponding to the
most recent reversal from imbibition to drainage, and ending at point
(AS§jAhc|j), corresponding to the reversal from drainage to imbibition on the
preceding imbibition curve, the interpolation formula is
/rfc) (
(8.12)
where dS*(*li& and dS*(*Md> are values of the dS*{h*) function given by
(8.10) with h* = *h$i and AflJd, respectively.
Application of S-P scaling relations
To clarify the use of (8.11) and (8.12), we consider the prediction of scaled
saturation-pressure head relations for a hypothetical case illustrated in Figure
8.1 for a path sequentially passing through reversals denoted by points 1-5.
Chord L-»2 lies on the main drainage branch described by (8.10) as well as by
(8.12) with (AS?d,Ah*d) = (1,0) and
-------
LI
X
QL
<
_l
_J
i— i
Q_
<
U
D
LU
_l
<
U
0.0 i.o
SCALED SATURATION
Figure 8.1. Hypothetical saturation path on S*(h*) functional (arbitrary head
units).
Chord 2-»3 represents a primary imbibition scanning curve described by (8.11)
with (AS$p*h§j) corresponding to point 2 and (AS*d»AMd) corresponding to
point 1. For primary imbibition scanning curves (^S'Jd^h^d) = diO). Note
that the location of (AS*d»*Md) w*^ De unaffected by any other
imbibition/drainage sequences that may have occurred en route to point 2
since those paths have achieved closure (on the main drainage curve)
whereupon they are assumed to have no further influence on system
behavior.
Chord 3-»4 is a segment on the drainage branch of closed loop 2-»3-»2 described
by (8.12) with (^d^Md) corresponding to point 3 and (As£iAh&)
corresponding to point 2. Finally, chord 4-»5 is a segment on the imbibing
branch of closed loop 3-*H3 which will be described by (8.11) with (*Scli,*hiJi)
corresponding to point 4 and (*S*d,A.h*d) to point 3. If imbibition were to
continue beyond point 3, then primary scanning curve 2-»3 would be resumed
and the system would retain no memory of loop 3-*4-»3.
89
-------
•tfndP»l"l.00' J"Pton»nting the foregoing method numerically, it is
(8'12> ^^ the 8ame f°rm W in an
and ( SJLi, Jif_i) represent the most recent and preceding
, respectively, regardless of path direction, and p is an index
a drainage or imbibition scanning curve is considered. If
-'uir« « »u • L li *S*'h*> = and define Ah*0 = *S*(Ah*0) = 0, then the
°dd « I m he previous example will be described by (8.13) with p = d for
°ccur« , P V for even '• Note that when closure of a scanning curve
*hen %* m"Si decremented as S*(h«) resumes its previous path. That is,
tl»en i t A * f°r even *
-------
generalized to arbitrary two fluid systems, we estimate effective residual
saturations corresponding to a given primary imbibition curve by
SjS = Jk (8.14a)
in which
—3 £• (8.14b)
where AS is the wetting fluid effective saturation at the reversal from the
main drainage curve to a primary imbibition scanning curve of the 5]^* ( h ,-fc )
functional and lSjfJ^ is the maximum effective residual saturation
corresponding to the main imbibition branch of SjfJ^(hjjf). Definitions of
variables pertaining to (8.14) are illustrated in Figure 8.2. Note that for the
main drainage curve, effective and apparent saturations are equal.
Entrapped nonwetting phase effective saturation, Sj^Jk, along a specified
primary imbibition path will vary from zero when SjjJ* = ^SfcJ* .(La., the
reversal point) to SjfJb when Sgjk s 1 (i.e., when hjifQ), where SjfJ^ • SjfJk
+ Sjfik is tne apparent wetting phase saturation (see eq. 8.6). To estimate
the amount of trapped nonwetting fluid at intermediate locations along a
primary imbibition path, we invoke the expedient assumption that all pore
claases entrap nonwetting fluid in proportion to their volumes and that
comPress*on °* entrapped fluids can be disregarded at pressures typical in
the vadose zone. As a result, Sj^ is predicted to vary linearly with SgJk
according to
c . _ ,?.
sjt ~ sjr
for SgJJt* *Sj(Jlt (8.15)
a
where Sjr'* ia tfiven b7 (8.14). Note that once SjgJ* = *SkJk during
Drainage event, that the saturation path resumes the main drainage branch
and "**.
assumed relationship between SjfJ* and SfcJ* has the advantage of
implicity and will facilitate subsequent analyses of relative permeability
functions. Mild nonlinearity in the actual relationship may be inferred from
fre form of (8.14), however, experimental studies are not available to directly
lest various alternative formulations. In the absence of specific information
concerning variations in fluid entrapment for higher order scanning curves,
will simply assume that (8.15) holds for arbitrary order scanning curves
with a specified imbibition path. Future laboratory investigations
we
be necessary to fully evaluate these assumptions.
91
-------
u
I
(E
Q.
U
0. 0
1. O
WETTING FLUID SATURATION CSHJk>
Figure 8.2. Hypthetical two phase wetting fluid saturation-capillary head
relationship showing main drainage and main imbibition curves,
primary scanning imbibition curve and relevant terminalogy
(arbitrary head units). sy
three Phase Pyatema
J, oil imbibes in the vadose zone, a three phase system ensues with air-water
interfaces supplanted by oil-water and air-oil interfaces. Air and oil
mtrapment or release in the three phase system are assumed to be controlled
V0*?10™ °f au-?U and o"-water interfaces, respectively, in the porous
«anH ^ a"Proach' ln *« limit, the behavior observed in two phase
and oil-water systems. Air entrapped in the two phase air-water
ysten, prior to a contamination event will remain initially entrapped by wa\er
l±i t * ^i 8! 78tem< Sub8«1ua""y. ^ may be entrapped due to
advancing air-oil interfaces, corresponding to increasing apparent total liquid
aturation. Partitioning of entrapped air between oil and ^ater phases w?l
-°n ^ P°««on of the air-oil interface, assuming entrapped fluTds
e w«± rthm.the »°™* medium« I" « -i^r manner, oil entrapment
the water phase is assumed to be controlled by the position of the
r
air
we c°nsider procedures for estimating trapped oil
three phase systems from two phase data.
92
-------
Oil entrapment by water. Oil entrapment is assumed to be controlled by the
position of the oil-water interface in the porous medium relative to its position
when free oil enters the system as tracked by apparent water saturation, S^r
s Sw 1- Satw + S0t. We define the minimum apparent water saturation Sj/mjn
as the lowest value of S^ that has occurred at the location under
consideration since the introduction of oil. The maximum amount of trapped oil
will occur if free oil is continuously present during water imbibition from Sff
_ Swmin io Sff = 1. We designate this maximum trapped oil saturation as the
residual oil saturation for the three phase system and estimate its value from
two phase oil-water system behavior by replacing ASWOW in (8.14) by Sjymjn to
obtain
, ., aia
1 - Sff
1 + Jf^l - S/2n)
(8.16)
where Sor is the effective residual oil saturation in the three phase system
and Row is given by (8.14b) with jk = ow. To estimate the effective trapped
oil saturation at a given apparent water saturation, we employ the analog of
(8.15). Recognizing that no more oil can become entrapped than is present in
the system leads to the relation
'or
( I- Sf
(8.17)
where Sot is the effective trapped oil saturation in the three phase system.
Noting that Sw = S^ + Satw -f Sofr (8.17) may be solved for Sot(Sw); however,
when implemented in a numerical model, the apparent saturation formulation
will be more convenient since apparent saturations are directly obtained from
the scaled functional. The second condition in (8.17) implies the absence of
free oil in the system, i.e., Sjy = Sj>. Strictly speaking, the validity of (8.17)
requires that if free oil saturation becomes zero and subsequently increases,
water saturation does not increase during periods in which free oil is absent.
Modifications to accomodate circumstances when such conditions are not met
may be imposed; however, since the frequency of such occurrences seems
likely to be low and the effects minor, we do not consider the eventuality.
Air entrapment by, oil and water. Analysis of air entrapment follows in a
similar fashion to that for oil, but is complicated by the fact that both oil and
water phases may entrap air. Subsequent to oil entering the system, air
entrapment is assumed to occur in response to an advancing air-oil interface,
corresponding to increasing apparent total liquid saturation. Release of
entrapped air follows the recession of the air-oil interface, or if free oil
ceases to occur, the air-water interface. Air initially trapped by one fluid is
assumed to transfer passively to the other as the position of the oil-water
interface, tracked by apparent water saturation, moves.
93
-------
Prior to the introduction of organic liquid into the system, air entrapment in
water ia described simply by (8.14) and (8.15). Let ua denote the effective
water saturation corresponding to last reversal from main drainage to primary
imbibition curves prior to introducing oil as *Swa**r, and the apparent water
saturation at the time oil enters the system as Syf*"*9 = Sw'"t3 + ^t*"*3 where
the latter represent corresponding effective water and trapped air saturations.
la the limit as oil saturation approaches zero* the effective saturation of air
trapped by water should correspond to that given by (8.14a) for a two phase
air-water system with reversal from ASwaw as
1 -
. (8.18)
1 +
where J^w is given by <8,14b) with jk r air. Similarly, in the limit as only oil
displaces air after .inception of the three phase system* the effective
saturation of air trapped by oil should correspond to that for the two phase
air-oil system with a reversal from S^*"*3 given by
1 -
(8.19)
where J^, is given by <8.14b) with jk = ao.
the distribution of air trapped in the three phase system will depend on
positions of oil-water and air-oil interfaces, aa tracked by apparent water and
total liquid saturations, respectively, and also on the minimum pore size into
nhich air-oil interfaces have receded since the inception of a three phase
system, as indexed by the minimum apparent total liquid saturation since oil
introduction, designated as Sy""0. To characterize air entrapment in the
ihrse phase system, it will be necessary to distinguish three cases: (1} $j»o>in
Sff*-*3, ttH *SwaHr * Sr**n * Sit?*-*1, and (III) Sj*i** * *SW*W. We will
'iscues these cases in turn.
lisa I.
"his case corresponds physically to circumstances in which air-oil interfaces
kve never receded into pores with radii smaller than those which were water
Sled upon inception of the three phase system. Hence, all air entrapment in
ine oil phase is controlled by processes associated with the air-oil interface.
to entrapment in water is due to that initially trapped in the two phase
w-water system and to transfer from the oil phase controlled by the position
$ the oil-water interface. Total entrapped air effective saturation in this
aaa is given by
94
-------
Sat =
>arw
>aro
sT- sw'
1-
(8.20)
Partitioning trapped air between water and oil phases requires consideration
of the current oil-water interface position relative to that of the oil-water
interface at the inception of the three phase system and of the air-water
interface when air first became entrapped by water in a two phase system.
Effective saturations of air trapped in water are calculated as
3aro
sw-
(8.21a)
Satw = Sgnt
(8.21b)
= 0
(8.21c)
Air entrapment in oil may be obtained from (8.20) and (8.21), noting that
Case II.
This case corresponds to circumstances in which air-oil interfaces have, at
some time since the introduction of oil, receded into smaller pores than those
corresponding to S^"*3 but larger than those corresponding to ASwaw. In
this event, some of the air initially trapped by water in the two phase system
will pass into the oil phase. If ST later increases, subsequent air entrapment
by oil will be controlled by processes at the air-oil interface. Total entrapped
air effective saturation will be given by
Sat =
*arw
aw
*aro
I - Sj
tun
(8.22)
95
-------
Effective saturations of air trapped in water are calculated as
Satw -
sT
a
_ A
1 - A;
aw
stf-
i - ST
'ia
Sf*1" (8.23a)
,S
(8.23b)
Satw
= 0
(8.23c)
Case
HI.
If air-oil interfaces have ever receded into pores smaller than those
corresponding to ASwaw, the total amount of air entrapped in water and oil
phases may be calculated as
Sat =
'aro
UJJ
i -
W1B
(8.24;
Entrapped air in water and oil phases can be partitioned dependent on the
current apparent water saturation path according to
Satw -
i -
UJJ
(8.25a)
Satw =' 0
Sf*
1"
(8.25b)
\B previousl/i air entrapment in oil may be calculated as the difference
between total entrapped air and water entrapped air.
96
-------
TWO PHASE PERMEABILITY MODEL DEVELOPMENT
Since nonaqueoua phase organics are not natural constituents in groundwater
systems, analysis of groundwater contamination problems involving NAPLs
requires consideration of the behavior of air-water as well as air-NAPL-water
systems and possible transitions between them associated with NAPL
introduction and subsequent removal through natural or induced means. In
this section, we consider the NAPL-free case and develop relations for air and
water relative permeabilites in two phase air-water systems. It may be noted
that while model complexities become apparent only in dealing with the vadose
zone in which nonzero air saturation occurs, the case of water-saturated zones
is accomodated naturally as an endpoint of the model and requires no special
treatment.
Water Relative Permeability
We will assume that hysteresis in relative permeability relations is a result of
nonwetting fluid entrapment in the mobile fluid phases. The presence of
immobile air bubbles occluded by the water phase will result in an obstruction
to water flow. Simultaneously, trapped air will displace water into larger
pores which will tend to increase the water relative permeability at a given
actual water saturation above that which would occur in the absence of air
entrapment. To account for these effects, we modify Mualem's (1976) model for
water relative permeability in two phase air-water systems as follows
an
aw
.1/2
-.1
2
(8.26)
where krwaw is the relative permeability to water, Swavr and Sw°w are
effective and apparent water saturations, Sataw is the effective trapped air
saturation (superscript aw on all variables denotes the two phase air-water
system), Se is the proportion of effective porosity in which flow can
potentially occur ( = SyP™ in the air-water system), and h(Se) is a surrogate
for the pore size distribution of the medium. We assume the latter to be
represented by the inverse of the scaled main drainage function dS*(h*J
defined by (8.10). It may be noted that any scaled or unsealed main branch
function would perform equivalently, since any constant coefficient multiplying
heads will cancel in numerator and denominator of (8.26). The term in square
braces on the right hand side of (8.26) represents the ratio of mean flow
channel radii for water conducting pores corresponding to the present
distribution of water and air in the system, to that for the completely water
97
-------
saturated system. The first terra in the numerator represents all pores with
hydraulic radii smaller than the largest water filled pores, including 'those
occupied by trapped air, and the second term is a correction for pore space
occuppied by trapped air. The leading term outside the brackets on the right
hand side of (8.26) is a tortuosity correction.
We will assume that the saturation-pressure model presented in a previous
ec ion, including the analysis of fluid entrapment, reasonably describes the
system behavior. Differentiating (8.15) for the air-water system, and noting
to betnerSeHr • ^ *??" "" **<*>"* «"» ta the numerator o?
to be expressed in terms of Se as
dSf
(8.27)
'w
where *S*,*w is the effective water saturation at the reversal from the main
drainage branch and Sar™ is the effective residual air saturation in the
""
tohe t K
into the second term in the numerator of (8.26) yields
1/2
.
rtt
'(O, V] -
aw
where
F(a.b) =
1 -
(8.28)
(8.29)
??n.Uchtei? (1*8<» and our analysis of Chapter 2, (8.28) may be
to obtain a closed form expression for krw** as
1-
aw
1-
>ar
aw
aw
>ar
aw
2
(8.30)
98
-------
Note that krwaw is expressed in (8.30) as a function of both effective and
apparent water saturations, as well as effective residual air saturation and
effective water saturation at reversal from the main drainage branch.
Apparent saturation for the air-water case may be written in terms of
effective water saturation (Swaw), reversal effective water saturation (ASwaw)
and maximum trapped air effective saturation (iSwaw) via (8.6a), (8.14) and
(8.15) and effective residual air saturation may be expressed in terms of
"SvP" and 1Svaw via (8.14) and (8.15). Thus, in terms of primary system
variables, we observe krwa" to be a function of SMra"r, ^Si/"*, 'Swawr and m.
With some algebraic manipulation, (8.30) may be modified to such a form;
however, the result is rather cumbersome and of no greater utility in practical
usage than (8.30) since intermediate variables which are eliminated are needed
in the course of evaluating the saturation-pressure relations in any case.
Hence, we take (8.30) as the final form for krwaw.
Air Relative Permeability
Effects of trapped air on air permeability are more straightforward than those
on water permeability, since entrapped air serves merely to reduce the
saturation of free air which can participate in convection. Thus, the only
changes required in Mualem'a equation, as adapted in Chapter 2 for gas phase
flow, are to utilize apparent rather than effective saturation in the pore size
integrals and free rather than total air saturation in the tortuosity term as
1/2
o
(8.31)
where Sa/a"r is the effective free air saturation defined by (8.5a). Integration
of (8.31) yields
(8.32)
which indicates that kraaw is also a function of Swaw, ASwawt *Swaw and ID
since SWaw and Sa/*w are functions of 8**", ASMraw, and iSwraw' via (8.6a),
(8.14) and (8.15) and the requirement that SaaMr = 1 - Swawr.
99
-------
THREE PHASE PERMEABILITY MODEL DEVELOPMENT
To calculate relative permeabilites in three phase air-oil-water systems, we
assume that fluid wettabilities follow the order: water > oil > air so that water
occupies the smallest flow channels and air the largest — trapped saturations
aside. In the three phase system, consideration will need to be given to
effects of trapped air in the free oil phase and trapped oil in water as well as
trapped air in water.
Water Relative Permeability
Modification of Mualem's equation to account for trapped air and oil in the
water phase yields the expression for water relative permeability
,2
- S l/2
~ &
(8.33)
where krw is the water relative permeability in the three phase system, and
saturations are as defined previously. The first integral term in the
numerator of (8.33) represents all pores with hydraulic radii smaller than the
largest water filled pores, including those occupied by trapped air and oil,
and the second and third terms are corrections for pore space occupied by
trapped oil and air, respectively.
To express dSot in terms of Sg, we differentiate (8.17) after making the
substitution Se = S w to obtain
dSot =
1 -
dSf
(8.34)
where S^min is the minimum apparent water saturation since introduction of
oil and Sor is the effective residual oil saturation as defined by (8.16).
The integral which corrects for Sot can now be written as
(8.35)
100
-------
To account for effects of air trapped in the water phase, the amount and
distribution of the entrapped air must be known, Air can either be entrapped
by advancing air-water interfaces in a two phase air-water system prior to an
oil contamination event, or it may be entrapped by advancing air-oil interfaces
in the three phase system. Air initially trapped by water may later become
entrapped in the oil phase or released into the free air phase; and air initially
trapped by oil may likewise later become entrapped in the aqueous phase or
released to free air.
We assume air entrapment in water to be described by (8.21), (8.23) and (8 25)
For our present purposes, it will be convenient to decompose the effective
saturation of air trapped in water into two components such that Satw =
Satww + Satwo where Satww is the effective saturation of air trapped in
**«"" ~-—"—• *• i» FT »r — —— ———»•»w»*w h#w*Mvt*.ai*&w*i ui ai.r urtippoQ in
water by advancing air-water interfaces prior to oil introduction and Satwo is
the effective saturation of air trapped by advancing air-oil interfaces in the
three phase system which has subsequently passed into the water phase.
We may calculate the effective saturation of air trapped in water by advancing
air-water interfaces by
n - Asv
1 -
atf
,n
,n
(8.36)
aw
where Sarw is the residual air effective saturation in water defined by (8.18)
*Swaw ia the effective water saturation at the last reversal from the main
drainage branch prior to oil introduction, and n is the lesser of Syf or
The amount of air trapped by advancing air-oil interfaces and subsequently
passed to the water phase is calculated by
Satwo -
- ST
I - Sf
tin
am
(8.37)
where Saro is the residual air effective saturation in oil defined by (8.19).
To express the integral which corrects for entrapped air in the aqueous phase
in the numerator of the right hand side of (8.33) in terms of S» when Sri»n *
*Swaw and SRT* ^a* we obtain the following
101
-------
0
,n
'aro
- x*™
S
n
(8.38a)
If
SW°" and S"" S^" "- S°""- - 0 and ,8.38., 8implifiM to
0
(8.38b)
or when S^r
and
tsatw
= 0
(8.38c)
Following substitution of (8.35)
carried out to obtain closed
water and total liquid saturation paths.
The result for S^n 4 ^s^aw and
,«
f **?*• inte«rati°n may be
* krw 8ubject to different
1/2
'a/w
L-AS-fc
or
(8.39a)
102
-------
For STmin ' *Swavr and
STmin we obtain
1/2
i-
\-Sff
St
I-ST
*aro
i -
,2
(8.39b)
Finally, for Sw * ST *' and Sp ' *$„ the result is
1/2
1-
in
°or /
=T U-'
l-Sif20 l
(8.39c)
Prom (8.39) and formulae for ^,r, Sa^w. S^, Satw, Saio and Sot, it is seen
that krw in the three phase system is a function of the current liquid
saturations: 5^ and So; the historical saturations: ASwaw, Stf0**1, S^3"*3, and
ir> an(j the fluid-porous medium parameters: *Saraw, ^S^0, JSorow and m.
Oil Relative Permeability
An analogous procedure can be used to predict oil relative permeability in
three phase systems. For ahis case, the integral expression for kro is
modified from that given in Chapter 2 to account for effects of air trapped
the oil phase as
in
f
/2
Sff
Sato
,2
0
h(Sc)
0
h(Se)
(8.40)
103
-------
The amount and distribution of entrapped air in the oil phase will be a
function of the current water and total liquid saturation paths relative to
AS^W an£l ST™"1. The integral in the numerator of (8.40) which corrects for
entrapped air can be expressed in terms of Se for different fluid saturation
paths in a manner similar to that followed for air entrapment in water. For
Sw * '
ST
i-sT
a
s*
(8.40a)
For
***** and Sff * ST
tin
0
where t is the greater of Sff or A
(8.40b)
Finally, for
0
saro
h(Se)
h(Se)
(8.40c)
0 is the greater of Sff or
nun
Closed form expressions for oil relative permeability may now be obtained upon
substituting (8.40) into (8.39) and integrating, The result for SRT * ~-"»'" -•-
ro
1/2
1-
'aro
1-
(8.41a)
104
-------
For S
and
1/2
1-
'aro
mm
*arw
aw
,2
'arc
(8.41b)
For Sf *• *SW the solution is
= sof
1/2
1-
'aro
- (1-
'aro
(8.41c)
Oil phase relative permeability is observed to be a function of the all of the
system variables and parameters which arose in the general expression for
krwi although the sensitivity to the parameters will clearly be different as will
be demonstated in a later section.
Air Relative Permeability
Evaluation of air phase relative permeability follows in a manner analogous to
that for the two phase case, only noting that the lower boundary for free air
in the system is defined by apparent total liquid saturation rather than
apparent water saturation. The integral expression takes the form
/2
h(Se)
(8.42)
105
-------
which yields upon integration
(8.43)
Although the form of (8.43) is quite concise, decomposition of fife/ a**d &p to
primary system variables indicates that Jrra is a function of the same variables
as krw anc' ^ro1
Example for Main Imbibition and Drainage Paths
For the special case in which both water and total liquid saturations are
monotonically draining along their respective main drainage branches, the
relative permeability equations given here reduce to the forms derived in
Chapter 2. In order to evaluate the possible significance of Jc-S hysteresis
arising due to fluid entrapment, we will apply the model to calculate relative
permeabiltites for various extreme cases involving saturation paths moving
along main drainage or imbibition branches. For these calculations, we assume
multiphase Jc-S model parameters to be given by: JSkraw = 0.2, *Sarao = 0.1,
teorow = °*15« n = 2*°« which are reasonable values for a well graded sandy
porous medium. Note that other model parameters (*«, cfa, £,>) are not
necessary to evalute k-S relations.
Figure 8.3 shows calculated water relative permeability as a function of
effective water saturation corresponding to water and total liquid saturation
paths that are both monotonically draining or imbibing. Four krw(Sw) curves
are given corresponding to: (1) the main water drainage curve in which no
fluid entrapment occurs, (2) the main water imbibition curve in a two phase
air-water system, (3) the main water imbibition curve in a three phase system
in which the effective water saturation at commencement of the three phase
system (Sw2"*3) is zero, (4) as for case 3 but SW2~** = 0.3.
In all cases, predicted imbibition relative permeabilities are larger than
drainage relative permeabilities at a given water content. This reflects the
displacement of water into larger pores on wetting saturation paths due to
nonwetting fluid entrapment. This effect is in agreement with theoretical
arguments presented by Bear (1972). Marie (1981) also reports the
phenomenon of higher imbibition permeabilites to be frequently observed
experimentally.
In general, differences between drainage and imbibition water relative
permeabilites are rather small in magnitude, although differences increase near
apparent water saturation. The relatively mild hysteretic effects on wetting
phase permeability are as observed experimentally (Saraf et al., 1982).
pjfferences in maximum water saturations, reflecting differences in nonwetting
fluid entrapment, are clearly evident in Figure 8.3 for the different saturation
paths. For Case 1 there is no fluid entrapment and complete water saturation
is achieved, while for Case 2 the limiting effective water saturation is l-'S *»
- 0.8. Since S** 3 = 0 in Case 3, air is trapped only by advancing ai^oil
interfaces resulting in a maximum amount of trapped air equal to iSAfao - 0 1
However, oil is concurrently trapped, by advancing oil-water interfaces in an
amount which reaches a maximum of *Sor°* r 0.15. Thus, the maximum water
106
-------
saturation attained is l-'Saj.^^,.0"^ 0.75. In Case 4 more air entrapment
occurs due to the assumed greater propensity for air-water interfaces to trap
air relative to air-oil interfaces (compare *SjiJk values) resulting in greater
total fluid entrapment and lower maximum water saturation.
Figure 8.4 shows predicted air relative permeabilities as a function of effective
total liquid saturation when water and total liquid saturation paths are on
their main branches corresponding to Cases 1-3 above. Differences between
main imbibition and drainage air relative permeabilities are greater than those
observed for water relative permeabilities. This occurs because only a
relatively small volume of wetting fluid is displaced into larger pores relative
to the volume of displaced nonwetting fluid. Hysteresis in relative
permeability has been experimentally observed to be substantially greater for
nonwetting than for wetting fluids (Corey, 1986; Honarpour et al., 1986).
Oil relative permeabilities as a function of effective oil saturation are shown in
Figure 8.5 for main drainage and imbibition branches corresponding to several
constant effective water saturations. It is evident that oil relative
permeabilities for a given effective oil saturation are strongly dependent upon
water saturation, since the latter governs the size range of flow channels oil
may occupy. Effects of hysteresis and nonwetting fluid entrapment on oil
relative permeability-saturation relations appear to increase with increasing
effective water saturations. Note that main drainage oil relative
permeability-saturation relations for Sw = 0 in Figure 8.5 are the same as for
main drainage water relative permeability-saturation relations in Figure 8.3 as
would be expected (Corey, 1986).
Figure 8.3.
0.0 0.2 0.4 0.6 0.8 l.O
EFFECTIVE WATER SATURATION
Water relative permeability versus effective water saturation for
hypothetical soil following main imbibition and drainage paths.
107
-------
>
l-
UJ
z
o:
ui
iii
UJ
o:
10
10
10
10
01
< ID'4
O. 0 0. 2 O. 4 0. 6 O. 8 1. C
TOTAL LIQUID SATURATION
Figure 8.4. Air relative permeability versus effective total liquid saturation
for hypothetical soil following main imbibition and drainage paths.
o
1U
>-
—— di~a 1 nage
'• Imbibition
0.0 O. 2 O. 4 0.6 0.8 1. O
EFFECTIVE OIL SATURATION
Figure 8.5. Oil relative permeability versus effective oil saturation at several
constant oil effective saturations for hypothetical medium
following main imbibition and drainage paths.
108
-------
MODEL CALIBRATION PROCEDURES
In this section we discuss two methods that may be employed to calibrate the
hysteretic k-S-h model presented above. The first method utilizes measured
S-h data for main drainage and primary imbibition scanning paths from each
of two phase air-water, air-oil and oil-water systems. The second method
requires much less experimental effort. Detailed S-h measurments of the main
drainage branch of the air-water system are utilized with interfacial tension
data to estimate scaling coefficients and relatively simple measurements to
evaluate residual saturations.
Method 1
The input data for this method involve wetting fluid saturation-capillary head
pairs (Sk-hjfr) for two phase air-water, air-oil and oil-water systems in a
given porous medium for a saturation path following the main drainage branch
and a primary imbibition scanning curve. (Recall that we designate actual
wetting fluid saturations as Sk, effective saturations as Su and apparent
saturations as Sjtf. Associated with each data pair are three additional
variables which, in combination, permit unique specification of the Su-h ,-j,
measurement point vis a vis the parametric k-S-h model. The variables are-
(I) an index specifying the particular two phase fluid system involved, (2) the
capillary head at reversal from main drainage to primary imbibition curve
(*hjk), and (3) the wetting fluid saturation at the reversal point (AS\r)
corresponding to *hjh. It the S^-bjk data pair reflects a measurement on the
mamo ?ofinag° branch tl»e value of 0.0 is used for *hjk and 1.0 for ASk (see
eq. 8.13).
X
The data are employed to fit model parameters via a nonlinear least squares
regression analysis with actual saturations as the dependent variable. An
outline of the steps involved in predicting wetting fluid saturations
corresponding to a given parameter vector in the regression analysis follows.
Step 1. Convert actual reversal point saturations to effective reversal point
saturatton8's ^ <8'4) **"*" ""* estimated value for ineducable
Step 2. Scale actual heads for all measurement points as tf = BiJrh^ from
the estimated scaling coefficients, ft ft. J J
Step 3. Calculate model predicted scaled saturations S* corresponding to all
data points from (8.13), (8.9) and (8.10) using effective reversa
saturations from Step 1, scaled heads from Step 2 and current
estimates of the parameters 'a, da and n>
Step 4. Calculate residual nonwetting fluid effective saturations, SiJk
corresponding to the primary scanning paths from £5, using
effective reversal saturations from Step 1 and current estimates of
maximum effective residual saturations, &/-/*. estimates of
Step 5. Calculate effective trapped nonwetting fluid saturations, S,W*
corresponding to each S* on an imbibition path from (8.15) using
ASkJK from Step 1, SjjJif from Step 4 and noting that SKJk = s* for
fluid system jk. K
109
-------
Step 6. Convert each predicted SgJk to the corresponding predicted actual
saturation via (8.6) and (8.4) employing SjtJ* from Step 5 and the
current estimate of Sm.
Parameters which are fitted in the regression analysis using Method 1 are: *a,
dat nt and Sm, which are regarded as porous medium dependent parameters;
ftao and Vow considered to be fluid t pair dependent scaling factors; and
maximum effective residual saturations iSm-**, ^S^ao and *Sorow, which will
in general depend on both fluid and porous medium characteristics. Since the
reference for scaling is arbitrary, we select ftaw • 1 by definition.
Method 2
In the second method, we attempt to reduce the experimental data
requirements to a minimum to facilitate practical application of the model in
circumstances when the detailed analyses employed in Method J are not
available. Detailed saturation-capillary head measurements are assumed to be
limited to main drainage curves in air-water systems which provide estimates
Of da, n and S^r Following results reported by Kopl and Parker (1987) for a
number of soils, we estimate *a assuming the ratio Wda * 2. The fluid pair
dependent scaling coefficients ftao and ftow are estimated from interfacial
tension data as described in Chapter 3.
The remaining model parameters are the maximum effective residual
saturations, JSjjJ*, which may be obtained from relatively simple displacement
experiments of various designs. The experimental approach we assume here
involves a two step drainage-imbibition test. For each fluick pair jk (=aw, ao,
ovf), a core is initially saturated with wetting fluid k, displaced by fluid j to a
relatively low effective wetting phase saturation, followed by fluid k imbibition
to near zero capillary head. If the final capillary head is exactly zero, *Si*Jk
may be calculated directly from (8.14) and the measured reversal saturation.
If the final capillary head is nonzero, then the computational procedure is
rather more complicated and may follow that described for Method 1 except
that the data set will be considerably reduced (only two data points for air-oil
and oil-water experiments) and a number .of parameters are assumed
independently known, i.e., 0aw, foo, /Jow and J'«/d«, leaving only do, n, &„,
aMr' Jsarao «"d lSorow to be estimated from the saturation-pressure data.
EXPERIMENTAL INVESTIGATIONS
Saturation-capillary head relations were measured, in triplicate on similarly
packed soil cores, for air-water, air-NAPL and NAPL-water two phase systems
in a sandy porous medium containing 0.92, 0.05 and 0.03 kg kg"1 sand, silt
and clay, respectively. The NAPL was p-cymene (4-isopropyltoluene) which is
essentially insoluble in water (Dean, 1979). Interfacial tensions were measured,
in quintuplicate, by the drop method (Adamson, 1967) at room temperature
110
-------
(22.5 * 1.5*0) using liquids passed through the sandy porous medium to obtain
values reflective of those in the porous medium. The mean and 95X confidence
limits of air-water, air-oil and oil-water interfacial tensions were 67.7 t 0.6,
34.3 * 1.1 and 31.8 * 2.6 mN m~l, -respectively. Water used in the air-water
measurements was uncontaminated by p-cymene.
Retention cell description and experimental methods used in the S-h
measurements for the main drainage branch followed procedures described in
Chapter 2. To measure S-h relations of a primary imbibition saturation path
after all measurements were completed for the main drainage path, the
pressure of the nonwetting phase was monotonically reduced resulting in
wetting fluid imbibition into the porous medium through an outlet at the
bottom of the retention cell connected to a reservoir of wetting fluid
maintained at atmospheric pressure. Inflow of wetting fluid was monitored
with a buret to determine when fluid saturations in the porous medium were
in equilibrium with applied wetting-nonwetting fluid pressure differentials.
Fluid saturations were then determined by gravimetric measurements of the
cells after equilibrium was obtained for the air-water and air-oil systems, and
by recording the water inflow volume for the oil-water system. Six pressure
increments were applied to all replicates for each main drainage wetting fluid
saturation path and five or six increments for the primary imbibition scanning
path. This yielded 105 S-h data pairs from the air-water, air-oil and oil-water
measurements which were all utilized in the regression analysis for Method 1.
For Method 2 the complete set of air-water drainage data was used along with
final points on the imbibition paths and their associated drainage reversal
points for air-water, air-oil and oil-water systems.
Results
Parameters of the hysteretic k-S-h model estimated by Methods 1 and 2 are
summarized in Table 8.1. Main drainage and primary imbibition saturation-
capillary head data employed for model calibration are given in Table 8.2.
The ratio 1et/dat assumed in Method 2 (= 2.0) deviates somewhat from the
best-fit ratio of Method 1 (= 1.63). This magnitude of discrepancy is within
the range reported by Kool and Parker (1987) from analyses of hysteretic
air-water saturation-capillary head relations. Scaling factors ^ and /Jow
calculated from interfacial tension data in Method 2 are within the 95%
confidence intervals of corresponding values determined by nonlinear
regression against observed saturation-capillary head data in Method 1. For
Method 1 96.5% of the variation in measured saturation-capillary head data was
accounted for by the best fit parameters (i.e., Ra = 0.965). When measured
Sk(^jk) relations were compared to predicted relations by using parameters
obtained by Method 2, 93.9% of the variation in measured relations were
accounted for by the model. Permitting the parameter n in (1) and (2) to
vary for the main branches (i.e., *n and dn), resulted in an increase in the Ra
for the Method I regression of only 0.001, indicating that the assumption of
constant n for main drainage and imbibition curves is a reasonable one (see
also Kool and Parker, 1987).
Main drainage and imbibition curves of S*(h») predicted using Method 1 and
Method 2 parameters are shown in Figures 8.6 and 8.7, respectively, along
with experimental data that have been transformed into equivalent main
111
-------
branches of S*(A*) using appropriate model parameters. Transformation of
main drainage data involves only scaling capillary heads, since main drainage
actual and apparent saturations are identical when Sm = 0, as in the present
case. Transformation of primary imbibition data involves first calculating
apparent saturations corresponding to each actual measured saturation using
(8.4), (8.5), (8.14) and (8.15). Prom the calculated apparent saturation of the
primary imbibition path and known reversal points, the apparent saturation of
the main imbibition branch can be calculated via (8.11). Transformation of the
data in this fashion is performed for graphical clarity. Since a total of nine
different primary imbibition scanning curves are associated with the three
replicates of the two phase systems, individual saturation paths would be
indistinguishable if all primary imbibition scanning curves were plotted.
Precision and reproducability in the saturation-capillary head measurments for
these experiments was generally good as evidenced by the low scatter in the
transformed data in Figures 8.6 and 8.7. Scatter in the data appear to be
about the same for the different two phase systems, except for rather large
variations in air-water drainage data at higher scaled capillary heads.
The transformed experimental data correspond very closely with main curves
predicted using parameters best fit to the entire set of two phase
saturation-capillary head data in Method 1 (Figure 8.6). Deviations between
observed and predicted results using the much less data intensive Method 2
parameters are only slightly greater than for the best fit case (Figure 8.7).
Note that although only a subset of measured Sfc-hjjf data were employed to
calibrate S*(h*) in Method 2, all of the data are shown in Figure 8.7 to enable
evaluation of the accuracy in predicting behavior for all of the two phase
systems.
Table 8.1. Parameters in hysteretic k-S-h model by Methods 1 and 2.
Parameter Method 1 Method 2
4* (on-i)
4 (on-l)
a
Pao1
Po*
JSar*"
isar*°
isorow
0.053
0.086
1.923
1.881
2.135
0.246
0.086
0.276
0.073
(0.145)a
1.745
(1.920)'
(2.080)3
0.241
0.074
0.221
1 Arbitrary assignment of 0aw=l to scale to air-water system.
3 Constrained by assuming W^x = 2.
3 Calculated from interfacial tension data.
112
-------
Table 8.2. Measured main drainage and primary imbibition S-h data for sandy
soil with p-cymene as NAPL.
Fluid system Replicate 1
/path haw
Air-water
/drainage
Air-water
/imbibition
10.3
22.5
49.7
78.6
97.4
149.6
99.4
63.5
41.3
28.9
13.0
1.6
hao
Air-NAPL
. /drainage
^
Air-NAPL
/imbibition
8.1
13.5
16.3
20.3
30.0
58.9
46.9
30.0
21.8
13.3
4.6
1.2
how
NAPL-water
/drainage
NAPL-water
/imbibition
4.7
9.2
15.0
23.9
35.4
62.7
40.2
22.2
11.0
6.7
2.4
Sw
0.78
0.58
0.34
0.16
0.13
0.10
0.11
0.15
0.24
0.34
0.51
0.69
So
0.81
0.70
0.54
0.45
0.31
0.23
0.24
0.27
0.33
0.54
0.79
0.89
Sw
0.87
0.73
0.56
0.36
0.28
0.24
0.25
0.29
0.44
0.55
0.68
Replicate 2 Replicate 3
haw
10.3
22.5
49.7
78.6
97.4
149.6
99.4
63.5
41.3
28.9
13.0
2.2
hao
7.2
13.6
15.7
20.3
29.7
60.9
45.8
29.2
22.5
12.1
3.3
0.8
how
5.5
9.2
14.7
22.6
35.2
62.0
40.5
21.8
10.9
4.6
1.8
Sw
0.84
0.64
0.39
0.28
0.25
0.23
0.24
0.28
0.34
0.43
0.61
0.80
So
0.83
0.70
0.60
0.48
0.32
0.22
0.24
0.28
0.33
0.59
0.82
0.91
Sw
0.85
0.72
0.56
0.38
0.30
0.26
0.27
0.30
0.44
0.63
0.72
haw &w
10.3
19.4
46.8
74.2
98.9
148.2
99.0
61.0
41.7
27.4
13.7
2.1
hao
7.2
13.6
15.7
20.3
29.7
60.9
45.8
29.2
22.5
12.1
3.3
0.8
how
5.5
9.2
14.7
22.6
35.2
62.0
40.5
21.8
10.9
4.6
1.8
0.82
0.67
0.39
0.25
0.21
0.18
0.20
0.24
0.31
0.42
0.57
0.73
So
0.78
0.69
0.54
0.41
0.31
0.21
0.22
0.26
0.30
0.52
0.77
0.89
Sw
0.88
0.75
0.56
0.35
0.28
0.24
0.25
0.29
0.44
0.64
0.73
113
-------
E
3 ISO
O
U
I
> 100
-I
**
Q.
U SO
a
u
u
W
•AIR-WATCR DATA
•OIL-WATER DATA
•AIR-OIL DATA
O.O O. Z Q.4 O. B O. 8 l.O
APPARENT SATURATION
figure 8.6. Comparison of Method 1 predicted main drainage and imbibition
S*(h«) curves to measured two phase main drainage data (closed
symbols) and transformed primary imbibition data (open symbols).
•AJR-VATER DATA
•OIL-WATER DATA
*AIR-QIL DATA
0.0 O. 2 O. 4 O. 6 O. 6 1.0
APPARENT SATURATION
Figure 8.7. Comparison of Method 2 predicted main drainage and imbibition
S*(J)S) curves to measured two phase main drainage data (closed
symbols) and transformed primary imbibition data (open symbols).
. 114
-------
PREDICTED Jc-S-h RELATIONS FOR A HYPOTHETICAL CONTAMINATION SCENARIO
To evaluate the significance of hysteresis and fluid entrapment on three phase
k-S-h relations under conditions reflective of a NAPL contamination event in
the vadose zone, we will apply the model parameters obtained in the preceding
section to a problem involving a hypothetical saturation history corresponding
to a realistic sequence of saturation paths. Under the assumptions of the
proposed hysteresis model, constitutive relations will follow single valued
functions so long as directions of change in water and total liquid saturations
remain constant. Thus, four generic saturation path types in three phase
systems may be distinguished in rigid porous media which may be denoted as
II3, DD3, ID3 and DI3, where the superscript designates a three phase system,
first and second letters refer to directions of change for water and total
liquid saturations, respectively, and I and D denote increasing and decreasing
saturations. Cases in which either water or total liquid saturations are
constant simply degenerate the classification scheme. In ground water
contamination scenarios, the native system is initially NAPL-free necessitating
consideration of two phase air-water system behavior. Saturation paths for
the latter may be characterized as Ia or Da where imbibition or drainage of
water is designated.
The hypothetical NAPL contamination scenario is assumed to involve the
following sequence of events and corresponding saturation paths:
(1) At time A the porous medium is assumed to be water saturated and
undergoes monotonic drainage to time B (D3 path)*
(2) Rainfall occurs which increases the water saturation during the interval
between times B and C (I3 path).
(3) NAPL contamination occurs at time C resulting in increasing total liquid
saturation as NAPL imbibes while water saturation concurrently begins
to gradually decrease (DI3 path).
(4) Redistribution of NAPL following the spill or leakage event results in
diminishing NAPL and total liquid saturations beginning at time D while
water continues to slowly drain (DD3 path).
(5) A second contamination event occurs with the same organic liquid
leading to increasing total liquid saturations beginning at time E while
water saturation continues to decrease (DI3).
(6) Between times F and G remedial measures are employed (e.g., free
product removal at the water table) which results in a reduction in
NAPL saturation and total liquid saturation at the point under
observation while rainfall induces an increase in water saturation (ID3
• path).
Specific fluid saturations which are assumed to occur at times A through G
and at intermediate 'times on the same path segment employed in subsequent
calculations are given in Table 8.3. Using the parameters determined by
Methods 1 and 2 in the preceding section, we will predict saturation-pressure
and relative permeability-saturation histories corresponding to the assumed
saturation paths. Saturation-pressure relations will be presented as water
saturation versus scaled air-water capillary head when two phase conditions
prevail or versus scaled oil-water capillary head for three phase conditions
115
-------
(Figure 8.8) and as total liquid saturation versus scaled air-water capillary
head for two phase conditions or versus scaled air-oil capillary head for three
phase conditions (Figure 8.9). Permeability-saturation relations calculated with
the hysteretic k-S-P model are given as water relative permeability versus
water saturation (Figure 8.10), oil relative permeability versus oil saturation
(Figure 8.11) and as air relative permeability versus total liquid saturation
(Figure 8.12). Tabulations of effective saturations of trapped oil (S0t), air
trapped in the water phase (Sajw) and air trapped in the oil phase (Sato)
during the simulation period are given in Table 8.3 corresponding to both
calibration methods.
Table 8.3. Assumed effective fluid saturations and corresponding calculated
effective trapped nonwetting fluid saturations for the hypothetical
NAPL contamination scenario.
Saturations Entrapped nonwetting saturations
Method 1 Method 2
A
B
C
D
E
f
G
Sw
1.00
0.60
0.40
0.20
0.30
0.40
0.38
0.35
0.33
0.30
0.26
0.22
0.20
0.26
0.30
S0
0.00
0.00
0.00
0.00
0.00
0.00
0.16
0.40
0.32
0.20
0.36
0.53
0.60
0.30
0.10
sa
0.00
0.40
0.00
0.80
0.70
0.60
0.46
0.25
0.35
0.50
0.38
0.25
0.20,
0.44
0.60
Sot
0
0
0
0
0
0
0
0
0
0
0
0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
0.0
0.
0.
050
083
Sattt
0.0
0.0
0.0
0.0
0.041
0.082
0.073
0.061
0.053
0.041
0.024
0.007
0.0
0.045
0.075
Sato
0.0
0.0
0.0
0.0
0.0
0.0
0.034
0.084
0.074
0.059
0.097
0.138
0.154
0.066
0.007
Sot
0
0
0
0
0
0
0
0
0
0
0
0
0
0.
0.
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
.0
035
058
sat*
0.0
0.0
0.0
0.0
0.040
0.079
0.071
0.059
0.051
0.040
0.024
0.007
0.0
0.037
0.062
Sato
0.0
0.0
6.0
0.0
0.0
0.0
0.030
0.074
0.066
0.055
0.089
0.126
0.141
0.066
0.017
116
-------
120
90 -
E
0
\/
Q
LJ
I
o:
< BO
Q.
U
O
Ut
_J
U
in
OC
O. 0 O. Z O. 4 0.8 0.8
WATER SATURATION
l.O
Figure 8.8. Water saturation as a function of scaled capillary head for path
described in Table 3. Solid and dashed lines denote Method 1
and 2 parameters, respectively.
120
E
U
\s
o 8°H
u
I
V
o:
< BO
0.
<
u
Q
U
u
(A
ao -
O. O O. 2 O. 4 0. B O. 8 1.0
TOTAL LIQUID SATURATION
Figure 8.9. Total liquid saturation as a function of scaled capillary head for
path described in Table 3. Solid and dashed lines denote Method
1 and 2 parameters, respectively.
117
-------
O. Z O. 4 O. 6 O, B
WATER SATURATION
i.o
Figure 8.10. Water relative permeability as a function of effective water
saturation for path described in Table 3. Solid and dashed lines
denote Method 1 and 2 parameters! respectively.
O. O O. 2 o. 4 O. 6 O. 8 1. O
OIL EFFECTIVE SATURATION
Figure 8.11. Oil relative permeability as a function of effective oil saturation
for path described in Table 3. Solid and dashed lines denote
Method 1 and 2 parameters, respectively.
118
-------
O. Q Q. 2 O. 4 O. 6 O. 8 1.0
TOTAL LIQUID SATURATION
Figure 8.12. Air relative permeability as a function of effective air saturation
for path described in Table 3. Solid and dashed lines denote
Method 1 and 2 parameters, respectively.
mon Srf Pif , ' ? ..**!? vpha8e air-water 8ysten> «»«ts which drains
montomcally from an initially water saturated condition to a final effective
water saturation of 0.2. During period B-C, water imbibition occurs and air
entrapment is predicted to occur via (8.16) and (8.15) using * V'= 0.2 and
Sufi" for either calibration method. Values of predicted air entrapment for
the two calibration methods differ very little (Table 8.3). Effective water
saturation versus scaled capillary head paths predicted using Method 1 and
Method 2 parameters diverge slightly (Figure 8.8). At time C, air entrapment
HrPreK-r« *K resultjin an approximately 30% higher water relative
permeability than occurred at the same saturation on path A-B (Figure 8.10).
Differences between water relative permeability predictions via (8.30) based on
the two calibration methods are considerably larger than air entrapment
Efflll*' 8,U*?e8tin* the !•««• •"«*• "*y be difficult to identify in practfce.
Effects of air entrapment on air relative permeability calculated using (8.32)
'5?,. * ^ ^ ^ *«"" 8aturations occurring befween B
At time C, oil enters the system and water and total liquid saturation paths
diverge. The water saturation path commenses a secondary drainage scanning
curve which reaches a saturation of 0.35 at time D (Figure 8.8), while the total
liquid saturation path follows a primary imbibition scanning curve reaching a
saturation of 0.75 at time D (Figure 8.9). As oil imbibes and the total liquid
saturation increases, air becomes trapped by advancing air-oil interfaces.
119
-------
Also, as water saturation decreases, receding oil-water interfaces transfer to
the oil phase some air that was originally trapped by water in the two phase
air-water system prior to introduction of NAPL. Total entrapped air effective
saturations are estimated using (8.20) and effective saturations of air trapped
in water (Satw) from (8.21b). Air entrapment in oil (Sato) is obtained by
difference. Water, oil and air permeabilities subsequent to initiation of a three
phase systems are calculated from (8.39a), (8.41b) and (8.43), respectively.
During interval D-E, the total liquid saturation path changes from a primary
imbibition to a secondary drainage scanning path (Figure 8.9), while the water
saturation path remains on a secondary drainage scanning path (Figure 8.8)
As effective water saturation decreases from 0.4 at time C to 0.3 at time E the
minimum apparent water saturation Swmui which is used in S-h and k-S
relations is continually changing. The minimum apparent total liquid
saturation since oil introduction (ST""**) and the apparent water saturation at
the time oil enters the system (S^-3), however, remain constant and equal to
each other, e.g., with Method 1 parameters Sp""11 r S^ +$,+&» r 0.4*0+0.082 -
0.482 (see Table 8.3 and eq. 8.3a).
A second contamination event is presumed to increase oil saturations between
g and F leading to a total liquid saturation path which follows a tertiary
imbibition scanning curve until the total liquid saturation reaches it previous
reversal saturation of 0.75 (point D, Figure 8.9). With closure of loop D-E-D
the path resumes along the primary imbibition curve originating at the time fl
reversal point on the main drainage curve. Note that the Sw(ht) and St(h*)
curves of Figures 8.8 and 8.9 will not close at h* = 0 as would S*(h*) since
effective rather than apparent saturations are plotted.
At time F, the water secondary imbibition curve closes with the main drainage
curve at the point of departure corresponding to time B. Upon closure of the
scanning drainage curve with the main drainage curve, the model predicts all
air originally trapped by air-water interfaces will have been transferred to
the. oil phase. At time F, the total amount of trapped air, which is all in the
oiJ phase, is calculated by (8.20). Water, oil and air relative permeabilities for
saturation path E-F are determined using (8.39a), (8.43) and (8.41b).
tfote that at point B for the two phase air-water system and at point F
corresponding to the three phase system, (8.39a) and (8.30) yield the same
water relative permeability (Figure 8.10). Furthermore, whereas hysteresis in
5_h relations for closed loop scanning curves is apparent (i.e., loop B-C-F of
figure 8.8 and loop E-D-E of Figure 8.9), there is no hysteresis in k-S
relations of i water and air relative permeabilities for those saturation paths
(F,gures 8.10 and 8.12). This occurs because we assume the amount of
nonwetting fluid trapped due to fluid interfaces advancing from pore size
class a to pore size class b is equal to the amount of nonwetting fluid
released as interfaces recede from pore size class b to pore size class a. Note
that effective saturations of water and of air trapped in water are identical
for the point intermediate between B and C in Table 8.3, which occurs on a
primary imbibition path in the two phase air-water system, and for point E
^rhich is on a secondary drying path of the three phase system. This reflects
tjie fact that air trapped by advancing air-water interfaces between B and C
120
-------
Is predicted to be reversibly released from the water phase as oil-water
interfaces recede. water
Between times F and G, we assume that oil is removed while water slowly
imbibes such that total liquid saturation decreases. As water imbibes? ofl
becomes entrapped in the continuous water phase by advancing oil-water
interfaces and as the total liquid saturation decreases air-oil interfaces recede
into smal er pores releasing entrapped air into the 'free' air phase. The F-G
water saturation path once again follows a primary imbibition curve (Figure
8.8), while the total liquid saturation path follows a secondary drying scanning
curve (Figure 8.9). Since the minimum total liquid saturation and the water
saturation at inception of the three phase system are equal (S,*to=Sw*-*)M
air entrapment in the water phase subsequent to time F arUTvia transfer of
air trapped in oil, which was either initially trapped by advancing air-water
interfaces prior to the NAPL contamination event and transfered into olPas
oil-water interfaces receded (times C to F) or trapped by advanc£g a?r-o*!
interfaces after initiation of a three phase system, to the water pha^as
water saturation increases. Advancing oil-water interfaces associated with
increasing water saturation are predicted to simultaneously traiToil and th«
effect of greater total fluid entrapment (i.e., air anTo*) in water on S-J
relations can be observed in Figure 8.8 by the displacement of the F-G water
saturation path to the left of the B-C path. The effect can also^bser^din
Figure 8.10 where greater nonwetting fluid entrapment has the effect of
increasing water relative permeability at a given water saturation (i.e°,
compare curves F-G and B-C). Fluid entrapment effects on air relative
permeability relations become very substantial at high total liquid saturations
(Figure 8.12) as apparent total liquid saturations approach unity.
Notice that values of Satw at the intermediate point on saturation path B-C
and at point E are different from Satw of point G. This phenomenon can be
explained by the way we model nonwetting fluid entrapment. Apparent
saturations, which reflect the pore size, class of the fluid interfaces trapping
or releasing trapped nonwetting fluid, are used to model nonwetting fluid
entrapment. The apparent water saturation at point G (SorrS*,+&,,+£., -0 4«?ft
for Method 1 parameters) is greater than the apparent water saturat^atlh!
other points (S^O.341 for Method 1) because STis aTsobe!ng trapped for
water saturation path F-G (Sot=0.083 for Method J), whereas for preceding
saturation paths Sot =0. Therefore, a greater amount of a^is pred'cUd to
become occluded by the water phase as oil-water interfaces advance into
**
SUMMARY AND CONCLUSIONS
We have presented a model for permeability-saturation-pressure relationships
in porous media systems confining up to three immiscible fluid phases whfch
accounts for hysteresis and fluid entrapment effects. The model ™«
formulated by proposing a single hysteretic scaled saturation-capillary held
functional which describes water saturation vs. air-water capillary head in two
phase systems and water saturation vs. oil-water capillary head or total
saturation vs. air-oil water capillary head in three phase systems and
121
-------
employed to predict fluid permeabilities via a parametric model. The scaling of
heads is obtained by a linear transformation and saturations are scaled by
introducing the concept of apparent fluid saturations which include entrapped
nonwetting fluid contents with those of the fluid under consideration.
By employing an adaptation of van Genuchten's retention function to
characterize the pore size distribution of the system in conjunction with
Mualem's relative permeability model, closed form expressions for relative
permeabilites are obtained. Calculations for a hypothetical problem indicate
that the proposed model predicts behavior that is in accordance with general
behavior that has been reported in the literature for such systems. It may
be noted that other parametric models could have been utilized to obtain
similar ends, e.g. that of Brooks and Corey (1966), however, the evaluation of
possible advantages in various formulations must be left to future
investigations. Detailed measurements of three phase permeability-saturation-
pressure relations will be necessary to rigorously validate the model.
Calibration of the proposed hysteretic k-S-h model was investigated for a
sandy porous medium with p-cymene as the NAPL. Two methods were studied:
Method 1 in which model parameters were fit a to combined set of two phase
air-water, air-oil and oil-water drainage and imbibition saturation-capillary
head data, and Method 2 utilizing air-water drainage path data, single
imbibition S-h points for air-water, air-oil and oil-water systems and
interfacial tension data. Predicted and observed drainage and imbibition data
for all three of the two phase systems were described very closely by Method
I parameters and only slightly less well by those using Method 2. This is a
very encouraging result since the use of Method 2 greatly simplifies model
calibration. The most tedious aspect of calibration in Method 2, and probably
the most crucial in practice, is the determination of maximum residual
nonwetting phase saturations. The development of simpler experimental
procedures or empirical estimation methods for the latter would be of
substantial value.
Model parameters estimated by the two methods were employed to predict
three phase k-S-h relations during a hypothetical NAPL contamination
scenario. Differences in paths predicted by the parameters obtained for the
iv/o calibration methods were generally small again providing encouragement
for use of the simpler procedure. Fluid entrapment effects were predicted to
load to greater hysteresis in three phase S-h relations than were evident in
the two phase data. Effects of fluid entrapment on water relative permeability
Y/aa not great for the saturation paths simulated, while effects on air relative
permeability were large at high total liquid saturations.
The results presented in this chapter indicate that hysteresis and nonwetting
fluid entrapment effects in k-S-h relations of three phase systems may be
quite substantial. However, since the proposed model is based on a number of
hypotheses regarding the extension of two phase behavior to three phase
systems, validation will require direct comparisions between model predictions
and three phase system behavior.
122
-------
CHAPTER 9
NUMERICAL MODEL FOR MULTIPHASE COMPONENT TRANSPORT
Among the most widespread and hazardous contaminants are organic liquids of
low water solubility introduced into the subsurface environment via surface
spills, leaks from underground storage facilities, or seepage from improperly
designed or managed landfills or land disposal operations. Numerous
contamination incidents have been documented resulting from petroleum
hydrocarbons, industrial solvents and other organic liquids (e.g., Abriola,
1984). Such contamination problems generally involve complex mixtures of
multiple organic constituents which may move in aqueous and nonaqueous
liquid phases and in the gas phase. Modeling of these systems generally will
require consideration of multiphase fluid flow and multicomponent transport.
Analyses of multiphase flow pertinent to water resources problems have been
discussed recently by a number of authors (Faust, 1985; Osborne and Sykes
1986; Parker et al., 1986; Kuppusamy et al., 1987; Chapter 5 of tMsr^ru'
Corapciaglu and Baehr (1987) and Baehr and Coraplioglu (1987) have presented
a model for multiphase, multicomponent organic transport under conditions in
MMS1 ^ndCp?VH rfKrVnnNoo^ J1™1 *M phMe8' Abriola and ^der
(1985) and Pmder and Abriola (1986) have discussed the analysis of coupled
fluid flow and component transport in three fluid phase systems with mobile
water and NAPL phases. In previous chapters of this report, we have
addressed the problem of modeling fluid flow in three phase porous media^
systems. In this chapter, we present a model for coupled multiphase flow and
component transport which follows in a vein similar to that of Abriola and
Pmder, while introducing a number of simplifications to ease the problem of
model calibration and numerical implementation.
MODEL DEVELOPMENT
Multiphase Fluid Flow
To describe flow of the water phase, we neglect variations in aqueous phase
density due to fluid compression or to changes in phase composition by
confining our attention to migration of organic contaminants having low water
solubility in near surface environments. Assuming also an incompressible
porous matrix, the water phase continuity equation may be written as
«t ' ~**<*+ -£, iriow - riws - riwa)
123
-------
mass transfer of organic component i between NAPL and water phases (+ for
NAPL to water) per porous medium volume, nwa is the corresponding quantity
for transfer between aqueous and solid phases, and r£wa for water to gas
phase. The notation Zj denotes summation from i=l,...,n for n "noninert"
organic components, where we use the term "inert" to refer to organic
components considered to have both zero solubility in the aqueous phase and
zero volatility to the gas phase. Assuming fluid wettability following the
order: water^NAPlAair implies no air-water interfaces occur when a continuous
"oil" (i.e., NAPL) phase exists. Thus, in three phase air-NAPL-water systems
with a continuous oil phase riwa=0, while in a two phase air-water system
riow =0-
A generalized form of Darcy's law is taken as the governing equation for fluid
flow which may be written for the water phase as
in which krw is the water relative permeability; K8W = k/>wg/7jw is the
saturated water conductivity tensor where k is an "intrinsic" permeability
tensor, p* is water phase density, g is gravitational acceleration and r^ is
water viscosity; hw=PVf/flwg is the water pressure head where Pw is gauge
water pressure; and z is elevation. The viscosity of water and other fluids
will be regarded as independent of fluid composition over the narrow ranges
imposed by the restriction of the model to low water solubility organic
compounds.
For flow of NAPL, we likewise consider the phase density to be time invariant,
yielding the "oil" phase continuity equation
rioa) (9.2a)
where So is oil phase saturation, q0 is the oil phase volumetric flux density,
Po is the oil phase density, rfoa is the rate of mass transfer of organic
component i between oil and air phases (+ for oil to air) per porous medium
volume, and other terms are as previously defined.
Darcy's law for the oil phase takes the form
) (9.2b)
where kro is the oil relative permeability; K^, r k/>wg/7»o is the saturated
water conductivity tensor with TJO the oil dynamic viscosity; h0=Po/pwg is the
water height-equivalent oil pressure head where Po is gauge oil pressure; and
is the oil phase specific gravity where PQ is the oil density.
124
-------
In subsurface hydrology, pressure gradients in the gas phase are commonly
assumed negligible due to the low viscosity of air, thus effectively reducing
the number of equations governing How in the system by one. Such an
assumption may be justifiable under certain conditions in the sense that
solution to the liquid phase flow equations are often little affected by the
resistance to gas counterflow. It should be noted that assuming negligible
gas flow resistance is not equivalent to assuming zero gas flow, since
counterflow clearly must occur to equalize gas pressures. Thus, volatile
organics will in general be subject to convective as well as diffusive transport
if liquid saturations vary in time. In the present study, we will assume for
simplicity that gas flow resistance is negligible and also that variations in
total liquid saturation are sufficiently low in magnitude and/or frequency that
component fluxes m the gas phase are dominated by diffusion, enabling gas
flow to be disregarded altogether. We emphasize that future studies should
directly address the analysis of gas convection effects on tranpsort.
To solve the liquid flow equations, functions relating krw, kro, Sy,, and So to
hw and ho must be known. The model which we will employ follows that
outlined in previous chapters, but will not consider hysteresis effects since
coding of modules for hysteretic k-S-P relations is in progress as of the time
of drafting this report (July87).
Component Transport Model
Governing equations. In a multiphase porous medium system, transport of
organic components in all phases must generally be considered. For transport
of organic component i in the aqueous phase, the component conservation and
flux equations, respectively, may be written as
(9.3a)
(9.3b)
where ciw is the mass of component i per volume of water phase, Jiw is the
mass flux density of component i in the water phase of the prous medium, Diw
is the dispersion tensor for water phase transport of component i, and other
terms are as previosly defined.
»
Combining (9.3a) and (9.3b), expanding time derivative and convective flux
divergence terms, and using (9. la) to eliminate saturation time derivatives
yields
- rjwa - rjwa) + riow - riws - riwa (9.4)
125
-------
*Umm*iion over a11 components j=l ..... n including the
For component transport in the organic liquid phase, the analogous
conservation and flux equations are analogous
_
~v'Jio ~ rioa - r^ (9.5a)
(9.5b)
where cio is the mass of component i per volume of oil phase, Jfc is the mass
flux density of component i in the oil phase per porous medium area, and Din
ifl the dispersion tensor for oil phase transport of component i. Elimination of
the saturation time derivative from (9.5a) and (9.5b) via (9.2a) yields
*cio
*So T
- rioa (9.6)
The corresponding component transport equations for the gas phase assuming
negligible convective fluxes are H«»»O ««Buraing
*ciaSa
* ,t -"Jia + rioa + riwa
Jia = -*
which yields via continuity for the gas phase
(g
vrhere cia is the mass of component i per volume of air phase, J^ is the mass
dtZ7- fCOIDp0nrt *> ih° "r Pha8e "« P°~U8 »edium" area, D™ "
dispersion tensor for air phase transport of component i, and J is the
phase density which has been assumed constant.
0rganoiCiHC°InKPOnent8 ?d8orbed ,°n the TO»d phase are assumed immobile, so that
the solid phase species mass balance equation is simply
Pb ~ *» (9.9)
126
-------
where «i is the solid phase concentration given as mass of component i per
mass of solid phase, and pb ia the porous medium bulk density.
We shall assume differential phase velocities to be sufficiently small to invoke
local equilibrium as the control of interphase mass transfer. Equilibrium
partition coefficients are introduced defined by
Tio = cio/ciw (9.10a)
I"ia = cia/ciw (9.10b)
I" is = «i/ciw (9.10c)
where it may be noted that Tio and Pfe are dimensionless and Tia has units of
L3 M~l.
The final form of the i-component transport equation can now be obtained by
adding phase equations (9.4), (9.6), (9.8) and (9.9) and making use of (9.10)
under the assumption of local equilibrium to eliminate all concentrations except
as
*i* ~TT = v"Di*wciw ~ li*'vciw ~ 7i* (9.11)
where
i* = +sw + *s0rio + +saria '+ Pbris
iv + +s0rioDio + *sariaDia
1 TJQ
7i* = ""I! EjCrjow-rjws-rjwa) ~ -~— Ij(rjow+rjoa)
m Po
Note that (9.11) has the form of the conventional water phase convection-
dispersion equation with modified coefficients accounting for oil and air phase
transport. In the absence of an oil or air phase, (9.11) degenerates naturally
to the corresponding two phase air-water or single phase water-saturated
cases.
Transport model parameters. In (9.11), porosity and bulk density are regarded
as experimentally determined constants. To determine Fja and Tio, we invoke
Henry's and Raoult's laws to obtain
127
-------
* (9.12a)
Tic = fi/ciw* (9.12b)
where starred quantities denote values at saturation, which are readily
obtainable experimentally or via various approximation formulae (Lyman et al.,
J982) and
in which Mj is the molecular weight of component j and qo is the activity
Coefficient of component i in the oil phase. Note that while Fig's are constant,
the Flo'8 are concentration dependent except for a single- component NAPL in
wnjch case «{ equals the NAPL phase density. The solid partition coefficient
which is identical to the commonly employed "distribution coefficient" '- '
be directly measured or estimated using correlations with soil organic
arbon content and water solubility or octanol-water partition coefficients
.» Lyman et al., 1982).
decribe gas, water and NAPL dispersion tensors, we employ the general
form of Bear (1972) with a first order correction for nondilute solution effects
&nd *itn tne Millington-Quirk model for tortuosity (Jury et al., 1983). The jk
* roponent of the Dip tensor, where p = a,o,w is a phase index, is given by
(XT-XL)qpJqpk/lqpl*Sp] (9.13a)
= *'/3 Sp»/3 (9.13b)
= 1 - cip/pp (9.13c)
is Kronecker's delta, Dip is the molecular diffusion coefficient of
o0iponent i in free p-phase, XT and XL are transverse and longitudinal
jiapersivitiea which we regard as porous medium constants, qpJ and qpk are
of qp in the j and k coordinate directions, lopl is the scalar
of qp, Tp is a porous medium tortuosity, factor, and g;p is a
ondilute solutions which causes DiJk to vanish as phase
-erection for nondilute solutions which causes Dip to vanish as phase
O0ipo8ition approaches pure component i. Note that if the nondilute solution
correction were ignored, large dispersive fluxes may be predicted for single
c0fllponent organic liquids when in fact no oil phase dispersion is feasible.
vlolecular diffusion coefficients for a given component in the various phases
be estimated via the Stokes-Binstein equation and its relatives, which
that diffusion coefficients vary inversely with molecular radius and
viscosity (Lyman et al., 1982). Since liquid-gas viscosity ratios are on
order of 10s, gas phase diffusion may be expected to be of considerable
importance for components of appreciable volatility.
128
-------
NUMERICAL APPROACH
In the mathematical model outlined in the preceding section, n+2 coupled
differential equations, corresponding to n organic component transport
equations and two flow equations, must be solved. Even for small n, the
solution of such a system of equations becomes quickly impracticable,
especially for higher spatial dimensions (Pinder and Abriola, 1986). However,
the solution efficiency can be improved by noting that coupling between the
flow and transport equations will be quite weak. Because we are interested in
organic chemicals of low water solubility, rates of mass transfer between
aqueous and nonaqueous liquid phases will be small, and since gas densities
are low compared to liquid, transfer between liquid and vapor phases will
necessarily also be low. Hence, interphase transfer terms in the flow
equations (9.1) and (9.2) will generally have a very small effect on the
solution. Furthermore, in the present formulation, coupling between the n
transport equations arises only due to interphase mass transfer terms, which
are small, and to the dependence of rio on oil phase composition, which will
change very slowly m time.
A general numerical solution approach to take advantage of the resulting weak
coupling is as follows:
i. solve the water and oil phase flow equations simultaneously using
interphase transfer rates from previous time step or global iteration
(i.e., Steps i-iii),
u. solve the n component transport equations sequentially using fluid flux
densities from Step i and interphase mass transfer and Fjo values from
the previous time step or global iteration,
Hi. update T{o and calculate corrected interphase transfer rates by
back-substitution of concentrations from Step j'i into discretized forms
of the component transport equations (9.4), (9.6), (9.8) and (9.9),
iv. if concentrations meet convergence criteria relative to previous global
iteration, increment time step, otherwise repeat Steps i - iv with
updated interphase transfer terms.
Note that because the flow equations are nonlinear, Step i will involve an
iterative solution procedure, while the individual transport equations solved in
Step u are linear (with Tio lagged) enabling direct solution.
We have implemented a finite element solution of the multiphase transport
model described above for the case of a single noninert organic component
(i.e., n=l). The finite element formulation for the coupled flow equations
without interphase transfer has been described in Chapter 5 for a two
dimensional spatial domain discretized using bilinear quadrilateral elements.
The finite element formulation for the phase summed transport equation follows
an analogous procedure. Initial numerical experiments with the algorithm
given above for updating interphase mass transfer terms indicated that when
phase partition coefficients pertinent for water immiscible organic liquids were
employed, iteration of Steps i-iii was unnecessary. That is, sequential solution
of the coupled flow equations, followed by the phase summed transport
equation using time lagged interphaae transfer rates without iterative
refinement, had no deleterious effect on the solution accuracy.
129
-------
Figure 9.1. Benzene mass in water phase per porous medium volume,
in kg m~3 at U1.01, 3.34, and 10.5 days. Vertical exageration 3x.
Figure 9.2. Benzene mass in air phase per porous medium volume, *§a
-------
HYPOTHETICAL PROBLEM
To illustrate application of the multiphase transport model, a hypothetical
problem will be investigated involving contaminant migration accompanying
seepage of organic liquid from an underground storage tank. The domain is a
10 m deep by 60 m long vertical slice through the vadose and groundwater
zones with a contaminant source located in the unsaturated zone (Figures 9.1
and 9.2) and discretized by 757 elements. The initial distribution of water
pressure was in vertical equilibrium with a water table having a constant 1:30
slope and the system was initially oil free. Upstream and downstream water
levels were maintained at 5.5 and 3.5 m, respectively, and zero water flux was
imposed across the bottom of the aquifer and on vadose zone boundaries.
For a period of 1 day, NAPL was introduced at the bottom of the storage tank
under constant head. Zero NAPL flux was prescribed at all other boundaries
during the injection event and subsequently on the entire domain. Boundary
conditions imposed for the transport problem were zero concentration
gradients on all boundaries except the tank bottom during the injection
period, corresponding to zero dispersive flux in air, oil and water phases.
During injection of a single-component NAPL, cVf=cyf*=po/To was imposed at the
tank bottom corresponding to saturation of all phases with the organic
component.
Porous medium properties were chosen to correspond to a well graded sand
with +=0.4, o=5 m'1, n=2, Kaw=4 m d~l, \L=0.5 m and X-r=0.1 m. Molecular
diffusion coefficients in water and gas phases were taken to be 9.0 x 10"5 and
6.3 x 10~l ma d~l, respectively. The first simulation considered assumed
benzene as the NAPL for which /JO/PW =0.899, »b/*)w=0«65, Fo=505 and ra=0.17
(Weast, 1985) and fao=2.2 and /Sbw=1.9 (Chapter 2).
Contours of benzene mass in the water phase per porous media volume, tSwCw,
and in the gas phase, IS^cfe, at various times after initiation of leakage are
shown in Figures 9.1 and 9.2, respectively. Over the short period of the
simulation, convection of the water phase plume in groundwater was small.
The influence of vapor phase diffusion on transport in the vadose zone is
pronounced, although the total masses in the vapor phase remain low. Space
integrals of the water and gas phase masses of benzene are shown in Figure
9.3 as a function of time. Larger proportionate masses in the water phase
reflect an air-water partition coefficient for benzene 'which is less than unity.
Total masses in water and gas phases at the end of the simulation were only
0.54 and 0.11 percent of that in the NAPL phase, illustrating the extremely
long term contamination potential posed by immiscible organic chemical
releases.
To evaluate effects of partitioning behavior on contaminant migration, a second
simulation was run using fluid properties for o-xylene, for which po/^Q.BW,
7)6/^=0.81, 1*0=5,900, ra=0.20, /%o=2.1 and fow=1.9. Comparison of benzene and
o-xylene indicates very similar properties, with the exception of fo which
reflects an order of magnitude lower aqueous solubility for o-xylene.
Predicted distributions of benzene and o-xylene were very similar in their
pattern, with concentrations of the latter scaled proportionately lower. A
131
-------
comparison of space integrated water and gas phase masses of benzene and
o-xylene versus time are shown in Figure 9.3, which clearly illustrates the
reductions in o-xylene concentrations associated with its lower water
solubility.
SUMMARY AND CONCLUSIONS
r»n * chaJlY' We4 *r? . Ascribed a model for multiphase containant
transport and demonstrated its numerical implementation for a simple problem
involving a single contaminant moving simultaneously in water, gas and
nonaqueous liquid phases. A computational scheme involving decoupling of the
solutions for flow and mass transport equations was implemented and found to
provide satisfactory results. For future extensions to multicomponent
transport, this procedure is expected to provide very substantial savings in
computer execution costs. Much work remains to be done pertaining to the
development of numerical solutions for multiphase transport to consider
SEES ȣ ^Pfe* . NAPL mixtures, constituent interactions (e.g., cosolvent
effects) and biochemical transformations, nonequilibrium phase partitioning.
gas phase convection, hysteresis and other effects, as well as to improve
devlop algorithms and computational procedures which lead to improved
computational efficiency. Experimental studies should be pursued concurrently
to further develop and refine constitutive models for multiphase flow and
transport and to validate numerical codes which implement these.
• 10
TIME COAYS)
ia
10
TIME
13
Figure 9.3.
Space integrated contaminant mass in water and air phases
versus time for benzene and o-xylene simulations.
132
-------
REFERENCES
mh«mH *fultiph,a8f ""Cation of organic compounds in a porous medium
n A o™!^ / } ?° * ' .Lectur? Notes in Engineering, Vol. 8, C. A. Brebbia and
D. A. Orszag (eds.), Sprmger-Verlag, New York, 232 p., 1984.
Abriola, L. M. and Finder, G. F,, A multiphase approach to the modeling of
porous media contamination by organic compounds, 1. Equation development,
Water Reaour. Res., 21, 11-18, 1985a.
Abriola, L. II. and Finder, G. P., A multiphase approach to the modeling of
compound" 2- Numericai ai
°f Surface8' 2 ed., Interscience Publishers,
p.; eoleum Reservoir simuiation'
* 9' MBd C?ra,plio5lu' M- YM Groundwater contamination by petroleum
, 2. Numerical solution, Water Resour. Res., 23:201-213, 1987.
*" ^^^ ****' American «««vier Publishing Co.,
Belytschko, T. and Liu, W. K., Computational method for analysis of transient
of Porou-
* SM A.8imple method '°r determining unsaturated hydraulic
from moisture retention data, Soil Sci., 17, 311-314, 1974.
uufl Tnd Greenspan.' D" Numerical simulation of miscible and immiscible
fluid flows in porous media, Soc. Petr. Eng. J., 22, 635-646, 1982.
' •*' J" *nd Brown8c°n*«. B- R-, Further developments in
. I«L Mining
Coats, K. H., In-situ combustion model, Soc. Petr. Eng. J., 20, 533-554, 1980.
133
-------
Coraplioglu, M. Y. and A. L. Baehr, Groundwater contamination by petroleum
products, 2. Numerical solution, Water Resour. Res., 23:191-200, 1987.
Corey, A, T., Rathjena, C. H., Henderson, J. H. and Wyllie, M. R. J.t
Three-phase relative permeability, Trans. Soc. Petrol. Eng. Amer. Inst. Mining
Eng., 207, 349-351, 1956.
Corey, A. T., Mechanics of Immiscible Fluids in Porous Media, Water Resour.
publ., Littleton, CO., 1986.
Corey, A. T., Mechanics of Immiscible Fluids in Porous Media, 2nd ed. Water
Resour. Publ., Fort Collins, CO., 1986.
Dumore, J. M. and Schols, R. S., Capillary pressure functions and the influence
Of connate water, Soc. Petr. Eng. J. 14:437-444., 1974.
pean, J- A. (ed.), Lange's Handbook of Chemistry, 12th ed., McGraw Hill Brook
Co., New York, 1979.
gckberg, D. and Sunada, D. K., Nonsteady three-phase immiscible fluid
distribution in porous media, Water Resour. Res., 20, 1891-1897, 1984.
Faust, C. R., Transport of immiscible fluids within and below the unsaturated
aone: A numerical model, Water Resour. Res., 21, 587-596, 1985.
payer, M. J. and Hillel, D., Air encapsulation: I. Measurement in a field soil,
Soil Sci. Soc. Am. J., 50, 568-572, 1986.
payers, F. J. and Ma thews, J. D., Evaluation of normalized Stone's methods for
estimating three phase relative permeabilities, Soc. Pet. Eng., Reservoir Eng.
j., 24, 230-242, 1984.
Faust, C. R., Transport of immiscible fluids within and below the unsaturated
zone — A numerical model, Water Resour. Res., 21, 587-596, 1985.
per rand, L. A., Milley P. C. D. and Finder, G. P., Dual-gamma attenuation for
the determination of porous medium saturation with respect to three fluids,
Water Resour. Res., 22, 1657-1663, 1986.
pritton, D, D., Resolving time, mass absorption coefficient and water content
with gamma-ray attenuation, Soil Sci. Soc. Am. Proc., 33, 651-655, 1969.
Gillham, R. W., A. Klute and D. F. Heermann, Hydraulic properties of a porous
medium: measurement and empirical representation, Soil Sci. Soc. Am. J., 40,
203-207, 1976.
Hiemenz, P. C., Principles of Colloid and Surface Chemistry, Marcel Dekker.
York. 516 p., 1977.
Hoa, H. T., Gaudu R. and Thirriot, C., Influence of the hysteresis effect on
transient flows in saturated-unsaturated porous media, Water Resour. Res., 13,
992-996, 1977.
134
-------
Hochmuth, D. P. and Sunada, D. K., Groundwater model of two-phase immiscible
flow in coarse material, Ground Water, 23, 617-626, 1985.
Honarpour, M., Koederitz L. and Harvey, A. H., Relative Permeability of
Petroleum Reservoirs, CRC Press, Inc., Boca Raton, Florida, 1986.
Hopmans, J. W., and Dane, J. H., Calibration and use of a dual-energy gamma
system, pressure transducers, and thermocouples to determine volumetric
water content, dry bulk density, soil water pressure head, and temperature in
soil columns, Agronomy and Soils Department Series No. 101, Alabama
Agricultural Experiment Station, Auburn University, Alabama, 1985.
Hopmans, J. W. and Dane, J. H., Calibration of a dual-energy gamma radiation
system for multiple point measurements in a soil, Water Resour. Res., 7,
1009-1114, 1986.
Huyakorn, P. S. and Pinder, G. P., New finite element technique for the
solution of two-phase flow through porous media, Adv. Water Res., 1, 285-298,
1978.
Huyakorn, P. S. and Pinder, G. F,, Computational Methods in Subsurface Plow,
Academic Press, New York, 1983.
Kaluarachchi, J. and Parker, J. C., Effects of hysteresis on water flow in the
unsaturated zone, Water Resour. Res., (in review).
Kazemi, H., Vestal C. R. and Shank, G. D., An efficient multicomponent numerical
simulator, Soc. Petr. Eng. J., 18, 255-268, 1978.
Killough, J. E., Reservoir simulation with history-dependent saturation
functions, Trans. AIME, 261, 37-48, 1976,
Kool, J. B. and Parker, J. C., Development and evaluation of closed-form
expressions for hysteretic soil hydraulic properties, Water Resour. Res., 23,
105-114, 1987.
Kool, J. B., Parker, J. C. and van Genuchten, M. Th., Parameter estimation for
unsaturated flow and transport models: A review J. Hydrol., (in press), 1987.
Land, C. S., Calculation of imbibition relative permeability for two- and
three-phase flow from rock properties, Trans. AIME, 243, 149-156, 1968.
Leverett, M. C., Capillary behavior in porous solids, Trans. AIME, 142, 152-169,
1941.
Leverett, M. C. and Lewis, W. B., Steady flow of gas-oil-water mixtures
through unconsolidated sands, Trans. AIME, 142, 107-116, 1941.
Leverett, M. C., Lewis, W. B. and True, M. E., 1942, Dimensional-model studies
of oil-field behavior, Petroleum Technology, Tech. Paper 1413, January,
175-193.
135
-------
Lewis, R. W., Morgan, K., and Johnson, K. H., A finite element study of
two-dimensional multiphase flow with particular reference to the five-spot
problem, Computer Methods in Applied Mechanics and Engineering, Elaevier
Science Publishing, North-Holland, 1984.
Low, P. F, Viscosity of interlayer water in montmorillonite, Soil Sci. Soc. Am.
J., 40:500-506, 1976.
Lyman, W. J., Reehl, W. P. and Rosenblatt, D. H., Handbook of Chemical
Property Estimation Methods, McGraw-Hill, 1982.
Marie, C. M., Multiphase flow in porous media, Gulf Publishing Co., Houston,
Texas, 1981.
Miller, E. E. and Miller, R. D., Physical theory for capillary flow phenomena, J.
Applied Physics, 27:324-332, 1956.
Morel-Seytoux, H. J., Two-phase flows in porous media, Advances in
Hydroscience, Vol. 9, Academic Press, 1973.
Mualem, Y., A conceptual model of hysteresis, Water Resour. Res., 187, 514-520,
1974.
Mualem, Y., A new model for predicting the hydraulic conductivity of
unsaturated porous media, Water Resour. Res., 12, 513-522, 1976.
Naar, J. and Wygal, R. J., Three phase imbibition relative permeability Soc
pet. Eng. J., 1, 254-258, 1961.
Nakornthap, K. and Evans, R. D., Temperature-dependent relative permeability
and its effect on il displacement by thermal methods, (Tech. Paper SPE 11217)
Soc. Petr. Eng., Reservoir Eng. J., 230-242, 1986.
Nofziger, D. L., and Swartzendruber, D., Material content of binary physical
mixtures as measured with a dual-energy beam of gamma rays, J. Appl. Phys..
45, 5443-5449, 1974,
Oak, M. J. and Ehrlich, R., A new X-ray absorption method for measurement of
three-phase relative permeability, Soc. of Petr, Eng., SPE 14420, 1985.
Ogden, M. W., Deactivation and preparation of fused silica open tubular
columns for gas and supercritical fluid chromatography, Ph.D. Thesis, Virginia
polytechnic Institute and State Univ., Blacksburg, VA, 1985.
Osborne, M. and Sykes, J., Numerical modeling of immiscible organic transport
at the Hyde Park landfill, Water Rerour. Res., 22, 25-33, 1986.
Parker, J. C., Kool J. B. and van Genuchten, M. Th., Determining soil hydraulic
properties from one-step outflow experiments by parameter estimation, II.
Experimental studies, Soil Sci. Soc. Amer. J., 49:1354-1360, 1985.
Peery, J. H. and Herron, Jr., E. H., Three-phase reservoir simulation, J. Petr.
Tech., 21, 211-220, 1969.
136
-------
Finder, G. F. and Abriola, L. M., On the simulation of on aqueous phase
organic compounds in the subsurface, Water Resour. Res., 22, 109S-199S, 1986.
Rao, P. S. C., Rao, P. V., and Davidson, J. M., Estimation of the spatial
variability of the soil-water flux, Soil Sci. Soc. Am. J., 41, 1208-1209, 1977.
Reid, S., The flow of three immiscible fluids in porous media, Ph.D., Thesis,
Chem. Eng. Dept., Univ. Birmingham, England, 1956.
Saraf, D. N. and Fatt, I., Three phase relative permeability measurement using
a NMR technique for estimating fluid saturation, Soc. Pet. Eng. J. 7, 235-242,
1967.
Saraf, D. N., Batycky, J. P., Jackson, C. H. and Fisher, D. B.f An experimental
investigation of three-phase flow of water-oil-gas mixtures through water-wet
sandstones, proc. California Regional Meeting, Society of Petroleum Engineers,
San Francisco, CA, March 24-26, SPE 10761, 1982.
Scheidegger, A. E., Physics of flow through porous media, University of
Toronto Press, Toronto, 1974.
Scott, P. S., Farquhar, G. J. and Kouwen, N., Hysteretic effects on net
infiltration, Proc. National Conf. on Advances in Infiltration, Chicago, Am. Soc.
Agric. Eng., 1983.
Sheffield, M., Three-phase fluid flow including gravitational, viscous, and
capillary forces, Trans. AIME Petr. Div., 246, 232-246, 1969.
Shutler, N. D., Numerical three-phase simulation of the linear steamflood
process, Soc. Petr. Eng. J., 9, 232-246, 1969.
Snell, R. W., Three-phase relative permeability in an unconsolidated sand,
Trans. Soc. Petrol. Eng. Amer. Inst. Mining Eng., 84, 80-88, 1962.
Stone, H. L., Probability model for estimating three-phase relative permeability,
Trans. Soc. Petrol. Eng. Amer. Inst. Mining Eng., 249, 214-218, 1970.
Stone, H. L., Estimation of three-phase relative permeability and residual oil
data, J. Canadian Petrol. Tech., 12, 53-61, 1973.
Su, C. and Brooks, R. H., Soil hydraulic properties from drainage tests,
Watershed Drainage Proceedings, Irrigation and Drainage Div., ASCE, Logan,
Utah, p. 516-542, August 11-13, 1975.
Su, C. and Brooks, R. H., Water retention measurement for soils, J. Irrig.
Drain. Div., ASCE, IR2, 106, 105-112, 1980.
Touma, J. and Vauclin, M., Experimental and numerical analysis of two-phase
infiltration in a partially saturated soil, transport in Porous Media, 1, 27-55,
1986.
van Genuchten, M. Th., A closed-form equation for predicting the hydraulic
conductivity of unsaturated soils, Soil Sci. Soc. Amer. J., 44, 892-898, 1980.
137
-------
van Genuchten, M. Th. and Nielsen, D. R., On describing and predicting the
hydraulic properties of unsaturated soils, Annales Geophysicae, 3, 615-628,
Weast, R. C. (ed.), Handbook of Chemistry and Physics, CRC Press, Boca Raton,
Florida, 1985.
Weiss, G. (ed.), Hazardous Chemicals Data Book, Noyes Data Corp., Park Ridge,
New Jersey, 1980.
Wilson, J. L. and Conrad, S. H., Is physical displacement of residual
hydrocarbons a realistic possibility in aquifer restoration?, Proc. Petroleum
iC °h°mic^a in Ground Water- Houston, Nat. Water Well
Zimkiewicz, 0. C., The Finite Element Method, 3rd Ed., McGraw-Hill, New York,
138
-------
PUBLICATIONS AND PRESENTATIONS ASSOCIATED WITH PROJECT
Refereed Journal Papers
1. Parker, J. C., Lenhard, R. J. and Kuppusamy.T., A parametric model for
constitutive properties governing; multiphase flow in porous media, Water
Resour. Res., 23, 618-624, 1987.
2. Kuppusamy, T., Sheng, J., Parker, J. C. and Lenhard, R. J., Finite element
23^626-632, "**1""* immi8Cible fl°W throu*h a°^> Water Resour. Res.,
3. Lenhard, R. J. and Parker, J. C., Measurement and prediction of
' in three phaae
4. Parker, J. C. and Lenhard, R. J., A model for hysteretic constitutive
ae
f°r kyfrtic constitutive
fl°W' 2' ^—ability-saturation relations,
6. Lenhard, R. J. and Parker, J. C., A model for hysteretic constitutive
flow- 3-
review
and simulation
media' Water Resour' *•••. in
rdtelto ,K KJ' °- ExPerimental validation of methods for
™«Hi. w^ B 6e Phaa° ^^ration-pressure relations in porous
media. Water Resour. Res., in review. porous
Published Proceedings
1. Parker, J. C., Lenhard, R. J. and Kuppusamy, T., Modeling multiphase
contaminant transport in groundwater and vadose zones, Proc. Conf. on
2. Parker, J. C., Lenhard, R. J. and Kuppusamy, T., Measurement and
estimation of permeability-aaturation-pressure relations in multiphase
porous media systems, Proc. Conf. on Design and Opimization of Processes
in Natural Porous Media, Nancy, Prance, p. 139-150, 1987.
3. Parker, J. C., Kuppusamy, T. and Lien, B.-H., Modeling multiphase organic
chemical transport m soils and groundwater, Proc. Groundwater
Contamination: Use of Models in Decision Making, Amsterdam, Martinus
Nllnoff. in nrnan.
Nijhoff, in press.
139
-------
Published Abstracts
1. Lenhard, R. J. and Parker, J. C., Extension of two-phase capillary
pressure relationships to three phase systems, Trans. Am. Geophys. Union.
66, 893, 1985.
2. Kuppuaamy, T, Sheng, J., Parker, J. C. and Lenhard, R. J., Finite element
formulation of three-phase immiscible flow through soils, Trans. Am.
Geophys. Union, 67, 279, 1986.
3. Lenhard, R. J., Parker, R. J. and Kuppusamy, T., Parametric model for
constitutive properties governing multiphase fluid conduction in
~err0oc drocarbon P°roua "edia systems, Trans. Am. Geophys. Union,
, 100D,
4. Lenhard, R. J., Parker, J. C. and Kuppusamy, T., Prediction of multiphase
flow phenomena from scaled two phase retention data, Soil Sci. Soc. Am.
Meeting Abst. p. 159, 1986.
5. Parker, J. C., Kuppusamy, T. and Lenhard, R. J., Modeling organic
chemical transport in three fluid phase porous media systems, Trans. Am.
Geophys. Union, 67, 945, 1986.
6. Lenhard, R. J. and Parker, J. C., Measurement of three-phase fluid
pressure-saturation relations in porous media, Trans. Am. Geophys. Union,
67, 945, 1986.
7. Huyakorn, P. S., Saleem, Z. A., Parker, J. C. and Lester, B. H.,
Alternative modeling approaches for subsurface transport of nonaqueous
phase organic wastes, Trans. Am. Geophys. Union, 68, 326, 1987.
8. Lenhard, R. J., Dane, J. H., Parker, J. C., and Kuppusamy, T.,
Measurement of oil and water flow through porous media, Trans. Am.
Geophys. Union, 68, 311, 1987.
Invited Presentations
1. Seminar at Battelle Pacific Northwest Laboratories, Richland WA, entitled
Modeling multiphase organic fluid transport," August 1985.
2. Seminar at University of Nevada Environmental Research Center, Las Vegas
NV, entitled "Modeling multiphase fluid flow in soils and groundwater,"
May 1986.
t
3. Presentation at Gordon Research Conference on Modeling Plow in Permeable
Media, Andover NH, entitled "Analysis of the inverse problem for flow in
porous media with multiple fluid phases," July 1986.
4. Paper at Society of Engineering Science Annual Meeting, Buffalo NY,
entitled "Modeling heterogeneous fluid conduction in soils," August 1986.
140
-------
5. Presentation at workshop on Modeling Oily Wastes, Washington DC, entitled
"Data requirements for modeling immiscible organic flow and transport,"
January 1987.
6. Seminar at Oak Ridge National Laboratory, Environmental Sciences Division,
Oak Ridge TN, entitled "Modeling multiphase flow and transport in soil and
groundwater," January 1987.
7. Seminar at Los Alamos National Laboratory, Environmental Sciences
Division, Los Alamos NM, entitled "Modeling multiphase flow and transport
in soil and groundwater," March 1987.
8. Seminar at Electric Power Research Institute, Environmental Sciences
Department, Palo Alto CA, entitled "Saturation-pressure relations for three
phase flow in porous media," March 1987.
9. Seminar at Stanford University, Civil Engineering Department, Stanford CA,
entitled "Modeling multiphase flow and transport in soil and groundwater,"
March 1987.
10. Seminar at the University of Stuttgart, Institute for Hydraulic
Engineering, Stuttgart, W. Germany, on "Modeling multiphase flow and
transport in geologic media," June 1987.
11. Seminar at the Swiss Federal Institute (ETH), Zurich, Switzerland, on
"Modeling multiphase flow and transport in geologic media, June 1987.
12. Seminar at^ the Technion Institute, Civil Engineering Department, Haifa,
Israel, on Modeling nonaqueous phase organic chemical transport," June
1987. *^
141
------- |