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

-------