-------
DISCLAIMER
The results reported in this document have been funded wholly or in part by the United
States Environmental Protection Agency under Cooperative Agreement No. CR-814320
to Virginia Polytechnic Institute and State University. It has been subjected to the
Agency's peer and adminstrative review, and it has been approved for publication as an
EPA document. Mention of trade names or commercial products does not constitute
endorsement or recommendation for use.
QUALITY ASSURANCE STATEMENT
All research projects making conclusions or recommendations based on environmentally
related measurements and funded by the Environmental Protection Agency are required
to participate in the Agency Quality Assurance Program. This project was conducted
under an approved Quality Assurance Project Plan. Information on the plan and
documentation of the quality assurance activities and results are available from the
Principal Investigator.
ACKNOWLEDGEMENTS
Jong Soo Cho with the U. S. Environmental Protection Agency, R. S. Kerr
Environmental Research Laboratory in Ada, Oklahoma and Lou G. Swaby with the
U.S. EPA Office of Research and Development in Washington, D.C., served as project
officers whose technical and administrative assistance throughout the course of the
study were much appreciated.
During the second half of the project, Bob Lenhard left VPI to work for Battelle Pacific
Northwest Laboratories, where he continued to contribute to the project with support
from the Ecological Research Division, Office of Health and Environmental Research,
U.S. DOE under contract DE-AC06-76RLO as part of OHER's Subsurface Transport
Program.
11
-------
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 implement 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 programs 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
characterizing 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 applicability and limitations
of using natural processes indigenous to the soil and subsurface environment, for the
protection of this resource.
This report describes the development of a computer model for multiphase flow and
transport in the subsurface environment and laboratory validation of the model. The
research assesses the behavior of immiscible phase flow and contaminant transport in
each phase.
Clinton W. Hall
Director
Robert S. Kerr Environmental
Research Laboratory
in
-------
ABSTRACT
Groundwater contamination due to surface spills or subsurface leakage of hydrocarbon
fuels, organic solvents and other immiscible organic liquids is a widespread problem
which poses a serious threat to groundwater resources. In order to model the movement
of such materials in the subsurface, it is necessary, in general, to consider flow of water,
nonaqueous phase liquid (NAPL) and air, and transport of individual chemical
components, which may move by convection and dispersion in each phase. A
mathematical model is developed for multiphase flow and multicomponent transport in
porous media with water, NAPL and air or any subset of these phases. Numerical
procedures for solving the system of coupled flow equations, based on various
formulations of the governing equations, were compared. An adaptive solution
procedure is developed that automatically eliminates or includes equations in
subdomain regions, to take advantage of the fact that flow may be near steady state in
large parts of a physical domain. A linear phase-summed component transport equation
is derived for the case of local equilibrium interphase mass transfer. A semi-decoupled
solution approach is employed with transport equations solved serially with flow
equations using time-lagged interphase mass transfer terms in the flow equations. For
the case of rate-limited interphase mass transfer, a phase-summed transport equation is
derived in terms of apparent partition coefficients, which are functions of the current
concentrations and interphase mass transfer rates, hence requiring an iterative solution
procedure.
Accurate representation of three-phase permeability-saturation-pressure (k-S-P)
relations is crucial to model multiphase fluid movement and accurate models for
interphase mass partitioning are critical to describe species transport. A detailed
physically-based model for hysteresis in three-phase k-S-P relations is described.
Simplified models, which consider effects of nonwetting fluid entrapment, are shown to
provide a reasonable compromise between accuracy, on the one hand, and efficiency and
robustness, on the other. Laboratory studies of light and dense NAPLs in a 1 x 1.5
meter sand tank, involving measurements of water and NAPL pressures and saturations
and component concentrations are described and were used to validate the
mathematical model for multiphase flow and transport. Column experiments performed
to study NAPL-water interphase mass transfer kinetics indicate that existing empirical
models for mass transfer coefficients are reasonably accurate.
IV
-------
CONTENTS
DISCLAIMER ii
QUALITY ASSURANCE ii
ACKNOWLEDGEMENTS ii
FOREWORD iii
ABSTRACT iv
CONTENTS v
1. Summary of project objectives and results 1
1.1 Project overview 1
1.2 Conclusions and recommendations 2
2. Numerical analysis of multiphase flow 4
2.1 Introduction 4
2.2 Development of implicit pressure-pressure formulation 6
2.2.1 Theory and mathematical formulation 6
2.2.2 Numerical investigations 17
2.3 Evaluation of alternative equation formulations 30
2.3.1 Formulation of governing equations 30
2.3.2 Solution procedure 34
2.3.3 Numerical results . 39
2.4 Development of an adaptive solution method 47
2.4.1 Description of ASD algorithm 49
2.4.2 Numerical results 50
-------
3. Hysteresis in three-phase flow relations 57
3.1 Evaluation of general hysteretic formulation 57
3.1.1 Model description 57
3.1.2 Numerical simulation of hysteretic flow 69
3.2 Development and testing of a simplified hysteretic model 80
3.2.1 Model description 81
3.2.2 Example problems 86
3.3 Experimental verification of hysteretic flow model 93
3.3.1 Static three-phase measurements 93
3.3.2 Two-phase dynamic measurements 106
3.3.3 Three-phase dynamic measurements 121
4. Equilibrium-controlled multiphase transport 126
4.1 Mathematical model for multiphase transport 126
4.1.1 Governing equations 126
4.1.2 Phase-summed equations for local equilibrium transport 130
4.1.3 Initial and boundary conditions 132
4.2 Numerical model description 134
4.2.1 General solution approach 134
4.2.2 Finite element formulation 136
4.2.3 Interphase mass transfer and density updating 138
4.3 Model applications 140
4.3.1 Example 1 140
4.3.2 Example 2 145
4.3.3 Example 3 148
VI
-------
5. Transport with nonequilibrium interphase mass transfer 152
5.1 Model for nonequilibrium interphase mass transport 152
5.1.1 Governing phase transport equations 152
5.1.2 Consideration of nonequilibrium interphase mass transfer. . . . 152
5.1.3 Solution approach 156
5.1.4 Numerical verification 157
5.2 Laboratory investigations 159
5.2.1 Experimental methods 159
5.2.2 Data analysis methods 163
5.2.3 Experimental results 168
6. Two dimensional laboratory studies of multiphase flow and transport .... 171
6.1 Experimental setup 171
6.1.1 Cell design 171
6.1.2 Gamma attenuation measurements 171
6.1.3 Pressure measurements 173
6.1.4 Sample collection and analysis 174
6.2 LNAPL experiment 174
6.2.1 Experimental methods 174
6.2.2 Model calibration and numerical analysis 177
6.2.3 Results 180
6.3 DNAPL experiment 193
6.3.1 Experimental methods 193
6.3.2 Model calibration and numerical analysis 195
6.3.3 Results 197
7. References ' 202
Appendix: Publications resulting from this project 206
vii
-------
1. SUMMARY OF PROJECT OBJECTIVES AND RESULTS
1.1 Project Overview
Spills and subsurface leaks of immiscible organic liquids are a frequent cause of
groundwater contamination. This project was undertaken to improve our understanding
of contaminant migration arising from such events and to develop improved
quantitative tools for their description. Attention was focussed on three fronts involving:
(1) the development of efficient and robust numerical methods for solving the difficult
problem of simulating simultaneous multiphase flow and transport, and (2) developing
and numerically implementing theoretical and empirical constitutive models governing
multiphase flow and transport processes, and (3) performing laboratory scale
experimental studies to validate the mathematical models developed in conjunction
with the first two objectives.
Chapter 2 of this report describes the development and testing of numerical methods for
solving multiphase flow problems. The basic governing equations for three-phase flow
are presented and various numerical formulations are derived and compared vis a vis
their efficiency, accuracy and robustness. A new algorithm is presented and discussed
which enables substantial reductions in computational effort by automatic elimination
and inclusion of elements in the global solution.
Chapter 3 presents results of numerical and experimental studies undertaken.to develop
and test constitutive models for permeability-saturation-capillary pressure (k-S-P)
relations governing three-phase flow. A rigorous, physically-based hysteretic k-S-P
model is described as well as various simplified models. Numerical comparisons of the
models are presented and results of static and dynamic laboratory column experiments
are compared with the models to evaluate their accuracy.
Chapter 4 describes the theoretical foundation for modeling coupled multiphase flow
and multispecies transport with equilibrium interphase mass transfer, leading to a
system of flow equations for each bulk phase and phase-summed species transport
equations for each species. A numerical formulation for solving the system of equations
is presented based on partial decoupling of flow and transport equations. Results of
several hypothetical numerical simulations are presented to verify the model and to
demonstrate its applicability to specific problems.
-------
In Chapter 5, the assumption of local equilibrium-controlled mass transfer is relaxed
and a formulation is derived which enables the form of the phase-summed equilibrium
model to be retained by introducing apparent partition coefficients that depend on the
mass transfer rates and first-order mass transfer coefficients, as well as the true
equilibrium partition coefficients. Hypothetical numerical simulations are presented to
verify the formulation and to demonstrate effects of nonequilibrium mass transfer.
Laboratory experimental studies are presented which indicate interphase mass transfer
can be described accurately at the laboratory scale by a first-order mass transfer
relation and that mass transfer coefficients may be estimated from empirical relations
previously reported in the literature.
Chapter 6 presents results of two-dimensional experiments involving simulated spills of
light and dense organic liquids. The design of the "sand tanks" and procedures for
measuring fluid saturations and pressures and component concentrations is described.
Experimental results are compared with numerical simulations to provide validation of
the coupled flow and transport model described in Chapter 4.
1.2 Conclusions and Recommendations
In this report, we present a mathematical model for multiphase flow and component
transport which is based on fundamental principles. Solution of the resulting set of
nonlinear coupled differential equations is computationally difficult. The numerical
results have indicated that the accuracy and robustness of multiphase flow problems is
quite sensitive to the method of formulating the numerical model. In particular, it has
been shown that chaining of saturation time derivative terms can lead to mass balance
errors and other difficulties that are greatly reduced by means of an alternative
formulation that treats saturation time derivatives directly as implicit functions of
pressure. Marked reductions in computational effort were achieved by an adaptive
solution procedure that takes advantage of the fact that flow may be near steady state
in large parts of the physical domain. Implementation of schemes to selectively
eliminate individual phase equations should prove even more effective in reducing
computational effort as well as improving solution robustness. Other facets affecting
numerical solution efficiency, accuracy and robustness, such as matrix solution
methodology, should be pursued particularly to facilitate practical extensions to
three-dimensional problems.
-------
Solution of coupled multiphase flow and transport problems introduces additional
numerical difficulties. The semi-decoupled phase-summed approach employed in this
study is an efficient formulation, although the efficiency undoubtedly comes at the cost
of accuracy. The most critical aspect of the decoupled approach is the back-calculation
of interphase mass transfer rates, since accumulated small errors can eventually
destabilize the solution. Such problems are greatly diminished by suppressing mass
transfer calculations during periods of highly transient oil flow. Since compositional
changes are generally small over short time periods, mass balance errors incurred by
this approach are very small. Future investigations to investigate alternative solution
formulations and to develop algorithms for automatic component mass balance controls
on the solution should be pursued (e.g., "re-inject ion" of mass balance errors or iteration
of transport and flow with corrected phase transfer terms).
Accurate representation of k-S-P relations is crucial in order to predict multiphase fluid
movement in the subsurface, and accurate models for mass interphase partitioning are
critical to describe species transport. We have developed a detailed model for hysteresis
in k-S-P relations which has proven to accurately describe static laboratory S-P data
and transient flow response. Simplified models which consider effects of nonwetting fluid
entrapment provide a reasonable compromise between accuracy, on the one hand, and
efficiency and robustness, on the other. Direct measurements of three-phase relative
permeabilities for non-monotonic saturation histories is a formidable task, but one that
would be highly desirable to undertake in the future.
Existing empirical relations for estimating column scale mass transfer rate constants
from basic physical properties of the system yielded surprising accuracy. However.
under common field conditions, it may be shown that mass transfer limitations should
be of minor importance. Field scale heterogeneity, however, may lead to preferred flow
paths that produce apparent nonequilibrium effects controlled by diffusive mass transfer
between "fast" and "slow" zones.
Since explicit treatment of heterogeneity is generally impractical, future studies should
address the feasibility of describing field scale behavior using effective large scale mass
transfer relations. Likewise, effective k-S-P relations at the field scale that implicitly
consider effects of heterogeneity may well differ from laboratory scale relations.
Substantial future efforts will be needed to more fully understand and model the
behavior of large scale heterogeneous systems.
-------
2. NUMERICAL ANALYSIS OF MULTIPHASE FLOW
2.1 Introduction
Groundwater contamination due to surface spills or subsurface leakage of hydrocarbon
fuels, organic solvents and other immiscible organic liquids is a widespread problem
which poses a serious threat to groundwater resources. In order to describe the
movement of such materials in the subsurface, it is necessary in general to account for
separate phase flow of nonaqueous phase liquid (NAPL) as well as water and air. A
substantial obstacle to modeling multiphase flow has been the difficulty in determining
constitutive relationships for three-phase flow. Recent work pertaining to the
parametric description and calibration of three-phase permeability-saturation-pressure
relations including hysteretic effects has diminished these difficulties to a significant
extent (e.g., Parker et al., 1987; Lenhard and Parker, 1988; Parker and Lenhard, 1987).
Numerical models for NAPL migration in the vadose zone and groundwater have been
presented recently by various researchers based on finite difference methods (Abriola
and Pinder, 1985b; Faust, 1985; Falta and Javandel, 1987) and on finite element
methods (Osborne and Sykes, 1986, Kuppusamy et al., 1987). Most of these efforts
have focused on basic aspects of numerical model formulation with little attention to
the effects of various techniques on computational efficiency and solution accuracy. A
variety of formulations of the governing equations for NAPL and water flow have been
taken as a starting point for numerical model development. Abriola and Pinder
(1985a,b), Osborne and Sykes (1986), Kuppusamy et al., (1987) and Kaluarachchi and
Parker (1989) used flow equations in terms of fluid pressures with mass storage terms
expressed as the product of time derivative of fluid pressure and saturation-pressure
derivatives. The formulations reported by Faust (1985) and Forsyth (1988) have oil
pressure and water saturation as the dependent variables. Forsyth and Sammon (1986)
presented a three-phase flow model with non-zero gas phase pressure gradients
employing oil pressure, oil saturation and water saturation as the primary unknowns
and Thomas and Thurnau (1983) reported a similar three-phase reservoir simulation
model in terms of oil pressure, water saturation and gas saturation.
With highly nonlinear single phase flow and transport problems, it has been shown that
upstream weighting techniques can improve nonlinear convergence behavior especially
-------
for problems with hyperbolic tendencies (Huyakorn and Nilkuha, 1979; Kaluarachchi
and Parker, 1989). Newton-Raphson iteration methods also offer greater robustness for
highly nonlinear problems, but at a cost of considerable additional computation
(Huyakorn and Finder, 1978; Huyakorn and Finder, 1983). The influence coefficient
method of evaluating element matrices has been employed in single phase finite element
simulations to improve computational efficiency (Huyakorn et al., 1984; Kaluarachchi
and Parker, 1987), but has not been applied to multiphase flow models.
Any numerical solution to flow and transport equations should have the capability of
conserving mass in the system under the applied boundary conditions. Most numerical
studies of multiphase flow simulations have not addressed this issue. Abriola (1984)
investigated mass balance accuracy for a finite difference simulator during continuous
NAPL infiltration only. Long term mass conservation behavior for pulse injection
scenarios, which are expected to be more troublesome, has not been investigated.
Difficulties associated with mass conservation require careful attention to the handling
of saturation-pressure derivatives in the problem formulation.
In most practical problems, large changes in fluid pressures and saturations do not occur
throughout the spatial domain at a given time step. Computational effort is thereby
inefficiently spent solving equations in areas where little activity occurs rather than
concentrating effort in the more active zones. However, since the locations of active
zones change with time, schemes designed to take advantage of these variations must be
capable of automatically adapting to the current conditions. One method for doing so
involves adaptive grid refinement in which the numerical mesh is automatically refined
in active zones and expanded in inactive zones. Another general approach which
employs a fixed grid varies the solution approach within the domain to gain efficiency.
Adaptive implicit methods have been described by Thomas and Thurnau (1983),
Forsyth and Sammon (1986) and Forsyth (1988) that involve solving implicit
formulations of the governing equations in parts of the domain where changes in the
primary variables are large, while in the rest of the domain a computationally less
intensive but less robust implicit pressure-explicit saturation (IMPES) formulation is
employed.
In this chapter, we will discuss various formulations and techniques for modeling
multiphase flow problems to evaluate their efficiency and accuracy.
-------
2.2 Development of Implicit Pressure- Pressure Formulation
2.2.1 Theory and mathematical formulation
Governing equations for multiphase flow. In the present analysis, we assume pressure
gradients in the gas phase to be negligible so that gas pressure remains effectively
constant at atmospheric pressure. Assuming also, an incompressible porous medium
and constant fluid properties, the flow equations for water (w) and organic liquid phase
(o) may be written in summation convention for a two-dimensional Cartesian domain as
h \\ dhw dh ,
dh dh
where i, are the Cartesian spatial coordinates (t=l, 2), Kwij and Koij are conductivity
tensors for water and oil (i.e., NAPL), respectively; hw and h0 are water height-
equivalent pressure heads of water and oil (Parker et. al., 1987) respectively; pro is the
ratio of oil to water density; Uj = dz/dij is a unit gravitational vector where z is
elevation, and t is the time. The terms Cow> C00, etc. are fluid capacities defined by
dSD
P>V=°>W (2-2)
where is the porosity of the medium, 5p is the saturation of phase p, and hg is the q
phase head.
Darcy velocities in water and oil phases may be written in the form
qwi = - Kwij + «, (2.3a)
and
,flk
(2.3b)
where qwi and qoi are velocity components along Cartesian coordinates i for water and
oil phases. Phase conductivites are assumed to be described by
-------
(2.4a)
o (2.4b)
where krw and kro are relative permeabilites of water and oil, respectively, rjro is the
viscosity ratio between oil and water and Ktw^ is the saturated conductivity tensor for
water. We shall assume here that the coordinate system is oriented with the
conductivity tensor, or otherwise that off-diagonal components may be disregarded, so
that Ktwij = 0 for t ? j.
Initial and boundary conditions for each phase p may be written as
hp(xi, 0) = hpl(xt) on R for t = 0 (2.5a)
hp(xi} t) = fcp2(z,., t) on S1 for t > 0 (2.5b)
9P, »,- = QPi(*i> *) °n 52 for t> 0 (2.5c)
where n, is the outward normal unit vector at the boundaries. Here (2.5a) denotes the
initial conditions described by hpl over the entire region R and (2.5b) denotes type-1
boundary conditions along the boundary segment 5j. Equation (2.5c) denotes flux-type
boundary conditions with Qpi(zi,t) representing normal boundary fluxes along boundary
segment 52 for phase p.
Upstream weighted finite element formulation. We will apply the upstream weighting
technique described by Huyakorn and Nilkuha (1979) in conjunction with Galerkin's
weighted residual method where terms involving spatial derivatives of the governing
equations are weighted using unsymmetric upstream weighting functions. Denote the
symmetric basis functions and unsymmetric upstream weighting functions as Nj and
Wj, respectively, which are given by Kaluarachchi and Parker (1989) for the case of
linear quadrilateral elements. The exact solutions for hw and h0 are assumed to be
approximated by
#/ (e, T,) hwl (t) (2.6a)
=1
and
-------
. h. (zt, t)=rNj (e, 77) hol (t) (2.6b)
where hwj and hoj are the values of hw and h0 at node 7, and e and 77 refer to local
coordinates within the element.
Using the weighted residual method based on Galerkin's principle, eq. (1) can be
written as
f W, £ fay ( ^ + u, )) dR - J AiC^ ^ tf -}
fl ' ^ j ' R R
(2.7a)
f w 9 ( K (dhw \\ , r oh0 , r r dhw . ,
J ' ~dx ( oij (~dz~ pro j)) J ' °° ~dt ~\ 'C°w ~dt dR = 0 v2-7b)
R ' ^ } ' R R
After further simplification using Green's theorem, (2.7 may be written in matrix form
as
w \ I r cx^i i o \ fc1^ /oo\
"dTJ + (-5TJ=^F-) (2-8a)
= {F0} (2.8b)
where
BWjdN^
v1- f TS dWrdN, JD
= e£ J ^ -g^- -g^ ^ (2-8c)
(2.8d)
(2.8e)
= E f CW^JV, iIA (2.8f)
C = 1
= E f C^/^J ^ (2.8g)
e = 1
-------
E°u = E f CowNjNj dR (2.8h)
e = 1 X
n I t AW r \
= EN Kwij ^ «,. d* - J W/fcl dS\ (2.8i)
* *" ^ \ ia n /
= £( J ** '" TE1"' dR' I tr"°' *)
e 1 \ D * r« /
(2.8J)
where n is the total number of elements. Here Se refers to the segment of the surface of
element e where flux-type boundary conditions prevail and Re is the element volume.
Evaluation of element matrices. An important consideration in the finite element
formulation is the handling of nonlinear terms in element integrals such that element
matrices may be computed in an efficient manner. In this work, we use the linear shape
functions described to evaluate soil properties within the elements. For example,
capacity terms are represented by
4
C/_ 4\ _ \~^ AT (~q (f) Q-\
pg \xi> l) ~ 2i B P9 \&-s&)
where
in which C|, denote Gauss point capacities and h*, and hs0 refer to Gauss point pressure
heads. In a typical approximation of nonlinear parameters by (2.9), local nodal heads
would be used. Our experience indicates that numerical instability and poor
convergence could occur with nodal heads especially with problems involving sharp
fronts or material nonhomogeneity. Present scheme is more appropriate as Gauss point
are somewhat optimal sampling points in the finite element method than the nodal
values (Zienkiewicz, 1986). Once this procedure is adopted with all other soil
properties, it is possible to evaluate the influence coefficient matrices for each integral in
(2.8c) to (2.8j). In the present analysis, we consider the formulation for linear
rectangular elements with sides parallel to the principle axes. Equation (2.8c) may be
modified as
-------
Equations (2.8d) - (2.8j) may be modified in a similar fashion. Element matrices in
(2.8a) and (2.8b) can now be written as
A», + 4 £ (0")S ^ (2.Ha)
9=1
(2.nd)
4 .=1
In (2.11g) and (2.11h), flux-type boundary conditions have been omitted and will be
treated separately. Also superscript g in each parameter corresponds to the values
evaluated at a given Gauss point g. Matrices [0*]£, [IP]', [V\eg and [G]eg in (2.11)
represent influence coefficient matrices which can be further simplified as
(2.i2a)
« (2.l2b)
10
-------
(2.12c)
where i and y represent coordinate directions, r denotes matrices involving regular
linear shape functions and u denotes those with contributions from the upstream
weights QJ and /?,-. Note that matrices [V\ are independent of upstream weights.
Influence coefficient matrices in (2.12a) to (2.12c) are described by Kaluarachchi and
Parker (1989).
Treatment of nonlinearities by Pi card method. The Picard method is one of the
commonly used schemes to solve a set of nonlinear ordinary differential equations by a
suitable finite difference approximation. With moderately nonlinear flow and transport
problems, the Picard scheme provides rapid convergence. A distinct advantage of the
Picard scheme is its simplicity and lower computational effort per iteration than more
sophisticated schemes.
The Picard iterative scheme involves approximating the resulting finite element
formulation by a finite difference approximation. First equations (2.8a) and (2.8b) can
be written as
where
(A-'} 0
0 [A°]
[!?} =
(2.13)
(2.14a)
(2.14b)
(2.14c)
Equation (2.13) can be written in finite difference form as
(2.14d)
(2.15)
11
-------
where matrices at time (k + 9} are evaluated using the heads hk + e with
(2.16)
in which k and k+1 denote the previous and current time levels, Aifc + i is the current
time step and 6 is a time weighting factor. For a fully stable numerical solution, it is
recommended to use 6=1 corresponding to the fully-implicit backward difference
scheme. Use of (2.15) combined with (2.16) to obtain the current time step head values
is straightforward.
Treatment of nonlinearities by Newton-Raphson method. The Newton-Raphson method
is often recommended for highly nonlinear problems for which the Picard scheme may
fail or provide slow convergence (Huyakorn and Finder, 1983). We may rewrite (2.8a)
and (2.8b) at matrix element level as
rwi = A0J ~ hoJ) ~ ~rwl (2.18a)
-oJ
wJ unoJ
s (ft i . Vo]] = -r'ol (2.18b)
un
where r and r + 1 are the iteration numbers at k + 1 time step and which may be
written as
££^1*31. M. (2.19a)
&tk + i ohwj dhwj v i
12
-------
l - hkoj)
(2.19b)
dh..,j
+
XV
wJ
+
E°
IJ
dhwj
(2.19c)
_
tloL
-gr -
dh
oj
(2.19d)
In (2.19a) to (2.19d), terms with superscript k refer to previous time step values where
the quantities are known. The subscript L is a dummy index denoting the appropriate
influence coefficient matrix to be selected for a particular set of / and J. Once the
recently updated values of heads are used to evaluate quantities in (2.19), the final set
of simultaneous equations to be solved by an iterative method can be given as
(2.20)
where
drwl ,r
(2.21a)
(2.21b)
13
-------
Mass lumping. For most nonlinear flow problems and especially for multiphase flow
problems, the use of a consistent mass matrix with integrals given by (2.8e) to (2.8h)
will cause instability in the solution that would produce poor convergence. To
overcome these difficulties and to obtain a more stable solution, a common approach is
to diagonalize the mass matrix -- a procedure referred to as mass lumping. The
procedure is given for (2.8e) as
NgNjdRC°WW for I=J
(2.22a)
g =
Bfi = 0 for 7 / J
(2.22b)
Similar procedures are followed for the remaining equations given by (2.8f) to (2.8h).
Flux-type boundary conditions. For flux-boundary conditions given by (2.5c), the
original finite element formulation given by (2.11g) and (2.11h) can be modified with
very little computational effort for both Picard and Newton-Raphson schemes. These
modifications in Equations (2.8i) and (2.8j) for both schemes can be written for any
given phase as
-I
2
-------
Cp'g = (1 - w) Cpkg + wCpkg+l P,q=o,w (2.24)
where C*g is the weighted capacity, superscripts k and k -f 1 denote previous and
current time levels, and w is a time- weighting coefficient varying between 0 and 1.
Mean head analytic scheme. This scheme was initially suggested and used by Osborne
and Sykes (1986) and is given as
Cp'q = Cpq (/C V) P,q=o,w (2.25)
where
with superscripts k and A -f 1 denoting previous time step and current iteration values of
pressure heads, respectively.
Modified Chord-Slope Scheme. According to the standard chord-slope scheme (Abriola,
1984), the capacity Cpq (p, q= o, w) is defined as
Due to very poor convergence under highly nonlinear conditions from preliminary
investigations with a standard chord-slope method, this method was abandoned in favor
of the modified scheme which introduces time-weighting in a manner given by (2.24)
with w = 0.5. The expression for this scheme is given by
-f Sp(hkp, &5 + >) - $(/£+', fcj) - Sf(hli, hkg] (2.27)
15
-------
Updating nodal pressure beads. During each iteration, nodal pressure heads need to be
updated for the next iteration. With the Picard method, a common procedure is to
simply employ (2.16). With highly nonlinear flow problems, schemes of this type may
not be fully effective as the heads are updated without due consideration to the
maximum convergence error for the entire mesh. The updating method introduced by
Cooley (1971) which introduces an optimal relaxation scheme, which accounts for the
maximum convergence error for the entire mesh, was used in conjunction with both the
Picard and Newton-Raphson schemes. This method is briefly described below. Let
where N is the number of nodes in the mesh
hrp
Step 1
=1 Jfc=0
Step 2
*, = (3 + Sp)/(3 + |5P|) 5P > 1
Step 3
ep+1
where eax is the maximum change in the hp allowed during any iteration.
The convergence criterion used in this study for a given phase p (= o, to) is as follows
where r+1 and r refer to current and previous iterations, v is the allowable relative
convergence error and cp is the allowable absolute convergence error.
16
-------
Automatic adjustment of time step is undertaken depending on the number of iterations
required for nonlinear convergence. If the number of iterations is less than /mtn, the
time step is increased by a factor 1 + F and if greater than /mox, the time step is
decreased by a factor 1/(1 + F). In the simulations in this paper, satisfactory results
were obtained using 7min=4, Imax=l2 and F < 0.04.
2.2.2 Numerical investigations
Several example problems will be analyzed to investigate various aspects of the
numerical model. The first example will focus attention on the effects of various
capacity computation techniques on mass balance accuracy. The second example will
evaluate the efficiency of Picard and Newton-Raphson iteration schemes combined with
upstream weighting in analyzing a highly nonhomogeneous soil profile. The last example
will address the effects of different fluid properties on the flow of NAPL in a two-
dimensional flow domain.
TABLE 1. Soil and Fluid Properties Used in Simulations
Parameter
n
a
sm
K.w
P.O
(iow
pn
nn
Example 1 Case
A B
3.25 3.25
0.05 0.05
0.0 0.0
50.0 50.0
1.8 3.0
2.25 2.5
0.8 0.8
2.0 2.0
0.40 0.40
Example 2 Layer
1 2 3
3.0 1.8 4.0
0.07 0.01 0.1
0.05 0.2 0.0
40.0 0.5 60.0
1.8 1.8 1.8
2.25 2.25 2.25
0.8 0.8 0.8
2.0 2.0 2.0
0.40 0.40 0.40
Example 3 Case
A B C D
2.1 2.1 2.1 2.1
0.007 0.007 0.007 0.007
0.02 0.02 0.02 0.02
2.1 2.1 2.1 2.1
1.8 1.8 1.8 1.8
2.25 2.25 2.25 2.25
1.2 1.2 0.8 0.8
2.0 0.5 2.0 0.5
0.43 0.43 0.43 0.43
All units are given in centimeters and hours. Isotropic conductivity is assumed with Km
17
-------
Example ii One-dimensional Infiltration and redistribution in uniform soil The
problem analyzed corresponds to a vertical soil column 100 cm long with an oil-free
initial condition in equilibrium with a water table located 75 cm below the top surface.
The simulation was performed in two stages. In stage 1, oil was allowed to infiltrate
into the system under a water equivalent oil head of 3 cm until a total of 5 cm3 cm ~2 of
oil had entered. Boundary conditions for the water phase were zero-flux at the top and
a constant head of 25 cm at the bottom during infiltration. A zero flux oil condition
was employed at the lower boundary. During stage 2, oil was permitted to redistribute
up to t=100 hours with zero-flux conditions for oil at both boundaries and conditions for
water the same as during infiltration. Soil and fluid properties used in the simulation
are given in Table 2.1 corresponding to nonhysteretic three-phase relations described in
Chapter 3.
E
U
Q.
01
D
t - 0. 09 h
Chordsi ope
Mean Head
100
0. 0 O. 2 0. 4 0. 6 O. 8 1.0
Saturat i on
Fig. 2.1 Initial water saturation distribution and total liquid saturation at the end of
infiltration for Example 1A.
18
-------
The initial condition of zero oil saturation was assigned by fixing the initial oil head at
each node to the critical oil head, hc0T, using the following expression which derives from
the form of the scaled constitutive model
(2.28)
where 0ao and 0OW are fluid-dependent scaling factors defined in Chapter 3. As pointed
out by Parker and Lenhard (1990), a jump condition in the water saturation vs. air-
water capillary head function, Sw(haw), will occur during the transition from a two-
phase air-water system to a three-phase air-oil-water system in circumstances for which
NAPL contamination diminishes the surface tension of water relative to the pristine
state. The latter condition corresponds to the criterion
I/flu, +
(2.29)
E
U
o_
01
D
20 -
40 -
60 -
80 -
100
MQon Hood
Chords1opQ
Equ ili br i um
0. O 0. 2 0. 4 0. 5 0. 8
Saturot. i on
1. 0
Fig. 2.2 Water and total liquid saturations at the end of redistribution period for
Example 1A and theoretical equilibrium distribution.
19
-------
To avoid numerical problems associated with this jump condition, a phase updating
scheme is adopted at the end of each time step to index whether the node is a two-
phase or three-phase system (i.e., oil is absent or present). Once a three-phase
condition occurs, reversion to a two-phase condition is not allowed. The criterion for a
node to change from a two- to three-phase system is that h0 > hcj. It is important to
note that in the part of the domain remaining as an air-water system, capacity and
conductivity terms related to the oil phase will become zero and the oil flow equation
solution reduces to the identity 0 = 0. To avoid this problem, we assign minimum
cutoff values to capacities and oil relative permeabilities. Minimum values used for
capacities Cwo, C00 and Cow are 10~6 cm"1 and the minimum kro is 10'6, which were
determined by trial to provide accurate results.
Table 2.2
Summary of Results for Example 1
Capacity
Computation Scheme
Time average
analytic, u> = 0
Time average
analytic, a/ = 0.5
Time average
analytic to = 1.0
Mean head analytic
Modified chord slope
Period T
I
R
I
R
I
R
I
R
I
R
CPU Time, ' s
Case A Case B
PMB*
PMB*
105
444
92
449
123
434
104
483
81
411
80
421
80
433
89
428
93
460
Oil Mass
Error
Case A
PMB*
PMBt
1.6
4.2
-0.8
-1.8
4.2
5.2
0.1
3.4
Balance
, %
Case B
-30.8
-30.6
11.6
-11.6
-20.8
-22.0
-9.2
9.8
-5.4
-6.4
*An IBM 3094 system was used for simulations.
tl, infiltration of 5 cm of oil; R, redistribution up to 100 hours.
tPMB, extremely poor mass balance (> 70%).
20
-------
To evaluate effects of the jump condition associated with transition from two to three-
phase systems, we consider two cases in this example from £!//? = 1 for Case A and
531//3 = 0.73 for Case B. Simulations were performed using the time-averaged analytic
capacity scheme with time-weighting of w = 0.0, 0.5 and 1.0, and with the mean head
analytic scheme and the modified chord-slope method. Upstream weights used in all
simulations were 0.2. The finite element mesh used in the simulation consists of 100
elements with a uniform spacing of 1 cm and the time step varied between 0.001 to 0.5
h with a time increment factor of 4%.
The duration of the infiltration stage was approximately 0.09 h for Case A and 0.1 h for
Case B. The initial water content distribution and liquid saturation distribution at the
end of the infiltration stage using the mean head analytic and modified chord-slope
schemes are illustrated in Figure 2.1 for Case A. It may be noted that the initial water
saturation at the upper surface was very low (Sw K 0.05). The water saturation
distribution at the end of infiltration was almost identical to the initial conditions.
Computed total liquid saturation distributions at the end of infiltration were almost
identical for the two schemes shown, as well as for the time-averaged analytic schemes
for w = 0.5 and 1.0, which are not shown. Although the infiltration fronts are very
sharp, no numerical difficulties were encountered owing to the effectiveness of the
upstream weighting procedure.
Saturation distributions at the end of 100 h of redistribution are illustrated in Figure 2.2
along with equilibrium distributions computed from hydrostatics for an oil volume of 5
cm subject to the imposed boundary conditions. Here, equilibrium condition is defined
as the fluid distributions at which the total head gradient of both phases with respect to
elevation approaches zero. Results for mean head and chord slope schemes are quite
comparable as were time-averaged analytic scheme results for w = 0.5 and 1.0.
Convergence toward the equilibrium results provides a verification check for the
numerical model. It is noted, however, that whereas the equilibrium distribution
indicates no oil above a depth of 43 cm the numerical model predicts an average oil
saturation remaining above this depth of 7-8% even though oil velocities at 100 h are
practically zero. This "residual oil" content predicted by the numerical model reflects
the equivalent of an oil "field capacity" due to rapidly diminishing drainage rates under
gravitational gradients alone. It is important to note that this oil is not immobile in
any strict sense and is predicted without invoking hysteresis, fluid entrapment or fluid-
dependent residual saturations. Due to ambiguity in the term residual oil, which is
21
-------
often used to imply. occluded or discontinuous oil, a more descriptive term for this
phenomenon may be nondrainable oil.
A comparison of execution times and mass balance errors for the simulations of
Example 1 are summarized in Table 2.2. The mass balance error was always higher
during redistribution than infiltration. For Case A, all of the capacity computation
schemes provided satisfactory mass balance results, except for the analytic scheme with
w = 0.0 which produced very poor results. Optimum time-weighting for the analytic
scheme appears to be between 0.5 and 1.0. Mass balance results for Case B with £!//?
^ 1 were consistently poorer than those for Case A. The modified chord-slope scheme
yielded clearly superior mass balance results with an increase in computational effort of
about 10% over the analytic schemes. The higher accuracy of the modified chord-slope
scheme when £!//? ^ 1 is attributable to the fact that it is the only method that
properly evaluates saturation derivatives during the jump condition associated with
transition from a two- to three-phase system.
Example 2l One-dimensional infiltration and redistribution jn nonhomogeneous soil. The
main objective of this example is to evaluate the efficiency of the Picard and Newton-
Raphson schemes in combination with different upstream weights under conditions of
strong nonlinearity and heterogeneity in material properties. The geometry of the
problem is the same as that of Example 1 but the soil profile is assumed to be
comprised of three layers. Layer properties are given in Table 2.1. The middle layer
represents a finer grained material than the other materials with a saturated hydraulic
conductivity contrast of two orders of magnitude. Fluid properties used are the same as
for Example 1A with £!//? = 1. Initial and boundary conditions are also identical to
those of the previous example except that redistribution time was continued to 150
hours. Simulations were carried out using both Picard and Newton-Raphson methods
each repeated for four sets of upstream weighting factors (0, 0.2, 0.5, 0.8). The mesh
and time steps used in this example is identical to that of Example 1.
Simulation results are illustrated in Figures 2.3 and 2.4 for Picard and Newton-Raphson
analyses with an upstream weighting factor of 0.5. The extreme variability with depth
in initial water contents (Figure 2.3) produces a numerically difficult problem on which
to test the capabilities of various numerical schemes. Simulated total liquid saturation
distributions at the end of infiltration exhibit very sharp fronts and are virtually
indistinguishable for the two iteration methods (Figure 2.3) as is the case at the end of
22
-------
the redistribution period (Figure 2.4a). At 150 h oil has penetrated to the water table
but liquid velocities have become so small that redistribution has virtually ceased. The
equilibrium liquid distribution shown in Figure 2.4b predicts no oil above a depth of 55
cm which is 5 cm below the clay layer. The simulations show fully half of the oil in the
system remains above this depth at 150 hours. This "nondrainable oil saturation"
averages about 12% for the upper two layers which is considerably higher than found for
the homogeneous soil in Example 1. This is attributable to the higher oil retention
within the fine layer itself and to its impeding effect on drainage from the overlying
layer.
A summary of mass balance errors and CPU times for Picard and Newton-Raphson
simulations with all upstream weighting factors are given in Table 2.3. The different
numerical schemes provide similar accuracy with respect to mass balance with larger
errors occurring during redistribution periods. For the case of zero weights, the Picard
scheme failed to converge after 12 time steps during the redistribution period. The
Newton-Raphson scheme converged, but required greater computational effort than
when upstream weighting was employed. The Picard scheme was able to converge only
after upstream weighting was increased to 0.2 and thereafter remained insensitive to the
weighting. Mass balance was affected little by upstream weights confirming that
upstream weighting serves mainly to improve convergence.
E
0
I
£
0.
01
a
20 -
40 -
60 -
BO -
100
t - 0. 122 h
Picard
° Nowton-Raphson
0.0 0. Z 0.4 0.6 0. B 1.0
Saturation
Fig. 2.3 Initial water saturation distribution and total liquid saturation at the end of
infiltration for Example 2.
23
-------
E
U
I
r
*j
CL
01
D
20 -
40 -
BO -
80 -
100
D NewtonRophson
0. O 0. 2 0. 4 0. 6 0. 8 1.0
- - - FEM
Equi1i brium
0. O O. 2 Q. 4 O. 6 0. 8 1. O
Saturat i on
Fig. 2.4 Water and total liquid saturations at the end of redistribution period for
Example 2. (Top) Comparison between Picard and Newton-Raphson
schemes. (Bottom) Comparison between numerical (150 hours) and predicted
equilibrium distributions.
24
-------
One important .aspect of these simulations was to investigate whether increasing the
weights would deteriorate the accuracy. The results in Table 2.3 do not indicate such
deterioration and show that a highly nonlinear problem of this type can be effectively
handled with a low weight of 0.2. It is evident however that computational cost can be
almost 1/3 more with the Newton-Raphson scheme while obtaining the same degree of
accuracy as the Picard method. The results indicate that the Newton-Raphson method
did not improve convergence substantially when upstream wieghting was used and took
nearly the same number of iterations per time step as the Picard scheme. The reason
for increased computational cost of the Newton-Raphson method is mainly due to
number of terms needed per node at each iteration. For the Picard scheme, 6
coefficients must be evaluated, namely krw, kro> Cww, C00, Cwo and Cow, whereas the
Newton-Raphson scheme requires 12 more terms corresponding to the partial
derivatives of the former terms with respect to both hw and h0. Computational costs
decreased slightly with increased upstream weights for the Newton-Raphson method but
increased somewhat for the Picard scheme beyond a weight of 0.2.
Table 2.3 Summary of Results for Example 2
Oil Mass
Balance Error, *
Upstream
Weighting
Factors
0.0
0.2
0.5
0.8
% CPU Time, s
Period t
1
R
1
R
1
R
1
R
P
-3.7
NC*
-2.8
-6.0
-1.6
-4.0
-1.6
-5.2
NR
-3.6
-6.5
-2.8
-4.4
-2.8
-6.4
1.6
-5.0
P
78
NC*
70
282
71
315
71
312
NR
99
485
88
413
89
414
89
405
*P, Picard scheme; NR, Newton-Raphson scheme.
tl, infiltration of 5 cm of oil; and R, redistribution up to 150 hours.
*NC, No convergence after 12 time steps during redistribution.
An IBM 3094 system was used in the simulation.
25
-------
A 5m B
E
in
Water table
E
n
23 m
Fig. 2.5 Geometry of flow domain for Example 3.
Example 3: Two-dimensional infiltration and redistribution for varying fluid properties.
The objectives of this example are to evaluate the model for a field scale problem and to
investigate the effects of fluid density and viscosity on the movement of NAPL. The
geometry of the flow domain which has an initial water table sloping from left to right
with a gradient of 2/23 is shown in Figure 2.5. Soil and fluid properties used in the
simulations are given in Table 2.1. Note that Cases A-D involve variations in the oil
density ratio (pTO = 0.8 and 1.2) and viscosity ratio (r)ro = 0.5 and 2.0). Oil was
assumed to infiltrate along segment AB in Figure 2.5 under a small head of 1.0 cm to
simulate an oil spill. The boundaries at the ends and below the water table for the
water phase correspond to type-1 conditions with water heads stipulated, while other
boundaries were zero-flux type for water. For the oil phase, all boundaries except AB
were zero-flux. For all cases a total of 7.5 m3 m"1 of oil was permitted to infiltrate,
following which redistribution was allowed up to a total simulation time of 200 days.
The finite element mesh used in the simulation consisted of 345 nodes and 308 elements.
Time step sizes vary between 0.01 to 2.5 days with a time increment factor of 4%.
Both Pi card and Newton-Raphson schemes were used with an upstream weight of 0.5 to
solve the set of nonlinear equations. Mass balance errors at the end of the simulations
for both schemes remained almost identical and varied between -2 to -6%.
Computational efficiencies however were much different between each scheme where the
26
-------
Newton-Raphson scheme required approximately 25 to 35% more computer time. The
Newton-Raphson scheme was slightly faster in convergence than the Picard scheme in
comparison to the average number of iterations taken per time step but as discussed
earlier, the former scheme requires many more computations per iteration due to larger
number of nonlinear parameters associated per node.
Cumulative oil infiltration with time during the infiltration stage is shown for the four
cases in Figure 2.6 where the infiltration times were 40.1, 10.1, 57.7 and 13.4 hours for
cases A-D, respectively. The effects of fluid density and viscosity on infiltration rate are
clearly evident. Analysis of the results indicates that the duration of the infiltration
stage is almost exactly proportional to the relative dynamic viscosity T)ro/pro, as
expected on theoretical grounds, since the oil plume does not encounter any material
inhomogeneities or the water table during the infiltration period. These results can also
be used as a verification exercise of the numerical model. Although the infiltration
times varied for the different fluids, at the end of the infiltration stage the distributions
of oil for all cases were almost identical with a maximum penetration depth about 50
cm above the water table.
10 2O 3D 40 so 60
Fig. 2.6 Cumulative infiltration versus time for Example 3A-D.
27
-------
ro
= 2-0
t = 200 d
Fig. 2.7 Oil saturation contours after 200 days for Example 3A.
t = 200 d
ro
ro
1. 2
O. 5
Fig. 2.8 Oil saturation contours after 200 days for Example 3B.
28
-------
r-o
= 0. 8
= 0. 5
t = 200 d
Fig. 2.9 Oil saturation contours after 200 days for Example 3C.
'r-o
0. 8
2. 0
0. 05
t = 200 d
Fig. 2.10 Oil saturation contours after 200 days for Example 3D.
29
-------
Figures 2.7 - 2.10 show oil saturation contours at the end of the simulations (t = 200
days). Cases A and B which involve dense NAPLs of high and low viscosity,
respectively, exhibit substantially different behavior (Figures 2.7 - 2.8). Whereas the
low viscosity NAPL has fully penetrated the aquifer and spread laterally about 14 m
along the bottom boundary, the more viscous NAPL does not reach the bottom of the
aquifer and remains confined to a soil volume 2/3 as large as that of the dense NAPL
with a center of gravity well above the original water table level. While the plume for
the low viscosity dense NAPL is skewed in the direction of water flow, the viscous
plume remains nearly symmetric and is evidently unaffected by the groundwater flow.
For the light NAPLs, limited penetration of the aquifer is observed as expected due to
the bouyancy effect (Figures 2.9-2.10). Effects of viscosity ratio are similar to those for
the dense NAPLs. The low viscosity plume has spread laterally to a much greater
extent and exhibits downgradient migration over the water table. By contrast, the
viscous plume remains much more compact, maintains a much shallower center of mass
and is virtually unaffected by the groundwater flow.
For all simulations, it was observed that NAPL velocities at the end of 200 days were
extremely small, indicating that in real field situations, very slow movement may be
observed over long periods of time. The reason for this slower movement is that as the
oil plume spreads over a larger area, oil saturations and hence oil permeabilities become
smaller.
2.3 Evaluation of Alternative Equation Formulations
In this section, we consider various alternative methods of formulating the governing
equations for flow to assess their relative merits in terms of efficiency, accuracy and
robustness. Five different formulation methods will be evaluated, which are described
in the following section.
2.S.I Formulation of governing equations
Method 1; Pressure formulation with saturation time derivatives. Assuming an
incompressible porous medium and incompressible and compositionally uniform liquids,
the governing equations for flow in a three fluid phase porous medium may be written
30
-------
as
<2-30a>
where ^ is porosity; Sp is the p-phase fluid saturation for p = o (oil) or w (water); t is
time; x, is the t-direction coordinate; Kpij is the p-phase conductivity tensor (i,j = 1,2,3)
for phase p; hp = Pp/pwg is the water height-equivalent pressure head of the p-phase
where Pp is the p-phase pressure, pw is the density of water, and g is gravitational
acceleration; prp is the specific gravity of the p-phase (prw = 1); and Uj = dz/dXj is the
unit gravitational vector where z is elevation. For brevity, we will refer to the
nonaqueous phase liquid simply as "oil." Phase conductivities are described by (2.4),
initial and boundary conditions are specified by (2.5), and Darcy fluxes are given by
(2.3).
The unknowns in Method 1 are taken to be water and oil head, hw and ft0, which will be
treated here using a finite element model with linear basis functions. Time derivatives
will be handled using a finite difference approximation. Water and oil relative
permeabilities, krw and fcro, and saturations, Sw and 50, are treated as nonlinear
functions of phase pressures which are updated iteratively for the current time step.
Details of the numerical methods will be described in a later section.
Method 2l Pressure formulation with split time derivatives. Method 2 also involves oil
and water pressure heads as unknowns, but in this case saturation time derivatives are
split into two terms each written as pressure time derivatives. Noting that Sp = f(hw,
h0) and expanding the time dervatives of (2.30) yields
- dh* _L n.... ^o _ A. I v \ dhw f tt A (2.31a)
where pg-phase fluid capacities are defined by
31
-------
. Cpg=t ^ p, q =o, w (2.31c)
9
Initial and boundary conditions are treated exactly as in Method 1 as given by (2.31).
In Method 2, water and oil relative permeabilities, krw and kro, and fluid capacities, CM,
are treated as nonlinear functions of phase pressures which are updated iteratively.
Method 2l Capillary pressure formulation. In Method 3, fluid pressure heads are
replaced by capillary heads defined for fluid pair pq by hpg = hp - hg. Substituting air-
oil and oil-water capillary heads into (2.30) and noting that in the present analysis
dhjdii = 0 by assumption yields
d
d
0 tf \ ., -oo 1 1 /0 ookl
~dt = H; K0ij 1 Wi ~ -fir J (2.32b)
In this method, air-oil and oil-water capillary heads are taken as the primary variables
and phase saturations and relative permeabilities are treated as nonlinear functions of
the capillary heads. Initial and boundary conditions for the model are defined by
hao(*i> °) = fc-oifo) on R at i = 0 (2.33a)
>W(*,, 0) = howl(Xi) on R at t = 0 (2.33b)
hao(z{, t) = ha^(xi}t) on Si for t> 0 (2.33c)
**.(«» *) = hortfaJ) on S2 for <> 0 (2.33d)
0 (2.33e)
s n
q0. n, = q0n on S4 for t> 0 (2.33f)
where (2.33a) and (2.33b) denote initial conditions, (2.33c) and (2.33d) indicate
prescribed capillary head boundary conditions, and (2.33e) and (2.33f) specify flux type
boundaries.
Method 4 Semi-diffusive oil pressure-water saturation formulation. Diffusive forms of
mass flow equations may be written by expanding pressure derivative terms into
saturation gradients and pressure-saturation derivatives yielding homogeneous equations
in terms of fluid saturation. Here, we consider a "semi-diffusive" two-phase model
which is formulated in terms of water saturation and oil pressure. The basic form of the
32
-------
governing equations in this formulation is the same as used in Method 1. However, the
unknowns are now taken to be water saturation, Sw, and oil pressure head, h0. Other
quantities in the equations are treated as nonlinear functions of the primary variables.
Water relative permeability is expressed as a function of water saturation, and oil
relative permeability is expressed as a function of water saturation and oil head via the
constitutive relations described in Chapter 3. Oil saturation is computed from the
equality: S0=St-Sw where total liquid saturation, 5,, is computed as a function of oil
pressure head. Initial and boundary conditions for Method 4 are prescribed by
Sw(xit 0) = Swl(xi) on R at t = 0 (2.34a)
h0(xi, 0) = hol(zi) on R at t = 0 (2.34b)
SJz,, t) = 5w2(x,-,t) on Si for t > 0 (2.34c)
h0(Xi, t) = hM(xi,t) on S2 for t> 0 (2.34d)
q«,. n,- = qu, on S3 for t> 0 (2.34e)
i n
q0. n, = q0 on S4 for t > 0 (2.34f)
i n
where (2.34a) and (2.34b) denote initial conditions, (2.34c) and (2.34d) indicate type-1
boundary conditions, and (2.34e) and (2.34f) specify flux-type boundaries.
Method & Semi-diffusive water pressure-liquid saturation formulation. Method 5 is a
semi-diffusive scheme in terms of total liquid saturation, St=S0 + Sw, and water pressure
head, hw. The governing equations for this method are identical to those in Method 1.
The transformation of the governing equations to obtain St and hw as the primary
variables is achieved by treating all other variables in the equations as nonlinear
functions of the primary variables. Oil pressure head, hol is expressed as a function of
St via the constitutive model described in Chapter 3; and Sw is expressed as a function
of hw and St. Water and the oil relative permeabilities, krw and kro, are expressed in
terms of hw and 5, also via the constitutive relations. Initial and boundary conditions
for Method 5 are specified as
Stfa, 0) = 5tl(z.) on R at t = 0 (2.35a)
MX,-, 0) = hwl(Xi) on R at t = 0 (2.35b)
St(xi} t) = St2(Xi,t) on Sj for t> 0 (2.35c)
hjii, t) = hrffat) on S2 for t> 0 (2.35d)
q«,. n, = q^ on S3 for t > 0 (2.35e)
q0. n,- = q0n on S4 for t > 0 (2.35f)
33
-------
where (2.35a) and (2.35b) denote initial conditions, (2.35c) and (2.35d) indicate type-1
boundary conditions, and (2.35e) and (2.35f) specify flux-type boundaries.
2.3.2 Solution procedure
The governing equations for each of the above formulations have been solved using a
finite element method. Spatial derivative terms in the flow equations are weighted by
unsymmetric upstream- weighting functions developed by Huyakorn and Nilkuha (1979),
while the remaining terms are handled using linear basis functions. In this paper, we
shall restrict ourselves to linear quadrilateral elements. In the interest of brevity,
details of the numerical methods will be given for Method 1 only. Computational
procedures employed for the other equation formulations were similar.
Finite element formulation. The exact solutions for the liquid pressure heads in (2.30)
over an element are approximated by trial functions of the form
4
(2-36a)
7=1
4
*>(***)= 2>Xf,*KX*) (2-36b)
7=1
where hwj and h0j are pressure heads for water and oil, respectively, at node 7; Nj is the
linear basis function for node 7; and £ and r; are the local coordinates. Applying Green's
theorem and the standard Galerkin principle to (2.30) with upstream-weighting for
spatial terms yields
R R
dR + lqnwNjdS=0 /, 7=1, 2, ...N (2.37a)
R ' S
34
-------
+ f NINr7T dR + Uno NrdS=Q I, J=l, 2, ... N (2.37b)
J at J
R S
where Wj is an upstream-weighted shape function described by Kaluarachchi and
Parker (1989) and N is the number of nodes in the domain. Since only diagonal terms
in the conductivity tensor are considered in the present analysis, the double subscript
notation is henceforth dropped.
As in Section 2.2, we expand the matrices in the first two integrals of (2.37) into four
submatrices by permitting the hydraulic conductivity within an element to vary
through the upstream-weighting function such that
4
KP.(*i, ') = E wfal) KP (2-38)
7=1
where Kp1 denotes the hydraulic conductivity attributed to node /. Previously, we
approximated Kp as values corresponding to 2x2 Gauss points nearest each node. In
the present analysis, we take Kp values corresponding exactly to node points.
Substituting (2.38) into (2.37) and integrating gives
JO
4./J h*J + BIJ -dr ~ F»I = ° [2-39a)
' Fol = 0 (2-39b)
where
(2.40a)
= E f f E [^]C/ < + 4 E IP?/ K1 }
e=l \ /=! g=l I
= E f f E I^]e/ ^x + 4 E [Uf, Ai ) (2.40b)
e=l \ 7=1 7=1 y /
35
-------
%=E ₯ (2-40c)
e=l 4
n , 4
FU>I= ~ 5Z ^ 52 \G\e]KL + boundary flux term (2.40d)
Si 27=i
n , 4
FoI=-J2 f ProU [£]7-Ko -I- boundary flux term (2.40e)
e=l ^ 7=1 tf
where m and d are the width and length of the rectangular element, respectively.
Influence coefficient matrices, [U*]f, [Uv]f, [G\je and coefficient vectors for boundary
flux terms are given by Kaluarachchi and Parker (1989). In the subsequent analysis,
the boundary flux terms in (2.40) have been omitted for simplicity. A mass lumping
procedure is employed to diagonalize the mass storage matrix in the flow equations
which is described earlier.
Treatment of nonlinearitv. Owing to the treatment of time derivative terms in Method
1, Picard iteration cannot be used, and a Newton-Raphson iteration procedure is used
for comparison for all five methods. In residual form, (2.40) may be written as
roj= AoIJhoJ+ BIJ -$ - FOI (2.41b)
Applying the Newton-Raphson procedure gives
^o/+1- hof) = - r./ (2.42a)
+*- h0f) = - r.f (2.42b)
where r and r-fl refer to the previous and current iterations, respectively. The
derivatives in (2.42) can be expressed as
36
-------
, 9AwIL dSu.j dFwI
= WL ~ + [ n ]/ ~ (2'43b)
, 9AoIL
=
IJ
Q O
n
where Ai =
-------
Table 2.4 Soil and fluid properties used in simulations.
Parameter
n
a, cm"1
5m
KSWI cm ^l"1
Pao
Pow
Pro
r)ro
*
1
3.00
0.05
0.00
20.00
3.20
1.455
O.SO
O.SO
0.40
Example
2
3.00
0.05
0.00
20.00
3.20
1.455
1.20
O.SO
0.40
3
3.00
0.05
0.00
50.00
3.20
1.455
O.SO
2.00
0.40
Table 2.5 Total time steps during infiltration.
Example 1
1 104
2 95
3 100
2
107
9S
111
Method
3
104
95
100
4
NC
NC
101
5
103
94
99
NC = no convergence
38
-------
2.3.3 Numerical results
Three example problems involving infiltration and redistribution of organic liquids in
partially water saturated soils will be investigated to compare the computational
efficiency, mass balance error and robustness of the different equation formulations.
Assumed fluid and soil parameters in the constitutive model for permeability-
saturations-capillary pressure relations are given in Table 2.4 for the three problems.
Example 1^. This problem involves vertical infiltration and redistribution of light oil in
a 200 cm long column initially free of organic liquid and in equilibrium with a water
table at a depth of 100 cm. For Methods 1, 2 and 4, oil infiltration was permitted
under zero oil pressure, while in Method 3, air-oil capillary pressure was fixed at zero.
In Method 5, the top boundary was simulated by maintaining the total liquid
saturation, St at unity. In all simulations, infiltration was continued until a total
volume of 3.95 cm3 cm"2 was absorbed. Following the infiltration period, fluid
redistribution was allowed with zero oil or water flux at the upper boundary.
At the lower boundary, zero oil flux and a fixed water head of -1-100 cm was assumed.
This was achieved directly for Methods 1, 2 and 5 by specifying zero oil flux and fixing
hw at 100 cm. For Methods 3 and 4, there is no way to directly specify the desired
boundary conditions so surrogate approaches were employed. For Method 3, the
capillary pressures were fixed at the lower boundary in such a way as to impose the
correct water head and to force zero oil saturation. The air-oil capillary head was
assigned a value hao ha - hcj where ha = 0 (atmospheric gas pressure) and hcj is the
oil head at which S0 -» 0, which for the soil and fluid properties in this problem is equal
to 31.25 cm. Thus, hao = -31.25 cm at the lower boundary and the oil-water capillary
head, how = h0 - hw, is 100 - 31.25 = 68.75 cm. For Method 4, with dependent
variables 5^ and A0, the lower boundary conditions were specified by setting Sw = 1 and
h0 = h". This effectively stipulates zero oil saturation at the bottom, but does not fully
specify the water head as will be discussed later.
Spatial discretization in all cases was achieved using 42 elements with 86 nodes.
Simulations were started with an initial time step of 0.0002 h at the beginning of
infiltration and and 0.0001 h at the beginning of redistribution using a time-step
acceleration factor, F, of 1.03 and a maximum time-step of 0.1 h.
39
-------
200.CO
150.00 -
E
o
-^ 100.00 -
'£ 3
50.00 ^
iquid
Water
O.CO ' !.
0.00 0.25 0.50 0.75 1.07-
Saturation
Fig. 2.11 Fluid saturation distribution at the end of infiltration for Example 1.
200.00
150.00 -
E
o
tf 100.00 H
CD
IE
50.00 -
0.00
N Joto! liquid
\
/
Water
0.00 0.25 0.50 0.75 1.03
Saturation
Fig. 2.12 Fluid saturation distribution at the end of redistribution for Example 1.
40
-------
Table 2.6 Toted iterations during infiltration.
Example
1
2
3
1
450
418
476
2
2S9
269
324
Method
3
456
423
476
4
NC
NC
298
5
299
282
324
NC = no convergence
Table 2.7 CPU seconds using IBM 3090 during infliltration.
Example
1
2
3
1
30
28
32
2
23
22
27
Method
3
30
29
33
4
NC
NC
25
5
26
25
29
NC = no convergence
41
-------
Table 2.8 Total time steps during redistribution.
Example
1
2
3
Method
12345
2,201 2,186 2,202 NC . MB
2,202 2,186 2,208 NC MB
2,188 NC 2,187 2,187 2,187
NC = no convergence achieved
MB = mass balance error exceeded 10% and the solution was halted
Table 2.9 Total iterations during redistribution.
Method
12345
1
2
3
3,218 2,922 3,483 NC MB
3,273 3,192 4,398 NC MB
4199 NC 4191 3619 4514
NC = no convergence achieved
MB = mass balance error exceeded 10% and the solution was halted
42
-------
Table 2.10 CPU seconds using IBM 3090 for redistribution period.
Example
1
2
3
Method
1 2 3
211 222 222
215 249 279
270 MB 274
4 5
NO 210
NC MB
261 331
NC = no convergence achieved during infiltration problem
MB = mass balance error exceeded 10% and solution halted
Table 2.11 Maximum percent mass balance error for oil.
Example
1
2
3
1
0.50
0.34
0.75
2
1.87
3.70
(d)
Method
3
0.47
(b)
0.71
4
NC
NC
0.10
5
(a)
(c)
0.07
NC = no convergence during infiltration
(a) Mass balance error of 0.6% prior to oil front reaching water table;
thereafter mass balance error became greater than 10%
(c) Mass balance error of 0.1% prior to oil front reaching water table;
thereafter mass balance error became greater than 10%
(b) Mass balance error of 0.34% prior to oil front reaching lower boundary;
thereafter oil mass lost by flow through boundary
(d) Mass balance error gradually increasing to greater than 10%
43
-------
Oil saturation distributions at the end of oil infiltration (0.095 h) and 200 h after the
start of redistribution are shown in Figures 2.11 and 2.12, respectively, for Method 1.
Results for other methods except Method 4 were visually identical. A comparison of
the number of time-steps and iterations and computational time during the infiltration
period are given in Tables 2.5 - 2.7 for all methods. Similar information, along with the
maximum mass balance error during the redistribution period, are compared in Tables
2.8-2.11.
Methods 1 and 3 exhibited similar convergence behavior, with an average of about 4.3
iterations/time-step during infiltration and 1.5 iterations/time-step during redistribution
(Tables 2.5 - 2.9), and yielded maximum mass balance errors of about 0.5% (Table
2.11). Method 5 exhibited more rapid convergence during infiltration (2.7
iterations/time-step) but similar behavior during redistribution to Methods 1 and 3. In
contrast, mass balance error was markedly higher at almost 2%. Total CPU time for
infiltration plus redistribution periods (Tables 2.7 and 2.10) within a few percent for all
three methods.
Method 5, the semi-diffusive scheme in terms of St and hw, performed well until the oil
front reached the water table, but produced large mass balance errors thereafter. To
permit comparison with other methods prior to this occurrence, Tables 2.8-2.11 indicate
results for the entire redistribution period of 200 h. The difficulty encountered by-
Method 2.5 may be due to nonlinearities caused by the boundary conditions. Fixing
water head (^=100 cm) and total liquid saturation (St=l) does not uniquely specify oil
pressure which may take on values h0 > 0. We conjecture that for all nodes with 5f=l
this lack of constraint on h0, which is not one of the dependent variables in the
formulation, leads to the observed numerical errors after the oil front reaches the water
table.
Method 4, the semi-diffusive scheme in terms of Sw and /i0, experienced serious
convergence difficulty during the infiltration period. The boundary conditions imposed
at the bottom nodes of the column for Sw and h0 does not guarantee that hw=WQ cm
will be imposed on these nodes during the current iteration. When hw > 0 and 5^=1,
the constitutive model gives hw > h0. We conjecture that this lack of constraint on hw,
which is effectively treated as a nonlinear coefficient of the primary variables, causes
instability at nodes below the water table.
44
-------
Example 2.. This case is the same as the previous example except that the organic
liquid is denser than water (Table 2.4). Initial and boundary conditions imposed on this
problem are taken identical to those of Example 1. Because of the high fluid density,
oil will penetrate through the capillary fringe and eventually reach the lower boundary
of the system. For Methods 1, 2 and 5, the zero oil-flux condition will force oil to
accumulate at the boundary while water outflow is permitted under constant water
head. However, the lower boundary conditions employed for Methods 3 and 4 have the
effect of forcing zero oil saturation at the lower boundary which will result in oil loss
from the lower boundary. This problem cannot be avoided as there is no way to
directly impose the desired boundary conditions for the latter methods. We will take
these problems into consideration in the following discussion.
200.00
\
Total liquid
150.00 -
i
j
100.00 -
CD
50.00 -
o.oo
Water
0.00 0.25 0.50 0.75 1.00
Saturation
Fig. 2.13 Fluid saturation distribution at the end of infiltration for Example 2.
45
-------
200.00 TTT
150.00 -
CJ
100.00 H
0
50.00 -^
0.00 -r
lotol liquid
/
Water
0.00 0.25 0.50 0.75 1.00
Saturation
Fig. 2.14 Fluid saturation distribution at the end of redistribution for Example 2.
Saturation distributions at the end of the infiltration period (0.07 h) and 200 h after the
start of redistribution are shown in Figures 2.13 and 2.14, respectively, for Method 1.
Oil was predicted to reach the lower boundary at about 120 h. The total number of
time steps and iterations, CPU time and the maximum mass balance error for all five
methods are given in Tables 2.5 - 2.11 for the infiltration and redistribution periods.
In Example 2 as in Example 1, Methods 1 and 2 both performed well, but the former
yielded much lower mass balance error (0.34% vs. 3.7%) and required 16% less
computer time. Method 3 yielded mass balance accuracy similar to Method 1 up until
the time oil reached the lower boundary, after which the imposed boundary conditions
permitted oil drainage from the lower surface. Method 4 failed to converge from the
start of the simulation for inferred reasons discussed for Example 1. Method 5 behaved
well up until the time the oil front reached the water table. Thereafter, the mass
balance error increased from only about 0.1% to in excess of 10% within 50 h.
46
-------
Example 2i In this example, we investigate a situation involving a domain in which
water pressure is everywhere negative. Column dimensions, spatial discretization and
time-integration parameters are the same as in Example 1. Fluid and porous media
properties are similar to those of Example 1, except for a higher organic liquid viscosity
and higher hydraulic conductivity (Table 2.4). Initial conditions correspond in this case
to equilibrium with a water table at the bottom of the column (200 cm depth) and the
water head is fixed throughout the simulation at the lower boundary at hw=0 for
Methods 1, 2 and 5. For Methods 3 and 4, the lower boundary conditions are as
described for Example 1 corresponding to the stipulation of zero oil saturation. The
results of Method 1 at the end of infiltration and at 200 h after the start of
redistribution are shown in Figures 2.15 and 2.16, respectively. The total number of
time-steps and iterations, CPU time and the maximum mass balance error are given in
Tables 2.5-2.11.
For the present example, Method 1 yielded results with somewhat greater but still
acceptable mass balance error (0.75%) compared with the previous problems reflecting a
higher degree of nonlinearity in this problem associated perhaps with the lower initial
water content and/or greater mobility contrast between oil and water. Method 2
succumbed entirely to these difficulties and eventually failed to converge after the mass
balance during redistribtion grew above 10%. Method 3 yielded nearly identical results
in terms of accuracy and computational effort to Method 1. Methods 4 and 5, which
behaved miserably for Examples 1 and 2, now exhibit remarkably good results with
mass balance errors 7-10 times lower than Method 1. Method 4 yields the most efficient
solution for this problem both during infiltration and redistribution, while Method 5 was
somewhat more computationally intensive.
2.4 Development of an Adaptive Solution Method
In most practical problems, large changes in fluid pressures and saturations do not occur
throughout the spatial domain at a given time step. Computational effort is thereby
inefficiently spent solving equations in areas where little activity occurs rather than
concentrating effort in the more active zones. We discuss in this section a new
"adaptive solution domain" (ASD) algorithm designed to reduce unnecessary
computational effort.
47
-------
200.00 -r
j Water
0.00 | I I : I : I I I I | I I i , I I I II I I M I I II I I | i I i , I I i I I
0.00 0.25 0.50 0.75 1.00
Saturation
Figure 2.15. Fluid saturation distributions at the end of infiltration foi Example 3.
200.00 -TT
T
\
\
150.00 H
E
o
100.00 -
CD
50.00 -
\ _ , , .. . .
" 'oial liquid
0.00 I I I I i I I I I | . I i i i I I I i i i . I | I I M I M I I
0.00 0.25 0.50 0.75 1.00
Saturation
Figure 2.16. Fluid saturation distributions at the end of redistribution for Example 3.
48
-------
2-4-1 Description of ASD algorithm
In the ASD method, elements are categorized as either "active" or "inactive" at a given
iteration. If the element is classified as active, then it is included in the global matrix
assembly, otherwise it is excluded. The criteria for changing a node from inactive to
active is based on the changes in water pressure head, hw, oil pressure head, hoi oil
saturation, S0, and water saturation, Sw at connected nodes from the last converged
time step to the current iteration. The criteria for switching an active element to an
inactive status is based on the difference in hw, h0 Sw and 50 at the connected nodes
between the converged values at the previous and the current time step. The specific
criteria for changing the classification of an element is as follows:
Inactive -» active if any nodes meet the conditions:
or
or
or
or
Ahw
Ah0 > Tj
ASW > r2
AS0 > T,
where Ahw, Ah0> ASW and A50 are the changes in hw, /i0, Sw and 50, respectively,
between the previous converged time step and the current iteration.
Active -» inactive if all nodes meet the conditions:
Ahw
and
Ah0
and
and
where Ahw, Ah0i ASW and A50 are the changes in hw, h0, Sw and 50 refer to the changes
between the converged values at the end of the current and the previous time step, and
£ hw is the cumulative change in hw from the last time the node met the criterion to
the present time.
To ensure an accurate solution, a bias is imposed towards maintaining active elements
by making T4 < rl and rs < r2. Preliminary simulations for a number of cases indicated
49
-------
good results using r3=2 cm and r4=r1/lQO subject to the constraint that
0.001 < r4 < 0.0001 cm and r5=r2/10 subject to the constraint that 0.0001 < rs < 0.00001.
The foregoing values will be employed in the present study. To avoid erratic switching
between active and inactive subdomains, an active element is labelled inactive only
after convergence conditions are satisfied for three consecutive time-steps. This allows a
smooth transition between active and inactive elements while achieving high
computational efficiency.
It may be noted that complete elimination of a node from the global matrix only occurs
if all elements associated with the node are classified as inactive. Equations for nodes
on the active-inactive frontier will have contributions only from active elements, and
the boundary will shift naturally as imbibition or drainage "fronts" pass through the
domain.
2.4-2 Numerical Results
Two examples will be described to illustrate the ASD method and to compare its
performance with a conventional full domain finite element approach. The first
example involves infiltration and redistribution of oil in a one-dimensional soil column
and the second example involves a two-dimensional domain.
Example 1^. This problem involves vertical oil infiltration and redistribution into a 225
cm long column initially free of hydrocarbon and in equilibrium with a water table at a
depth of 200 cm. A dense oil was added at the top of the column under a constant head
of ha = 0 over a period of 0.14 h until the cumulative oil infiltration reached 6.92 cm.
At the end of the infiltration period, the upper boundary condition for oil was specified
as zero flux. Fluid redistribution was permitted for a period of 250 h after the end of
the infiltration stage. A zero flux condition for water was imposed at the upper
boundary at all times. At the lower boundary, water pressure was maintained
continuously at hw = 25 cm, and zero oil flux was imposed. Lateral boundaries were
zero flux for both fluids. Spatial discretization for the problem was achieved using 45
rectangular elements each 5 cm wide and spaced 5 cm apart in the vertical direction
yielding a total of 92 nodes. All simulations were started with an initial time-step of
0.0002 h, and a time acceleration factor of F = 1.03 was used to increase the time-step
size up to a maximum time-step of 0.5 h. Soil and fluid properties are given in Table
2.12 which correspond to the constitutive model described in Chapter 3.
50
-------
Table 2.12 Soil and fluid properties for example problems.
Parameter
4>
n
Q, cm
Sm
Ksw, cm h~
Pao
ft ow
Pro
r]ro
1
0.4
3.0
0.05
0.0
20.0
3.2
1.455
1.2
O.S
Example
2
0.4
2.7
0.02
0.0
2.1
3.2
1.455
O.S
2.0
Table 2.13 CPU time (seconds) and mass balance error on IBM 3090 for Example 1.
*
TI TI CPU time (s) MBE
(cm) (-) infiltration redistribution (%)
Conventional FEM 30 127
0.002 0.00005 9 92
0.010 0.00010 6 77
0.050 0.00050 6 76
0.100 0.00050 6 73
0.33
0.66
0.37
0.43
0.46
51
-------
200.00 -
ASD o
CFtM
Total liquid
o.oo
o.oo
0.25 0.50
Saturation
i.oo
Figure 2.17 Water and total liquid saturation distributions at the end of redistribution
for CFEM and ASD method with ^=0.01 cm for Example 1.
Simulations were carried out using the ASD method with various values of r^ and r2
and using a conventional full-domain finite element solution (CFEM) based on the same
numerical formulation employed for the ASD solution but with no element exclusion.
Simulations were performed in two stages corresponding to the period of oil infiltration
and the subsequent redistribution period. Table 2.13 shows the CPU time required for
infiltration and redistribution periods using different ASD tolerances and for the CFEM.
The maximum percentage mass balance error for the oil phase (MBE) during the
redistribution period is also given in Table 2.13. Water and total liquid saturation
distributions predicted by the ASD method for the r^= 0.01 cm case and the CFEM are
virtually identical (Figure 2.17) and the other ASD results were visually indiscernible.
Small differences among the various solutions were evident in the mass balance error
(Table 2.12). The ASD cases with rl > 0.01 cm yielded nearly the same mass balance
error as the CFEM, while the ASD case with rl = 0.002 cm exhibited substantially
higher, although still tolerable, mass balance error. We attribute the poorer behavior of
52
-------
the latter case to erratic changes in the number of active elements for consecutive
iterations of the same time-step which result when imposed ASD tolerances undercut
the precision in the numerical solution. At this point, the ASD criteria are simply
driven by noise in the solution.
Marked savings in CPU time are obtained with the ASD method (Table 2.13). For the
TI > 0.01 cm cases, the ASD method yields an overall speed-up factor (CPU time for
CFEM/CPU time for ASD method) of 5 relative to the CFEM for the infiltration
period. For the entire redistribution period, the overall speed-up factor diminishes to
1.7. The large speed-up factor for the infiltration period reflects the fact that only a
small fraction of the total elements are active early in the simulation. As redistribution
proceeds, the number of active elements increases as the oil front advances down the
column.
ci I
II-
6.3S:
E
0.00
I I
2.55
7.5
X, m
Figure 2.18 Geometry of Example 2.
53
-------
Example 2^ A two-dimensional domain is considered in this example involving a
vertical slice through an unconfmed aquifer which is 7.50 m wide and 6.38 m in height
(Figure 2.18). The system initially contains no oil and is in equilibrium with a water
table at an elevation of 2.55 m from the bottom boundary. A low density oil is added on
a 1.1 m strip (A-B) at a constant pressure of h0 -0.1 m for a period of 3.43 d yielding
a cumulative oil infiltration of 0.32 m3 m"1. Thereafter, A-B is a zero-flux boundary for
oil. All other boundaries are zero-flux for oil at all times. Redistribution was permitted
out for a period of 45 days following the end of infiltration. Water pressure is fixed at
the initial condition on C-D and all other boundaries are zero-flux for water at all times.
A total of 294 rectangular elements was used to discretize the domain yielding 330
nodes. All simulations were started with an initial time step of 0.002 d, and a factor of
F=1.03 was used to increase the time step size up to a maximum time-step of 0.2 d.
Information on soil and fluid properties for the problem is given in Table 2.12.
Simulations were carried out with the ASD method using different tolerances and using
CFEM. Comparisons of water and total liquid saturation distributions on vertical
sections at x = 0 and 1.2 m are shown in Figures 2.19 and 2.20. It is observed that at z
= 0, the CFEM results slightly undershoot the ASD solution and slightly overshoot at x
= 1.2 m. This reflects a slight oscillation in the CFEM solution which could be
observed in the y-direction. Such behavior was not evident in the ASD results. The
superior stability of the ASD method may also be evidenced by the more rapid
convergence behavior of the nonlinear iterations which averaged 4.0 iterations per time-
step for the CFEM and 3.5 for the ASD method for the TI = 0.1 cm case. The greater
accuracy of the ASD method is confirmed by the consistently lower mass balance error
for all ASD solutions (Table 2.14). As in Example 1, the lowest mass balance error for
the ASD method was obtained for the case with rl = 0.1 cm with smaller and higher
tolerances resulting in somewhat greater error.
Total CPU times for infiltration and redistribution periods showed even more marked
contrasts between CFEM and ASD results than were observed for the one-dimensional
problem. Overall speed-up factors for the present case were 6 for the infiltration period
and 2.5 for the redistribution period. This greater speed-up reflects the larger
proportion of the domain which can on average be excluded from the active domain as
dimensionality increases. As in the first example, as redistribution proceeds, the fraction
of active elements gradually increases. At sufficiently long times, the active elements
will diminish as steady state is approached.
54
-------
Table 2.14 CPU time (seconds) and mass balance error on IBM 3090 for Example 2.
T\
(cm)
T2
(-)
Conventional FEM
0.05
0.10
0.15
0.0001
0.0001
0.0001
CPU
infiltration
511
86
84
84
time (s)
redistribution
871
371
349
323
MBE
(%)
0.68
0.39
0.12
0.14
600.00 -E
500.00 i
400.00 i
J= 300.00 ^
200.00 I
100.00 -E
CD
0.00
Water
ASD
CFEM
Total liauid
x=O.Cb
0.00
0.20
0.40
0.60
0.80
1.00
Saturation
Figure 2.19 Water and total liquid saturation distributions at the end of redistribution
at z=0 m for CFEM and ASD method with Ta=0.1 cm for Example 2.
55
-------
CD
600.00 -E
500.00 '-.
400.00 \
300.00 \
200.00 \
100.00 \
o.oo
Water
ASD o
C'EX
o.oo' '0.20
otal liquid
0.40 0.60 0.80
Saturation
1.00
Figure 2.20 Water and total liquid saturation distributions at the end of redistribution
at i=1.2 m for CFEM and ASD method with Tj=0.1 cm for Example 2.
56
-------
3. HYSTERESIS IN THREE-PHASE FLOW RELATIONS
3.1 Evaluation of General Hysteretic Formulation
S.I.I Model description
In order to model three-phase flow, functional relationships between fluid pressures,
saturations and permeabilities must be known. In most cases, the NAPL is less wetting
to the solid than water, so that water will occupy the pore space in immediate contact
with the solid, although oil-wet or mixed oil-water-wet systems may occur for
example, in soils with high organic matter contents or in cases where mineral surfaces
exhibit natural organic coatings. We will assume in the following that water (w) is the
wetting phase, organic liquid (o) is the intermediate wetting phase, and air (a) is the
nonwetting phase. Hysteresis will occur in three-phase systems due to nonwetting fluid
entrapment, contact angle hysteresis, ink bottle effects and other phenomena.
Fluid entrapment. To model hysteresis in three-phase air-NAPL-water systems, the
concept of apparent saturation may be used. Apparent water saturation in a two-phase
air-water system is defined by
IT = sr + 5 (3.1)
where 5^"" is the effective water saturation and Sataw is the effective trapped air
saturation of a two-phase air-water system. Effective water saturation is defined by
Sw = (Sw STW}I(\ 5ru;), and effective nonwetting phase saturations (i.e., air) have the
form 5n = Sn/(l Srw) where 5n is the actual nonwetting phase saturation.
For three fluid phase systems, apparent water and total liquid saturations are defined.
Apparent water saturation is defined as the sum of the effective water saturation and
effective trapped air and oil saturations, and apparent total liquid saturation is defined
as the sum of the effective water and oil saturations and the effective trapped air
saturation. The apparent water saturation can be used to determine the flow channel
radii where NAPL-water interfaces reside. The apparent total liquid saturation can be
used to determine the flow channel radii where air-NAPL interfaces reside. Apparent
water and total liquid saturations are given as
57
-------
Sw = Sw + Sot + Salu, (3.2a)
St = Sw + S0 + Sat (3.2b)
where Sgi is the effective trapped oil saturation (i.e., in water), Saiw is the effective
saturation of air occluded by water, and 5B< is the effective total trapped air saturation
(i.e., air that is occluded either by oil or water). The effective trapped saturations are
defined analogous to that for the two-phase case. Using apparent saturations in lieu of
effective or actual saturations to model relationships among fluid pressures, saturations,
and permeabilities has physical merit. Apparent saturations can be used to determine
fluid flow channel dimensions that separate continuous fluid phases that are immiscible
with each other. This modeling approach will yield a more accurate description of fluid
distributions in porous media for arbitrary saturation paths. Correlating fluid
distributions to flow channel radii is important when predicting relative permeabilities
from S-P relations via a pore size distribution model. Note that the model assumes the
entrapped fluids are discontinuous and immobile.
Fluids that wet other fluid surfaces can become trapped by the fluid phase they wet.
That is, fluids with lesser wettability can become occluded by fluids with greater
wettability. For example, air may become trapped by oil or water, and oil can become
trapped by water. The model assumes that water cannot become trapped by oil or air,
and that oil cannot be trapped by air. Entrapment of fluid occurs as fluid-fluid
interfaces advance into larger flow channels, which implies that one fluid phase is
displacing another.
Trapping of a fluid phase can occur only in flow channels larger than those
corresponding to the smallest or narrowest flow channels where the fluid-fluid interfaces
resided. For example, air can be trapped by advancing air-water interfaces in a two-
phase fluid system in flow channels larger than those corresponding to the historic
minimum apparent water saturation. Historic minimum saturations are the smallest
fluid contents for a given saturation path history. For two-phase systems, this value
would be associated with the lowest saturation obtained on the main water drainage S-P
branch. Note that apparent and effective saturations are equal for the main drainage 5-
P branch because this path has no trapped fluids. For all other saturation paths, the
58
-------
model assumes some trapped fluid exists and the apparent water or total liquid
saturations are greater than effective water or total liquid saturations for the same
actual water or total liquid content, respectively.
The amount of trapped air in a two-phase system, therefore, is zero at the historic
minimum apparent water saturation and is at a maximum when the porous medium is
apparently water saturated, which is at an air-water capillary pressure of zero on an
imbibition path. An algorithm proposed by Land (1968) is employed to estimate
trapped fluid saturations of two-phase fluid systems that occur at a zero capillary
pressure on imbibition saturation paths. Land's algorithm predicts the amount of
nonwetting fluid that becomes trapped in pore spaces when the wetting fluid saturation
increases from some point along the main drainage branch to an apparently saturated
condition (i.e., saturated with respect to wetting fluid and entrapped nonwetting fluid)
indicated by a capillary pressure of zero. The effective wetting fluid content at the
saturation path reversal from the main drainage branch to an imbibition scanning path
is given by the symbol Sv. Land's empirical relationship may be written as
V = 1 + R (1 J\^} (3.3a)
in which
Rv = 7TU - l (3-3b)
tr
where 5t>" and ISiji are the maximum effective trapped nonwetting fluid saturations or
the effective residual fluid saturations corresponding to the imbibition scanning path
beginning from ^S» and to the main imbibition branch, respectively. This algorithm
predicts that as A5';»0, the effective residual fluid saturation for that scanning path
(Sjji) approaches the maximum effective entrapped nonwetting fluid saturation
associated with the main imbibition path (75-r''). In other words, as more pore space
becomes occupied by a continuous nonwetting phase as wetting fluid drains, the amount
of nonwetting fluid that potentially can become trapped will be greater because of the
larger volume of flow channels that can trap nonwetting fluid. Note that S^ = ISij>
when A5^ = 0 and, by definition, A5;" = A5;-*'. Note that we use the term irreducible
saturation to refer to the minimum wetting fluid content associated with wetting fluid
drainage, while residual saturation refers to the nonwetting fluid saturation at zero
capillary pressure. The irreducible saturation thus represents a continuous fluid phase,
59
-------
albeit that only thin films may be present, and the residual saturation represents
trapped, discontinuous nonwetting fluid.
The algorithm proposed by Land (1968) places bounds on the amount of nonwetting
fluid that can become trapped in a porous medium with two immiscible fluids. That is,
at a wetting fluid saturation of A5-^, there is no trapped nonwetting fluid, and at an
apparent wetting fluid saturation of 1 (i.e., at a capillary pressure of zero), the effective
entrapped nonwetting fluid saturation is equal to SlV«. To interpolate between these two
end-points, it is assumed that all pore size classes will entrap nonwetting fluid in
proportion to their volumes. Accordingly, the amount of entrapped nonwetting fluid is
predicted to vary linearly as
_ / 5..; _ 45 .«; \
* j (3'4)
where 5,^ is given by (3.3). Note that for apparent saturations equal to A5^, the main
drainage path is assumed, and for apparent saturations greater than A5-^', scanning
saturation paths are assumed. When modeling S-P behavior, A5.v needs to be updated
every time the main drainage path is followed. Implicit in the modeling scheme is that
as interfaces between continuous nonwetting and wetting fluids advance into larger-
sized flow channels, nonwetting fluid becomes trapped (i.e., discontinuous) by the
wetting fluid as it displaces nonwetting fluid into larger pore cavities. The process is
assumed to be reversible during drainage.
Algorithms to predict nonwetting fluid entrapment in two-phase systems are employed
to estimate fluid entrapment in three-phase systems. This is because in three-phase fluid
systems the entrapment of nonwetting fluids are essentially two-phase fluid processes.
Air in a three-phase system will be trapped by advancing air-oil interfaces and oil will
be trapped by advancing oil- water interfaces. However, because the fluid system is
initially an air-water system, air entrapment by advancing air-water interfaces must
also be considered.
Distinguishing between air trapped by air-water interfaces and air trapped by air-oil
interfaces is done by tracking the historic minimum apparent total liquid saturation
(Stmin) and *SW. If 5tnin > *3W, then air is trapped by air-water interfaces in flow
channels larger than those corresponding to *SW*W but smaller than those corresponding
60
-------
to 5/71"1, and air is trapped by air-oil interfaces in flow channels larger than those
corresponding to 3tmin. The historic minimum apparent total liquid saturation (3t""B) is
used to determine the smallest flow channel radii where air-oil interfaces resided. If
§min < *Swaw, then all air is entrapped by advancing air-oil interfaces. Once a two-
phase system becomes a three-phase system, the value of *Swtw becomes fixed (i.e.,
invariant). When the interfaces between the continuous air and oil phases recede into
flow channels smaller than what air-water interfaces historically occupied (i.e.,
51T7"n
-------
and the effective entrapped air saturation occluded by water for a given saturation path
history can be predicted from the following relations:
? aw
Jw
for 5
_ _ / C mi'n Ac
C C I °< ^u
_
-1
for 5 5,, > *Swau'
"rt ~C> I *^1H ^~ *^ll
1 A C au'
I Ou
_ C tntn
" ar° V 1 _ C miti
\ J. O<
for 5^ > Stmin
/ C _ C mi
C C ( ''ui "t
^S/fim ^~ ^/I^A 1 ^^^^^^^^^^^^^^"
for 5^ > 5,"11'"
_ _ / ^min _ A^ ou; \ _ / ^ _ ^ min \
- (3-6b)
for 5, <*SW
Stttw = 0 (3.6d)
for Stmin <
for 5W < Sn
S*t* = 0 (3.6g)
62
-------
where Sat is the effective total entrapped air saturation, Saiw is the effective entrapped
= s=
air saturation in the aqueous phase, and 5, and Sw are the current apparent total liquid
and water saturations, respectively. Note that the effective saturation of air trapped in
oil (5a<0) can be obtained from Saio = Sai 5alu, for any saturation path history.
In a three-phase fluid system, displacement of oil by advancing oil-water interfaces will
result in discontinuous blobs or ganglia of oil being formed in flow channels larger than
those corresponding to the minimum apparent water saturation (Swm>n). Within the
trapped oil phase, there may be entrapped air that resulted from air-water or air-oil
interfaces. To account for concomitant entrapment of air and oil by oil-water interfaces.
an effective total trapped oil saturation (Sgii) is defined as
Sott = Sot + Sotw + Soio (3.7)
where Soi is the effective trapped oil saturation, Sotw is the effective trapped air
saturation contained within the trapped oil that resulted from air-water interfaces, and
Soto is the effective entrapped air saturation contained within the trapped oil that
resulted from air-oil interfaces.
The maximum amount of trapped oil for a given saturation path history is estimated by
- / 1 - ? "»«'» \
<>
s ~
i D fl _ C miiA / v '
T "ouA.1 &w ) /
where 5or is the effective residual total oil saturation, Swmin is the historic minimum
apparent three-phase water saturation, and Rgw is a calibration term defined in (3.3).
5wTOI'n represents the smallest flow channel radii where oil-water interfaces have
occurred.
Effective total trapped oil saturations between 5^ = Swmin and 5^, = 1 are estimated
by
_ / C _ C tnin
_ C f Jw &w (n Q\
- \A fe . I (3.9)
\
_ W|.n
63
-------
where 5^ references the current location of the continuous oil-water interfaces in a
porous medium. The amount of trapped air contained in Sott (i.e., Soio and 50 S *SU
~ SaroSor(5w- $r»)
= a - S^KI - 5,m'n)
for §,*» > *SW
5~ C IC min C min^
5arw Jor \J1 Jw I /O inu\
olw = s ' A^ (O.1UDJ
for I,,,"1"1
mtn ^
'U) i
for 5tt > 51rai'n and Stmin < ^Sw
5C / C mtn Ac nu'^
Soru1 Jor \Jt Jw ) I*
^~ i i i^» .1-^^^^^^^^^^-. I
0TW /^ 7r ' \/* A7r \ ^
fnr ^ ^ C min aT1H C min ^ A C ou/
lOl Jw v O^ ollU »jj «^ iJ^
50<0 = 0 (3.10d)
for Swm>n > A5u,au'
^otw = ^C°minVl _ Ao aw\ (3.10e)
for Swmtn < Sv/
-------
5otu, = 0 (3.10h)
for 5,,, < 5,""'" and Stmi'n < A5u,au'
S,to = 0 (
Sot* = 0 (3.10J)
Concomitant entrapment of air and oil needs to be considered for mass balance purposes
as oil-water interfaces move through a porous medium. The historic saturations need to
be continually updated so that current apparent saturations are not less than their
historic values.
Scaled saturation-pressure function. Three-phase capillary pressure curves are estimated
from two-phase relation using an extension of the scaling proposed by Leverett (1941) as
3£WOUAJ = S*(fc*) (3.Ha)
where 5 £,J/ is the apparent water saturation in the three fluid phase system, S \IJ is the
total liquid saturation in the three-phase system, S*(h*} is the S-P function for a
reference two fluid phase system, and flao and (3OW are fluid-dependent scaling factors.
Capillary heads are defined by
haw = ha - hw (3.12a)
hao = ha- h, (3.12b)
h,w = hc - hw. (3.12c)
If we take a pristine two-phase air-water system as reference such that
5*(/i*) = 5 (haw), then the scaling factors are given by
&,o = TT* (3.13a)
v ao x '
Pow = -^ (3.13b)
65
-------
where 0 for the three-phase system, it may be inferred from
(3.11) that the air-water capillary pressure relations for a contaminated system is
described by
3 I"(l3'aw haj = S*(h*) (3.14a)
ft* = % (3.14b)
where a'aw is the surface tension of water with dissolved components of the oil present.
Organic contaminants generally decrease the surface tension of water so that 0'aw > 1,
which indicates that the capillary pressure in an air-water system with dissolved organic
contaminants will be less than that in the pristine system. It may also be shown that
ao
If /3'aw ^ 1, then Sw(hatv) will exhibit a discontinuity when a transition is made from an
oil-free system to one with oil present. The oil pressure at which this transition occurs
may be found by letting St>Sw, indicating the necessary condition to classify a location
as oil-free is
h0 < V (3.16a)
h? = ^A + JoA, (3 16b)
POW ' "ao
where hcj is a critical pressure head at the transition from a two-phase air-water system
to a three-phase air-oil-water system.
66
-------
Apparent saturation-capillary pressure hysteresis is modeled by describing main
drainage and imbibition S-P relations using the van Genuchten function, while arbitrary
scanning paths are described by scaling the main S-P branches through appropriate
reversal points. The main drainage and imbibition branches are parameterized by
(A*) = (l 4- (D«haw)n)' (3.17a)
*(O = l 4- (7a/Un~m (3.17b)
where ^a, Ja, and n are curve shape parameters with m = 1 - 1/n, and superscript D
and / denote main drainage and imbibition curves, respectively. To predict S*(h*)
scanning paths, the main drainage and imbibition branches are scaled to pass through
the appropriate saturation path reversal points such that closure of scanning 5*(/i*)
paths is enforced. The main imbibition branch is scaled to pass through reversal points
to give the current imbibition scanning path as
+ IDS* (3.18)
where h* and 5* are the current scaled capillary head (0^) and apparent saturation
(Swaw, Sw or 5,), (DIS*,DIh*) is the most recent drainage to imbibition reversal point and
(IDS*,IDh*) is the most recent imbibition to drainage reversal, and 75*(/i*), 75*(£l//i*)
and IS*(IDh*) are apparent saturations of the main imbibition branch that correspond
to the specified scaled capillary heads. Note that when using this function to predict S-P
relations of a primary imbibition scanning path (i.e., an imbibition path that originates
from the main drainage branch), IDS* and 7S*(7I>h*) are equal to 1.
For drainage S*(h*) scanning paths, an analogous procedure is followed to scale the
main drainage branch to pass through the reversal points. The interpolation formula for
an arbitrary drying S*(h*) scanning path is
- DiA
67
-------
35*(7ZV) - rS*(c/fc*)j
(3.19)
where (IDS*,IDh*) is the most recent imbibition to drainage reversal and (DIS*,DIh*) is
the most recent drainage to imbibition reversal, and other symbols follow the notation
of the imbibition formula.
Relative permeability relations. Three-phase relative permeability functions may be
derived using the theory of Mualem (1976) after making suitable corrections to account
for effects of trapped fluids. For the case of the water phase relative permeability,
trapped oil and air trapped within the water phase will occupy some of the pore space
which would otherwise be filled with water, thus displacing water into slightly larger
pores. For the air and oil phases, dispacement of the fluid interfaces within the pore
space will affect the hydraulic radius of the phase, while the portion of the phase which
is trapped is assumed to have no contribution to the phase permeability. Correcting
Mualem's integral to account for these effects yields for the water, oil and air phase
permeabilities
k -
"Vu>
7
J 0
»dSw
h*(S*)
"ot
J 0 h*(S*)
f1
-------
1
0
dS*
(3.20c)
Using the van Genuchten model for S*(h*) and the expressions for apparent and
trapped saturations described above, closed form expressions for phase permeabilities
may be derived. For the case of no trapped fluids, the expressions reduce to the simple
form derived earlier by Parker et al. (1987)
krw = sj/2 [ i - (i-
kro = (5t - 5
]2 (3.21a)
)m - (l-5t1/rn)m ]2 (3.21b)
2m. (3.21c)
Results for the general case in which fluids are trapped following arbitrary saturation
path histories are given by Lenhard and Parker (1987).
3.1.2 Numerical simulation of hysteretic flow
Problem description. The hysteretic three-phase hysteretic k-S-P model described above
was incorporated into an upstream-weighted finite element multiphase flow code, and
the numerical model was employed to simulate a pulse injection of NAPL into the
vadose zone under conditions of a fluctuating water table to assess the effects of
hysteresis.
Parameters used in the simulations (Table 3.1) correspond to a sandy porous medium
and a NAPL less dense and more viscous than water. A one dimensional vertical soil
section is modeled which has a depth of 150 cm discretized in 1 cm intervals. Initially,
the system is assumed to be free of oil and in static equilibrium with a water table at
the lower boundary. Initial fluid saturations are taken to correspond to main drainage
branches of the S-P relations. Water phase pressure at the lower boundary as a function
of time is shown in Figure 3.1. During the first 10 hours, the water table rises as the
69
-------
water pressure head at the lower boundary increases from 0 to 50 cm. At 10 hours, a
pulse of NAPL is applied under a constant oil pressure head of 3 cm water-equivalent
height at the upper surface. NAPL is applied until 7.5 cm3 cm'2 infiltrates into the
porous medium after which a zero flux upper boundary condition is imposed for NAPL.
The upper boundary condition for water and the lower boundary condition for NAPL
are at all times zero flux conditions (with the water head of 50 cm at the bottom, 5-P
relations preclude NAPL from reaching the lower boundary). After NAPL infiltration,
a period of redistribution follows during which water head at the lower boundary is
fixed at 50 cm (A-B of Figure 3.1). Subsequently, the lower boundary water head again
increases (B-C), decreases (D-E), increases (F-G) and finally returns to the 50 cm level
(H-I) with constant pressure periods between period of changing head (Figure 3.1). The
total simulation time is 500 hours.
Three different simulations of the above physical scenario were carried out using
variations of the constitutive model. Hysteretic analyses were carried out using the full
hysteretic model described above which accounts for reversible hysteresis effects and
nonwetting fluid entrapment. Additionally, analyses were conducted using a simplified
model which considers fluid entrapment effects only. Since the scaled saturation
function becomes single-valued for the latter case, storage requirements axe reduced
(scaled reversal saturations are not relevant) and computational requirements are
slightly diminished. Finally, simulations were also carried out with hysteresis entirely
disregarded, for which case the present model reduces to the model described by Parker
et al. (1987).
Simulation results. To assess the system behavior for full hysteretic, trapped fluid only
and nonhysteretic simulation cases, we first consider the temporal system response at a
depth 50 cm below the soil surface. Water, total NAPL and free NAPL saturation
histories at the 50 cm depth are illustrated for the three cases in Figures 3.2 - 3.4.
Water saturations predicted for the trapped fluid only case are similar to those of the
full hysteresis analysis, while the nonhysteretic analysis indicates substantially higher
water saturations since air entrapment is disregarded (Figure 3.2).
70
-------
Table 3.1 Parameters Employed in Numerical Simulations
Parameter
Value
Porous Mciliiini Clmrticierimk's
'o 0.10 cm"1
Ja 0.05 cm"'
S,,, 0
n 2.0
S"; 0.25
'5';' 0.10
's;" 0.25
* 0.4
Saturated water phase 50 cm h"'
conductivity
Fluid Properties
A,,, -25
A.,. '-80
NAPL specific gravity 0.8
(P.,'PJ
NAPL to water viscosity 2.0
ratio (TJ./T;,,.)
100 ZOO 300 400 500
TIME Chours)
Fig. 3.1 Variation in water head at lower boundary for simulations.
71
-------
0. 9
100 200 300 400
TIME Choui-s)
500
Fig. 3.2 Water saturation histories at 50 cm depth.
100 200 300 400
TIME Chour-s)
soo
Fig. 3.3 Total NAPL saturation histories at 50 cm depth.
72
-------
100 200 300 400
TIME
-------
Q_
U
25 -
50 -
75 -
100 -
125 -
t. - 100 n
FULL
HYSTERESIS
NDNHYSTERETJC
/ FLUID
'/ ENTRAPMENT
ONLY
150
0. 0 O. 2 0. 4 0. 6 0. 8 1.0
FREE NAPL SATURATION
Fig. 3.6 Free NAPL saturation distributions at 100 hours.
25 -
SO -
75 -
D.
O 1DD
125 -
150
t - 100 h
FLUID
ENTRAPMENT
ONLY
O. O O. 2 0. 4 0. 6 0. 6 1.0
WATER SATURATION
Fig. 3.7 Water saturation distributions at 100 hours.
74
-------
Q.
LU
D
300 -
125
150
0. 0 20. 0 40. O 60. O 80. 0
NAPL-WATER CAPILLARY HEAD
Fig. 3.8 NAPL-water capillary head distributions at 100 hours.
so -
75 -
£L
D 100
150
. \ NONHYSTERETIC
FULL
HYSTERESIS
FLUID
ENTRAPMENT
ONLY
0. O O. 2 0. 4 0. 6 0. 6 1.0
WATER SATURATION
Fig. 3.9 Water saturation distributions at 200 hours.
75
-------
so -
75 -
Q.
D 1DO
125 -
150
FULL HYSTERESIS
FLUJD ENTRAPMENT ONLY
NONHYSTERETIC
t - 200 h
D. D 0. 2 0. 4 0. 6 0. 6 1.0
FREE NAPL SATURATION
Fig. 3.10 Free NAPL saturation distributions at 200 hours.
E
0
Q.
UJ
D
100 -
125 -
150
PULL
HYSTERESIS
FLUID
ENTRAPMENT
ONLY
AT ZOO HRS
0. 0 0. 2 0. 4 0. 6 0. 6 1.0
TRAPPED NAPL SATURATION
Fig. 3.11 Trapped NAPL saturation distributions at 200 hours.
76
-------
Table 3.2
Free and trapped fluid saturations at a depth of 50 cm which correspond to
the time in Figure 3.1 for different simulations
Time
A
B
C
D
E
F
G
H
1
J
A
li
C
D
E
F
C
H
1
J
A
B
C
D
E
F
G
H
I
J
Su-
0.199
0.299
0.663
0.667
0.493
0.472
0.544
0.582
0.385
0.342
0.200
0.215
0.548
0.579
0.355
0.332
0.405
0.425
0.300
0.266
0.199
0.249
0.549
0.574
0.385
0.347
0.408
0.458
0.320
0.274
s.,,
Nonhvxterrlic
o'.o
0.141
0.337
0.333
0.123
0.076
0.403
0.3S9
0.085
0.057
Full Hvxierexix
0.0
0.143
0.191
0.146
0.112
0.086
0.126
0.159
0.093
0.074
Fluid Entrapment Onl\
0.0
0.136
0.182
0.146
0.102
0.082
0.257
0.252
0.094
0.071
s.,,
0.0
0.004
0.144
0.158
0.063
0.054
0.0X4
0.092
0.040
0.026
0.0
0.014
0.143
0.154
0.072
0.056
0.082
0.104
0.044
0.025
s«,
0.002
0.038
0.117
0.117
0.064
0.055
0.077
0.086
0.049
0.039
0.001
0.050
0.126
0.126
0.076
0.064
0.105
0.116
0.060
0.046
For imbibing water saturation paths corresponding to periods of rising water tables,
large differences occur between total and free NAPL saturations predicted by the
hysteretic and nonhysteretic analyses (Figures 3.3 and 3.4). On the other hand, for
conditions in which water saturation paths do not induce NAPL entrapment, good
agreement occurs between NAPL saturations for the different simulations.
The history of apparent water saturation versus capillary head at the 50 cm depth for
the full hysteretic simulation is shown in Figure 3.5. Imposed boundary conditions
generated two closed scanning loops. The first loop is barely distinguishable from the
primary imbibition curve (Figure 3.5a). Greater detail is shown in Figure 3.5b. Since
initial conditions reflect main drainage branch 5-P relations, a .primary imbibition path
ensued when the water table was raised at time A. Letter symbols in Figure 3.5
correspond to the times designated in Figure 3.1.
77
-------
NAPL injection was initiated immediately following time A. As NAPL moved
downward, water was displaced resulting in a water front preceeding the NAPL front.
The hollow square symbol in Figure 3.5a designates when the fluid system changed from
an air-water to an air-NAPL-water system which is also the origin of the first scanning
loop.
As NAPL displaces water, a secondary drying path ensues. A drying water saturation
path continues until the NAPL front moves beyond 50 cm at which time the NAPL-
water capillary heads adjust to the NAPL front passage causing water saturation to
increase and inducing a saturation path reversal. At time B, which is after the passage
of the NAPL front, the water saturation path is on a tertiary wetting scanning path.
Subsequent rising of the water table to time C produced higher water saturations and
an imbibition saturation path continued until time D. The internal scanning loop
(Figure 3.5b) closed when the water saturation increased beyond the water saturation at
the inception of the three-phase fluid system.
In interval D-F, the water table is lowered initiating a secondary drying scanning path
until time F. Later, the water table elevation increases during period F-H producing a
saturation path reversal to a wetting path. During the final simulation period H-J, the
water table is lowered causing scanning loop F-H-F to close whereupon the secondary
drying path D-J is resumed.. The second scanning loop F-H-F is more distinguishable
than the first because the range of apparent saturations encountered in F-H-F is
significantly greater than that encountered in the first scanning loop. Free and trapped
fluid saturations corresponding to the times in Figure 3.1 are listed in Table 3.2. For
nonhysteretic and fluid entrapment only simulations, scanning S-P loops do not exist
because S-P relations are modeled to follow the main drainage branch.
To investigate the spatial distribution of fluids for the various cases in greater detail, we
consider saturation distributions at t = 100 h (time B) and at 200 h (time D). At 100
h, injected NAPL has redistributed for a period of approximately 90 h while boundary
conditions remained constant leading to nearly static conditions. Predicted free NAPL
distributions at 100 h are shown for the three simulation cases in Figure 3.6. Nearly
identical free NAPL distributions are predicted for the three cases reflecting the fact
that water saturation was nearly constant over the period of NAPL- redistribution so
78
-------
that no NAPL entrapment is predicted. In the upper part of the profile where total
liquid saturation is draining, all cases effectively use the same S-P functions (i.e.,
drainage), so that results are indistinguishable. Slightly greater deviations arise lower in
the profile where saturation paths are more complicated, but it is clear that effects of
reversible hysteresis phenomena are of secondary importance at this time.
Finally concerning the NAPL predictions, it is worth noting that all simulations predict
a nearly constant NAPL saturation in the upper zone of about 17% which for practical
purposes is irreducible under the influence of gravity. This nondrainable residual
saturation should be clearly distinguished from the residual oil saturation introduced in
the constitutive model which represents the maximum amount of trapped NAPL. The
nondrainable saturation evidenced in the present simulation results is not trapped and it
is not parametrically defined by the constutitive model. It is simply an outcome of the
relative permeability model which indicates that as NAPL saturation decreases, NAPL
permeability eventually diminishes to the point that flow virtually ceases and precludes
attainment of true equilibrium conditions.
Water saturation distributions versus depth at 100 h are shown for the three simulation
cases in Figure 3.7. Full hysteretic and fluid entrapment results are quite similar;
however, the nonhysteretic case overpredicts water saturation since air entrapment
associated with the rising water table is unaccounted for. Comparison of NAPL-water
capillary head distributions at 100 h shows very little difference among the three cases
(Figure 3.8) indicating that predicted apparent water saturations are nearly identical.
At 200 h the water table has been held at its highest elevation for an extended period.
Deviations between hysteretic and nonhysteretic model predictions of water saturation
exhibit similar trends to those observed at 100 h, but are now more pronounced due to
increases in air entrapment associated with the rising water table (Figure 3.9). Full
hysteresis and fluid entrapment only simulations again exhibit rather close
correspondence except near the upper surface where reversible hysteresis effects
apparently become more significant due to increasing separation between scaled
scanning curves.
In contrast to the results at 100 h, free NAPL saturation distributions at 200 h exhibit
marked differences for the hysteretic and nonhysteretic simulation cases (Figure 3.10).
79
-------
Full hysteretic and fluid entrapment only analyses indicate free NAPL volumes which
are only about 1/3 of that predicted by the nonhysteretic case. The discrepancy reflects
the substantial amount of trapped NAPL for the hysteretic cases (Figure 3.11) which is
by assumption nonexistent in the nonhysteretic analysis.
Mass balance calculations for NAPL were performed during all three simulations by
integrating NAPL contents over the spatial domain and normalizing to the injected
NAPL volume. The mass balance errors ranged from 1.7 to 4.9% using the pressure-
pressure formulation described in Section 2.2.1, with the lowest error for the full
hysteresis model. Apparently, for this particular problem, inclusion of hysteresis had a
stabilizing effect on the numerical solution, not only reducing the mass balance error
but slightly reducing the average number of iterations per time-step for nonlinear
convergence. The latter effect partially offset the additional computational burden for
the hysteresis model.
3.2 Development and Testing of a Simplied Hysteretic Model
The model described in Section 3.1 distinguishes two sources of hysteresis: (I) hysteresis'
in apparent water saturation (actual water saturation trapped fluid saturation in
water) versus oil-water capillary pressure relations and apparent liquid saturation
(actual water saturation -f actual oil saturation trapped fluid saturations in water
and oil) versus air-oil capillary pressure relations due to contact angle hysteresis and
differences in pore-scale fluid entry and drainage pressures, and (II) hysteresis in
saturation-capillary pressure relations and in relative permeability relation-saturation
relations due to nonwetting fluid entrapment associated with pore bypassing during
wetting fluid imbibition. The results of numerical simulations discussed in the preceding
section indicated that a simplified model that considers only Type-II hysteresis yielded
nearly identical results to the full model with both Type-I and Type-II hysteresis.
Simplification of the model to consider only Type-II hysteresis is an attractive option as
this will markedly reduce the computational effort and storage requirements of the
model. Furthermore, by invoking this simplification it becomes feasible to treat various
aspects of fluid entrapment with greater rigor than is possible for the full model. In
particular, although the detailed model describes maximum trapped fluid saturation for
80
-------
a given saturation path as a hyperbolic function of wetting fluid saturation at the
drainage-imbibition reversal point, actual fluid entrapment at a given wetting fluid
saturation is taken as a linear function of wetting fluid saturation. The latter
assumption is in fact inconsistent with the former, but is invoked for the sake of
expediency.
In this section, we describe a simple model for three-phase k-S-P relations that accounts
for effects of nonwetting liquid entrapment in a stringent fashion while avoiding
complications arising from less important factors. A number of one- and two-
dimensional simulations involving low and high density hydrocarbons will be performed
to demonstrate the behavior of the model.
8.2.1 Model description
Saturation-capillary pressure relations. The proposed constitutive model is based on the
nonhysteretic three-phase k-S-P relationships of Parker et al. (1987) with extensions to
consider hysteresis due to nonwetting liquid entrapment. Explicit consideration will not
be given to gas entrapment during periods of increasing total liquid saturation, since gas
entrapment has a relatively minor effect on liquid phase relative permeabilities. Effects
of trapped gas may be accomodated implicitly by employing the maximum liquid-filled
porosity in lieu of actual porosity and the hydraulic conductivity at "field saturation" in
lieu of the fully saturated conductivity (Kool and Parker, 1987). It will be assumed in
the following that these definitions of porosity and saturated conductivity have been
employed where relevant.
Adopting the notation of the preceding section, effective water and total liquid
saturations are defined by
~SW= S~m (3-22a)
5'= ^2 (3-22b)
81
-------
where St = Sw -\- -S0 is the total liquid saturation and 5m is the irreducible water
saturation. Effective oil saturation is defined such that S0 = St - Sw, hence
~S0 = -y&jr- (3.23)
Furthermore, oil saturation is divided into two fractions so that 50 = Soi + Sot where
subscripts / and t denote "free" and "trapped" oil, respectively. In terms of effective
free and trapped oil saturations, defined by analog to (3.23), we have
~S0 = ~Sot + ~Soi. (3.24)
Finally, we define a quantity referred to as apparent water saturation, S w, which
corresponds to the sum of effective saturation of water and trapped oil, i.e.,
1 , = ~SW + 50, (3.25)
Following Parker and Lenhard (1987), apparent saturation will be assumed a simple
function of capillary pressure.
Prior to the occurrence of oil at a given location, the porous medium is treated as a two-
phase air-water system described by the van Genuchten (1980) function
5, = [ 1 + (a haw)« J- (3.26)
where a and n are porous medium parameters and m = 1-1/n, and haw is the air-water
capillary head. Following the occurrence of oil at a given location, the system is
described by the three-phase relations
1 w = [ 1 + (a $ow how)» ]'m (3.27a)
3, = ( 1 + (a 0ao hao)n ]-» (3.27b)
where how and hao are oil-water and air-oil capillary heads, and /?ao and 0OW are scaling
factors. To avoid numerical oscillations associated with changes from a two-phase air-
water system to a three-phase air-water-oil system, once a node is classified as a three-
82
-------
phase node, it is not be allowed to revert to a two-phase air-water system. The
necessary condition to classify a location as a three-phase node is
0 + 0
"ow ' ~a
(3.28)
To compute water saturation in the three-phase system, an expression for Sot is
required. Since gas entrapment is disregarded in the simplified model, we assume oil
entrapment in a three-phase air-oil-water system can be inferred directly from that in a
two-phase oil-water system. Consider an oil-water system with main drainage and main
imbibition curves described by (3.27a). On any given scanning curve, the effective
water saturation at how = 0 is less than 1 due to oil entrapment given by 1 - Sor where
Sor is the effective residual oil saturation which may be estimated using an empirical
relation given by Land (1968) as
_ t _ c
ST =
OT
R l-
1 (3.29b)
where 5J^'n is the minimum effective water saturation corresponding to the reversal
from water drainage to imbibition and S?" is the effective residual oil saturation for
the main imbibition curve. Residual oil saturation corresponds to the trapped oil
saturation at zero capillary pressure.
In Section 3.1, the trapped saturation at nonzero capillary pressure was estimated by
linear interpolation with ~Sot = 0 at f , = Sj7"" and ~Sot = 507" at f , = 0. This
method has the advantage of simplicity, but is inconsistent with (3.29). Linear
interpolation can also lead to numerical problems at capillary heads near zero due to
excessive changes in trapped oil saturation with small head changes. A procedure which
is fully consistent with (3.29) involves estimation of trapped oil saturation as the
difference between residual oil saturation for the actual scanning curve and that for a
curve with a reversal point equal to the saturation on the actual path.
83
-------
To clarify, consider Point 2 on Figure 3.12 with apparent water saturation 5 w on a
path starting from Point 1 where 5,,, = 5J""1. The effective trapped oil saturation at
Point 2 is computed as the difference between the residual saturation on the path
starting from Point 1 and that starting from Point 3. From (3.29) we have
_ { 1 - S» 1-5,_ \
"11+ R(l-S) i + tf(l-SJ /
(3.30)
Equation (3.30) is valid if 5^ > S ,""", otherwise 50, = 0. For the three-phase air-oil-
water system, an additional constraint which must be enforced is that 5ot at a given
node is the maximum of Sot computed by (3.30) and 50 at the node that is, the
system cannot trap more oil than is present at the node. Additionally, for the three-
phase system, 5Jn>n takes on the operational meaning of being the minimum effective
water saturation since the first occurrence of oil at the node. The value of SJ711" is
updated at the end of each time step after convergence and stored for each node.
Relative permeability relations. At a given water saturation, oil entrapment affects
water permeability by displacing water into larger pores. The net effect can be shown to
be small and experimental studies often report no observable hysteresis in wetting phase
permeability vs. wetting fluid saturation relations (e.g., Saraf et al., 1982). Therefore, in
the present analysis, we will assume water relative permeability to be a simple single-
valued function of water saturation as described by Parker et al. (1987). Since gas
entrapment is not considered in the present analysis, we will assume that gas relative
permeability can be described as a unique function of gas saturation also as given by
Parker et al. (1987). To describe oil relative permeability, it is assumed that only free
oil is a continuous phase subject to convection while trapped oil cannot move. Parker
et al. (1987) derived an expression for oil relative permeability in the absence of fluid
entrapment which is a function of effective water saturation and effective total liquid
saturation. In the presence of trapped oil, the oil-water interface in the porous medium
will be displaced to larger pores identified with the apparent water saturation.
Replacing effective water saturation by apparent saturation in the previous model
should therefore yield a good representation of oil relative permeability.
84
-------
Phase relative permeabilities in the present model are described by the following
relations
* = ~SJ'* [ 1 - (1-5^/T ]2 (3.31a)
kro = (5,-f J1/2[(l-1^/m)m-(l-5t1/rn)m]2 (3.31b)
kra = (1 - 5t)1/2 [ 1 - St1/m ]2m (3.31c)
where fcruj, kro and kra denote relative permeabilities to water, oil and air, TO = 1-1/n is
the van Genuchten parameter and saturations are as previously defined.
Table 3.3. Soil and fluid properties used in simulations.
Property
*
sm
n
a [L-1]
K,W[L T'1]
OV71OT
Pao
"ow
Pro
Iro
Example 1
0.4
0
2.5
0.05
40
0-0.4
2.67
1.6
0.8
2
Example 2
0.4
0
3
0.8
3
0.25
1.8
2.25
0.8
2
Example 3
0.4
0
3
0.8
3
0.25
1.8
2.25
1.2
2
Units for Example 1 in cm and hours and for Examples 2 and 3
in meters and days. Conductivity assumed isotropic.
85
-------
3.2.2 Example problems
Example li One-dimensional infiltration and redistribution of low density NAPL. The
purpose of this problem is to investigate model sensitivity to the maximum residual oil
saturation, SJ£ai. The flow domain is taken to be a 75 cm long soil column which is
initially free of oil and at hydrostatic equilibrium with a water table located 10 cm from
the lower surface. The soil hydraulic and fluid properties used in the simulation are
given in Table 3.3. A total of seven simulations were performed with S?0ax values
ranging from 0 to 0.5 which is the extreme range which is likely to occur in geologic
media. Note that 5J^OZ = 0 corresponds to no oil entrapment in which case the model
reduces to that described by Parker et al. (1987). Boundary conditions for each of the
simulations consist of following stages:
Stage 1: A slug of oil is allowed to infiltrate from the upper surface at a water
equivalent head of 1 cm until the total accumulation is 5 cm3cm~2 while a zero-flux
condition was stipulated for the water phase. At the lower boundary, the water head is
maintained at its initial value of 10 cm and the oil phase is zero-flux.
Stage 2: Oil is allowed to redistribute for 25 h under zero-flux conditions for both water
and oil at the upper surface, while the lower boundary conditions remain the same as in
Stage 1.
Stage 3: The water head at the lower boundary is raised from 10 cm to 35 cm over a 1
h period and maintained at this value for 175 h. A zero-flux condition is imposed for
water at the upper boundary and for oil at both boundaries.
Simulations were carried out with a finite element mesh consisting of 75 elements using
a time-step size varying from 0.001 to 1 hour. Predicted oil saturations versus elevation
at the end of Stage 3 are shown in Figure 3.14 for the cases S%?x = 0 and 0.25. During
Stage 3, water saturation increases as the water table rises resulting in oil entrapment
when this is considered by the model. As a result, vertical displacement of oil is much
less efficient and the zone of oil saturation is much more dispersed when oil entrapment
is considered. At the end of the simulation nearly 40% of the oil in the system is
trapped for the case with 5^QI = 0.25.
86
-------
Variations in the percentage of added oil which has become trapped at the end of Stage
3 in response to variations in 5J£ai are shown in Figure 3.13. The results indicate that
the fraction of the 5 cm3cm~2 spill volume which becomes trapped following a 25 cm
increase in water table elevation increases in a markedly nonlinear manner with
increasing S%?x. The entire spill becomes immobilized under these conditions when S?0ax
« 0.4. For a given maximum residual saturation, immobilization of the spill would be
greater if the total spill volume were smaller or the water table rise greater.
0
01
r
x
L
0
Q_
D
U
Effective water saturation
-------
100
0. 0
D. 1
0. 2
C. mcx
or
0. 3
Fig. 3.13 Variation of percentage trapped oil volume in domain with maximum
residual oil saturation for Example 1.
75
Total oil
-------
Example 2i Two-dimensional infiltration and redistribution of low density NAPL. The
objective of this example is to demonstrate the effects of fluid entrapment under field
conditions in which lateral migration of hydrocarbon is significant. The flow domain is
taken to be a 4 m deep x 17 m wide vertical slice through the unsaturated and
saturated zones which is initially free of NAPL (Figure 3.15) and in equilibrium with a
water table along plane BI with a gradient of 1.6%. Soil hydraulic and fluid properties
used in the simulation are given in Table 3.3. Simulations were performed with S^x =
0 and 0.25 using a finite element mesh consisting of 232 elements with the time-step
varying from 0.0001 to 0.4 days. Boundary conditions were changed in three stages as
follows:
Stage 1: A slug of oil is allowed to infiltrate into the system under a water equivalent
head of 1 cm along a 0.5 m long strip (EF) until the total accumulation is 0.8 n^m"1.
All other boundaries are zero-flux for oil. Water head is fixed at the initial conditions
on sides AD and GJ and other boundaries are zero-flux for water.
Stage 2: Oil is allowed to redistribute under ambient hydraulic conditions to t = 8
days with zero-flux for oil on all boundaries and water boundary conditions the same as
in Stage 1.
Stage 3: Water heads on the segments AD and GJ are increased by 0.7 m raising the
water table to points C and H on the boundaries over a 0.5 day period and other
conditions are maintained as in Stage 2 for another 12 days bringing the total
simulation time to 20 days.
The duration of Stage 1 was determined to be 0.45 days irrespective of the value of
S^". Figure 3.16 shows the variation in percentage of the oil volume in the domain
which is trapped over the course of the simulations for the case with S?*x = 0.25. The
results indicate that negligible oil entrapment occurs during Stages 1 and 2 since water
saturation paths are predominantly draining. As the water table rises during Stage 3,
the percentage of trapped oil increases sharply to almost 55% of the spill volume.
Distributions of total oil saturation (free + trapped) at the end of Stage 3 are shown in
Figure 3.17 for simulations with S? = 0 and 0.25. Trapped oil saturation
distributions are shown for the case of S?*x 0.25 (for S" = 0 no oil is trapped).
89
-------
The results show that fluid entrapment markedly reduces the lateral migration of the
liquid plume since much of the oil volume is trapped and thus immobile. When fluid
entrapment is disregarded, the upward displacement of the oil by the rising water table
is more efficient and lateral migration occurs more readily.
3.5
E F
o
B
17
Fig. 3.15 Two-dimensional flow domain used in Examples 2 and 3.
90
-------
100
eo
a
a.
o
<
20 H
Water-table rise
Low density
High density
5 10 IS
Time (days)
Fig. 3.16 Variation of percentage trapped oil volume with simulation time for
Examples 2 and 3.
Example 3j. Two-dimensional infiltration and redistribution of high density NAPL. The
objective of this simulation is to investigate the effects of fluid entrapment on the
migration of dense NAPL plumes. The simulation will be identical to that described in
Example II except that the fluid density ratio, pro, will be taken as 1.2 instead of 0.8
(Table 3.3). The flow domain, initial and boundary conditions and numerical
discretization is the same as in the previous example. The first stage of the simulation
predicted an infiltration time of 0.36 days for a total oil accumulation of 0.8 nAir1.
The predicted variation of the percentage of oil volume in the domain which is trapped
with simulation time is shown in Figure 3.16. The dense NAPL simulation predicts
almost 10% of the oil to be trapped at the end of Stage 2 whereas no trapped oil is
predicted for the light hydrocarbon at the end of Stage 2. The reason for the greater
entrapment for the dense fluid is the more complex saturation history of the medium
associated with the penetration of dense NAPL into the saturated zone. As the front
initially moves downward, high NAPL saturations develop in the capillary fringe and
upper saturated zone which gradually diminish as gravity drainage continues. As water
saturation rebounds, oil entrapment occurs. As the water table is raised, additional
entrapment occurs in the capillary fringe zone.
91
-------
Predicted oil saturation contours are shown in Figure 3.18 for simulations with and
without fluid entrapment effects. Marked differences in oil saturation distributions occur
for the two cases. A much larger fraction of the spill volume is predicted to be retained
in the unsaturated zone and capillary fringe when fluid entrapment is considered and
the NAPL plume migrates laterally along the aquifer lower boundary to a far lesser
extent. The time required for NAPL to first reach the lower boundary was also much
longer when fluid entrapment was considered.
Total oil
S ox-0. 25
Trapped oil
Sromox-0. 25
Figure 3.17 Predicted oil saturation distribution for Example 2 at end of Stage 3.
92
-------
"Propped oil
i-ox-0. 25
Fig. 3.18 Predicted oil saturation distribution for Example 3 at end of Stage 3.
3.3 Experimental Verification of Hysteretic Flow Model
S.S.I Static three-phase measurements
Experimental methods and materials. Measurements of fluid contents and pressures as a
function of saturation path history were taken in two different porous media to evaluate
the hysteretic S-P model. The experimental apparatus used to conduct the two- and
three-phase S-P measurements consists of alternating plexiglass sleeves containing
treated and untreated porous ceramic rings on their inner surfaces (see Figure 3.19).
There are two treated and two untreated ceramics rings composing the retention cell.
The treated ceramics were immersed in chlorotrimethysilane to render them
hydrophobic. This is done so that NAPL could be continuous from NAPL-filled flow
channels within a porous medium to a pressure transducer located outside the retention
cell. In this manner, the NAPL pressure within a porous medium could be measured
with a pressure transducer mounted outside the experimental cell. Pressure transducers
93
-------
were manufactured by Soil Measurement Systems in Tucson, Arizona, and had a
sensitivity of 1 mbar. The untreated ceramics maintained a continuous water phase
from water-filled flow channels within a porous medium to the pressure transducer for
the aqueous phase also located outside the retention cell. Therefore, NAPL and water
pressures could be measured simultaneously for a three-phase fluid system. The gas
phase within the porous media was always in contact with the ambient atmosphere and
its pressure was assumed to be atmospheric (i.e., zero gage pressure) when liquid
movement within the retention cell ceased.
SCHEMATIC OF THREE-PHASE FLUID
PRESSURE-SATURATION APPARATUS
POROUS MEDIA
FITTINGS
TO TRANSDUCER.
BURET AND
PRESSURE/VACUUM
REGULATOR *-
DETAIL OF CELL SEGMENT
FITTING
CERAMIC RING
TO TRANSDUCER.
BURET AND
PRESSURE/VACUUM
REGULATOR *
* TO TRANSDUCER.
BURET AND
PRESSURE/VACUUM
REGULATOR
PLEXIGLASS CAVITY
Figure 3.19. Experimental saturation-pressure apparatus.
94
-------
Two vacuum-pressure regulators, one for water and the other for NAPL, controlled
movement of liquids into and out of the cell. Each regulator was connected to a buret
that was connected to the retention cell via tubing. Fluid contents were determined by
recording the volume of fluid that drained from the porous medium or that imbibed
into the porous medium. The apparatus also was used to measure two-phase air-water.
air-NAPL, and NAPL-water S-P relations. For a more detailed description of the
experimental apparatus, see Lenhard and Parker (1988).
The NAPL was Soltrol 220, which is manufactured by Phillips Petroleum Company,
Bartlesville, Oklahoma. Soltrol 220 is a mixture of branched alkanes having low
solubility in water and a liquid density of 0.80 g cm"3. The aqueous phase was distilled
water. Air-NAPL and NAPL-water interfacial tensions, as measured by the ring
method, were 0.026 and 0.036 N m" , respectively. The NAPL-water interfacial tension
reflects prolonged contact between NAPL and water.
The porous media employed in the experiments were unconsolidated. Porous medium 1
was a glass bead mixture with 99% of the beads between 0.125 and 0.25 mm, 0.06%
between 0.106 and 0.125 mm, and 0.04% less than 0.106 mm. Porous medium 2 was a
sand consisting of approximately 97%, 1%, and 2% sand, silt and clay sized particles
(USDA), respectively. A measured amount of air-dried porous media was packed into
the retention cell to achieve a constant bulk density. Preliminary investigations were
undertaken to determine what air-dry weight would yield a mass density where porous
media grains would not move significantly during liquid wetting and drying cycles (i.e.,
a constant pore volume). The resulting bulk densities of the glass bead and sand
packings were 1.59 and 1.72 g cm"3, respectively.
Selection of the two porous media used in the experiments was based on satisfying a
major constraint of the model that it be rigid porous media (i.e., no swelling clay
minerals that may induce changes in the pore-size distribution with fluid content
changes). Although both porous media can be classified with respect to grain-size
distributions as sands (USDA), flow channels that result from packing glass beads will
have a much different geometry than those that result from packing natural angular
porous material. The glass beads are significantly smoother and more rounded than are
the sand grains (porous medium 2). This variation is expected to produce differences in
the amount of nonwetting fluid that becomes entrapped and in contact angle changes
95
-------
with saturation path reversals. Porous media were added in several stages to the
retention cell containing de-aired water. After each stage, the porous media were mixed
to homogenize the packing and to dislodge any air bubbles that might have become
trapped. The initial condition for the S-P measurements, therefore, was a completely
water-saturated porous medium.
Two- and three-phase saturation path histories were initiated by desorbing wetting fluid
(i.e., water for two-phase air-water, two-phase NAPL-water, and three-phase systems,
and NAPL for a two-phase air-NAPL system) from the porous media via adjustments in
the vacuum-pressure regulators. As the wetting fluid content decreased, a continuous
nonwetting phase (i.e., air for two-phase air-water, two-phase air-NAPL, and three-
phase systems, and NAPL for a two-phase NAPL-water system) entered the void
spaces. Measurements along this saturation path corresponded to the main wetting fluid
drainage branch. After the wetting fluid content decreased to some low arbitrary level,
wetting fluid was permitted to imbibe back into the porous media by adjusting the
vacuum-pressure regulators. Measurements along this second saturation path
corresponded to a wetting fluid imbibition scanning path.
For two-phase fluid systems, only the main wetting fluid drainage and imbibition
scanning paths were measured. In addition, measurements for the imbibition scanning
path were conducted until the porous media appeared to be apparently saturated (i.e.,
containing wetting fluid and entrapped nonwetting fluid at a zero capillary pressure).
This yielded a measurement of the residual nonwetting fluid content (5,r'J) from which
the maximum amount of entrapped nonwetting fluid (i.e., residual nonwetting fluid
content) for the main imbibition branch (/5,>^) could be calculated by inverting Land's
(1968) algorithm. The residual nonwetting fluid contents associated with the main
wetting fluid imbibition branch are parameters in the hysteretic S-P model.
For three-phase saturation path histories, NAPL was allowed to imbibe into the porous
media at a point along the imbibition scanning path. When NAPL imbibed, the fluid
system changed from a two-phase to a three-phase system. For subsequent three-phase
S-P measurements, NAPL and/or water were permitted to imbibe into the porous
media or were extracted from the porous media by adjusting the vacuum-pressure
regulators. Both NAPL and water pressures and saturations were recorded for each
three-phase measurement. This yielded simultaneous measurements of water and total
96
-------
liquid saturation path histories that can be compared to model predictions. The total
liquid path refers to water saturations during air-water S-P measurements and to the
sum of water and NAPL saturations during three-phase S-P measurements. The
saturation path history was different for each porous media packing. During the three-
phase S-P measurements, cyclic changes in fluid contents (i.e., wetting and drying) were
imposed to investigate scanning path behavior.
Experimental results. Figure 3.20 shows the water saturation path history for porous
medium 1. The open square symbols are two-phase air-water measurements, the closed
diamond symbols are three-phase measurements, and the dashed lines connect
successive measurements. The first series of measurements involved draining water from
an initially water-saturated porous medium to Swaw = 0.42 (i.e., ASu,au/ = 0.42). As
water drained, air entered the porous medium. This saturation path corresponds to the
main drainage branch. An imbibition scanning path was then initiated from Swaw =
0.42 to Swaw = 0.55. During this path, air should theoretically become trapped and
occluded by water as air-water interfaces advance into larger flow channels.
NAPL was then allowed to imbibe into the glass beads creating a three-phase system
while the water saturation was held constant at 0.55. Thereafter, water was drained
forming a drying scanning path until the main drainage path was rejoined at Sw =
0.42. Note that when scaled capillary heads (i.e., /?y\) are employed, a closed 5-P loop
is formed when the main drainage branch is rejoined. Air trapped by air-water
interfaces and occluded by water during the imbibition segment of the 5-P loop should
become occluded by NAPL during the drainage segment. This follows from the
assumption that entrapped fluid ganglia axe immobile as fluid-fluid interfaces move
within a porous medium. Upon rejoining the main drainage path, air initially trapped
by air-water interfaces should be occluded by NAPL or be released into the continuous
gas phase. Whether trapped air remains occluded by NAPL or released into the gas
phase depends on the corresponding total liquid saturation path history.
97
-------
I
E 60-
U
O 50 -
^
UJ
1 40 -
>
tr
< 30 -
'
'
£20-
U
10 -
o
UJ
_) n -
\^
\\ \
\\v
t ^ SP-
\ \ h \ ^°"
^ ^ B^ ^^""B-^
\ ^^VV,^ "%
\ \ Nj. S
V v 1
~* \
\
1
1
WATER SATURATION
PATH HISTORY I
1
U
0.0 0.2 0.4 0.6 0.8
WATER SATURATION
1 .0
Fig. 3.20 Experimental two- (open squares) and three-phase water (closed diamonds)
saturations measured for the water saturation history in porous medium 1.
IU J\J
X
E
JJ 25
D
< 20 -
UJ
I
<
M 10 "
OL
U 5-
D
UJ
-J o -
l *
\ ^-^
\ \ ^
\\\ \
^X^il \
^ Ri a
Vt»l T
]L\ i i
1 1
-I1 '
TOTAL LIOUIO * { '
SATURATION PATH \ '
HISTOBY
U
tn
0.0 0.2 0.4 0.6 0.6 1.0
TOTAL LIQUID SATURATION
Fig. 3.21 Experimental two-phase water (open squares) and three-phase total liquid
(closed diamonds) saturations measured for the total liquid saturation path
history in porous medium 2.
98
-------
After rejoining the main drainage path at Sw = 0.42, water contents were further
reduced, yielding a continuation of the main drainage branch. As water drained, NAPL
filled the water-vacated flow channels. Note that the three scaled main drainage S-P
points for the three-phase fluid system appear to be an extension of the air-water main
drainage data. At 5^ = 0.24, which was the lowest water saturation attained (i.e.,
Swmtn _ 0.24), another imbibition scanning path was produced by allowing water to
imbibe into the porous medium. During this imbibition path, however, initially only
NAPL is trapped as NAPL-water interfaces advance into larger flow channels. This is
because air has never entered these flow channels. The lowest total liquid saturation for
this experiment was 0.37 (i.e., 5
-------
saturations in the three-phase fluid system, and the dashed lines connect successive S-P
measurements. As in Figure 3.20 and for all S-P experiments, the main drainage path
was initially followed. Water was drained from porous medium 2 until Swaw = 0.40,
then an imbibition scanning path ensued.
At Swaw = 0.60 on the first imbibition scanning path, NAPL was permitted to imbibe
into the sand creating a three-phase fluid system. For the next three measurements, the
total liquid saturation continued to increase. These measurements should lie on the
same imbibition scanning path that was formed in a two-phase fluid system. No
saturation path reversals occurred that would generate a different saturation path. Note
that when scaled capillary pressures are employed, these three experimental three-phase
S-P points appear to be an extension of the imbibition path formed in the two-phase
system. Scaling the capillary pressures by best-fitting essentially removes effects that
are due to differences in interfacial tensions and contact angles of the different fluid
systems (Lenhard and Parker, 1987b). Theoretically, scaled S-P main branches should
be unique for a given porous medium. Scaled S-P scanning paths, however, are
dependent on the saturation path history.
At St = 0.78 on the first imbibition scanning path, the total liquid content was
decreased, resulting in a drying scanning path. This path was followed until St = 0.53.
At this point, the total liquid content within the sand was increased again forming
another imbibition scanning path, which was internal to the first imbibition scanning
path. Note that the second imbibition scanning path did not close with the first at
St = 0.78. This is because the saturations that index the historic fluid-fluid interface
locations changed. However, if apparent saturation scanning paths were shown in lieu of
actual saturation paths in Figure 3.21, they would close at the scaled capillary pressure
Lead corresponding to the beginning of the drying scanning path (i.e., St = 0.78).
The explanation why scanning paths would close when apparent saturations are used
but not when actual saturations are employed follows from tracking the movement of
5tro". On the first imbibition scanning path, air was trapped by air-water interfaces
from 5| = 0.40, which corresponded to *SW, to St = 0.60, which was the saturation
where the fluid system changed from a two to a three-phase system (i.e., 5tmi'w = 0.60).
For the three-phase segment of the first imbibition scanning path, air was trapped by
advancing air-NAPL interfaces from St = 0.60 to St = 0.78. On the following drying
100
-------
scanning path, air was released to the continuous gas phase as air-NAPL interfaces
receded into smaller flow channels. As air-NAPL interfaces receded into flow channels
smaller than those corresponding to 5, = 0.78, but larger than 5, = 0.60, air trapped
by air-NAPL interfaces was released into the continuous gas phase. As air-NAPL
interfaces receded into flow channels smaller than those corresponding to St = 0.60, but
larger than St = 0.40 (note that ^Swaw = 0.40), air trapped by air-water interfaces
was released into the continuous gas phase. Therefore, as 5, became less than 0.60, Stmtn
changed. Stmin at the termination of the drying scanning path was 0.53 (see Figure 3.21).
For the ensuing second imbibition scanning path that originated at St = 0.53, all air is
trapped by air-NAPL interfaces as they move into larger flow channels. The total
amount of air trapped on the second imbibition scanning path, however, will be less
than the amount of air released on the preceding drying scanning path because ISarao is
less than ISa,raw. Hence, for the same scaled capillary pressure head corresponding to the
initiation of the drying scanning path, there should be a higher water saturation for the
second imbibition scanning path than for the first imbibition scanning path. This is
because less air is entrapped on the second imbibition scanning path, which agrees with
what was experimentally observed in Figure 3.21. When apparent saturations are
modeled, these differences are taken into account by definition. Apparent saturations
track the movement of fluid-fluid interfaces that separate immiscible fluid phases.
As the total liquid content was increased beyond 5, = 0.78, the measurement point at
St = 0.84 appears to lie on an extension of the first imbibition path, which began in the
two-phase system. It is assumed that once a drying or imbibition scanning path has
closed, the previous drying or imbibition saturation path is respectively followed. The
experimental measurements in Figures 3.20 and 3.21 confirm this behavior. The final
total liquid saturation path for the sand was another drying scanning path. However,
the latter drying scanning path originated from a saturation reversal of 5^, = 0.84, and
the former drying scanning path was at a reversal of Sw 0.78. Note how the latter
drying scanning path does not cross over the former drying scanning path and the
second imbibition scanning path does not cross over the first imbibition scanning path.
This suggests that higher order imbibition and drying scanning paths are internal to
their corresponding lower order paths. From the data shown in Figures 3.20 and 3.21,
there appears to be orderly progression of saturation paths as the saturation path
history develops. Hence, an accurate representation of history dependent S-P relations
may be obtained with a model that incorporates the major observed features.
101
-------
Modeling results. Model parameters were obtained by best-fitting the model, via a
nonlinear regression algorithm, to experimental water and total liquid saturation path
histories. Results of these regression analyses are given in Table 3.4. Also shown are the
residual nonwetting fluid saturations of the main imbibition branch, *5,>';'.
Figure 3.22a shows a comparison between predicted and experimental water saturation
path histories for porous medium 1. The dashed lines connect successive S-P
measurements and the solid line is the model predictions. There is close agreement
between predicted and experimental water saturation path histories. Significant
discrepancies between predicted and experimental S-P relations are found only for the
larger scanning loop. Parameters used for predicting the water saturation path history
in Figure 3.22a were obtained by best-fitting the model to the experimental two- and
three-phase water saturation path history.
An alternative method for calibrating the model is to best-fit the van Genuchten (1980)
retention function to only hysteretic air-water S-P data to obtain the parameters, ca, ;a,
n and Sm. Ratios of air-water to NAPL-water and air-NAPL interfacial tensions can be
used to predict the scaling parameters /3OW and /?ao, respectively (Lenhard and Parker,
1987). The residual nonwetting fluid saturations of the main imbibition S-P branch for
two-phase fluid systems (/5irt^), however, will have to be experimentally determined by
imbibing wetting fluid into the porous medium that is initially saturated with
nonwetting fluid until the porous medium is apparently saturated. Using this approach
to obtain model parameters, Figure 3.22b was constructed. Note that this alternative
calibration method, which can be used in the absence of three-phase S-P data, predicted
the measured two and three-phase water saturation path history in porous medium 1
just as accurately as did the best-fit calibration method.
Figure 3.23 shows the total liquid saturation path history measured in porous medium
1. Again, there is relatively close agreement between predicted and experimental total
liquid saturation path histories. Model parameters were obtained from best-fitting to the
experimental data (see Table 3.4). Figure 3.23 shows two closed scanning paths;
however, they are difficult to discern. The coefficient of rmiltiple determination (R2) in
the regression analysis for the total liquid saturation path history was 0.96.
102
-------
If the alternative method of calibrating the model to a porous medium is used for
predicting the total liquid saturation path history in porous medium 1, the results are
not as accurate as those for the water saturation path history. This can be deduced by
noting that the ratio of air-water to NAPL-water interfacial tensions (i.e., 0OW) required
for predicting water saturation paths is 2.0. The best-fit Pow was 2.05 (see Table 3.4)
when the model was regressed to the experimental data. These values are very close to
each other; hence, model predictions by both calibration methods should be similar, as
shown in Figures 3.22a and 3.22b. However, when predicting the total liquid saturation
path history by using the alternative calibration method, the ratio of air-water to air-
NAPL interfacial tensions is 2.77, whereas the best-fit /?ao was 2.39.
A noteworthy observation is the agreement between predicted and measured total liquid
saturations for the imbibition scanning path at a scaled capillary pressure of zero. For
this point, the porous medium was apparently saturated with NAPL and water. That is,
the sum of the NAPL, water and entrapped air contents equaled the porosity or flow
channel volume. At this 5-P point, air is predicted to be trapped by advancing air-
water interfaces in a two-phase fluid system and by advancing air-NAPL interfaces in a
three-phase fluid system. The entrapped air saturation that resulted from air-water
interfaces was predicted to be 0.015, and the entrapped air saturation that resulted from
air-NAPL interfaces was predicted to be 0.137. The sum of these amounts differed from
the measured amount of trapped air in the glass beads by a saturation of 0.002.
A similar water saturation path history, as measured in porous medium 1, is shown for
porous medium 2 in Figure 3.24. A closed scanning path that developed from the main
drainage 5-P branch was measured. Agreement between experimental data and model
predictions is very close, particularly for the scanning paths. Figure 3.25 shows the
corresponding experimental and predicted total liquid saturation path histories. Again,
very close agreement was found between predicted and experimental total liquid
saturation path histories. Model predictions in Figures 3.24 and 3.25 were generated by
best-fitting the model to the experimental data.
Although the same three-phase fluid system was used for porous medium 1 and 2, the
best-fit scaling factors are different. These differences may reflect contact angle effects
or experimental uncertainty. Theoretically, contact angles in dissimilar porous media
could be different. The observed differences in the best-fit scaling factors, 0OW and /?
00'
103
-------
for the two porous media suggest that estimating the scaling factors via ratios of
interfacial tensions may not always yield satisfactory results. Previous investigations by
Lenhard and Parker (1987, 1988) for drainage saturation paths, however, have shown
that interfacial tension ratios provides a good estimate of the best-fit scaling factors.
Table 3.4. Model parameters used to predict experimental water and total liquid
saturation path histories for porous media 1 and 2.
Parameters Porous Medium 1 Porous Medium 2
Da 0.022 cirf; 0
*a 0.032 cm-' 0
n 10.70 5
5m 0.22 0
P.. 2.39 1
/*, 2.05 2
/5arai" 0.06 0
/5arao 0.16 0
75erei" 0.15 0
o o _
» /\,
I
S60~
D 50 -
^
UJ
1 40 -
EC
< 30 -
_J
M 20 -
a
CJ
10 -
D
UJ
_J n
«
N
*\ »t
tV^t*^
VV L^^*^*^^
V \ »4c\ ^'^''^k-.
N. \ V''Sv\* ^^XB
>^ l^i ^^n
>«»4fc '
WATER SATURATION
PATH HISTORY
I
j>-
D 50 -
^
LU
1 40 -
-^
EC
< 30 -
_)
fT 20 -
CL
i U
10 -
O
> Lb
i 1 0 -
.042 cirf;
.072 an1
.86
.10
.93
.30
.19
.14
.09
\
\x\
\\w^_
\ NV v\r^^c^^
\ ^ ^A ^t&^afc^.
\/. "SLv >1>:^*
X. \V
-------
o
(u 50
E
U 50
40 -
LJ
I
£30-
_J 20 -
a
u 10 -|
D '
111
d 0
TOTAL LIQUID
SATURATION PATH
HISTORY
CJ
en
0.0 0.2 0.4
0.6 O.B
1.0
TOTAL LIQUID SATURATION
Fig. 3.23 Comparison of experimental (dashed lines) and predicted (solid line) total
liquid saturation path histories for porous medium 1 when the model was
calibrated by best-fitting the model to the experimental data.
cv 4B
I
E
U 40
LU
I
32 -
24 -
16
D.
CJ
D
UJ
_l
U
a -
WATER SATURATION
PATH HISTORY
0.0 0.2 0.4 0.6 O.B 1.0
WATER SATURATION
Fig. 3.24 Comparison of experimental (dashed lines) and predicted (solid line) water
saturation path histories for porous medium 2 when the model was calibrated
by best-fitting the model to the experimental data.
105
-------
a- 30
I
E
U 25
20 -
tu
I
I1""
_J
-J 10 -
M
D.
0 5 H
Q
ID
_J
TOTAL LIQUID
SATURATION PATH
HISTORY
U
in
0.0 0.2 0.<
0.6 O.B
1.0
TOTAL LIQUID SATURATION'
Fig. 3.25 Comparison of experimental (dashed lines) and predicted (solid line) total
liquid saturation path histories for porous medium 2 when the model was
calibrated by best-fitting the model to the experimental data.
8.8.2 Two-phase dynamic measurements
Experimental methods. A one dimensional air-water flow experiment was conducted to
measure water content and pressure distributions over time that resulted from a
fluctuating water table scenario. The objective of this exercise was to compare measured
hysteretic fluid distributions that included effects of nonwetting fluid entrapment to
results of numerical simulations of the experiments that incorporated the hysteretic k-S-
P model for the two-phase case.
The porous medium employed in the experiment was an unconsolidated sandy material
comprising approximately 97.5, 0.8 and 1.7 percent sand, silt and clay sized particles
(USDA), respectively. The sand was packed in a glass flow cell with a cross sectional
area of approximately 39 cm2 (6 cm x 6.5 cm) and 1 m in height. A 10-cm-thick gravel
layer underlaid 72 cm of packed sand. Porous ceramic tensiometers connected to
pressure transducers were inserted into the packed column through access ports located
106
-------
10 cm apart along the flow cell. Ceramic tensiometers were cylindrical with
approximate dimensions of 3 cm long by 5 mm in diameter. Pressure transducers,
manufactured by Omega Engineering (model 240), were constructed with silicon
diaphragms and possessed a range of ± 17 kPa or ± 175 cm of water equivalent height.
Transducers were accurate to within 1.5% and the response time was 1 ms. Seven access
ports were outfitted with tensiometers. Drainage and imbibition of water from and into
the flow column were accommodated by an outlet at the cell base. The sand-packed
flow column was flushed with C02 before saturating with water to minimize gas
entrapment on the inital water wetting of the dry sand. Initial condition for the
experiment, therefore, was assumed to be a water saturated porous medium.
A gamma radiation attenuation system was employed to nondestructively measure
water saturations. Gamma radiation measurements were taken within 1 to 2 cm of the
*
ceramic tensiometers. This was done to avoid gamma radiation attenuation by the
ceramic tensiometers and to avoid conducting measurements in areas that might be
disturbed (i.e., possessing a different soil particle arrangement) as a result of inserting
the tensiometer tips into the packed column. The radiation collimator was circular with
a 6.6-mm diameter. By moving two interconnected parallel platforms, one supporting
the radiation sources and the other supporting the gamma detector, to desired locations
with stepper motors under computer control, the gamma radiation beam could be
returned repeatedly to any elevation. The seven measurement elevations were 70, 60,
50, 40, 30, 20 and 10 cm, where 0 cm corresponds to the base of the packed sand.
The gamma system was calibrated with a procedure reported by Lenhard et al. (1988).
A sectioned small Plexiglas cell and repeated measurements of 10-min durations were
used to determine water and soil attenuation coefficients. In all calibration
measurements, 10 min counting times were employed. Any gamma attenuation caused
by air was neglected. Flow cell path lengths for each measurement position were
determined via empty and water filled gamma counts of the flow cell and by direct
measurements with calipers. Close agreement was found between path lengths
determined from the gamma attenuation measurements and those that were physically
measured with calipers.
Water contents were determined from the relationship
107
-------
nwpwewx) (3.32a)
where
V = I0exP(-n,P,x) (3.32b)
in which //, and nw are soil and water gamma attenuation coefficients, respectively, pt
and pw are bulk soil and aqueous phase densities, 6^ is the volumetric water content, i
is the path length, and 70 is the average empty cell gamma count of each measurement
location. Bulk densities employed in (3.32) for calculating water contents were
physically measured by sectioning the soil column at the end of the experiment.
Differences between soil bulk densities as determined by gamma counts of the dry soil
column prior to the experiment and those physically measured after the experiment
were less than 0.04 g cm'3, except for the uppermost two positions where differences
were larger. The average soil bulk density of the packed column was 1.70 g cm"3.
Because of the sandy nature of the porous medium, changes in bulk densities resulting
from water saturation path reversals were not expected (i.e., porous medium was
assumed to be rigid). Gamma counting periods during the transient experiments were
90 s and all gamma counts were corrected for resolving time and Compton scattering
following methods reported by Stillwater and Klute (1988) and Lenhard et al. (1988).
After initially wetting the dry soil column, the water table was positioned at 72 cm,
which corresponded with the soil surface, for approximately 12 h. At the onset of the
experiment (i.e., t = Oh), the water table was lowered 5 cm, and every 10 min
thereafter it was lowered an additional 5 cm until the water table reached an elevation
of 7 cm at t = 2 h, where it remained stationary for an hour. At t = 3 h, the water
table was raised 5 cm every 10 min to a final elevation of 42 cm and remained
stationary there for an hour until t = 5 h. All changes in water table elevations
proceeded at the rate of 5 cm every 10 min. The first lowering and raising of the water
table should produce the main water drainage and primary water imbibition scanning
paths, respectively. An internal water drying scanning path was then generated by
lowering the water table at the prescribed rate from 42 cm at t = 5 h to 17 cm where
again the water table remained stationary for an hour until t = 6.67 h. The final
saturation path, a water imbibition scanning path, was produced by raising the water
table elevation from 17 cm back to the soil surface at 72 cm, which was attained at t
= 8.33 h. In this last path as the water table was raised past 42 cm elevation, all
internal scanning S-P loops should have closed and, thereafter, primary water
108
-------
imbibition scanning paths should have been followed. A distinct S-P path history should
have developed for each measurement elevation. The rate of lowering and raising the
water table was arbitrary and was imposed, in lieu of more abrupt elevational changes,
to minimize potential numerical oscillations in the simulations that could be induced by
large, instantaneous changes in boundary conditions. Water content and pressure
measurements were conducted until t = 10 h.
The upper water boundary condition as well as the lower air phase boundary condition
were no-flow. Evaporation from the upper boundary was minimized by covering the
flow cell with plastic that contained numerous small holes to ensure that the air
pressure within the flow cell would be atmospheric. The upper air phase boundary
condition, therefore, was constant at atmospheric pressure. The lower water boundary
conditions were constant pressures that corresponded to the changing water table
elevation scenario described above.
The experiment was simulated using the numerical model described in Chapter 2. A one
dimensional 72 cm vertical homogeneous soil section was discretized in 1 cm intervals to
model numerically. Although the multiphase flow code was developed to simulate fluid
flow in air-oil-water systems, two-phase air-water flow can be simulated by regulating
the oil phase boundary conditions such that the system remains oil free. The code
assumes that the air phase, when present, is at atmospheric pressure.
To calibrate the k-S-P model, the drainage van Genuchten parameters, do, n and Sm,
were fit via a nonlinear regression algorithm to transient S-P measurements taken
during the first drainage sequence (i.e., main drainage). Two sets of parameters were
generated and five different simulations of the experiment were conducted. Parameter
Set 1 was obtained by fitting da, n and Sm to the main S-P drainage data and
Parameter Set 2 was obtained by fitting only da and n to the main drainage data
assuming Sm = 0. In Simulation 1, the full hysteretic model was employed using
Parameter Set 1. In Simulation 2, a simplified hysteretic model was employed that
accounts only for nonwetting fluid entrapment effects in k-S-P relations (i.e., no
hysteresis in apparent saturation-capillary pressure relations). We refer to the latter
analyses as "fluid entrapment only" simulations. In Simulation 3, hysteresis was
disregarded altogether (i.e., no apparent saturation hysteresis and no fluid entrapment
effects). Simulation 4 was a variant of Simulation 1 in that Parameter Set 2 was used
109
-------
and Simluation 5 was identical to Simulation 4, except the tortuosity term in Mualem's
(1976) relative permeability integrals was changed from Sw°-5 to Sw2, which conforms to
experimental work conducted by Burdine (1953) and Corey (1954).
Figure 3.27 shows the resulting best fit main drainage van Genuchten function with
experimental data of the measurement elevations. The solid line is for Set 1 parameters
and the broken line is for Set 2 parameters. Although there is scatter in the data, the
soil column was assumed to be homogeneous and isotropic for the simulations. The
residual air saturation of the main water imbibition branch used in simulating the
experiment was estimated to be 0.25, which was based on results observed in earlier
studies of similar porous media (Parker and Lenhard, 1987). The saturated water
conductivity was measured to be 119 cm h"1 in the packed flow column by the falling
pressure head method before initiating the experiment. The porosity (i.e., = 0.36)
corresponded to the average bulk density of all measurement elevations. Parameters
used in the 5 simulations are listed in Table 3.5. The hysteretic k-S-P routines can be
employed to simulate nonhysteretic flow (i.e., Simulation 3) by equating the a for main
imbibition S-P relations with the a for main drainage relations ('a = da) and using a
negligible quantity for '£,. (i.e., 10'^). A nonzero number is needed for *Sm or (4) will be
undefined.
Table 3.5. Model parameters for simulations 1-5.
Parameters
"a
n
sn
X
*
Kftu
Parameter Set 1.
0.042 cm'1
5.25
0.17
0.25
0.36
119 cm h*1
Parameter Set 2
0.039 cm'1
4.02
0.0
0.25
0.36
119 cm h'1
110
-------
ID
cr
D
cn
w
UJ
a:
a
a
<
u
o i
APPARENT SATURATION
Fig. 3.26 Apparent saturation-capillary pressure relations showing main drainage and
imbibition branches and a scanning path scenario.
Ill
-------
E
u
UJ
I
LJ
o:
D
C/)
tn
UJ
DC
a.
oc
LU
H-
-10 -
0.0 0.2 0.4 O.B 0.8 1.0
WATER SATURATION
Fig. 3.27 Comparison of measured water saturations (symbols) as a function of air-
water capillary head and best fit main drainage S-P relations for 5m > 0
(solid line) and for 5m = 0 (broken line) where the closed diamond, closed
square, closed triangle, open inverted triangle, open circle, open triangle, and
open squares are the symbols used for S-P relations measured at the 70, 60,
50, 40, 30, 20 and 10 cm measurement locations, respectively.
112
-------
1 .0
1 .0
z
2 o.e H
I-
in 0.4 -
DC
0.0
4 6
TIME (h)
i .0
z
'
o . e H
<
en 0.4 -
cr
UJ
0.2 -
0.0
10
1 .0
z
2 O.B -j
50.
<
U5 0 . 4 -
o:
0.0
1.0
z
2 O.B 1
0.4 -
« -, J
0.2 -
0.0
TIME (h)
4 B
TIME (h)
4 6
TIME (h)
10
10
:o
Fig. 3.28 Comparison of measured (solid symbols) and simulated full hysteretic (solid
lines), fluid entrapment only (broken lines with equal sized segments) and
nonhysteretic (broken lines with alternating long and short segments) water
saturation distributions at measurement locations of a) 70 cm, b) 60 cm, c) 50
cm, d) 40 cm and e) 30 cm.
113
-------
-55
E
u
-45 -
HI
X
HI
cr
D
in
en
UJ
tr
D.
oc
UJ
i-
-35 i
-25 -
-15 -
-5 -
MEASUREMENT
ELEVATION
70 em
0.0 0.2 0.4 0.6 O.B 1.0
WATER SATURATION
E -45
U
HI
-30 -
UJ
EC
D -15
in
in
UJ
cc
a o
EC
UJ
15
MEASUREMENT
ELEVATION
50 cm
0.0 0.2 0.4 O.B O.B 1.0
WATER SATURATION
-40
E
U
-30 -
UJ -20
ID
CC -10
in
in
UJ o
cc
a
cc
UJ
10 -
20
MEASUREMENT
ELEVATION
50 cm
-40
0.0 0.2 0.4 0.6 O.B
WATER SATURATION
i .0
E
u
D -20
<
UJ
I
UJ
DC
in
in
UJ
cr
a.
cc
UJ
0 --
20 -
40
MEASUREMENT
ELEVATION
40 cm
0.0 0.2 0.4 0.5 O.B i.D
WATER SATURATION
Fig. 3.29 Comparison of measured (solid symbols) and simulated hysteretic (solid lines)
S-P histories at measurement locations of a) 70 cm, b) 60 cm, c) 50 cm and d)
40 cm.
114
-------
«JU -
"E
u
-30
D
UJ -20 -
X
UJ
H -10 -
D
cn
c/i
UJ o .
DC °
a
oc 10 -
UJ
i-
S pn -
.
*V \ >x
^^t N
*N 1
^ ^
t I
*
I *
t
0.0 0.2 0.4 0.6 0.6 1.0
WATER SATURATION
Fig. 3.30 Measured 5-P history of the 50 cm measurement location showing a closed
internal scanning loop where broken lines are employed to connect successive
5-P measurements.
z
2 o.e -
i-
U) 0.4 -
cc
t- °-2 "
2
0
*L
V MEASUREMENT
I ELEVATION
\ 70 en
\
V
w .
2 4 6 E
TIME (h)
['£"*'"?
*
*
) 1C
1.0
go.aH
U) 0.4 -
UJ A _ ,
,_ 0.2 -
0.0
MEASUREMENT
ELEVATION
40 cm
4 6
TIME (h)
1C
Fig. 3.31 Comparison of measured (solid symbols) and simulated full hysteretic water
saturation distributions using Parameter Set 1 (solid line), Parameter Set 2
(broken line with equal sized segments) and Parameter Set 2 with a
tortuosity modification in the k-S relations (broken line with alternating long
and short segments) for measurement locations of a) 70 cm and b) 40 cm.
115
-------
Results and discussion. Measured and predicted water saturations are compared in
Figure 3.28 for the different measurement locations. Experimental water saturations are
shown as diamond symbols. Full hysteretic, fluid entrapment only and nonhysteretic
simulations are shown. Results of Simulation 1, the full hysteresis model, are shown as
solid lines, results of Simulation 2, the fluid entrapment only option (FT), are shown as
equal spaced broken lines, and results of Simulation 3, the nonhysteretic simulation
(NH) that employs only k-S-P relations of the main drainage branch, are shown as
broken lines with alternating long and short segments. In these three simulations,
Parameter Set 1 was employed. Very close agreement was found between experimental
and predicted water saturations at all measurement locations when the full hysteretic
constitutive relations were employed (i.e., solid lines - Simulation 1). Predicted water
saturations using the fluid entrapment only and nonhysteretic k-S-P relations, however,
deviated significantly from measured water saturations of the scanning paths.
Note in Figure 3.28 how poorly nonhysteretic predictions compare with water saturation
measurements of the scanning paths. At the 70- and 60-cm positions (Figures 3.28a and
b), measured water pressure heads decreased approximately 10 and 20 cm, respectively,
during the first raising of the water table to an elevation of 42 cm, yet measured and
predicted full hysteretic water saturations increased only slightly. The reason for this
phenomena will be explained later. Predicted saturations in the fluid entrapment only
and nonhysteretic simulations, however, increased significantly during the same period.
For measurement locations closer to the water table, there were only relatively minor
differences among experimental water saturations and predicted full hysteretic and fluid
entrapment only simulations. The nonhysteretic simulations, however, poorly matched
the measured water saturations at all elevations. Furthermore, the early increase in
water saturations predicted by the nonhysteretic k-S-P relations for the final wetting
scanning path (i.e., t = 6.7 to 8.3 h) in Figures 3.28a-c may suggest that effects of
nonwetting fluid entrapment on wetting phase relative permeabilities may have
significant consequences. For the remaining two measurement locations that are not
shown in Figure 3.28 (i.e., at elevations of 10 and 20 cm), differences between measured
and simulated water saturations are generally less than 0.04. At these positions the soil
was nearly water saturated at all times.
116
-------
Finally, note that for the 40-cm-measurement elevation (Figure 3.28d) experimental and
predicted water saturations are constant when the water table is positioned above this
location. During these periods, constant water saturations infer constant entrapped air
saturations. This occurred twice since the maximum water table elevations attained in
the first and second water imbibition paths were 42 and 72 cm, respectively, and the
minimum water table elevations attained for the first and second drainage sequences
were 7 and 17 cm, respectively. All measurement locations between 17 and 42 cm
should, therefore, exhibit plateaus in water and entrapped air saturations during periods
in which the water table is positioned above those elevations. Unfortunately, however,
only the measurement elevation at 40 cm (Figure 3.28d) drained sufficiently during the
second drainage sequence to exhibit recurring plateaus. The pore size distribution of the
porous medium (i.e., height of water saturated fringe above the water table) was such
that no water drained during the second drainage sequence at the 20-cm-measurement
location and only an insignificant amount drained at the 30-cm location. These distinct
and equal plateaus observed at the 40-cm elevation suggest that the entrapped air
saturation corresponding to an apparent saturation of unity may be a constant for a
given saturation path history, and can be reproduced following cycles of wetting and
drying provided the saturation path reversal from the main water drainage branch
remains unchanged. Similar results have been reported recently by Stonestrom and
Rubin (1989) where they observed a reproducible functional dependence of entrapped air
content on water saturations.
The constant water saturations measured during periods when the water pressure is
greater than atmospheric infers that there is negligible movement of trapped air out of
the porous medium for conditions imposed during the experiment. This is to be
expected considering that water was not continually flowing and was probably saturated
with respect to air at atmospheric pressure. Water-occluded air bubbles in porous media
may move in response to increasing the viscous forces that tend to push the air bubble
through the porous medium (i.e., water pressure gradient), to decreasing the interfacial
forces that hold the air bubbles stationary, and to dissolution of the air in the aqueous
phase and subsequent transport by convective and dispersive mechanisms. Researchers
in the petroleum industry investigating various techniques to enhance oil recovery in
petroleum reservoirs have reported a correlation between the amount of entrapped
nonwetting fluid that remains in porous media and a dimensionless group called a
capillary number, which is defined as a ratio of viscous forces to interfacial forces, and
117
-------
can be expressed by
Nc = fy^f (3.33)
where qw is the Darcian water flux, rjw is the water absolute viscosity, and 10~5, the amount of mobilized
entrapped nonwetting fluid increased rapidly with increasing -Afc. For air-water fluid
systems in porous media, the Darcian water velocity needs to exceed approximately 2.6
m h , which is a substantial flow rate, in order to mobilize trapped air bubbles
according to (3.33). For boundary conditions imposed during our experiment such a
large flow rate would not be attained.
Dissolution of air bubbles in water has been studied by Bloomsburg and Corey (1964)
and others based on theoretical considerations. The time required for a nonwetting fluid
ganglion to dissolve in water is a function of the pressure difference between the
nonwetting fluid and water and the nonwetting fluid solute concentration in the aqueous
phase. For small ganglia, the dissolution time is orders of magnitude faster than for
significantly larger ganglia, however, as shown by Bloomsburg and Corey (1964) the
time period for significant amounts of entrapped air to dissolve in the aqueous phase
may be on the order of weeks to months. Considering the time required for entrapped
air bubbles to dissolve in water and the Darcy water velocities required to dislodge
entrapped bubbles, the modeling assumption that entrapped nonwetting ganglia (i.e.,
entrapped air bubbles) are immobile is reasonable, particularly for a fluctuating water
table scenario where the duration between fluctuations is less than several weeks.
Incorporating dissolution of entrapped air in k-S-P constitutive models may yield more
accurate fluid distribution predictions. However, markedly increased computational
effort will be required since a transport equation for dissolved gas must be solved.
Measured and simulated water S-P histories are compared in Figure 3.29. Experimental
S-P relations are shown as diamond symbols and simulated S-P relations employing the
full hysteretic constitutive model (i.e., Simulation 1) are shown as solid lines. There are
several experimental features that are worthy of attention in Figure 3.29. First, note
that scanning S-P loops formed by the first rising of the water table and subsequent
water table lowering (i.e., t = 3 to 6.67 h) become more pronounced as the reversal
118
-------
water saturation from the main drainage branch increases. At the 70-cm-measurement
location (Figure 3.29a), the lowest water pressure head measured was approximately -45
cm, which occurred during the first lowering of the water table (i.e., main drainage
branch). During the first rising of the water table from 7 to 42 cm, the water pressure
heads increased to approximately -34 cm; however, the slope of the S-P scanning curve
is quite steep in this region and a negligible increase in water saturation is predicted,
which matched the observed behavior. This is why the predicted full hysteretic water
saturation distributions in Figures 3.29a and b show little response to the first table
elevation increase at t = 3 to 5 h. Slopes of water imbibition 5-P curves will always be
steeper than the drainage curves at lower water saturations (i.e., fo > *a). As the water
table was lowered for a second time from 42 to 17 cm, the measured water pressure
heads approached -45 cm and water contents very closely followed the saturation path
as for the initial wetting. During the final rising of the water table from 17 cm to the
soil surface, the water saturation path mirrored the path measured during the first
increase in water table elevations suggesting that the scanning S-P loops close.
Forming closed 5-P loops for cyclic changes in boundary conditions can be better
observed in Figure 3.30 where the experimental 5-P history is plotted for the 50-cm-
measurement location. Successive experimental measurements are connected by broken
lines in Figure 3.30. Note that a closed internal 5-P loop is formed during the cyclic
raising and lowering of the water table. Also note that during the second water table
rise (i.e., to the soil surface), the same 5-P path is followed as that during the first
increase in water table elevations. This closed loop behavior can be observed for all
measurement locations shown in Figure 3.29.
Predicted 5-P relations employing the full hysteretic model and Parameter Set 1
matched fairly well the experimental 5-P paths. Main features of the 5-P histories were
captured by the hysteretic constitutive relations. This is very encourging in that the
hysteretic model was calibrated only from main water drainage 5-P measurements and
an estimate of the maximum amount of entrapped air corresponding to the main water
imbibition 5-P branch. Water imbibition 5-P relations were predicted by employing an
a (i.e., *a) of the van Genuchten retention function to be twice that of the drainage o
(i.e., a) and using the same parameter n for imbibition as for drainage. This approach
greatly simplifies calibrating the hysteretic flow model and appears to perform
satisfactorily, at least for the experimental conditions employed in this study.
119
-------
Two additional simulations were conducted to evaluate consequences of assuming an
irreducible wetting fluid saturation, 5m, of zero. For Simulation 4, the full hysteretic
constitutive model was employed. However, it was calibrated to the main water
drainage S-P experimental data assuming STO = 0 (i.e., Parameter Set 2). For
Simulation 5, the full hysteretic constitutive model and Parameter Set 2 were
employed. However, the tortuosity term in Mualem's (1976) relative permeability
integrals was changed from Sw°-s to 5^*. That is, in Simulation 5 we used a slightly
different k-S model. Water distribution results from these two additional simulations are
shown in Figure 3.31 for the 70 and 40 cm measurement locations. Simulation 4 is
shown as broken lines with equal sized segments and Simulation 5 is shown as broken
lines with alternating long and short segments. For comparison purposes, Simulation 1,
which is the full hysteretic model with Parameter Set 1 (i.e., Sm > 0), is also shown as
solid lines. The parameter sets employed in the simulations are listed in Table 3.5.
In Figure 3.31a, it can be observed that Simulation 4 predicts lower water saturations
during drier periods than were measured or predicted by Simulation 1. This is primarily
attributed to imposed k-S relations. In Simulation 1, as water saturations approached
0.17 (i.e., Sm = 0.17), the water relative permeabilities were negligible. Hence, water
drainage essentially ceased as predicted water saturations approached 0.17. In
Simulation 4, k-S relations were such that water relative permeabilities would not
approach a negligible value until the water saturations were close to zero, which
permitted water to continually drain to lower water contents. It is clear from comparing
Simulation 1 and 4 results to the experimental data that an irreducible water
saturation, or a residual water saturation as is commonly refered to in soil physics and
hydrology literature, may be required to accurately model subsurface water behavior or
that the k-S relations of the van Genuchten-Mualem model may not be appropriate for
this soil. To further explore this supposition, we conducted Simulation 5, which used the
same parameters as Simulation 4, but with modified k-S relations. We changed the
tortuosity term in Mualem's (1976) relative permeability integrals to conform with
experimental work completed by Burdine (1953) and Corey (1954). The results of
Simulation 5 are shown in Figure 3.31a for the 70 cm measurement location where it
can be observed that the predicted water distribution compares favorably to the
experimental measurements. Furthermore, hysteretic S-P histories predicted by
120
-------
Simulation 5 compared favorably to the experimental data shown in Figure 3.29. At
some measurement locations, Simulation 5 predictions yielded better agreement with
the experimental data than did Simulation 1 predictions. However, at other locations
the reverse was true.
These results seem to suggest that employing an irreducible water saturation in
numerical modeling may not be necessary provided k-S relations axe described
adequately. Considering that Mualem (1976) obtained the tortuosity term in his relative
permeability integrals by best fitting to data from 45 soils in which recorded S-P and k-
S relations were available, it is not unreasonable to assume that this tortuosity factor
may depend on flow channel geometry and may be empirically correlated to other
retention parameters. Mualem (1976) noted that the tortuosity term may be porous
medium specific.
In Figure 3.31b there is better agreement among Simulations 1, 4 and 5. This is
primarily because at the 40 cm measurement location, measured and simulated water
saturations are significantly greater than 5m of Parameter Set 1 (i.e., 0.17) and predicted
k-S relations are similar for the different simulations. For higher water contents, i.e., 30,
20 and 10 cm measurement locations, Simulations 1, 4 and 5 predict almost identical
water saturation paths. Finally, we note that the coefficient of determination (R2 value)
obtained from Parameter Set 1 was 0.966 and that obtained for Parameter Set 2 was
0.963 (see Figure 3.27). Although this difference in R2 values is small, substantially
different predicted water saturation paths may be obtained as observed in Figure 3.31a.
It is not sufficient to accurately describe S-P relations; k-S relations must also be
accurately described.
3.8.8 Three-phase dynamic measurements
Experimental methods. A three-phase column experiment was conducted to evaluate
the hysteretic constitutive model with fluid entrapment effects. The experimetal setup
consisted of a 71 cm long soil column (6 cm x 6.5 cm) initially free of oil and the water
table located at the upper surface. Fluid saturations along the column were measured
using a dual-energy gamma radiation apparatus and liquid phase pressure measurements
were taken using ceramic tensiometers connected to pressure transducers. The
121
-------
experiment was conducted with time dependent boundary conditions that can be
summarized as follows:
Stage 1: At t = 0, the water table was lowered at a rate of 1 cm per minute until it
reached an elevation of 7 cm at t = 1 h. The flow system therafter was allowed to
redistribute for a period of 12 h.
Stage 2: At the end of stage 1, a slug of oil containing 90% Soltrol and 10% 1-
iodoheptane was allowed to infiltrate from the upper surface under a zero water
equivalent oil head until the accumulation was 250 cm3. At the end of the infiltration
period, the system was allowed to redistribute under natural hydraulic conditions until
t = 17.45 h.
Stage 5: The lower surface water table at the end of stage 1 was varied as shown in
Figure 3.33 and the experiment was conducted until t = 22.1 h.
Soil hydraulic properties of the medium and the fluid properties of the oil mixture are
given in Table 3.6. The experimental procedure described above was simulated using
the simplified hysteresis model described in Section 3.2. For the simulation, oil head at
the upper and lower surfaces was fixed at zero flux condition except during infiltration.
The water boundary conditions at the upper surface water head was fixed at zero flux
condition during the entire simulation.
122
-------
Table 3.6. Soil and fluid properties
Parameter
Value
n
Pro
Iro
POO
&
OW
ax
0.4
2.48
0
0.05 cm'1
121 cm h'1
0.83
3.62
2.85
1.7
0.25
75. O
6O. O
T3 45. O
0
0
3O. O
L
0
0. O
10 IS 20
Time Chrs)
Fig. 3.33 Water table elevation during column experiment.
123
-------
' U
D. a -
C
0
; O.B-,
0
L
D 0. 4 -
JJ
0
t/1
0. 2 -
0. O -
o
\
\ w"°" «"'-«'>"
j Oopth - 14 cm
\
\ r
" \ /* ! "
\ I*/
" * * i * /
m \ IB g ' *^ *
IW^L-^JM " I ^** * '. * S
~* ^ ^^f \
v J
i i 1 1
5 10 IS 20 25
Tima
i . U
o. e -
C
0
~ 0. 6 -
0
L
3 C. 4 -
4J
0
l/l
0. 2 -
0. 0 -
\
\ Vota,- ...u,nr ,
\ °t0" "*"-«">«
\ DoptM - 34 cm __
\ /'""^ '»"""
\ 1 I /
,_\ ^ /"""''.'"
>_ \«_ _ . /^ ^ ^
"
" * " J*
* «r "
*
10
15
Fig. 3.34 Predicted and observed water saturations for the column experiment. Solid
line and symbols indicate predicted and observed saturations, respectively.
1. U -
0. B
C
0
" 0. 6 -
0
L
3 0. 4 -
0
in
D. 2 -
0. 0 -
Oil otvi-otlon
C>pth - 14 cm
1
\ m
if fv
v r\. .
"^\^ v^-
"TT"" " . '
*
1 1
J. U -
o. e -
c
0
- 0.6-
0
L
D 0. 4 -
*J
0
in
o.z-
v.o-
Ol 1 ct.wi-otlon
Ooptn - 34 cm
»
\ , *
'K "
\^_^^^ /^"
" . «^~~^
i
i i
1O :S 2O 25 10 15 20 2!
Timo Chrs) Timo
Fig. 3.35 Predicted and observed oil saturations for the column experiment. Solid line
and symbols indicate predicted and observed saturations, respectively.
124
-------
Predicted water and oil saturations at 14 and 34 cm depths are shown in Figures 3.34
and 3.35, respectively. Predicted water saturations at the 14 cm depth indicate a drop
in saturation at approximately 12 h due to the advancing oil front. During the same
period, water saturation at the 34 cm depth increased due to the downward movement
of water associated with the advancing oil front and decreased thereafter for a short
period as the oil front reached this depth. Fluctuations in the water table caused smaller
changes in the water saturation at the 34 cm depth compared to the 14 cm depth since
during the most period of the simulation, the 34 cm depth remained liquid saturated.
Predicted oil saturations after 17 h at the 14 cm depth indicate less oil entrapment
compared to the 34 cm depth.
The results show good agreement between the observed and predicted fluid saturations.
The more pronounced peak in oil saturation as oil moves downward relative to the
observed measurements probably reflects the simplifications involved in treating the
saturation-capillary pressure relations. In particular, disregarding hysteresis in apparent
water saturation produces slightly sharper fronts than are observed. However, the
discrepencies diminish markedly as the distance from the source increases, indicating
that ignoring such effects will be of little consequence if predictions very near the source
are not of primary concern.
125
-------
4. EQUILIBRIUM-CONTROLLED MULTIPHASE TRANSPORT
4.1 Mathematical Model for Multiphase Transport
4-1.1 Governing equations
The mass conservation equations for water (w), organic liquid (o) and air (a), assuming
an incompressible porous medium, incompressible liquid phases and compressible gas
phase, may be written in summation convention for a two-dimensional Cartesian
domain as
dpaSa_ >
* §r = + a ( J
where tf is porosity, 5p is the p-phase saturation, x, (and !_,) are Cartesian spatial
coordinates (i,j = 1, 2), qp. is the Darcy velocity of phase p in the i-direction, pp is the
density of phase p, Rp is the net mass transfer per unit porous media volume into ( + )
or out of (-) phase p, and t is the time. Darcy velocities in the p-phase are defined by
<4-2>
»
where Kp.. is the p-phase conductivity tensor, hp = Pp/
-------
Boundary conditions may be stipulated as type-1 or type-2 as follows
hp(xt> t) = hpi(Xi, t) on Sl for t > 0 (4.4a)
-n = . *, * on $ f°r t> 0
Relationships between phase permeabilities, saturations and pressures and numerical
solution procedures for the flow equations have been discussed in Chapter 3.
To model component transport, continuity and mass flux equations for each
partitionable component in each phase must -be specified. Mass conservation of species
a in the p-phase requires that
dcnnsn dJOD.
+ Rap + 7OP (4.5)
where Cap is the concentration of the noninert a component in p-phase expressed as the
mass of a per phase volume [M L"3], Jop, is the mass flux density of a in p-phase per
porous media cross section in the i-direction [M L"2 T"1], Rap is the net mass transfer
rate per porous medium volume of species o into (-I-) or out of (-) the p-phase [M L"3
T"1], and 7op is the net production (-f) or decay (-) of a within phase p per porous
medium volume due to reactions within the p-phase [M L"3 T'1] described by
fop = ~ "op Cap (4-6)
where nap is an apparent first-order decay coefficient.
The mass flux density of component o in phase p due to convection, diffusion and
mechanical dispersion is described by
dC
Jopi = Cop qpi - t Sp Dapij -gg (4.7)
where Dopij is a dispersion tensor given by
127
-------
(4-8)
in which D^ is the molecular diffusion coefficient of o in the p-phase of the porous
medium, Dptyd is a mechanical dispersion coefficient, and gap is a nondilute solution
correction factor. For the case of transport of low solubility organic components, small
volume fractions of organic components in water and gaseous phases will occur and gQW
= gaa = 1 is assumed. For nonaqueous phase liquids, the phase composition may
reach 100% (e.g., single component organic liquid) at which point dispersive transport
becomes nonexistent. We accommodate nondilute solution diffusion-dispersion by
approximating the oil phase nonideal solution factor by
900 = 1 - C00/.Po (4.9)
Employing the tortuosity model of Millington and Quirk (1959) yields the expression for
*^op
J~> d\1 _ jl/3 C 7/3 r>o 1 A in\
Vap ~ ' Jp ISap (4-lUJ
where D^p is the diffusion coefficient of o in bulk p-phase. The mechanical dispersion
coefficient is shown by Bear (1972) to have the form
(4.n)
where AL and AT are longitudinal and transverse dispersivities [L], qpi and qpj are p-
phase Darcy velocities in the t and j directions, ~qp = /E^,2/1/2 is the absolute
magnitude of the p-phase velocity, and £,- is Kronecker's delta.
Combining the phase continuity equation and the mass flux equation for transport of
component o in the p-phase yields
Expanding the first and third terms in (4.12), employing the bulk p-phase continuity
equation (4.1), and assuming density derivative terms to be of second order importance
128
-------
within a given time-step yields
C dC<>p ,
~ 9" ~~
op
where the total phase mass transfer rate, Rp, it may be noted, is related to the
individual component mass transfer rates by
RP = £ Rap (4-14)
0 = 1
where n, denotes the number of "noninert" or partitionable species. In the present
context, we will use the term "inert" to refer to components of the NAPL phase which
are for practical purposes insoluble and nonvolatile. For the three fluid phase system,
(4.13) constitutes a system of three equations which we may write for the water phase
(p = w), the organic liquid phase (p = o) and the gas phase (p = a) as
(4.15b)
+ Raa - ( + ) Coa (4.15c)
To accommodate adsorption of a by the solid phase, an additional continuity equation
is required which may be written as
dCg, _ p r u -\R\
-0T ~ R°' ~ "« C<" (4-16)
where Cat is the solid phase concentration expressed as mass of adsorbed component a
per porous medium volume [M L"3], nat is the first-order decay coefficient (T"1) in the
solid phase and Rai is the mass transfer rate per porous medium volume to (+ ) or from
(-) the solid phase [M I/3 T1].
129
-------
4-1.2 Phase-summed equations for local equilibrium transport
Coupling between the phase transport equations arises due to the interphase transfer
terms. Explicit consideration of interphase transfer kinetics may often be justifiably
avoided by assuming phase transfer to be equilibrium controlled. We consider first the
case of linear partitioning and introduce the approximate thermodynamic relations
Cao = T00 Cow (4.17a)
Caa = roa Cow (4.17b)
Ca. = ra. Cow (4.i7c)
where roo is the equilibrium partition coefficient for species o between water and
organic liquid (Raoult's constant), Toa is the equilibrium partition coefficient between
water and gas (Henry's constant), and Tat is a dimensionless equilibrium partition
coefficient between water and solid phase. Using the equilibrium relations, we may
rewrite the phase transport equations in terms of a single-phase concentration. For
water-wet systems, it is logical to retain the water phase concentration, since water will
always be present in the system. Using (4.17) to eliminate oil, gas and solid phase
concentrations from (4.15) and summing the equations noting that
Raw + Rao + Ra. + Raa = 0 (4-18)
leads to the a-component phase-summed transport equation in terms of water phase
concentrations
*:
where
*« = *$«, + * SOTOO + * 50 roa + ra. (4.i9b)
Dafj = Sw Dawi; + 4 S0 Dooij roo + * 5a Doaij Taa (4.19c)
9«* = 9«i + 9oi roo + qai roa (4.19d)
D /? r /? r
* 1 T» I T» I T* I W i ^O CrO I fl QQ / A i r\ \
it ~ u + u 1 -T- ii 1 -T- 11 1 -I- U- .-"v -L "u (A 1 Qo i
'^^k ~^t*» ' ~f*n nn i *v*i rt/i i ~f*» * /v« i /i I /) f /) 1~*X^C f
130
-------
Note that (4.19) has the same form as the simple single-phase transport equation.
However, the coefficients represent pooled effects of transport in all phases. Note also
that interphase mass transfer terms occur in the phase summed equation only as a sum
over all components.
An alternative form of the phase-summed transport equation may be written in terms
of the oil phase species concentrations. The a-component phase-summed transport
equation in terms of oil phase concentrations has the form
**
** r> (A t)n.\
" "« GO° (4.20a)
where
i V'-'w i v-'o x oo i * as i A onv^
+ -P -(- f -I- TT (4.20b)
oo oo oo
i T^a^aaij a a i/, nr\\
-I- f (4.20c)
^00
**
(4.20d)
(4.20e)
The form of the oil-based equation is identical to that of the water-based phase-summed
equation. The water-based equation has the conceptual advantage that Caw is
physically meaningful at all locations in the porous medium, whereas Cao only has
physical meaning at locations where S0>0. The latter constraint, however, imposes no
mathematical difficulty since (4.17) can in any event be employed as a definition of Cao
as well as Caa whether or not the phase exists. Since most of the species mass
commonly resides in the oil phase for organic liquid spill simulations, the oil-based
equation has the advantage that mass balance accuracy of the solution may be more
readily controlled.
131
-------
When simulating two-phase liquid flow consisting of water and oil with constant gas
pressure, only diffusive transport is considered in the gas phase i.e., qa = 0 is
assumed for this case. It should be noted that even if gas pressure gradients have
negligible effect on liquid phase flow, this does not strictly imply zero gas flow. If fluid
saturations change in time, gas flow will occur. Furthermore, phase partitioning may
result in density gradients in the gas phase which can induce significant flow.
4-1.3 Initial and boundary conditions
Since the transport equation is written in the phase-summed form, it is necessary to
specify initial and boundary conditions in terms of a single concentration. In the
following, as well as in the program itself, we will specify initial conditions and type-1
boundary conditions in terms of water phase concentrations. If the oil concentration
form of the phase-summed transport equation is employed, oil phase concentrations are
internally computed from the specified water phase values and the equilibrium
relations. Type-3 or flux-type boundary conditions are defined in terms of phase-
summed fluxes. Initial conditions for component a are given as
Cow(Xi, t = 0} = COWo(x,) on R for t = 0 (4.21a)
where Caw (i,) represents the initial water phase concentration of component a at
O
location i,. Note that the concentrations of o in all other phases are implicit in (4.21a)
via the phase partition relations (4.17). Boundary conditions for the phase-summed
equation are stipulated by
CUM) = Cawi(xi,t) on 5: for t > 0 (4.21b)
= 0 on 52 for t> 0 (4.21c)
/ on S3 for O 0 (4.21d)
3
Equation (4.21b) describes a type-1 boundary condition on boundary region 5j with
specified o-species concentration Cawi(xitt) within the porous medium at the boundary.
Note again that stipulation of the water phase concentration simultaneously fixes the
organic liquid, gaseous and solid phase concentrations via the partition relations. Type-
1 boundary conditions for transport will be used rarely since direct control over porous
132
-------
media concentration seldom occurs. One case which may be relevant would be the
occurrence of a source of known concentration where diffusion is the dominant transport
mechanism (a type-3 condition would be used if significant convection occurred at the
boundary).
Equation (4.21c) describes type-2 boundaries with zero normal concentration gradient
specified on region S2. The zero gradient condition indicates zero normal dispersive flux
at the boundary. If normal fluid velocities are nonzero, transport across the boundary
will be permitted by convection. This boundary condition is commonly applied at
boundaries which have fluid efflux to permit constituent transport across the
boundaries. In the event that fluid fluxes are zero normal to the boundary, condition
(4.21c) stipulates zero component flux. Caution needs to be exercised in applying the
type-2 condition on boundaries with inward fluid velocities of uncontaminated fluid. If
the concentration at the boundary is nonzero due to presence of chemical near the
surface, the type-2 condition will stipulate an inward mass flux density equal to the
water phase concentration at the boundary times the phase-summed velocity, g*, which
is incorrect. In such a case, a type-3 boundary condition with zero influent
concentration should be specified to force the desired zero-flux condition. The type-2
zero gradient condition is the default condition.
Equation (4.21d) is a type-3 boundary condition for use on boundary regions, 53, which
have inward NAPL phase fluxes predicated on the assumption that contaminant enters
the system solely via inflow of nonaqueous phase liquid. Here, qon is the normal Darcy
velocity of the oil phase at the boundary and C£{,fe is the equivalent water phase
concentration of species a in the influent oil phase defined by C£^e = C^/T^ where
Cof/ is the actual concentration of a in the influent oil. Note that applying the
condition C^e = 0 corresponds to imposing zero mass flux of component a on the
boundary segment. In applying the type-3 boundary condition during periods of oil
infiltration, it is important to ensure consistency in the specified oil phase
concentrations such that
(4.22)
where Cao (a = 1, . . . , n,) is the oil phase concentration expressed as mass of
component per phase volume for partitionable component a and Cj is the oil phase
133
-------
concentration of "inert" or nonpartitionable components. Equivalent!}-, in terms of
mass fractions the condition may be expressed as
// + £ /00 = 1 (4-23)
0 = 1
where fao and fj are mass fractions of partitionable and inert components, respectively,
in the oil phase. Mass fractions and phase concentrations are related by
Cap = fap Pp (4.24)
where pp is the bulk phase density. Simple mixing theory permits component mass
fractions and volume fractions to be related as
Cr n»
77 + £
* (4-25)
0=1
where pa and pj are the densities of pure noninert components and inert oil,
respectively. The oil phase density may be related to phase composition by
(4-26)
fl/Pj + £
0 = 1
In applying the type-3 oil infiltration boundary condition, influent oil phase
concentrations, C£^e must be specified such that (4.22) - (4.26) are obeyed.
4.2 Numerical Model Description
4-2.1 Central solution approach
Fluid flow equations are highly coupled with each other due to interdependence of fluid
permeabilities, saturations and pressures mandating a simultaneous solution approach as
discussed in the previous chapter. The flow equations are also coupled with the
transport equations via the occurrence of interphase mass transfer terms and through
the dependence of various coefficients in the flow equations on fluid composition e.0.,
phase density, viscosity, scaling factors, etc. Since mass transfer rates are necessarily
134
-------
small for the case of low solubility organic fluids and changes in fluid properties over
short time periods will accordingly also be small, the dependence of the flow equations
on transport will be very mild over short time spans permitting computational
decoupling. Over extended periods of time, cumulative phase mass transfer may have
significant effects on flow as dissolution and volatilization deplete the nonaqueous phase
liquid requiring some means of updating phase transfer rates as well as compositionally
sensitive fluid properties. We will consider fluid compositional effects on phase density
but disregards effects on other coefficients. Interphase mass transfer rates and phase
densities are time-lagged in the flow solution and updated at the end of each time-step
after solution of the transport equations.
The transport equations are highly dependent on the solution of the flow equations due
to the occurrence of fluid velocity and phase saturation terms directly in the transport
equations and in the functional forms for the dispersion coefficients. Thus, the flow
equations must be solved concurrently with or prior to evaluating transport. Due to the
weak back-coupling, the most opportunistic approach is to solve the transport equations
serially with the flow equations. Furthermore, since the individual component transport
equations axe weakly coupled with each other in the present model due to the
assumption that components do not interact, the transport equations for different
components may be solved serially in arbitrary sequence.
The solution approach for solving the coupled multiphase flow and multicomponent
transport problem with equilibrium mass transfer is as follows:
1. Solve the fluid flow equations simultaneously for the current time-step using time-
lagged phase densities and interphase mass transfer rates.
2. Solve the phase-summed transport equation using the previous time-step phase
densities and interphase mass transfer rates.
3. Back-calculate new interphase mass transfer rates and update phase densities for
the current time-step.
4. Proceed to the next time-step.
135
-------
4-2.2 Finite, dement formulation
The phase-summed transport equations given by (4.19) or (4.20) are solved using an
upstream-weighted Galerkin finite element method with linear rectangular elements.
We will present the solution for (4.19). The formulation in terms of (4.20) is
fundamentally identical. Employing linear shape functions, Nj, the concentration
C*OU)(ij,<) can be approximated by
4
Caw(xt,t) = £ Nfa) CowJ(t) (4.27)
where CQWj denotes the concentration at node J and c and rj are local coordinates.
Detailed information pertaining to the upstream weighting functions have been
discussed by Kaluarachchi and Parker (1989). The Galerkin finite element
approximation of (4.19) over the entire domain R can be written as
J NjL(Caw)dR = 0 (4.28)
R
which leads to the set of equations given by
[Pi iCkw} + M { *} = {R} (4.29)
Using the type-3 boundary condition given by (4.21d) and disregarding cross derivative
terms, the matrices [P], [M\ and {R} are defined by
fj N
£ [NI»ZNJdR - £ <9*> / WjN
= lJZe c = 1 S
+ NI»NJdR - <9> WjNjdS (4.30a)
MIJ = E ^/ *Z NJ dR (4-30b)
« = 1 Re
N> r
RI = - E rao C^i/e / W7 dS (4.30c)
e = l
136
-------
where N is the total number of elements, Nt is the total number of boundary elements,
WT is the asymmetric upstream weighting function, is the average normal oil
flux along the surface 5 where contaminant enters the system and is the average
normal phase-summed velocity on surface 5. If an evaporative boundary is present, the
term in (4.30c) will have an additional term (-.D^1/ roo) to account for the
atmospheric boundary layer and the terms containing C^J,' in (4.30c) will vanish.
Treatment of the integrals in (4.30) employs the method of influence coefficients
described in the flow analysis and will not be repeated here. Handling of the convective
term in (4.30a) using the influence coefficient method has not been previously described
and will be discussed below. Phase-summed velocities within an element are represented
by
4
fcf (**-,«) = E *M foSW (4-31)
9=1
where qQf*p is the phase-summed velocity at the evaluation points in the element, which
may be at the nodes, at Gauss points or elsewhere as controlled by the parameter GP
(see discussion in flow analysis). Using equation (4.31) on the convective term of
(4.30a) and using linear upstream weighting shape functions (Kaluarachchi and Parker,
1989), it is possible to show that
dR =
J l "" "*!
e = lR*
}
/
+ [Hv9}a} (4-32)
where m and d are the width and height of a rectangular element, x and y axe the
coordinate axes and matrices [H\ with superscripts xg and yg refer to the contributions
from standard linear shape functions while matrices with superscripts xu and yu refer to
contributions from upstream weighting coefficients. Details of matrices [H\ are
described in detail by Kaluarachchi and Parker (1989) and will not be repeated here.
The surface integrals in (4.30a) and (4.30c) can be described by
137
-------
= (4±3w) if I=J (4.33a)
= ^ (2 ± 3w) if 7?fc 7
and
dS = (1 ± w) (4.33b)
where £ is the length of the side and w is the upstream weighting coefficient. The
appropriate sign of the integrals given in (4.33) will be determined by the side of the
element exposed to the boundary condition.
Time integration of (4.29) is performed using a standard finite difference approximation.
The resulting final set of equations can be written as
(4<34)
where k and k + 1 refer to previous and current time-steps and 6 is a time weighting
factor. In the present work, we use a Crank- Nicolson time-stepping scheme
corresponding to 6 = 0.5 to provide an unconditionally stable second-order accurate
temporal approximation. Time-step size is varied under program control by a time
incremental factor which increases or decreases time-step size depending on the number
of iterations required for solution of the flow problem. The time-step size for solution of
the transport problem is taken as identical to that for the flow problem. If equilibrium
mass transfer is assumed the transport equation is linear and no iteration of the solution
is required. For the case of nonequilibrium mass transfer, nonlinearities are handled by
a simple Picard iteration scheme with apparent partition coefficients, mass transfer
rates and phase densities updated at each iteration.
4-2.S Interphase mass transfer and density updating
After solution of the phase-summed transport equations, the interphase mass transfer
terms Rap must be updated. This is performed after each iteration of the transport
equations in the case of nonequilibrium mass transfer. For the equilibrium case, the
138
-------
transport model is linear and no iteration is required. Hence, for the latter case,
updating is performed only once at the end of each time-step. To compute the mass
transfer rates, the concentrations of each component in the oil, gas and solid phases are
first calculated from the partition coefficients and from the control phase concentrations
(Caw for the case of eq. 4.19 or Cao for eq. 4.20). The mass transfer rate terms, Rap.
are computed by back-substitution into the component transport equations for each
phase (eqs. 4.15 and 4.16) using a finite difference approximation. Summing over the
components yields the desired Rp terms for each phase for use in the subsequent
solution of the flow and phase-summed transport equations.
For equilibrium mass transfer problems, effects of mass transfer terms on the flow and
transport equation solutions at any given time-step will be quite small. However,
cumulative effects can be very important since these terms will control the long term
removal of NAPL associated with dissolution and/or volatilization. The accuracy of
interphase mass transfer calculations will generally be most troublesome during periods
of highly transient NAPL flow. During such periods, it may be preferable to suppress
mass transfer updating to avoid numerical instability as well as to reduce computational
effort and improve accuracy.
With increasing simulation time, the density of each phase will change depending on the
net mass of each component leaving or entering the phase. Although density
derivatives were neglected in the equation development, on the assumption that these
terms are small over a given time-step, cumulative changes in density cannot be ignored
since over long times these can accumulate and become quite large in magnitude. To
accommodate these effects, liquid phase densities are updated at the end of each time-
step as follows
nt n n>
O f 1 ^"^ ^QtU 1 i V* /t /1 n v \
Pw - Pw I * - 2^ ~f~ J + Z- C (4.35a)
a=l a=l
n, Q n,
P0 = Pj [ 1 ~ ^_- -p22 ] "I" 2_s ^oo (4.35b)
o=l o=l
where nt denotes the number of "noninert" or partitionable species, p°w is the density of
uncontaminated water, pa is the density of pure a-component, pr is the density of
"inert" oil components, and Ma is the molecular weight of component a. Density ratios,
139
-------
pru), pro and pra, are computed as the ratio of pw, p0 and pa to p°.
For the gas phase, density depends on both pressure and composition. Pressure
dependence is treated as a nonlinear effect in the gas flow equation according to
Pa = *ha + Pa" + Pa' (4'36)
where A is the gas compressibility taken to be 1.17 x 10~6 g cm'4, pa° is the density of
native soil air taken to be 1.12 x 10~3 g cm"3, and pac is the density of contaminants in
the gas phase computed as
Pac= £ Caa (4.37)
0 = 1
where Coa is the gas phase concentration of species a and na is the number of noninert
organic species.
4.3 Model Applications
4-3.1 Example 1. 1-D Spill of Two-Component Hydrocarbon
The first example problem involves a spill of a mixture of toluene and o-xylene at the
soil surface followed by gradual leaching by rainfall. The domain of the problem is a
200 cm vertical section with a water table at a depth of 150 cm described by a single
column of forty 5 x 5 cm elements using 82 nodes. Properties of the porous medium are
given in Table 4.1. Note that since the problem is one dimensional, Kfw alone controls
the flow. However, since the solution employs two-dimensional elements Ktw is still
employed. To minimize the possibility of oscillations in 1-D vertical problems, it is
advisable to input Ktw > Klw . Here, we use Ktw 2 Ktw . The nonaqueous liquid
entering the system is assumed to be composed of 50% mass fractions of toluene and o-
xylene. The oil specific gravity, relative viscosity and scaling factors of the mixture
estimated as mass fraction-weighted averages of the individual component values are
given in Table 4.1. Component values of air, oil and water phase diffusion coefficients,
and oil-water and air-water partition coefficients are also summarized in Table 4.1.
Solid-water partition coefficients and decay coefficients were assumed to be zero for
both components. The problem is simulated in three stages involving one initial run
140
-------
and two restart runs as follows:
Stage 1: Constant head oil infiltration
Stage 2: Redistribution with transport for 25 days
Stage 3: Constant flux water infiltration for 100 days
In Stage 1, oil infiltration is considered into soil which is initially oil-free and in
equilibrium with a water table at a depth of 150 cm. Oil is added under a constant
head of h0 0 cm at the top two nodes until the cumulative oil infiltration is 21 cm3
( = 4.05 cm3cm"2). Other boundaries for oil are no flow. The bottom nodes are
maintained at a constant water head of hw = 50 cm corresponding to the initial
conditions throughout Stage 1 and other boundaries are no flow for water. Transport is
not considered during Stage 1 since the total duration of the first stage is only 0.0074 d.
At the end of infiltration, the oil front is at a depth of approximately 20 cm (Figure
4.1).
Stage 2 involves a 25 day period during which oil redistribution occurs with no flow
conditions at the upper boundary. Initial conditions for flow are read from the restart
file. Transport is considered during Stage 2 with initial conditions computed internally.
Initial oil phase concentrations of toluene and o-xylene corresponding to 50% mass
fractions are 436.7 mg cm"3. Equilibrium water phase concentrations for toluene and
xylene are 0.2583 and 0.0762 mg cm"3, respectively. Stage 2 is carried out for a period
of 25 days with no flow boundary conditions for oil and with boundary conditions for
water as in Stage 1. Type-3 boundary conditions for transport are employed with the
influent concentration, %£> equal to zero corresponding to zero component mass flux
on all boundaries. The water and total liquid saturation profiles at the end of Stage 2
are shown in Figure 4.1. Note that a residual oil saturation of about 7% occurs in the
upper unsaturated zone due to a limitation on drainage under gravity when the oil
permeability is low. Oil permeability increases near the water table as the water
saturation increases resulting in faster oil drainage near the water table.
141
-------
Table 4.1 Input parameters for Example 1.
Soil properties: Bulk fluid properties:
KIWi = 800 cm cf1 /3ao = 2.1 pow = 1.83
K.w'z = 400 cm fl r,ro = 0.695 pro = 0.873
= 0.4 Component properties:
Sm = 0.05 o-Xylene Toluene
5or = 0.2 Da°w = 0.620 0.821
o = 0.05 cm"1 Da°0 = 0.70 0.987 cm2^1
n = 2.5 D0°a = 6099. 6765. cm2^1
AL = 2.0 cm roo = 5729. 1683.
AT =0.2 cm roo = 0.22 0.28
pa. = 880. 862. mg cm'3
142
-------
JOO.O -i
Stoje 1
150.0
'
-g- :
u, :
JS 100.0 -
.Q :
«
\
v
X.
^^
»«) St
*-»-*-*-*s*
50.0 -
0.0'
0.03
'olio'
0.4O 0.60
Seturction
aoo
200.0 -,
150.0 -
£100.0
0.
50.0 -
0.0
0.00
0.20
0.4O 0.60
Soturetien
Stc«e 2
0.60
1.03
200.0 -i
150.0
5*100.0
I
50.0
0.0
Stoje 3
o.oo 0.20 0.40 0.60 blab "i!6o
Saturation
Figure 4.1 Water and total liquid saturation vs depth at end of Stages 1, 2 and 3 for Example 1.
143
-------
At the beginning of Stage 3, water infiltration is initiated at the upper surface at a
constant flux of qw = 10 cm d~l via a type-2 flow boundary condition. The water head
was maintained as in the previous stages at 50 cm at the bottom of the column. Zero
oil flux was assumed on all boundaries. Initial conditions for flow and transport were
read from a restart file from the end of Stage 2. Assuming infiltration of clean water, a
type-3 transport boundary condition with C^w' 0 was imposed at the upper surface.
At the lower surface and elsewhere a zero concentration gradient (type-2 condition) was
employed to allow convective mass loss in the water phase from the bottom of the
column.
Concentrations of xylene and toluene at the lower boundary are shown as a function of
time in Figure 4.2. Due to its higher solubility, toluene is lost from the system more
rapidly than xylene leading to an increase in the mass fraction of xylene in the oil
phase, and hence an increase in aqueous phase concentration at short times, followed by
a reduction as extraction in the water phase continues.
0.3
Q 25
E 0.2-
: 0.15-
a as-
Toluene
Xylene
ZQOO 40.00 eaao earn taaao 120.00
Time
Figure 4.2 Aqueous xylene and toluene concentrations at outflow boundary for Example 1.
144
-------
4.3.2 Example 2. 2-D Planar Spill with Sloping Water Table
This problem involves a light hydrocarbon spill in a two-dimensional Cartesian section
through saturated and unsaturated zones with a water table gradient. The problem is
simulated in two stages involving one initial run and one restart run as follows:
Stage 1: Constant head oil infiltration on strip
Stage 2: Redistribution with transport for 25 days
The domain is 8 m in the vertical and 11 m in the horizontal with a water table initially
located 4 m from the bottom on the left boundary and 3.5 m on the right. The lower
surface is an aquitard and the top is the soil surface. An oil leak occurs on a strip 2 m
in width centered 5 m from the left boundary along the top. Spatial discretization is
achieved with 88 elements of 1 x 1 m spacing for a total of 108 nodes. Initially the
water table is assumed to vary linearly from the left to right boundaries and water
pressures are hydrostatic with no oil in the system. Boundary conditions for the water
phase in Stage 1 are no flow on the top and bottom boundaries and prescribed head on
the side corresponding to the initial conditions allowing water flow in the aquifer from
the left to the right. Boundary conditions for oil are no flow except on the strip source.
(Note: this presumes the physical domain is large enough that free oil never reaches the
boundaries.) At the beginning of Stage 1, oil infiltration commences on the strip source
under a constant head of h0 = -0.1 m and is continued until the total accumulation of
oil is 1.0 m3 (per m in the third dimension). The organic liquid is assumed to be 89.5%
inert oil and 10.5% benzene. Soil properties and bulk fluid properties are given in Table
4.2. Transport is not simulated during Stage 1 which lasts for only 4.17 d.
During Stage 2, oil redistribution and benzene transport is simulated subject to the
same boundary conditions for the water phase as in Stage 1 and with no flow for oil on
all boundaries. A zero concentration gradient was imposed on the downstream
boundary to allow mass loss with the water phase by convection. The zero gradient
condition was employed on the upper and lower surface which have zero fluid flux
resulting in the equivalent of a zero mass flux condition for transport. On the upstream
boundary, a type-3 condition with zero influent concentration was imposed. Stage 2
was continued for a period of 25 days during which time the solution indicated a mass
loss of benzene due to dissolution and transport of about 2%. Contours of oil saturation
and water and gas concentrations at the end of Stage 2 are shown in Figures 4.3 - 4.5.
145
-------
Table 4.2 Input parameters for Example 2.
Soil properties:
Ktw = 10.0 m d'1
Ktw = 5.0 m fl
$ 0.35
5m = 0.05
Sor = 0.2
o = 5.0 m"1
n = 2.8
AL = 0.2 m
AT = 0.04 m
Bulk fluid properties:
Pao = 2.69 /3OW
Iro = 2-0 Pro
Component properties
D0°w = 0.94x10-" m'
D°0 = 1.13xlO'4 m'
^^ O f\ 7CO MVk* J"l
^/ =^ U. /OO T7Z O
roo = 493.
Toa = 0.24
Pa = 0.87 xlO6 5
= 1.59
= 0.832
(benzene):
!(Ta
lfl
m-3
1.0 2.0 3.0 4.0
5.0
6.0
7.0
B.O
9.0
10.0
11.0
o
o
11.£P
1.0
2.0
3.0
4.0
S.O 6.0
X(m)
7.0
6.0
9.0 10.0
Figure 4.3 Oil saturation contours at the end of Stage 2 for Example 2.
146
-------
1.0 20 3.0 4.0 5.0 6.0 7.0 6.0 9.0 10.0 11.0
2.0 3.0 4.0 5.0 6.0 7.0 6.0 9.0 10.0 11.0
0.0
0.0
Figure 4.4 Water phase concentration contours (g m~3) at the end of Stage 2 for Example 2.
00 1.0 2.0 3.0 4.0 S.O 6.0 7.0 6.0 9.0 10.0 11.0
3.0 4.0 S.O 6.0 7.0 6.0 9.0 10.0 11.0
0.0
Figure 4.5 Gas phase concentration contours (g m"3) at the end of Stage 2 for Example 2.
147
-------
4-8.8 Example S. 2-D Radial Dense Solvent Spill with Vacuum Extraction
This problem involves a spill of tetrachloroethylene (PER) in a radial domain and
remediation using vacuum extraction. The simulation is performed in three stages with
two restarts as follows:
Stage 1: Oil infiltration event
Stage 2: Redistribution and transport under natural gradients
Stage 3: Remediation using gas vacuum extraction
The first stage of the problem involves infiltration of PER on a circular area with a
radius of 2.0 m at the soil surface. A total of 3 m'3 of NAPL is assumed to infiltrate
under a water-equivalent oil head of -0.25 m. A water table occurs at a depth of 3.2 m
and an impermeable layer occurs 5 m below the surface. The soil is uniform and has
properties given in Table 4.3. Fluid properties for PER are also given in the table. The
problem is analyzed as a 2-D radial section with an inner radius of 0.1 m (in order to
facilitate subsequent analysis of pumping from a well of the same radius) and an outer
radius of 8 m. A mesh with 12 nodes in the vertical direction and 16 nodes in the
horizontal direction is employed.
Initial conditions for the water phase are assumed to correspond to equilibrium with the
water table with no oil. Boundary conditions for the water phase involve a type-1
condition on the right side below the water table corresponding to vertical hydrostatic
conditions (i.e., same as the initial conditions). All other boundaries are no flow for the
water phase. Boundary conditions for the oil phase are a type-1 condition with h0 = -
0.25 TO on the 5 nodes on the upper surface between r = 0 and r = 2.0 m. All other
nodes are no flow for oil. Gas flow was not considered during infiltration. Termination
of Stage 1 occurred at t = 6.48 d after infiltration of 3.128 m3 of oil.
Stage 2 of the problem involves a continuation of the Stage 1 using the restart option.
Redistribution of the PER is permitted for a period of 25 d under no flow boundary
conditions for the oil phase on all boundaries. Boundary conditions for the water phase
are maintained as in the infiltration stage. Gas phase flow is again disregarded.
Transport is considered with the initial aqueous concentration of PER at nodes with oil
phase set to the solubility of 150 g m'3. Boundary conditions for transport are type-2
(zero dispersive flux) on all boundaries which is equivalent to zero flux since there is
148
-------
essentially no flow on the boundaries during Stage 2. By the end of the redistribution
period, the NAPL plume has reached the lower aquifer boundary. However, a
substantial volume of the spill is still retained in the zone above the water table as
"residual saturation" retained by capillary forces. The oil saturation distribution at the
end of Stage 2 is shown in Figure 4.6.
In Stage 3, remediation is simulated with vacuum pumping in the unsaturated zone. A
vacuum well with a 0.1 m radius is assumed to be placed at the left boundary screened
over a 2.5 m interval from a depth of 1.25 m to 3.25 m and regulated at a pressure head
of ha = -1.5 m. The gas boundary at the well is treated as a type-1 boundary condition.
Gas inflow is permitted along the upper surface from r = 3.5 m to the outer perimeter
under a type-1 condition for gas with a constant pressure of ha = 0. The inner 3.5 m of
the surface, which is assumed covered, and other boundaries are treated as no flow for
the gas phase. Water head is prescribed on the right boundary below the water table as
in the earlier simulations and all other boundaries are treated as no flow for water. No
flow conditions are imposed for the oil phase on all boundaries. Boundary conditions for
transport are type-3 with zero influent concentration on the top and right side
boundaries and zero dispersive flux elsewhere. Gas pumping results in water table
upwelling as shown in Figure 4.7. The gas flow rate stabilizes at 319 m3 d~l after 2 days
and the PER recovery rate reaches a corresponding steady rate of 16.77 kg d'1. PER
mass in the soil vs venting time (Figure 4.8) shows a small increase ( ~ 1.5%) during the
first day due to numerical error before the solution stabilizes and mass removal rate
becomes nearly constant.
Table 4.3 Input parameters for Example 3.
Soil properties:
K»w^
K.w*
* *
sm
SOT
a
n
AL
A f
=
=
=
=
=
=
=
zz
:=
5.0 m a
2.0 m
-------
£00
.10 100 190 260 370
403 H
450
550 640
750
10 100 190 280 370 460 550 640
R(cm)
503
400
300
200
100
713
Figure 4.6 Oil saturation contours at the end of Stage 2 for Example 3.
10.0 100.0 190.0 260.0 370.0 460.0 550.0 640.0 7300
500.0
403.0
300.0
N
E
N 200.0
100.0
1C
II 1 1 1 1 1 1
°"s
" \>^r °35 -ass
^^^~ 0.55 0.55
0^ ^-0.95
^ , ^~^"-93
.0 100.0 190.0 260.0 370.0 460.0 SSO.O 640.0 730.0
5OO.O
H 400.0
-\ 300.0
^200.0
4 100.0
0.0
R(cm)
Figure 4.7 Water saturation contours at the end of Stage 3 for Example 3.
150
-------
4.75E+C06 =
4.70E4006 !
4.65E+006 \
4.60E-f006
0.00
* I I I [ I tI I I I I i
5.00 10.00 15.00
Time(days)
Figure 4.8 PER mass in system vs time during Stage 3 of Example 3.
151
-------
5. TRANSPORT WITH NONEQUILIBRIUM INTERPHASE MASS TRANSFER
5.1 Model for nonequilibrium intetphase mass transport
5.1.1 Governing phase transport equations
The transport equations for a component developed in Chapter 4 can be expressed for
the water phase (w), the organic liquid phase (o) and the gas phase (a) as
To accommodate adsorption of Q by the solid phase, an additional continuity equation
is
= R0, - vas Cat (5.2)
where Cat is the solid phase concentration expressed as mass of adsorbed component Q
per porous medium volume [M L"3], »ot is the first-order decay coefficient [T"1] in the
solid phase and Ras is the mass transfer rate per porous medium volume to (+) or from
(-) the solid phase [M I/3 I"1].
5.1.2 Consideration of nonequilibrium interphase mass transfer
The phase-summed formulation of the transport equations may be generalized to
account for nonequilibrium phase partitioning by introducing the concept of apparent
partition coefficients. Under transient field conditions, at a given location in time and
152
-------
space, actual phase concentration ratios may differ from the true equilibrium ratios
defined by [4.17]. With this in mind, we define apparent partition coefficients - which
will vary in time and space by analogs of [4.17] as
Cao = T0'0 Cow (5.3a)
Caa = Ta'0 Cow (5.3b)
C0. = C CQW (5.3c)
For any two phases which are in physical contact, the rate of mass transfer will be
described by first order mass transfer functions of the form
*al2 = *.12 (Cal ~ Co\) (5.4a)
= -kol,(Ca2 - C0\) (5.4b)
where Ran is the rate of mass transfer of a per porous medium volume from phase 1 to
phase 2, (7ol is the actual concentration of o in phase 1, C^i is the concentration of a
which would occur in phase 1 if it were in equilibrium with phase 2, Co2 is the actual
concentration of a in phase 2, C^2 is the concentration of a that would occur in phase 2
if it were in equilibrium with phase 1, and A;ol2 is a mass transfer rate coefficient [T"1].
Consider first, the case of mass transfer between oil and water phases when both phases
exist at a point in time and space. Employing [5.4a] along with [5.3a] for the actual oil
phase concentration and [4.17a] for the equilibrium concentration gives
= *
which may be solved in terms of ra'0 as
D
Ft r _L "pom (f £\
ao - Fao T L /nr (5.6)
which indicates that the apparent partition coefficient may be expressed in terms of the
actual partition coefficient plus a correction term, which depends on the sign and
magnitude of the actual mass transfer rate and on the concentration in the water phase.
153
-------
If free oil is present in the system at a given point in time and space, so that air and oil
are physically in contact, mass transfer between oil and gas phases may be described in
a similar fashion. Employing in this case [5.4b] we obtain
oa
(ro
which yields
"000
-
In the absence of free oil in the porous medium at a given location, mass transfer will
occur between water and gas phases rather than between oil and gas phases. Proceeding
in the same manner as for oil-gas mass transfer yields
n
r ' r aau>a
1 *
era oa L
"'awa
Note that either [5.8] or [5.9] will be relevant at a given location and time, depending on
whether free oil occurs. Finally, mass transfer between water and solid phase may be
considered by employing [5.4b] to obtain
r.f. = r0, - ."!, (5.10)
.
atvf
In order to describe transport with nonequilibrium phase partitioning, the individual
phase transport equations, [5.1a-c] and [5.2], may be summed as described in section
3.1.2 but using apparent rather than equilibrium partition coefficients. Thus [4.19] will
hold for the nonequilibrium case with apparent partition coefficients used in lieu of the
equilibrium coefficients. Nonlinearity is introduced into the resulting phase-summed
transport equation due to the dependence of apparent partition coefficients on
concentration as well as mass transfer rates. Thus, an iterative solution procedure will
be required for the nonequilibrium problem as will be discussed later.
If the transport equations are written in terms of oil phase concentrations, similar
expressions to those given above may be derived to define apparent partition coefficients
as functions of mass transfer rates and oil phase concentrations.
154
-------
Using the nonequilibrium relations, we may rewrite the phase transport equations in
terms of a oil phase concentration. Using [5.3] to eliminate oil, gas and solid phase
concentrations from [5.1] and summing the equations noting that
leads to the a-component phase-summed transport equation in terms of water phase
concentrations
* dCaw _
where
*Z = 4SW + * S0T'ao + 4>Sa roo + T'a, (5.12b)
£Q* =Sw Dawij + S0 DQoij T'ao + 4Sa Doaij T'aa (5.12c)
««,* = Qwi + 9oi r^0 + 9ai rL« (5.12d)
p P r1 7? r'
* i rii i r' _L r' _L w _L o oo _L a QQ /c 10Q^
^o - "aw ^ ''ao rao + "aa ^0 ^ "o« T<» + ~P^ + p^~ + pa (5.12eJ
Note that [5.12] has the same form as the simple single phase transport equation.
However, the coefficients represent pooled effects of transport in all phases. Note also
that interphase mass transfer terms occur in the phase summed equation only as a sum
over all components.
An alternative form of the phase-summed transport equation may be written in terms
of the oil phase species concentrations. The a-component phase-summed transport
equation in terms of oil phase concentrations has the form
** 00 d f r) ** >o 1 ., ** oo **
ri /K 10 \
°° ( ^
where
, Q, /K1^k^
+ p- (5.13b)
1 00 * 00 * 00
155
-------
**
W W
1 oo * ao
/T
** _ i "mi , "oi
i - 9o. + p- + p
1 oo *
(5.13e)
'a1 oo
The form of the oil-based equation is identical to that of the water-based phase-summed
equation. The water-based equation has the conceptual advantage that Caw is physically
meaningful at all locations in the porous medium, whereas Cao only has physical
meaning at locations where S0>0. The latter constraint, however, imposes no
mathematical difficulty since [5.3] can in any event be employed as a definition of CQO
as well as Coa whether or not the phase exists. Since most of the component mass
commonly resides in the oil phase, the oil-based equation has the advantage that mass
balance accuracy of the solution may be more readily controlled.
When simulating two-phase liquid flow consisting of water and oil with constant gas
pressure, only diffusive transport is considered in the gas phase i.e., qa = 0 is
assumed for this case. It should be noted that even if gas pressure gradients have a
negligible effect on liquid phase flow, this does not strictly imply zero gas flow. If fluid
saturations change in time, gas flow will occur. Furthermore, phase partitioning may
result in density gradients in the gas phase which can induce significant flow.
5.1.S Solution approach
The solution approach for solving the coupled multiphase flow and multicomponent
transport problem with nonequilibrium mass transfer is as follows
1. Solve the fluid flow equations simultaneously for the current time-step using time-
lagged phase densities and interphase mass transfer rates.
2. Solve the phase-summed transport equation using current values of apparent
partition coefficients, interphase mass transfer rates and phase densities.
156
-------
3. Back-calculate interphase mass transfer rates, update phase densities and
apparent partition coefficients and repeat (2) until transport solution converges.
4. Proceed to the next time-step.
5.1.4 Numerical verification
Example 1^. 1-D Spill of Two-Component Hydrocarbon. Details of the input data and
the results for this example using equilibrium partition coefficients were presented in
Chapter 4 (Section 4.3.1). In this section, stage 3 of example 1 has been repeated using
nonequilibrium interphase partition coefficients. Several simulations of stage 3 were
performed using nonequilibrium interphase mass transfer analyses and results were
compared with those for the equilibrium model described in Chapter 4. To provide
greater mass balance stability for the nonequilibrium analyses, the oil-based phase-
summed transport formulation was employed. In addition to an equilibrium simulation,
analyses were performed using rate constants for all phase pairs of 1.0, 0.01 and 0.005
d "1. Concentrations of xylene and toluene at the lower boundary are shown as a
function of time for the various runs in Figure 5.1 and 5.2, respectively. Due to its
higher solubility, toluene is lost from the system more rapidly than xylene leading to an
increase in the mass fraction of xylene in the oil phase, and hence an increase in
aqueous phase concentration at short times, followed by a reduction as extraction in the
water phase continues. As the phase transfer rate constants increase to 1.0 d"1, the
nonequilibrium model is seen to give results identical to the equilibrium model. This
degeneration to the equilibrium case provides verification of the nonequilibrium model
formulation. As the rate coefficients decrease, the exit concentrations, and rate of
leaching of contaminants from the oil phase, both diminish as expected.
157
-------
0.08
0.00 20.00 40.00 60.00 80.00 100.00 120.00
Time (Days)
Fig. 5.1 Xylene effluent concentration in water vs time for Stage 3 of Example 1.
0.3-r
E
to
e
^"
c
o
0.2H
a
i 0.1 5-\
C
u
u
C
U 0.1-
V
C
jp o.osH
K=0.005
4-
K=o.or
»K
K=T.O
EQUILIBRIUM
0.00 20.00 40.00 60.00 80.00 100.00 120.CO
Time (Days)
Fig. 5.2 Toluene effluent concentration in water vs time for Stage 3 of Example 1.
158
-------
5.2 Laboratory Investigations
A series of laboratory column experiments were conducted to investigate the
nonequilibrium interface mass transfer of soluble constituents from a multicomponent
nonaqueous phase liquid (NAPL) at residual saturations during steady-state water flow.
Breakthrough curves obtained from these experiments were inverted to determine the
mass transfer rate coefficients. The objectives of the laboratory investigations were: (i)
to determine the effects of aqueous phase pore velocity, NAPL saturation, and NAPL
composition on oil-water mass transfer rate coefficients, (ii) to test the experimentally
derived rate coefficients against that generated by existing models for predicting mass
transfer rate coefficients, and (iii) to develop a new empirical model for nonequilibrium
mass transfer if results from (ii) are not satisfactory.
5.2.1 Experimental methods
NAPL and porous medium properties. The NAPLs used in experiments were composed
of soluble and insoluble (inert) components. The inert component of NAPL was allowed
to comprise more than 50% of the total NAPL volume in order to minimize volume
changes during the course of experiments, and to keep the NAPL-aqueous phase
interfacial area nearly constant. Two sets of experiments were conducted using two
different inert components, namely soltrol and hexadecane. For purposes of distinction,
the NAPLs with soltrol and hexadecane as their inert components will be referred to as
NAPL-1 and NAPL-2, respectively. The soluble constituents of NAPL were benzene,
toluene, ethylbenzene and xylene in varying amounts. All constituents used with the
exception of soltrol are well defined compounds. Soltrol is a light hydrocarbon composed
of straight and branched chained alkanes with carbon numbers ranging from 10 to 16,
hence it has a low aqueous phase solubility. Table 5.1 lists the relevant properties of the
NAPL constituents.
The porous medium consisted of a well graded sand with fractions of 3.4% very fine
sand, 13.2% fine sand, 38.7% medium sand, 29.4% coarse sand and 14.6% very coarse
sand.
159
-------
Table 5.1. Properties of NAPL constituents used in the experiments.
Constituents Mol. Weight (g) Solubility (mg/1) Spec, gravity Diff. Coef (m2/d)
Benzene
Toluene
Ethylbenzene
Xylene
Hexadecane
Soltrol
78.11
92.10
106.14
106.17
226.44
160.00
1780
515
152
192
0.0009
0
0.878
0.867
0.867
0.880
0.775
0.785
9.42 x 10'5
8.21 x 10'5
6.21 x lO'5
6.21 x 10'5
-
-
Experimental Setup 1-D Column Studies
Soil/Oil Mixture
Constant head
Sampling Port
Sand gravel mixture
Fig. 5.3 A schematic of the soil column used for nonequlibrium mass transfer studies.
160
-------
Experimental procedure. Experiments were conducted in a 45 cm long glass column
with a cross section of 7.7 by 7.7 cm (Figure 5.3). The column was furnished with a
Marriotte bottle attachment to achieve steady state flow conditions. A 5 cm thick layer
of pea sized gravels and a wire screen were placed at the bottom of a dry, clean column.
Clean sand was added to the gravel layer leaving the sampling port just above this sand
mixed gravel layer (Figure 5.3). In a separate container, 1500 g of well graded sand was
mixed uniformly with a known volume of NAPL to produce a desired NAPL saturation.
The column was packed uniformly using NAPL-sand mixture as quickly as possible.
During the mixing and packing processes, care was undertaken to minimize
volatilization losses. The top of the column was sealed with a plastic wrap to reduce
further volatilization. The column was slowly saturated from the bottom, to prevent air
entrapment. The duration of the saturation period ranged from 12 to 15 minutes. Upon
completion of saturation, the water source was moved to the top of the column to
maintain a constant head at the top of the column with the Marriotte bottle
attachment. Flow was initiated through the column adjusting the elevation of the flow
outlet until the desired flow rate was achieved. Outflow samples were collected for a
duration of 20 s periodically at 3 to 10 minute intervals during the course of the
experiment to analyze for aqueous phase concentrations. Flow rates were constantly
monitored and adjusted as required. Experiments were conducted at discharge rates of
50 and 100 ml/minute for a duration of 3 and 4 hours, respectively.
The experimental variables were NAPL saturation of the porous medium, NAPL
composition and aqueous flux. Table 5.2 summarizes the experiments performed for
NAPL- 1 and NAPL-2, respectively. Each experiment was replicated at least twice to
ascertain the reproducibility of observed results.
Sample collection and handling. At each sampling interval, a syringe was used to
extract approximately 5 cm3 of solution after the dead volume of the sampling tube had
been withdrawn and discarded. The sample was directly filtered into a 3.7 cm3 vial and
tightly capped, ensuring no head space was created in the process. Samples were stored
at 8 "Cand were analyzed within 24 hours.
Samples from NAPL-1 were analyzed using high pressure liquid chromatography
(HPLC), while those from NAPL-2 were analyzed using a packed-column type gas
chromatography (GC), utilizing standard solutions prepared in methylene chloride.
161
-------
Table 5.2 Specifics of the nonequilibrium mass transfer column experiments (NAPL composition, pore
water velocity, NAPL content), estimated mass transfer rate coefficients, and W values for
estimations.
Exp.
*
1
2
3
4
5
6
7
6
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
Chem.
Comp.
Index
1
1
1
1
1
1
1
1
1
1
1
2
2
2
2
2
2
2
2
2
2
2
2
2
3
3
3
4
4
5
5
5
V (cm/s)
0.0303
0.0321
0.0308
0.0324
0.0336
0.0665
0.0653
0.0659
0.0642
0.0631
0.0642
0.030
0.033
0.035
0.033
0.034
0.036
0.033
0.033
0.032
0.069
0.065
0.067
0.063
0.034
0.034
0.033
0.061
0.061
0.061
0.062
0.060
©oil
0.0193
0.0395
0.0397
0.0287
0.0283
0.0388
0.0373
0.0194
0.0284
0.0276
0.0284
0.029
0.030
0.029
0.030
0.029
0.029
0.029
0.020
0.019
0.029
0.030
0.0198
0.0197
0.029
0.0297
0.019
0.0294
0.0294
0.0291
0.0290
0.0296
K (1/sec)
0.3005
0.3120
0.2402
0.4442
0.4533
0.4305
0.5872
0.4233
0.5569
0.3672
0.3791
0.0067
0.2662
0.2475
0.1537
0.3553
0.2487
0.2587
0.2943
0.2735
0.3701
0.3143
0.4234
0.3713
0.2407
0.2373
0.1594
0.3394
0.3618
0.2580
0.2590
0.2575
R'
0.945
0.957
0.941
0.965
0.880
0.985
0.983
0.970
0.902
0.955
0.959
0.830
0.781
0.732
0.742
0.676
0.970
0.964
0.874
0.860
0.956
0.980
0.972
0.984
0.984
0.978
0.902
0.107
0.975
0.990
0.993
0.991
1: 10% Benzene + 10% Toluene + 10% Ethylbenzene + 10% Xylene + 60%
Hexadecane (NAPL-1)
2: 10% Benzene + 10% Toluene + 80% Soltrot (NAPL-2)
3: 10% Benzene + 3% Toluene + 87% Soltrol (NAPL-2)
4: 10% Toluene + 90% Soltrol (NAPL-2)
5: 10% Benzene + 90% Soltrol (NAPL-2)
162
-------
5.2.2 Data analysis methods
Aqueous transport of constituents through the experimental soil column, considering
nonequilibrium mass transfer from a multi- component NAPL at a residual saturation,
can be described by one- dimensional ad vective- dispersive- reactive equation
(5.14)
... n
subject to the following initial and boundary conditions
C(x,0) = C0' (5.15a)
C(0,t) =0 (5.15b)
dC(L,t) n
dx
(5.15c)
where Rld is the retardation factor for species t; D is the dispersion coefficient [cm2 s"1];
Vis the aqueous phase pore velocity [cm s"1]; C" is the concentration of species i in
the aqueous phase [mg L"1]; K is the mass transfer rate coefficient [t~l\; Ce'q is the
equilibrium solubility of component i [mg L~1]; C0' is the initial concentration of species
in the aqueous phase [mg L~1}; x is the spatial coordinate [cm]; t is time [s]; and n is the
number of soluble components in NAPL. Rj were estimated for the porous media from
batch experiments outlined by Abdul and Gibson (1986). C0' was taken as the
concentration of species determined just before flow was initiated. Assumptions
associated with [5.14] and [5.15] are (i) first order nonequilibrium partitioning between
NAPL and aqueous phase, (ii) linear equilibrium partitioning between aqueous and solid
phases, (iii) no partitioning between NAPL and solid phase, and (iv) single rate
coefficient describing nonequilibrium mass transport for all species.
The mass transfer rate coefficient (K) is defined as
K = k a (5.16)
where k is the mass transfer coefficient [cm s'1] and a is the aqueous phase-NAPL
163
-------
interfacial area per unit porous media volume [cm"1]. It is more convenient to determine
K rather than k since there is no rational and independent method for estimating a.
Since Ce'g is dependent on the concentration of all other species in the mixture, [5.14]
was solved numerically using a Crank-Nicholson finite difference scheme to simulate the
experimental breakthrough curves. The dependence of Ceg on the concentration of other
species at the current time step was eliminated by using C e'q calculated from the
previous time step. Ce'g at a given node was calculated as a linear function of the
component mole fraction at the same node and its aqueous phase solubility limit from
the previous time step as
t) = C^x, t-st) = jrir- C* (5.17)
where Ceq(x,t) is the equilibrium solubility at the current time step; Ceq(x, t 6t) is the
equilibrium solubility at the previous time step; Cel is the aqueous phase solubility of
component i; Af is the number of moles of component «; n+1 is the number of soluble
plus inert components. The denominator is the total number of moles present in the
NAPL.
Soltrol used as the inert component of the NAPL is not a pure compound but a mixture
of several compounds. Hence, the determination of an effective molecular weight for
soltrol, which is needed to determine the number of moles of soltrol in NAPL-1,
required special attention. A batch experiment was conducted, in which excess of a
mixture of 80% soltrol, 10% benzene and 10% toluene was allowed to equilibrate with
water. The equilibrium solubilities of benzene and toluene in the mixture determined by
combining the excess of the mixture with water in a volumetric flask and gently stirring
over night with a magnetic stirrer, maintaining the interfacial contact. Aqueous samples
were analyzed by the GC for benzene and toluene equilibrium concentrations. The
number of moles of soltrol in the mixture was calculated from [5.17] using the measured
values of equilibrium solubilities for benzene and toluene. Subsequently, the effective
molecular weight of soltrol was estimated from the relationship
" m' (5-18)
164
-------
where M, is the number of moles of soltrol ; mt is the mass of soltrol [g] used in the
experiment; and Wt is the effective molecular weight of soltrol [g]. The estimated values
of Wt corresponding to experimentally determined equilibrium solubilities of benzene
and toluene was 156.0 and 163.8 g, respectively. An average molecular weight of 160 g
was adopted for the analysis, which corresponds to that of hexane. Since soltrol is a
mixture of c!0-c!6 alkane, the estimated molecular weight for soltrol appears to be
reasonable.
The finite difference solution of [5.14] was incorporated in a nonlinear least squares
routine to estimate the parameters Rld, D and K through a numerical inversion
procedure. The nonlinear least squares routine was implemented using the Levenberg-
Marquardt algorithm to minimize the objective function <£
/ - cy (5.19)
where m is the number of measured data points, each corresponding to a specified time;
n is the number of soluble constituents; Cj' is the predicted concentration of component
A .
i for observation number ;; and Cj is the measured concentration of component i for
observation number j. Hence, the estimated values of parameters Rj, D and K
correspond to the best-fit between the predicted and observed breakthrough curves.
Figures 5.4 through 5.7 show some of these fits between predicted and observed
breakthrough curves. The estimation of parameter Rld is included in the least square
routine in order to account for the adsorption of the dissolved NAPL components by a
piece of cotton at the sampling port, which is mistakenly used during the column
experiments to prevent the sampling port from clogging. The values of Rld estimated
through optimization serve for the purpose of explaining the apparent retardation of
observed peak concentrations, in turn, improving the fits between the observed and
predicted breakthrough curves. These values of R'd are different from those used for the
soil in that they are limited to the last node (representing the sampling port) where no
mass transfer is allowed.
165
-------
400
R squore - 0.993
Simulated (Benzene)
Obnrved (Benzen*)
5000 10000
Time (seconds)
150CC
Fig. 5.4
Comparison of observed and predicted breakthrough curves for experiment 31
(10% benzene with pore water velocity of 0.062 cm/sec and NAPL content of
0.029)
300-q
R square 0.984
Simulated (Benzene)
Ob9«rv«d (Benzene)
Simulated (Toluene)
Observed (Tolutne)
5000 10000
Time (seconds)
15000
Fig. 5.5 Comparison of observed and predicted breakthrough curves for experiment 25
(10% benzene, 3% toluene with pore water velocity of 0.034 cm/sec and
NAPL content of 0.029)
166
-------
400'
300-
I 200-
c
V
u
o
O
100-
R «quore -' 0.980
Simuloted (Benjene)
Obterved (Beniene)
Simuloled (Toluene)
Observed (Toluene)
5000 10000
Time (seconds)
15COO
Fig. 5.6 Comparison of observed and predicted breakthrough curves for experiment 22
(10% benzene, 10% toluene with pore water velocity of 0.065 cm/sec and
NAPL content of 0.030)
300
R squore 0.674
Simulated (Benzene)
Objerved (Benzene)
Simulated (Toluene)
Observed (Toluene)
,1 I I
5000 10000
Time (seconds)
15000
Fig. 5.7 Comparison of observed and predicted breakthrough curves for experiment 19
(10% benzene, 10% toluene with pore water velocity of 0.033 cm/sec and
NAPL content of 0.020)
167
-------
5.2.3 Experimental Results
Mass transfer rate coefficients K estimated using the numerical inversion model are
given in Table 5.2 for all experiments consisting of combinations of different NAPL
composition, NAPL saturation, and aqueous phase pore velocity. Increasing the aqueous
pore velocity generally increased the mass transfer rate coefficient. However, no
relationship could be discerned between the rate coefficient and NAPL saturation. The
species concentration in effluent increased with aqueous pore velocity (Figures 5.4 and
5.7), indicating that increases in mass transfer rate coefficient contributed more than
compensating for the effects of reduction in detention time. With only one exception, all
the estimated K values are within the same order of magnitude. The majority of
experiments produced high values for the coefficient of determination R2 indicating a
satisfactory model performance. Coefficients of variation CV of estimated K range from
1.5% to 92% with a mean CV of 18% and a standard deviation of 21%. These CV values
correspond to a reasonably narrow confidence interval for the vast majority of estimated
K values emphasizing the reliability of estimated K values. In addition, reasonably low
correlation coefficients obtained between the estimated K and D indicate that estimates
of K are insensitive to the estimated D values.
To check the consistency of results presented in Table 5.2 with those reported in the
literature, a comparison was performed between experimental K values and K values
obtained from an empirical relationship reported by Miller et al. (1990). This empirical
relationship is of the form
Sh = /30Re ^0/zSc1/2 (5.20)
where /30 = 12, /^ = 0.75, /32 = 0.6, Sh is a Sherwood number, Re is a Reynolds
number, 5c is a Schmidt number and 0n is the NAPL volume fraction. The
dimensionless Sh, Re and Sc numbers are defined as
.. _ j
(5.21a)
(5.21b)
(5.21c)
168
-------
where v is the mean pore water velocity, p is density of solute, d is the mean diameter
of the porous media grains, /z is dynamic viscosity of solute, K is aqueous-NAPL mass
transfer rate coefficient and D0 is the molecular diffusion of the component in the
aqueous phase.
Based on grain size distribution data for the sand material used as the porous media in
column experiments, d10 was determined to be 0.012 cm and dsp was determined to be
0.029 cm. Using p = 878 mg/cm3, p = 6.028 mg/cm.stc, D0 = 1.0903xlO"5 cm2/sec for
benzene, and d10 = 0.012 cm with the experimental v, and 0n values in the above
empirical relationship, corresponding K values were calculated. Results are presented in
Table 5.3 along with the experimental K values. Empirical K values are very close to
experimental K values. The same calculations were repeated using dso = 0.029 cm. The
results, which are also presented in Table 5.3, produced K values that are consistently
lower than experimental K values. These results indicate that although empirical K
values are sensitive to particle diameter, use of d10 in [5.20] seems to be more
appropriate than d50 as the representative porous media particle diameter to estimate K
values consistent with the experimental results.
169
-------
Table 5.3 Comparison of experimental (K *) and empirical (K } mass transfer rate
coefficients for effective porous media particle diameter corresponding
to d10 and dso.
d10 = 0.012 cm.
«.
0.019
0.029
0.039
tP=0.033
K*
0.3005
0.4442
0.3120
v=0.033
K
0.2430
0.3131
0.3740
v=0.066
K*
0.4233
0.3791
0.4305
v=0.066
K
0.4086
0.5266
0.6291
v=0.033
Sh
3.3278
4.2888
5.1231
v=0.066
Sh
5.5967
7.2128
8.6159
= 0.029 cm.
0.019
0.029
0.039
0.3005
0.4442
0.3120
0.0826
0.1064
0.1272
0.4233
0.3791
0.4305
0.1389
0.1790
0.2138
6.3707
8.2104
9.8076
10.7142
13.8081
16.4943
170
-------
6. TWO-DIMENSIONAL LABORATORY STUDIES OF
MULTIPHASE FLOW AND TRANSPORT
Two experiments were conducted to simulate flow and transport of a light nonaqueous
phase liquid (LNAPL) and a dense nonaqueous phase liquid (DNAPL) in an unconfmed
aquifer with two-dimensional planar symmetry. Liquid pressures were monitored using
pressure transducers and a dual gamma attenuation system was employed to measure
fluid saturations. Aqueous phase samples were collected at specific locations over time
to monitor aqueous phase mass transport of soluble components. The experimental
apparatus and procedures common to both experiments are described below first,
followed by a separate discussion of the LNAPL and DNAPL experiments.
6.1 Experimental Setup
6.1.1 Cell design
A 2-D steel reinforced tank (Figure 6.1) with amber transparent plastic sides (Ultem)
resistant to the organic chemicals used in this study was used to perform these
experiments. The flume was 1.01 m tall x 1.5 m long x 0.085 m thick. Gamma readings
were taken at 66 locations in an incomplete rectangular matrix. At 16 of these locations,
water and oil pressures were measured using hydrophilic and hydrophobic ceramic cups,
respectively, instrumented with pressure transducers tied to a data acquisition system.
Twenty-four Teflon tube inserts served as aqueous phase sampling ports in the lower
sections of the cell. Aqueous phase samples were analyzed on a gas chromatograph using
a packed column.
6.1.2 Gamma attenuation measurements
Beer's law can be used to describe the attenuation of a radioactive beam passing
through a porous medium with three fluid phases (water, NAPL and air).
/ = I0 ezp (-/v,z - nwpjwx - »npnOnx - fgPg8gx) (6.1)
171
-------
Oil Sou re*
, « - , oSurtte9
\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\
Fnttun Trtntduotr*
1.1
tn
outflow
1*0
Figure 6.1 Experimental cell layout.
Table 6.1 Properties of NAPL constituents used in experiments.
Constituent
Benzene
Toluene
Soltrol
PER
1-Iodoheptane
Specific Water Viscosity
Gravity Solubility Ratio
(ppm)
0.877
0.862
0.785
1.60
1.375
1,780
515
0
150
0
0.65
0.58
4.32
0.90
N.A
Surface
Tension
(dyne cm~1)
27.9
29.8
26.0
31.3
N.A
Interfacial
Tension with
Water (dyne cm~l)
35.0
36.1
43.3
44.4
N.A
172
-------
where / is the exciting radiation intensity in counts per second (cps), I0 is empty cell
count (counts s"1), p,- is the mass attenuation coefficient (cm2 g~l) for phase t, pi mass
density of phase t (g cm"3), and x is the path length (cm) over which attenuation takes
place. Assuming the gas attenuation to be negligible and the porous media bulk density
is invariant in time, (6.1) reduces to
/ = T0 exp (»wpjwx - »npnenx) (6.2a)
where
/. = I0 exp (-ntptx] (6.2b)
Equation (6.2) was used to determine water and NAPL saturations (Lenhard et al.
1988) from dual gamma measurements from a radiation source which utilizes low energy
Am-241 and high energy Cs-137. To accomplish this, ^,, p{ and x at every measurement
location, needs to be determined for Am-241 and Cs-137. The product of /j, and p{ !>,?,,)
for oil and water, pa for soil and i at each measurement location were determined for
Am-241 and Cs-137 according to the protocols of Lenhard et al. (1988), for each
experiment. During the calibration exercise, all counting times were 300 seconds in
duration, however, during the experiment due to the transient nature of flow processes
and the number of nodes to be scanned, the counting times were reduced to 90 seconds.
Counts were also adjusted for the high energy Cs-137 detected in the Am-241 window
according to the method of Nofziger and Swartzendruber (1974) and for the resolving
time of both sources (Fritton, 1969).
6.1.3 Pressure measurements
Hydrophobic and hydrophilic ceramic tips connected to pressure transducers were used
to measure NAPL and water pressures. Hydrophobic tips were obtained by treating
ceramic tips with chlorotrimethylsilane as outlined by Lenhard et. al. (1988).
Calibration of individual transducers involved determining voltage responses to known
but varying pressure heads from which simple linear regression coefficients were
determined. In all cases, the square of the coefficient of determination (R2) was greater
than 0.99.
173
-------
6.1.4 Sample collection and analysis
Aqueous solution samples of about 5.0 m3 were collected 2 or 3 times a day using a
syringe. Sample were directly filtered into a 3.7 cm3 vials and tightly capped, ensuring
no head space was created in the process. Samples were stored in a walk in refrigerator
at 8 °C and analyzed within 24 hours.
6.2 LNAPL Experiment
6.2.1 Experimental methods
The LNAPL used in this study consisted of a mixture of soltrol, benzene, toluene and 1-
iodoheptane in a 7:1:1:1 ratio by volume. Soltrol is a mixture of straight and branched
chained alkane with carbon numbers ranging from 10 - 16, hence, relatively insoluble in
water. Benzene and toluene are much more soluble than the other constituents and
their transport in the aqueous phase was monitored during the experiment. lodoheptane
was added to the mixture to increase the gamma attenuation of the oil phase enabling
better resolution of fluid saturations from gamma measurements. Relevant properties of
the constituents of the LNAPL are given in Table 6.1. The specific gravity of the
LNAPL mixture was 0.826. The soil used in the experiment consisted of 3% very fine
sand, 13% fine sand, 39% medium sand 29% coarse sand and 16% very coarse sand.
After filling the experimental cell with soil, bulk density measurements were performed
with the gamma system at each observation location for the air dry soil and the water
saturated soil. At the end of the experiment, bulk density measurements were
determined again via direct measurements of core samples taken from the disassembled
cell. The bulk density measurements determined from the saturated cell were used in
calculations as these yielded the lowest overall variance and were deemed most
accurate. The average bulk density of the packed cell was 1.50 g cm"3, yielding a
porosity of 0.43 assuming a particle density of 1.50 g cm"3.
Table 6.2 lists the values of /i,p, and pt determined for both energy sources. Variation
in x values between energy sources were less than 1.0% for both sources at all
measurement locations, which will yield a precision in computed NAPL and water
saturations of <1%.
174
-------
Table 6.2 Gamma attenuation parameters of porous medium and liquid for the
LNAPL Experiment
Coefficient
^ (cm'1)
*Vn (cm~1)
nt (gm.cm )
Am-241
0.196189
0.755835
0.253297
Cs-137
0.083462
0.071227
0.07534
The LNAPL experiment was conducted over a period of 3 weeks following calibration of
the gamma system. The experiment was performed in "stages" with time varying
boundary conditions as described below.
Stage 1 ^ Water Drainage. A two-phase air-water experiment was performed to
calibrate the k-S-P (permeability-saturation-pressure) model. The soil in the cell was
initially water saturated and was allowed to drain for 9 hours after the water table was
lowered to 36.5 cm and maintained at this level. In this experiment, the origin for the z
(horizontal) and z (vertical) coordinates was the lower left corner. The boundary
conditions employed during this stage were as follows.
Top and bottom boundaries:
water and oil flux = 0
Left Vertical Boundary:
water and oil flux = 0
Right Vertical Boundary:
hw -r z= 36.5 cm
water flux = 0
oil flux = 0
everywhere
everywhere
(0 < z < 36.5 cm)
(z > 36.5 cm)
everywhere
175
-------
During this stage of the experiment, water saturations and pressures as well as water
outflow were monitored. The total outflow for the period was 15,000 cm3. After nine
hours, a horizontal gradient of 1.35% was established by imposing a constant head of
36.5 cm on the left-hand side of the cell (Figure 6.1) while the water table was dropped
to 34.5 cm on the right boundary. At steady state, the aqueous flux density was
measured and was used to compute the saturated hydraulic conductivity of the soil.
Stage 2 i LNAPL Infiltration. Over a period of 11.0 hours 1,500 cm3 of NAPL was
added as strip source 50 cm from the top left-hand corner (average rate = 136 cm3 h~l).
Initial conditions for this stage are the final conditions for Stage 1 of the experiment.
Using the coordinate system as previously defined, the boundary conditions for Stage 1
are as follows.
Top Horizontal Boundary:
water flux = 0
oil flux = 136 cm3 h~l
oil flux = 0
Bottom Horizontal Boundary:
oil and water flux = 0
Left-Hand Vertical Boundary:
hw + z = 36.5 cm
water flux = 0
oil flux = 0
Right-Hand Vertical Boundary:
hw + z = 34.5 cm
water flux = 0
oil flux = 0
everywhere
at x = 50 cm
elsewhere
everywhere
0 < z < 36.5 cm
z > 36.5 cm
everywhere
0 < z < 34.5 cm
z > 34.5 cm
everywhere
Stage 21 LNAPL Redistribution. At the end of LNAPL infiltration, oil was allowed to
redistribute for 170 hours. Boundary conditions during this phase of the experiment
were the same as those during Stage 2 except that the oil flux is zero everywhere on the
top boundary.
176
-------
Stage 4 ;. Water Flushing. To simulate the effects of water flushing on LNAPL
dissolution and distribution, nine intermittent pulses of water at constant flux were
applied to the soil surface over a period of 76 hours. Each pulse has a volume loading
rate of 2.84 cm3/s per unit width for a duration of 78 seconds. Pulses were applied at 0,
8, 14, 22, 30, 39, 48, 56, 63 hours, respectively. All other boundary conditions were the
same as Stage 3.
Stage £ - LNAPL Entrapment. To simulate LNAPL entrapment by the aqueous phase,
the water table elevation was instantaneously raised at inlet and outlet by 11.43 cm and
this constant head condition was maintained for 60 hours. The boundary conditions on
the top and bottom of the cell were the same as during Stage 3. However, on the
vertical sides, the conditions for water phase were:
Left-Hand Vertical Boundary:
hw + z = 47.8 cm 0 < z < 47.8 cm
water flux =0 47.8cm < z
Right-Hand Vertical Boundary:
hw + z - 45.8 cm 0 < z < 45.8 cm
water flux = 0 45.8 cm < z
6.2.2 Model calibration and numerical analysis
The experimental results in the LNAPL experiment were compared with simulations
performed using the multiphase flow and transport model described in Chapter 4. The
length (148 cm) and height (101 cm) of the experimental cell were discretized
numerically by a mesh with 16 columns and 21 rows resulting in 336 nodes and 300
elements. The cell width (8.5 cm) was taken as an axis of symmetry (i.e., 2-D Cartesian
analysis). Oil-water viscosity ratio, oil specific gravity, oil surface tension and oil-water
interfacial tension for the NAPL mixture were calculated as volume-weighted averages
of the component values, using the data in Table 6.1. Fluid scaling factors were
estimated as
177
-------
O ^ U Q i V OW
Pow ~ O
where ff0 and aow are the volume fraction-weighted average surface tension and oil-
water interfacial tensions, respectively. For purposes of modeling, soltrol and
iodoheptane were considered as a single pseudo-species of negligible solubility referred to
as "inert." Tabulated pure component densities and diffusion coefficients were obtained
from the literature. Equilibrium oil-water partition coefficients were computed as the
ratio of water to oil phase concentrations determined from batch experiments on the
NAPL mixture, air-water partition coefficients were obtained from tabulated Henry's
constants for the compounds, and solid-water partition coefficients were assumed
negligible since the organic content of the sand used in the experiments was very low.
The true (average) density for the inert components was specified to enable correct
calculation of the bulk oil phase density by the code which is corrected for
compositional changes, and the oil-water partition coefficient was assigned a large value
to prevent partitioning. Other parameters for the inert components are arbitrary.
Table 6.3 Soil and fluid properties used for LNAPL simulation.
Soil properties: Bulk fluid properties:
KIW
£U/«
^n
sm
a
n
AL
AT
= 231 cm h-1
r = 231 cm h'1
' = 0.42
= 0.00
= 0.05 cm"1
= 3.50
= 10.0 cm
= 0.01 cm
0ao = 2.84
/?ao= 1.543
r,ro = 4.620
Pro = 0.826
Component properties:
cm2 /T1
D~/a, cm2 h'1
p0, g cm'3
Benzene
0.035
317.0
260.8
0.24
0
0.877
Toluene
0.027
281.0
1058.5
0.28
0
0.862
"Inert"
0.0001
0.001
106
0.2
0
0.785
178
-------
The saturated hydraulic conductivity of the soil to water was computed from horizontal
flow rates through the cell under a fixed gradient and the soil was assumed to be
hydraulically isotropic. Soil porosity was determined from gamma attenuation
measurements and van Genuchten parameters (a, n, Sm) were estimated by matching
measured and simulated water outflow and vertical water saturation distributions at the
end of the Stage 1 drainage period. Dispersivities were assigned small enough values to
avoid oscillations in the solution with the selected mesh. A summary of the soil and the
fluid properties used in the simulation of the LNAPL experiment are given in Table 6.3.
Stage 2 involved infiltration of oil in the domain. During this stage a flux of 7 cm/hr
was applied on the top of element 291. Simulation was terminated when the cumulative
oil infiltration was 165 cm3/cm. Initially, the water table was allowed to vary linearly
from the left boundary to the right boundary. Water pressures were specified on the left
and right boundaries as hydrostatic through Stage 2 corresponding to a water table
elevation of 36.5 and 34.5 cm at the left and right sides, respectively. The rest of the
domain had no flow boundary conditions for all fluid phases.
Stage 3 involved redistribution for 31 hours after the end of the infiltration. The
boundary conditions for water phase during this stage were unchanged, while a zero flux
was imposed for the oil phase on all boundaries. The fluxes in the oil phase during the
early part of the redistribution period were relatively high. It was observed during some
earlier simulations that a more stable transport solution is obtained when liquid phase
fluxes are small. Therefore, simulation of transport was delayed until the next stage.
Nine cycles of water application and drainage were imposed during Stage 4. Each cycle
consisted of an application of 2,000 cm3 of water applied at the top surface of the cell in
78 seconds, followed by 8 hours of redistribution, while all other flow boundary
conditions were same as in Stage 3. During Stage 4, simultaneous flow and transport of
the mixture comprised of benzene, toluene and "inert" components was simulated.
Initial oil phase concentrations for benzene, toluene and "inert" species were assigned to
be 0.070, 0.075 and 0.68 g/cm3, respectively. Initial water and gas phase concentrations
were assumed to be in equilibrium with the oil phase and to be zero at locations where
no oil phase occurred. Since transport prior to Stage 4 occurred primarily due to
convection of the oil phase, this assumption for initial conditions should lead to little
error. For transport, a type-2 boundary condition (zero concentration gradient) was
applied at the right vertical boundary to allow mass efflux, and a type-3 boundary
179
-------
condition with zero influent concentration was applied elsewhere.The initial time step
during the simulation was 0.0001 hours and the time incremental factor and the
maximum time step were 1.03 and 0.05 hours, respectively.
During Stage 5, there was no flow of water through the top boundary. The water table
was raised in the cell by 11.4 cm on the left and right sides equally, while all other
boundary conditions for all phases remained unchanged. Flow and transport was
simulated.
6.2.8 Results
Stage 1 - Water Drainage. A comparison of observed and simulated water saturations at
the end of Stage 1 is shown in Figure 6.2. Good correspondence is observed. The van
Genuchten parameters, which describe the air-water capillary pressure curve, were
estimated by matching the drainage data.
Stage 2 i LNAPL Infiltration. Water and oil saturation distributions at the end of Stage
2 are shown in Figures 6.3 and 6.4. The oil plume at this time is largely confined to the
unsaturated zone and capillary fringe, with some vertical penetration into the aquifer.
Correspondence between observed and predicted water and oil saturations is good.
Stage 2 i LNAPL Redistribution. Water and oil saturation distributions at the end of
Stage 3 are shown in Figures 6.5 and 6.6. Water saturation contours differ only slightly
from those of Stage 2. Oil has drained from the unsaturated zone and spread laterally
on the water table. Measured oil saturations in the unsaturated zone exhibit much
greater variability than predicted values, perhaps reflecting some heterogeneity in the
sand packing. In general the degree of oil drainage, vertical penetration and lateral
spreading for the simulations is greater than was observed. Our suspicion is that this
may reflect simplifications invoked in the handling of terms in the hysteresis model, but
we have yet to confirm this by performing more rigorous model simulations. Measured
and predicted aqueous concentrations of benzene and toluene at the end of Stage 3 are
shown in Figures 6.7 and 6.8, respectively. Since aqueous samples could only be taken in
the saturated zone, the number of measured locations is rather small. The simulations
show much more vertical penetration of dissolved species into the aquifer than observed.
This reflects the greater predicted oil phase penetration as well as possible numerical
dispersion effects. Overshoot in the transport solution occurred at the downstream face.
180
-------
Stage 4 i Water Flushing. Water and oil saturation distributions at the end of Stage 4
are shown in Figures 6.9 and 6.10. Trends observed at the end of Stage 3 towards too
rapid oil movement are continued and accentuated at this time. Deviations between
measured and predicted aqueous concentrations of benzene and toluene at the end of
Stage 4 (Figures 6.11 and 6.12) are likewise similar to those observed for Stage 3.
Stage £ - LNAPL Entrapment. At the end of Stage 5, the water saturation contours
have been clearly displaced upwards (Figure 6.13). Observed and simulated oil
saturations (Figure 6.14) exhibit small upward displacements, reflecting the fact that oil
has already distributed to close to a "residual" state at which its permeability is very
low so that movement is impeded. Additional aqueous sampling ports at higher
elevations show high concentrations corroborated by the model, but the problem of
excessive vertical penetration of the simulated aqueous plume remains (Figures 6.15 and
6.16).
Saturation-time curves at selected locations. To evaluate the model predictions in
greater detail, several locations in the cell were selected to compare water and oil
saturation histories over the course of the experiment. The 6 locations are shown in
Figure 6.17. Measured and predicted water and oil saturations versus time for the
selected locations are shown in Figures 6.18 - 6.23. The results show generally good
agreement between observed and predicted data, with a tendency to underpredict the
oil saturation within the oil plume as time passes due to vertical penetration and
horizontal spreading which is somewhat overpredicted by the model.
181
-------
100 IB
loo F r
N
-04-
-OA-
-OJ-
-0.4-
-0.1
*M-
OJ
10
X (cm)
100 120
140
Figure 6.2 Water saturation at the end of Stage 1 for LNAPL experiment.
100 IB
=1 100
140
Figure 6.3 Water saturation at the end of Stage 2 for LNAPL experiment.
182
-------
1» HO
140
Figure 6.4 Oil saturation at the end of Stage 2 for LNAPL experiment.
140
Figure 6.5 Water saturation at the end of Stage 3 for LNAPL experiment.
183
-------
4i
M
140
Figure 6.6 Oil saturation at the end of Stage 3 for LNAPL experiment.
IIrnn 1 1 1 n Q
20
Figure 6.7 Benzene concentration at the end of Stage 3 for LNAPL experiment.
184
-------
too m
Figure 6.8 Toluene concentration at the end of Stage 3 for LNAPL experiment.
140
ioo p
U
s^
N
i I I I i
I ^ 100
-fcj-
-u-
-OJ-
-M.
40
*> to
X (cm)
100 1» 140
Figure 6.9 Water saturation at the end of Stage 4 for LNAPL experiment.
185
-------
- M
- M
140
Figure 6.10 Oil saturation at the end of Stage 4 for LNAPL experiment.
tO *) 100 120 140
- to
- 40
- 10
140
X (cm)
Figure 6.11 Benzene concentration at the end of Stage 4 for LNAPL experiment.
186
-------
100 IB) 140
Figure 6.12 Toluene concentration at the end of Stage 4 for LNAPL experiment.
too F
- »
Figure 6.13 Water saturation at the end of Stage 5 for LNAPL experiment.
187
-------
100 IS
- 10
140
Figure 6.14 Oil saturation at the end of Stage 5 for LNAPL experiment.
100 12D 140
U
M
10-
100
40 « « 100
X (cm)
140
Figure 6.15 Benzene concentration at the end of Stage 5 for LNAPL experiment.
188
-------
100 1JD 140
Figure 6.16 Toluene concentration at the end of Stage 5 ior LNAPL experiment.
OII Sourct
0.6 m
\\\\\\\\\\\\
Sell Surface
\\x\\\\\\\\\\\\\\\\\\\\\\\\\\\
x 1.
x 2. x 3.
x 6.
X 4.
M44.44)
(114,)
m
t.em-
outflow
Figure 6.17 Measurement positions for plots of concentration versus time.
189
-------
1.(
0.80 -
0.60 ^
c
o
"p
3
-------
irations
"o
1.00 -
-
0.80 -_
~
:
0.60 '-
0.40 -_
0.20 '-
.00 -
o.c
I
I .
Oil Infiltration
I
i£T ^~> *A *
i A 1"*' * AA>* " *
i
i
i
i
i Oil Redistribution
1
1
1
1
L, g _
J <*-\v*5vT v ^jy-V* "
,i ' - - - - .
1
^r i i i i i i i i | i i i i i i r~
)0 100.00
Time (hr
.
4
.
X1
^ * ' .
r* v.*/0 100.00 200.00 " 300.00
jkAA^A Water saturation (observed
Oil saturation (observed)
Water saturation (simulated)
Oil saturation (simulated)
Time (hrs)
Figure 6.21 Water and oil saturation versus time for position
191
-------
1 .UW
0.80 -_
-
g 0.60 -_
.2 -
"o ;
-4_< -
(8 0.40 -_
0.20 i
o no
>
!
'
i
\» . * * *
* "* / "** *
Oil Infiltration
Oil Redistribution
«* -n"" «VrV * Xv*v*> "
i i 1 1 1 1 > i i i 1 1 1 1 1 1
* V ^^
» * * *
****** *
*
A A
Water
Infiltration
"W^VVrf*
1 1 1 I 1 1 I
4t
*' "V^»
Oil
Entrapment
^^* -BB_.
1 1 1 | 1
AAAAA Water saturation (observed
Oil saturation (observed)
Water saturation (simuloU "
Oil saturation (simulated)
0.00 100.00 200.00
Time (hrs)
300.00
Figure 6.22 Water and oil saturation versus time for position
1.00
0.80 -
0.50 ^
o
'*;
o
0.40 ^
0.20 -
0.00
Oil Infiltration
Oil Redistribution
Water
Infiltration
Oil
Entrapment
IAHA Water saturation (observed
immmm oil saturation (observed)
Water saturation (simulated)
Oil saturation (simulated)
'^^^*"| I I 1 I I I I I 1 I I I I I I I I I I I I 11 1] 1"
0.00 100.00 200.00 300.00
Time (hrs)
Figure 6.23 Water and oil saturation versus time for position
19:
-------
6.3 DNAPL Experiment
6.S.I Experimental methods
The second experiment was conducted to study the movement of dense nonaqueous
phase liquid (DNAPL). The experimental setup used in this study is the same as that
employed for the LNAPL experiment. The DNAPL consisted of a mixture of
tetrachloroethylene (also known as tetrachloroethene, perchloroethylene or PER),
soltrol, toluene and 1-iodoheptane, mixed in a 4:4:1:1 ratio by volume, respectively,
yielding an overall specific gravity of 1.16. Properties of PER, soltrol, toluene and
iodoheptane are given in Table 6.1.
As in the LNAPL experiment, iodoheptane was included to provide a larger gamma
attenuation contrast between the NAPL and water to provide better resolution of fluid
saturations. The material used in this experiment was a medium sand, consisting of a
mixture of 2% very fine sand, 17% fine sand, 39% medium sand, 28% coarse sand and
14% very coarse sand. The gamma attenuation coefficients of various components are
shown in Table 6.4.
Table 6.4 Gamma attenuation coefficients of porous medium and liquids
for DNAPL experiment.
Coefficient
WW(-1)
iv.r-1;
PI (sm-cm^)
Am-241
0.19577
0.94179
0.24712
Cs-137
0.08361
0.09118
0.07466
The cell was filled with soil to a depth of 99 cm and saturated with water from the
bottom. The average bulk density, evaluated under saturated conditions, was 1.46 g
cm'3. The experiment was performed in three stages as described below.
193
-------
Stage 1 - Water Drainage. A two-phase air-water experiment was performed to
calibrate the k-S-P model. The soil in the cell was initially water saturated and was
allowed to drain for 9 hours after the water table was lowered to 36.5 cm and
maintained at this level. In this experiment the origin for the x (horizontal) and z
(vertical) coordinates was the lower left corner. The boundary conditions employed
during this stage were as follows.
Top and bottom boundaries:
qw = 0 everywhere
Left Vertical Boundary:
qw = 0 everywhere
Right Vertical Boundary:
hw + z = 36.5 cm (0 < z < 36.5 cm)
qw = 0 (z > 36.5 cm)
Water saturations and pressures as well as water outflow were monitored. Total outflow
over a 24 hour period was 15,500 cm3. After the drainage period, a horizontal gradient
of 1.35% was established by imposing a constant head of 36.5 cm on the left-hand side
of the cell (Figure 6.1) while the water table was dropped to 34.5 cm on the right
boundary. At steady state, the aqueous flux density was measured and was used to
compute the saturated hydraulic conductivity of the soil.
Stage 2 r DNAPL Infiltration. Over a period of 3 hours, 500 cm3 of NAPL was added on
a strip 50 cm from top left-hand corner (Figure 6.1) yielding an average application rate
of 167 cm3 h~l. Using the coordinate system previously defined, the boundary conditions
can be expressed as follows.
Top Horizontal Boundary:
water flux = 0 everywhere
oil flux = 167 cm3 h~l at i = 50 cm
oil flux = 0 elsewhere
Bottom Horizontal Boundary:
water and oil flux = 0 everywhere
194
-------
Left-Hand Vertical Boundary:
hw + z = 36.5 cm 0 < z < 36.5 cm
water flux = 0 z > 36.5 cm
oil flux = 0 everywhere
Right-Hand Vertical Boundary:
hw + z = 34.5 cm 0 < z < 34.5 cm
water flux = 0 z > 34.5 cm
oil flux = 0 everywhere
Stage 2. 2. DNAPL Redistribution. At the end of NAPL infiltration, oil was allowed to
redistribute for 122 hours. Boundary conditions during this stage of the experiment were
the same as those during Stage 2 except that the oil flux was zero everywhere on the
top boundary.
6.S.2 Model calibration and numerical analysis
Experimental results in the DNAPL experiment were compared with simulations
performed using the multiphase flow and transport model described in Chapter 4. The
length and height of the experimental cell were discretized numerically by a mesh with
19 columns and 19 rows resulting in 361 nodes and 324 elements. The cell width (8.5
cm) was taken as an axis of symmetry (i.e., 2-D Cartesian analysis). Oil-water viscosity
ratio, oil specific gravity, oil surface tension and oil-water interfacial tension for the
NAPL mixture were calculated as volume-weighted averages of the component values,
using the data in Table 6.1. Fluid scaling factors were estimated as
000 =
V =
'ow
where a0 and aow are the volume fraction-weighted average surface tension and oil-
water interfacial tensions, respectively. For purposes of modeling, soltrol and
iodoheptane were considered as a single pseudo-species of negligible solubility referred to
as "inert." Tabulated pure component densities and diffusion coefficients were obtained
from the literature. Equilibrium oil-water partition coefficients were computed as the
195
-------
ratio of water to oil phase concentrations determined from batch experiments on the
NAPL mixture, air-water partition coefficients were obtained from tabulated Henry's
constants for the compounds, and solid-water partition coefficients were assumed
negligible since the organic content of the sand used in the experiments was very low.
The true (average) density for the inert components was specified to enable correct
calculation of the bulk oil phase density by the code which is corrected for
compositional changes, and the oil-water partition coefficient was assigned a large value
to prevent partitioning. Other parameters for the inert components are arbitrary.
The saturated hydraulic conductivity of the soil to water was computed from horizontal
flow rates through the cell under a fixed gradient and the soil was assumed to be
hydraulically isotropic. Soil porosity was determined from gamma attenuation
measurements and van Genuchten parameters (a, n, Sm) were estimated by matching
measured and simulated water outflow and vertical water saturation distributions at the
end of the Stage 1 drainage period. Dispersivities were assigned values small enough to
avoid numerical oscillations with the selected mesh. A summary of the soil and the fluid
properties used in the simulation of the DNAPL experiment are given in Table 6.5.
Table 6.5 Soil and fluid properties used for DNAPL simulation.
Soil properties: Bulk fluid properties:
K,w =226 cm /f1 /? = 2.587
K '= 226 cm /f1 0ao= 2.587
^ w = 0.43 f,ro = 2.700
5m = 0.01 /jro = 1.158
a =0.15 cm'1
n =2.0
AL = 5.0 cm
AT = 0.2 cm
Component properties:
PER
D0°w, cm2 h'1 0.022
Da°a, cm2 h~l 266.0
Toluene "Inert"
0.027 0.00001
_aa, _ 283.0 0.001
roo 6920. 1110. 106
raa 0.35 0.28 0.2
rM o o o
pa,gcm-3 1.60 0.862 0.785
196
-------
Stage 2 involved infiltration of oil in the domain on a strip source 8.5 cm wide, centered
50 cm from the top left corner of the cell. Stage 2 simulation was terminated after 3
hours when the cumulative oil infiltration reached 500 cm3. Initially, the water table
was allowed to vary linearly from the left boundary to the right boundary. The water
pressure distribution was hydrostatic and no oil was present initially in the system.
Water pressures were specified on the left and right boundaries as hydrostatic through
Stage 2 corresponding to a water table elevation of 36.5 and 34.5 cm at the left and
right sides, respectively. The rest of the domain had no flow boundary conditions for all
fluid phases. Transport was not simulated during Stage 2 since the duration was only 3
hours and transport was dominated by oil phase convection.
Stage 3 involved redistribution for a period of 122 hours after the end of Stage 2. The
boundary conditions for water phase during this stage were unchanged, while a zero flux
was imposed for the oil phase on all boundaries. Initial oil phase concentrations for
PER, toluene and "inert" species were assigned the values of the initial DNAPL
mixture. Initial water and gas phase concentrations were assumed to be in equilibrium
with the oil phase and to be zero at locations where no oil phase occurred. Boundary
conditions for transport were a type-2 boundary condition (zero concentration gradient)
on the right vertical boundary to allow mass efflux, and a type-3 boundary condition
with zero influent concentration on the upstream boundary to account for inflow of
clean water.
6.S. S Results
A comparison of observed and simulated water saturations at the end of Stage 1, used
to estimate van Genuchten air-water capillary pressure curve parameter, is shown in
Figure 6.24. Observed water saturations below the water table ranged from 82-95%,
indicating significant quantities of trapped air. Since the constitutive relations used for
the simulations did not consider air entrapment, the model predicts full saturation of
water below the water table.
Observed and predicted water and oil saturation distributions at the end of the
redistribution period (Stage 3) are compared in Figures 6.25 and 6.26. Water saturations
are generally slightly underestimated by the model, except below the water table due to
trapped air. The model predicts an oil phase plume which is largely confined to a
deformed spheroidal region above the water table, although some oil is predicted to
197
-------
have sunk through the aquifer and spread over the tank bottom. Between the water
table and the bottom pool, oil saturations after drainage are predicted to be less than
1%. Measurements could not be made near the tank bottom, so confirmation of oil
spreading along the tank bottom could not be definitely confirmed. However, nonzero
oil saturations at measurement locations 15 cm from the tank bottom suggest oil
penetration below the water table. The predicted and measured oil distributions in the
unsaturated zone agree with one interesting exception. The observed data show a small
amount of oil has apparently spread laterally downgradient along the top of the
capillary fringe at oil saturations in the vicinity of 1%. The simulation does not predict
such lateral migration along the water table. This is an interesting and unexpected
phenomena. However, whether this is due to physical effects disregarded by the model
(e.g., nonzero spreading pressures at oil-water interfaces) or to subtle heterogeneities
due to packing anomalies cannot be ascertained.
Observed and predicted aqueous concentrations of PER and toluene at the end of Stage
3 are compared in Figures 6.27 and 6.28. Substantial concentrations along the capillary
fringe downgradient of the source confirm the lateral migration phenomena noted based
on oil saturation measurements. Sampling locations under the water table and
immediately beneath the source indicate no dissolved phase contaminants, although
gamma measurements indicate DNAPL at or below these depths. This suggests
downward migration of DNAPL through the saturated zone occurs through isolated
fingers rather than as a uniform front.
198
-------
0 10 20 30 40 50 60 70 80 90 100 110 120 130 140
1 1 1 I 1 1
90 - S *?
.? / «^
8U 0.2 U2 i
*
70 - , S
*
60 T,-, /-, . J
u.j t,j s
«
40 -0.0 t7.U » .
«?«?«?
30 - **
20 -
10 -
Q 1 1 1 1 1 1
0 10 20 30 40 50 60
i i
*
» *
*
4*J£-
IL ft?
.' *«
* *
.* J* ^ ^ - 7°
* *
S 3 - tfj "* - 60
Jli « bU
V « ri*^ 4^ ~
"" » "" * J" U.ti -
£ / .»" .^
* *
.* .< - 30
- 20
f .*
- 10
1 1 1 1 1 1 Q
90 100 110 120 130 140
Figure 6.24 Water saturation distributions at the end of Stage 1 for DNAPL experiment.
0.00 18.50 37.00 55.50 74.00 92.50 111.00 129.50 148.00
33. UU
86.63
74.25
61.88
49.50
37.13
24.75
12.38
n nn
i i i i i i i i i i i i i i i i i i i i i i i
A $ N* & N1* £*
o> o- *' * O' O-
nt- as "*' -'*' ' '" -*" n-
O' O O* O O' O
~
°-55 "' ~f S? rf" °-5^«-~# ^r5 ^# 0.55
ft 05 £ ._, £f 03^' J _ J" f J 03^
. . . °, *.
i i I i t i I i i i i i i i i i I i i l i i i
86.63
74.25
61.88
49.50
37.13
24.75
12.38
nnn
0.00 18.50 37.00 55.50 74.00 92.50 111.00 129.50 148.00
Figure 6.25 Water saturation at the end of Stage 3 for DNAPL experiment.
199
-------
0.00
20.00 40.00
60.00
80.00 100.00 120.00 140.00
80.00
60.00
40.00
20.00
0.00
I I I I
0.00
20.00
40.00 60.00 80.00 100.00 120.00 140.00
80.00
60.00
40.00
20.00
0.00
Figure 6.26 Oil saturation at the end of Stage 3 for DNAPL experiment.
0.00 18.50 37.00 55.50 74.00 92.50 111.00 129.50 148.00
0.00
0.00
0.00
18.50 37.00 55.50 74.00 92.50 111.00 129.50 148.00
Figure 6.27 PER concentration at the end of Stage 3 for DNAPL experiment.
200
-------
0 10 20 30 40 50 60 70 80 90 100 110 120 130 HO
90
80
70
60
50
40
30
2C
10
0
\\r
o
i i i i i i r
s
o
10 20 30 40 50 60 70 80 90 100 110 120 130 140
90
80
70
60
50
40
30
20
10
0
Figure 6.28 Toluene concentration at the end of Stage 3 for DNAPL experiment.
201
-------
REFERENCES CITED
Abdul, A. S. and T. L. Gibson, Equilibrium batch experiments with six polycylic
aromatic hydrocarbons and two aquifer materials, Hazardous Waste & Hazardous
Materials, 3, 125-137, 1986.
Abriola, L. M., Multiphase migration of organic compounds in a porous medium,
Lecture Notes in Engineering, Brebbia and Orszag (ed.) No. 8, Springer-Verlag, Berlin,
1984.
Abriola, L. M. and G. F. Finder, A multiphase approach to the modeling of porous
media contamination by organic compounds, 1. Equation development, Water Resour.
Res., 21, 11-18, 1985a.
Abriola, L. M. and G. F. Finder, A multiphase approach to the modeling of porous
media contamination by organic compounds, 2. Numerical simulation, Water Resour.
Res., 21, 19-26, 1985b.
Aienkiewicz, 0. C., The finite element method, McGraw-Hill, NY, 1986.
Bear, J., Dynamics of fluids in porous media, 764p., Elsevier NY, 1972.
Bloomsburg, G. L. and A. T. Corey, Diffusion of entrapped air from porous media,
Hydrology Papers No. 5, Colorado State University, 27 p., 1964.
Burdine, N. T., Relative permeability calculations from pore size distribution data,
Trans. AIME, 198, 71-77,1953.
Cooley, R. L., A finite difference method for unsteady flow in variably saturated porous
media: Applications to a single pumping well, Water Resour. Res., 7, 1607-1625, 1971.
Corey, A. T., The interrelation between gas and oil relative permeabilities, Producer's
Monthly, 19, 38-41,1954.
202
-------
Falta, R. W. Jr. and I. Javandel, A numerical method for multiphase multicomponent
contaminant transport in groundwater systems, Trans. Amer. Geophys. Union, 68, 1284,
1987.
Faust, C. R., Transport of immiscible fluids within and below the unsaturated zone: A
numerical model, Water Resour. Res., 21, 587-596, 1985.
Forsyth, P. A., Simulation of nonaqueous phase groundwater contamination, Adv.
Water Resour., 11, 74-83, 1988.
Forsyth, P. A. and P. H. Sammon, Practical considerations for adaptive implicit
methods in reservoir simulation, J. Comp. Phys., 62, 265-281, 1986.
Fritton, D. D., Resolving time, mass absorbtion coefficient and water content with
gamma-ray attenuation, Soil Sci. Soc. Am. J.,33,651-655, 1969.
Huyakorn, P. S. and K. Nilkuha, Solution of transient transport equation using an
upstream finite element scheme, Appl. Math. Modeling, 3, 7-17, 1979.
Huyakorn, P. S. and G. F. Pinder, Computational method in subsurface flow, Academic
Press, Inc., NY, 1983.
Huyakorn, P. S., S. D. Thomas and B. M. Thompson, Techniques for making finite
elements competitive in modeling flow in variably saturated porous media, Water Res.
Research, 20, 1099-1115, 1984.
Kaluarachchi, J. J. and J. C. Parker, Effects of hysteresis on water flow in the
unsaturated zone, Water Resour. Res., 23, 1967-1976, 1987.
Kaluarachchi, J. J., and J. C. Parker, An efficient finite element method for modeling
multiphase flow, Water Resour. Res., 25, 43-54, 1989.
Kool, J. B. and J. C. Parker, Development and evaluation of closed-form expressions for
hysteretic soil hydraulic properties, Water Resour. Res., 23, 105-114, 1987.
203
-------
Kuppusamy, T., J. Sheng, J. C. Parker and R. J. Lenhard, Finite element analysis of
multiphase immiscible flow through soils, Water Resour. Res., 23, 625-632, 1987.
Land, C. S., Calculation of imbibition relative permeability for two- and three-phase
flow from rock properties, Soc. Petro. Engr. J., 6, 149-156, 1968.
Lenhard, R. J., J. H. Dane, J. C. Parker and J. J. Kaluarachchi, Measurement and
simulation of transient three-phase flow in porous media, Water Resour. Res., 24, 853-
863, 1988.
Lenhard, R. J. and J. C. Parker, Measurement and prediction of saturation-pressure
relationships in three-phase porous media systems, J. Contain. Hydrol., 1, 407-424,
1987.
Lenhard, R. J. and J. C. Parker, Experimental validation of the theory of extending
two-phase saturation-pressure relations to three fluid phase systems for monotonic
saturation paths, Water Resour. Res., 24, 373-380, 1988.
Leverett, M. C., Capillary behavior in porous solids. Trans. AIME, 142, 152-169, 1941.
Miller, C. T. M. M. Poirier-McNeill and A. S. Mayer, Dissolution of trapped
nonaqueous phase liquids: mass trasfer characteristics, Water Resour. Res, 26, 2783-
2796, 1990.
Millington, R. J. and J. P. Quirk, Permeability of porous media, Nature (London) 183,
387-388, 1959.
Mualem, Y., A new model for predicting the hydraulic conductivity of unsaturated
porous media, Water Resour. Res., 12, 513-522, 1976.
Nofziger, D. L., and D. Swartzendruber, Material content of binary physical mixtures as
measured with a dual-energy beam of gamma rays, J. Appl. Phys. 45, 5443-5449, 1974.
Osborne M. and J. Sykes, Numerical modeling of immiscible organics transport at the
Hyde Park landfill, Water Resour. Res., 22, 25-33, 1986.
204
-------
Parker, J. C. and R. J. Lenhard, A model for hysteretic constitutive relations governing
multiphase flow, 1. Saturation-pressure relations, Water Resour. Res., 23, 2187-2196,
1987.
Parker, J. C., R. J. Lenhard and T. Kuppusamy, A parametric model for constitutive
properties governing multiphase flow in porous media, Water Resour. Res., 23, 618-624,
1987.
Saraf, D. N., J. P. Batycky, C. H. Jackson and D. B. Fisher, An experimental
investigation of three-phase flow of water-oil-gas mixtures through water-wet
sandstones, paper presented at Califronia Regional Meeting, Soc. Pet. Eng., San
Francisco, Calif., SPE paper #10761, March 24-26, 1982.
Stillwater, R. and A. Klute, Improved methodology for a collinear dual-energy gamma
radiation system, Water Resour. Res., 24, 1411-1422, 1988.
Stonestrom, D. A. and J. Rubin, Water content dependence of trapped air in two soils,
Water Resour. Res., 25, 1947-1958, 1989.
Thomas, G. W. and D. H. Thurnau, Reservoir simulation using an adaptive implicit
method, Soc. Pet. Eng. J., 23, 759-768, 1983.
van Genuchten, M. T., A closed form equation for predicting the hydraulic conductivity
of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892-898, 1980.
Willhite, G. P., Waterflooding, SPE Textbook Series Vol. 3, Society of Petroleum
Engineers, Richardson, TX, 326 p., 1986.
205
-------
APPENDIX: PUBLICATIONS RESULTING FROM THIS PROJECT
Anderson, M. A. and J. C. Parker, Sensitivity of organic contaminant transport and
persistence models to Henry's law constants: case of polychlorinated biphenyls, Water,
Air and Soil Pollution, 50, 1-18, 1990.
Kaluarachchi, J. J. and J. C. Parker. An efficient finite element method of modeling
multiphase flow in porous media. Water Resour. Res., 25, 43-54, 1989.
Kaluarachchi, J. J. and J. C. Parker. Modeling multicomponent organic transport in
three fluid phase porous media, J. Contaminant Hydrol., 5, 349-375, 1990.
Kaluarachchi. J. J., J. C. Parker and R. J. Lenhard, A numerical model for areal
migration of water and light hydrocarbon in unconfined aquifers, Adv. Water Resour.
Res., 13, 29-40, 1990.
Kaluarachchi, J. J., J. C. Parker and R. J. Lenhard, Modeling flow in three fluid phase
porous media with non wetting fluid entrapment, Proc. Conf, on Subsurface
Contamination by Immiscible Fluids. Calgary. Canada, April 17-20, 1990.
Kaluarachchi, J. J., J. C. Parker. Multiphase flow with a simplified model for oil
entrapment, Transport in Porous Media, (in press).
Katyal, A. K., J. J. Kaluarachchi and J. C. Parker, Evaluation of methods for
improving the efficieny and robustness of multiphase flow models, Proc. Conf. on
Subsurface Contamination by Immiscible Fluids, Calgary, Canada, April 17-20, 1990.
Lenhard, R. J., J. H. Dane, J. C. Parker and J. J. Kaluarachchi, Measurement and
simulation of transient three-phase flow in porous media, Water Resour. Res., 24, 853-
863, 1988.
Lenhard, R. J. and J. C. Parker, A model for hysteretic constitutive relations governing
multiphase fluid flow, 2. Permeability-saturation relations, Water Resour. Res., 23,
2197-2206, 1987.
206
-------
Lenhard, R. J. and J. C. Parker, Experimental validation of the theory of extending
two-phase saturation-pressure relations to three-fluid phase systems for monotonic
drainage, Water Resour. Res., 24, 373-380, 19SS.
Lenhard, R. J., J. C. Parker and S. Mishra, On the correspondence between Brooks-
Corey and van Genuchten models, J. Irrig. and Drainage, ASCE, 115, 744-751, 1989.
Lenhard, R. J. and J. C. Parker, A model for hysteretic constitutive relations governing
multiphase fluid flow, 3. Refinements and numerical simulations, Water Resour. Res.,
25, 1727-1736, 1989.
Lenhard. R. J. and J. C. Parker. Modeling and laboratory validation of multiphase fluid
hysteresis, Proc. Workshop in Indirect Methods for Estimating Hydraulic Properties of
Unsaturated Soils, Riverside CA, 1989.
Lenhard, R. J. and J. C. Parker, Modeling multiphase fluid hysteresis and comparing
results to laboratory investigations, Proc. Indirect Methods for Estimating Hydraulic
Properties of Unsaturated Soils, Riverside, CA, Oct. 11-13, 1989.
Lenhard, R. J. and J. C. Parker, Estimation of free hydrocarbon volume from fluid
levels in observation wells, Ground Water, 28, 57-67, 1990.
Parker, J. C., Multiphase flow and transport in porous media, Reviews of Geophysics,
27, 311-328, 1989.
Parker, J. C. and J. J. Kaluarachchi, A numerical model for design of free product
recovery systems at hydrocarbon spill sites, Proc. 4th Int. Conf. on Solving Ground
Water Problems with Models, Indianapolis, IGWMC/NWWA, Feb 7-9, 1989.
Parker, J. C. and J. J. Kaluarachchi, Modeling multicomponent organic chemical
transport, Proc. Conf. on Subsurface Contamination by Immiscible Fluids, Calgary,
Canada, April 17-20, 1990.
Parker, J. C., J. J. Kaluarachchi and A. K. Katyal. Areal simulation of free product
recovery from a gasoline storage tank leak site. Proc. Petroleum Hydrocarbons and
Organic Chemical in Groundwater. National Water Well Assoc. Houston, 1988.
207
-------
Parker, J. C., J. J. Kaluarachchi, V. J. Kremesec and E. L. Hockman, Modeling free
product recovery at hydrocarbon spill sites, Proc. Petroleum Hydrocarbons and Organic
Chemical in Groundwater, NWWA, Houston, 1990.
Parker, J. C., J. J. Kaluarachchi and R. J. Lenhard, Experimental and numerical
investigations of constitutive relations governing multiphase flow, Proc. Conf. on
Validation of Flow and Transport Models for the Unsaturated Zone, Ruidoso, New
Mexico, 1988.
Parker, J. C., A. K. Katyal, J. L. Zhu and S. Mishra. Estimation of spill volume from
monitoring well networks, In Proc. 4th Nat. Outdoor Action Conf., Las Vegas, May
1990.
Parker, J. C., T. Kuppusamy, and B. H. Lien, Modeling multiphase organic chemical
transport in soils and groundwater, Proc. Groundwater Contamination: Use of Models in
Decision Making, Amsterdam, Martinus Nijhoff, October 1987.
Parker, J. C. and R. J. Lenhard, A model for hysteretic constitutive relations governing
multiphase fluid flow, 1. Saturation-pressure relations, Water Resour. Res., 23, 2187-
2196, 1987.
Parker, J. C. and R. J. Lenhard, Vertical integration of three-phase flow equations for
analysis of light hydrocarbon plume movement, Transport in Porous Media, 5:187-206,
1989.
Parker, J. C. and R. J. Lenhard, Determining three-phase permeability-saturation-
pressure relations from two-phase system measurements, J. Petrol. Sci. Eng., 4, 57-65,
1990.
208
-------