United States National Risk Management EPA/600/R-97/102
Environmental Protection Research Laboratory October 1997
Agency Ada, OK 74820
vvEPA NAPL: Simulator
Documentation
-------
NAPL: Simulator Documentation
by
Joseph Guarnaccia and George Finder
Research Center for Groundwater Remediation Design
The University of Vermont
Burlington, VT 05401
and
Mikhail Fishman
U. S. Environmental Protection Agency
Robert S. Kerr Environmental Research Center
Ada, Oklahoma 74820
Cooperative Agreement CR-820499
Project Officer
Thomas E. Short
U. S. Environmental Protection Agency
Robert S. Kerr Environmental Research Center
Subsurface Protection and Remediation Division
Ada, Oklahoma 74820
National Risk Management Research Laboratory
Office of Research and Development
U.S. Environmental Protection Agency
Cincinnati, OH 45268
-------
-------
DISCLAIMER
The U. S. Environmental Protection Agency through its Office of Research and Development
partially funded and collaborated in the research described here under assistance agreement number
CR-820499 to The University of Vermont. It has been subjected to the Agency's peer and
administrative review and has been approved for publication as an EPA document. Mention of trade
names or commercial products does not constitute endorsement or recommendation for use.
When available, the software described in this document is supplied on "as-is" basis without
guarantee or warranty of any kind, express or implied. Neither the United States Government (United
States Environmental Protection Agency, Robert S. Kerr Environmental Research Center), The
University of Vermont, nor any of the authors accept any liability resulting from use of this
software.
11
-------
FOREWORD
The U.S. Environmental Protection Agency is charged by Congress with protecting the Nation's
land, air, and water resources. Under a mandate of national environmental laws, the Agency strives
to formulate and implement actions leading to a compatible balance between human activities and
the ability of natural systems to support and nurture life. To meet these mandates, EPA's research
program is providing data and technical support for solving environmental problems today and
building a science knowledge base necessary to manage our ecological resources wisely, understand
how pollutants affect our health, and prevent or reduce environmental risks in the future.
The National Risk Management Research Laboratory is the Agency's center for investigation of
technological and management approaches for reducing risks from threats to human health and the
environment. The focus of the Laboratory's research program is on methods for the prevention and
control of pollution to air, land, water, and subsurface resources; protection of water quality in public
water systems; remediation of contaminated sites and ground water, and prevention and control of
indoor air pollution. The goal of this research effort is to catalyze development and implementation
of innovative, cost-effective environmental technologies; develop scientific and engineering
information needed by EPA to support regulatory and policy decisions; and provide technical support
and information transfer to ensure effective implementation of environmental regulations and
strategies.
This report focuses on the simulation of the contamination of soils and aquifers which results from
the release of organic liquids commonly referred to as Non-Aqueous Phase Liquids (NAPLs). The
approach used in this simulation is applicable to three interrelated zones: a vadose zone which is in
contact with the atmosphere, a capillary zone, and a water-table aquifer zone. The simulator
accommodates three mobile phases: water, NAPL, and gas. The numerical solution algorithm is
based on a Hermite collocation finite element discretization. The simulator provides an accurate
solution of a coupled set of non-linear partial differential equations that are generated by combining
fundamental balance equations with constitutive and thermodynamic relationships.
Clinton W. Hall, Director
Subsurface Protection and Remediation Division
National Risk Management Research Laboratory
-------
IV
-------
Abstract
A mathematical and numerical model is developed to simulate the transport and
fate of NAPLs in near-surface granular soils. The resulting three-dimensional,
three phase simulator is called NAPL. The simulator accommodates three mo-
bile phases: water, NAPL and gas, as well as water- and gas-phase transport of
NAPL contaminants. The numerical solution algorithm is based on a Hermite
collocation finite element discretization. Particular attention has been paid to the
development of a sub-model that describes three-phase hysteretic permeability-
saturation-pressure (k-S-P) relationships, and that considers the potential entrap-
ment of any fluid when it is displaced. In addition rate-limited dissolution and
volatilization mass transfer models have been included. The overall model has
been tested for self-consistency using mass balance and temporal and spatial con-
vergence analysis. The hysteretic k-S-P and mass exchange models have been
tested against experimental results. Several example data sets are provided, in-
cluding a setup of the artificial aquifer experiments being conducted at the EPA's
Subsurface Protection and Remediation Division of the National Risk Manage-
ment Research Laboratory in Ada, OK (formerly RSKERL) at this writing.
-------
CONTENTS
1 INTRODUCTION 11
1.1 PROBLEM 11
1.2 FOCUS 12
1.3 APPROACH 13
1.4 MODEL CAPABILITIES 14
1.5 ORGANIZATION 17
2 CONCLUSIONS 19
3 RECOMMENDATIONS 20
4 THEORETICAL DEVELOPMENT 21
4.1 OVERVIEW 21
4.2 MASS BALANCE EQUATIONS 22
4.3 PRIMARY VARIABLES 25
4.4 FLUID PROPERTIES 27
4.4.1 Density 27
4.4.2 Viscosity 28
4.4.3 Interfacial Tension 28
4.5 PHASE ADVECTION 29
4.6 DISPERSION COEFFICIENT 30
4.7 NATURAL DEGRADATION 30
4.8 MASS TRANSFER 31
4.8.1 Liquid-Liquid Mass Transfer 35
4.8.2 Liquid-Solid Mass Transfer 37
4.9 SUMMARY 38
5 HYSTERETIC k-S-P MODEL 40
5.1 CONCEPTUAL MODEL OVERVIEW 44
-------
5.2 TWO-PHASE MODEL 45
5.2.1 Entrapment and Release Sub-Model 46
5.2.2 Saturation-Pressure Sub-Model 54
5.2.3 Relative Permeability-Saturation Sub-Model 61
5.3 THREE-PHASE MODEL 64
5.3.1 Entrapment and Release Sub-Model 64
5.3.2 Saturation-Pressure Sub-Model 65
5.3.3 Relative Permeability-Saturation Sub-Model 67
5.4 MODEL IMPLEMENTATION 70
5.4.1 Phase Entrapment and Release 71
5.4.2 S-P Curve Pressure Scale Transition 71
5.4.3 S-P Curve Restriction Parameters 72
5.4.4 Mass Balance and Consistency 74
5.5 CAPILLARY PRESSURE SCALING 75
6 NUMERICAL MODEL DEVELOPMENT 78
6.1 FINAL FORM OF THE BALANCE EQUATIONS 78
6.1.1 Water-Phase Flow and Contaminant Transport 79
6.1.2 NAPL-Phase Flow 81
6.1.3 Gas-Phase Flow and Contaminant Transport 81
6.2 SEQUENTIAL, ITERATIVE SOLUTION 82
6.2.1 Time-Discrete/Linearized Form of the Flow Equations . . 82
6.2.2 Time-Discrete/Linearized Form of the Transport Equations 88
6.3 SPATIAL APPROXIMATION 90
6.4 IMPOSED CONDITIONS 92
6.4.1 Initial Conditions 92
6.4.2 Boundary Conditions 93
6.4.3 External Flux Conditions 98
6.5 DIAGNOSTIC TOOLS 98
6.5.1 Peclet Constraint 99
6.5.2 Time Step Control 104
6.5.3 Phase Discontinuities 106
6.6 SOLVING THE SYSTEM OF LINEAR EQUATIONS 107
6.7 SUMMARY OF SEQUENTIAL ITERATION 108
-------
7 SIMULATOR DOCUMENTATION 109
7.1 DATA INPUT 109
7.1.1 Set up for the include, f file 113
7.1.2 Data Input Driver 114
7.2 INPUT- AND OUTPUT-FILE DESCRIPTION 138
7.2.1 Input Files 138
7.2.2 Restart Files 140
7.2.3 Compilation Files 141
7.2.4 Screen output 141
7.2.5 Output files 143
8 MODEL TESTING AND EXAMPLE PROBLEMS 146
8.1 CONVERGENCE AND MASS BALANCE 146
8.1.1 Compatibility of the grid and the flow model 146
8.1.2 Analysis of the three-phase hysteretic k-S-P model 147
8.1.3 Analysis of the mass transfer model 150
8.2 COMPARISONS TO EXPERIMENTAL RESULTS 153
8.2.1 LNAPL Spill 153
8.2.2 DNAPL Spill 158
8.2.3 DNAPL Dissolution 166
8.2.4 DNAPL Vapor Transport 166
8.3 SOFTWARE 167
A PARAMETER LIST 177
B PARTICULARS OF HERMITE COLLOCATION 183
B.I Nodal Degrees of Freedom 183
B.2 Basis Function Definition 183
B.3 Hermite Interpolation of Capillary Pressure 184
B.4 Boundary Condition Specification 184
C INITIALIZE TRAPPING PARAMETERS 187
D PECLET CONSTRAINT 190
E SOURCE FILE DESCRIPTION 192
-------
F NAPL PROJECT 198
F.I Purpose 198
F.2 Scope 198
F.3 Experimental Setup and Data Base 199
F.4 TASK 1 202
F.4.1 Purpose 202
F.4.2 Procedure 204
F.4.3 Results 205
F.5 TASK 2 205
F.5.1 Purpose 205
F.5.2 Procedure 205
F.5.3 Results 207
F.6 TASK 3 208
F.6.1 Purpose 208
F.6.2 Procedure 209
F.6.3 Results 213
F.7 TASK 4 215
F.7.1 Purpose 215
F.7.2 Procedure 215
F.7.3 Results 216
-------
LIST OF FIGURES
1.1 Definition illustration of NAPL contamination in near-surface soils
due to an intermittant surface release 12
4.1 Summary of recent studies of mass exchange processes in NAPL
contaminated soils 32
4.2 Continued: Summary of recent studies of mass exchange processes
in NAPL contaminated soils 33
4.3 Continued: Summary of recent studies of mass exchange processes
in NAPL contaminated soils 34
5.1 Summary of conceptual and numerical models describing three-
phase flow in granular soils 42
5.2 Continued: Summary of conceptual and numerical models describ-
ing three-phase flow in granular soils 43
5.3 Definition plot of the hysteretic relationship between saturation
and capillary pressure employing the empirical model used in the
simulator. Curve position and shape is governed by the mobility
status and the magnitude of the phase saturations when the curve
is spawned and whether the displacement process is drainage (D)
or imbibition (/) with respect to the wetting phase. Primary (P)
and Main (M) curves are spawned when only one phase is mobile
[curve numbers 1 and 2]. Scanning (S) curves are spawned when
both phases are mobile [curve numbers 3 to 6]. The curve-type
numbering scheme is set such that odd numbers are aligned with
drainage curves [1 = PDC or MDC, 3 and 5 = SDC] and even
numbers are aligned with imbibition curves [2 = PIC or MIC, 4
and 6 = SIC] 47
-------
5.4 The resulting hysteretic k — S functionals for the wetting phase
generated from the empirical model used in the simulator and the
data denning the S — P relationship. Note that the MDC and
MIC', shown as dashed lines, are practically coincident, and that
the scanning curves, shown as dashed lines, are group-labled be-
cause by model definition they are coincident 48
5.5 The resulting hysteretic k — S functionals for the nonwetting phase
using the data defining the S — P relationship and the empirical
model used in the simulator. Note that the curve labeled MIC*
is obtained upon reversal from a PDC where no nonwetting phase
was previously trapped. Subsequent reversals follow the MDC
and MIC which are practically coincident. Also note that the
scanning curves are group-labled because by model definition they
are coincident 49
5.6 An illustration of the effect that the entrapment model blending
parameter e has on the shape of the k — S — P functionals. For the
solid curves, e = 1, and for the dashed curves e = 0 72
5.7 The effect of using the blending rule described in Table 5.2 to define
-------
8.1 Analysis of appropriate grid spacing to compute capillary rise for
different soil-types. Parts (a) and (b) are for a relatively fine sand,
and parts (c) and (d) are for a relatively coarse sand 148
8.2 Results of a one-dimensional, three-phase, DNAPL injection and
redistribution simulation, highlighting spatial convergence and mass
balance 151
8.3 Computational analysis of the dissolution model. Parts (a) and
(b) illustrate the effect that the rate constant (ex in the figure, ex-
pressed in units of I/day) has on the solution. As the dissolution
front sharpens, oscillations appear indicating that a finer grid spac-
ing is required. Parts (c) and (d) illustrate spatial convergence for
ex = 24/d. For the parameters chosen a grid spacing of approxi-
mately 5 cm is appropriate 152
8.4 Plot of the primary and main S-P functionals defined by the current
model for the LNAPL spill simulation, where the drainage curves
are represented by the thick lines and the imbibition curves are
represented by the thin lines. Here the fitting parameters, assumed
to be valid for a water-gas system, have been scaled to represent
the water-NAPL and NAPL-gas systems 154
8.5 Definition sketch for the LNAPL spill simulation, showing spatial
scale and boundary conditions 154
8.6 Comparison of results from the physical experiment, the current
model, and the model used by Van Geel and Sykes (VGS). Part (a)
shows the vertical distribution of water pressure head and Part (b)
shows the vertical distribution of water saturation 155
8.7 Comparison of results for the LNAPL spill problem. The plots
on the left show the NAPL saturation contours as computed by
the current model at the times indicated from the initiation of the
LNAPL spill. The plots on the right compare results taken along
the instrumented vertical section (the vertical dotted line in the
plots on the left) 156
8.8 A comparisson of the cumulative LNAPL mass which has entered
the domain as a function of time. At time = 1120 s the LNAPL
source was removed. The solid line is the computed cumulative
mass which has crossed the boundary. The dashed line is the change
in mass in the domain. The dash-dot line is the experimental data. 157
-------
8.9 A comparison plot of the experimental and model moisture profile
used as the initial condition for the DNAPL flood 160
8.10 A comparison plot of the experimental and model data quantifying
the cumulative volume of PCE infiltrated as a function of time.
Specific to the experiment, it took 143 seconds for 200 cm3 of PCE
to infiltrate 161
8.11 Comparison between experiment and model results at time = 143
seconds (the time when the DNAPL source was removed) 162
8.12 Comparison between experiment and model results at time = 283
seconds after the DNAPL source was first applied 163
8.13 Comparison between experiment and model results at time = 1195
seconds after the DNAPL source was first applied 164
8.14 Comparison between experiment and model results at time = 3595
seconds after the DNAPL source was first applied 165
8.15 An illustration of the DNAPL vapor transport model problem do-
main including initial and boundary conditions 167
8.16 A comparison plot of the experimental (solid line) and model (dashed
line) results at time = 12 hours 168
F.I Photo of the experimental apparatus just after the PCE was re-
moved, showing dimensions and vertical constant head boundaries.
The PCE (dark grey) has been dyed red to maximize contrast. . . 202
F.2 Idealization of the experimental setup superimposed on a video
image of the box (5 seconds after the PCE source was applied). . 203
F.3 Plot of the initial static saturation profile. The computed curve is
fit to the experimental data by altering the S-P model curve fitting
parameters 204
F.4 Illustration of the one-dimensional water-gas displacement experi-
ment showing boundary and initial conditions and dimensions. The
mesh has one 2 cm element in the horizontal direction, and the ver-
tical direction is discretized in an appropritate manner 206
F.5 Illustration of the 2-D model setup for the three sequential simula-
tions, showing boundary and initial conditions and the time frame
for each experiment 210
F.6 Illustration of the uniform distribution of forcing conditions about
a node. For this example, NAPL is ponded with head hN over an
area equal to 2.5 cm2 212
-------
F.7 Video image at time = 143 seconds just after the DNAPL source
was removed (i. e., at the end of Experiment 2), showing dimen-
sions, where the superimposed grid is for reference purposes (ele-
ments are 5.5cm by 5.5 cm) 218
F.8 Video images at time = 283 and 683 seconds after the DNAPL spill
began (i. e., take Time = 0 as the initial condition for Experiment
2) 219
F.9 Video images at Time = 1195 and 1795 seconds after the DNAPL
spill began 220
F.10 Video image at Time = 3595 seconds after the DNAPL spill began. 221
-------
LIST OF TABLES
5.1 TWO-PHASE k-S-P CURVE-TYPE DEFINITION 50
5.2 DEFINITION OF THE SCALING PARAMETER 'a' FOR THE
TWO-PHASE HYSTERETIC S-P MODEL 57
5.3 DEFINITION OF THE SCALING PARAMETER 'Sr' FOR THE
TWO-PHASE HYSTERETIC S-P MODEL 59
5.4 DEFINITION OF THE SCALING PARAMETER 'Ss' FOR THE
TWO-PHASE HYSTERETIC S-P MODEL 60
5.5 THE RELATION BETWEEN TWO- AND THREE-PHASE S-P
MODEL PARAMETERS 67
6.1 SUMMARY OF THE SYSTEM USED TO DEFINE BOUNDARY
CONDITIONS FOR THE FLOW VARIABLES 95
7.1 TIME INDEPENDENT BUT SPATIALLY VARYING INPUT PA-
RAMETERS AND THEIT ASSOCIATED INPUT FILES .... 110
8.1 PARAMETERS USED TO MODEL THREE-PHASE PCE MI-
GRATION IN OTTAWA SAND 149
8.2 PARAMETERS USED IN THE LNAPL SPILL PROBLEM ... 153
8.3 PARAMETERS USED IN THE DNAPL SPILL PROBLEM ... 158
8.4 PARAMETERS USED IN THE DNAPL VAPOR TRANSPORT
EXPERIMENT (note, parameters with an asterisk are estimated) 168
F. 1 EXPERIMENTAL DATA - MOISTURE CONTENT AS A FUNC-
TION OF DEPTH 201
F.2 EXPERIMENTAL DATA - PCE VOLUME ENTERING THE DO-
MAIN AS A FUNCTION OF TIME 201
10
-------
1. INTRODUCTION
1.1. PROBLEM
The physical problem which is addressed herein is the contamination of a pristine
porous medium as the result of releases of organic liquids, commonly referred to
as Non-Aqueous Phase Liquids (NAPLs), in near-surface heterogeneous granular
soils. The organic liquids can be either lighter than water (identified as LNAPLs,
i. e., petroleum hydrocarbon-based) or heavier than water (identified as DNAPLs,
i. e., chlorinated hydrocarbon-based). By near-surface we mean that the scale
which characterizes fluid pressure is on the order of atmospheric pressure. In
addition the soil domain can be idealized as consisting of three interrelated zones:
a vadose zone which is in contact with the atmosphere, a capillary fringe zone,
and a water-table aquifer zone. A particular problem of interest may include all
three zones or a subset thereof. By granular soils we mean those soils which are
stable (non-deforming) and relatively chemically inert (the soil particles do not
interact with the soil fluids). Therefore, the soil is idealized as containing a high
percentage of quartz particles and only a minor percentage of clay particles and
organic matter.
A conceptual illustration of surface-release-generated NAPL migration in the
vadose, capillary fringe and aquifer zones is provided in Figure 1.1. There are
three fundamental mechanisms for NAPL migration. First, the NAPL infiltrates
into the soil and migrates both vertically and laterally under the influence of grav-
itational and capillary forces. The distribution of the NAPL liquid is a function
of fluid properties (density, viscosity, interfacial tension, wetting potential and
variable chemical composition), soil properties (grain size distribution, mineral
content, moisture content, porosity, hydraulic conductivity and spatial hetero-
geneity), and system forcing history. If the source is periodic in nature, then
during drying periods, not all the NAPL will drain from the pore space, leaving
behind an immobile residual, held in place by capillary forces. If the NAPL is
more dense than water, it will migrate through the capillary fringe and continue
11
-------
NAPL source area
1
ts
soil
moisture
profile
^
J
vado
i
se zone
_i
|
capillary fringe zone 1
water content
aquifer zone
t
residual NAPL
dissolved NAPL
vaporized NAPL
mobile NAPL
_water_table
roundwater flow
Figure 1.1: Definition illustration of NAPL contamination in near-surface soils
due to an intermittant surface release.
its vertical migration until either the mobility becomes zero (all the NAPL liquid
is at the immobile residual state) or the NAPL front encounters an impenetrable
geologic horizon
The second contaminant transport mechanism is dissolution and consequent
advection in the downward-flowing water-phase, with precipitation providing the
water source in the vadose zone. In the case of a DNAPL, flowing groundwater
picks up dissolved NAPL constituents.
The third transport mechanism is transport as a vapor NAPL constituent in
the soil gas, where the increased gas-phase density induces downward movement.
Partitioning between the gas- and water-phase contaminants further enhances the
migratory potential of the NAPL constituents.
1.2. FOCUS
The focus of this investigation is to develop a physically complete subsurface flow
and transport mathematical model (also referred to herein as the simulator) to
study the movement and fate of NAPL contaminants in near-surface granular
soils. Specifically, three fundamental, interrelated, physical processes have been
identified: multiphase flow, interphase mass transfer, and constituent mass trans-
port. The multiphase flow process defines the time-dependent volumetric extent
of the mobile and immobile components of the water, NAPL and gas phases. The
12
-------
interphase mass-transfer process defines how the NAPL contaminants partition
between phases. The constituent mass transport defines the temporal and spatial
distribution of the NAPL contaminants within a given phase. The three processes
are related in a nonlinear way since the phase velocity and mass exchange terms
are both functions of, among other things, phase volumetric content, and the
phase fluid properties (density, viscosity, and interfacial tension) are in general a
function of the chemical composition of the phase.
The model presented herein is developed with a focus on better quantifying
these three fundamental physical processes as they exist in the natural environ-
ment. Particular attention has been paid to quantifying the following processes:
1. fluid entrapment and release;
2. hysteresis in the relative permeability-saturation-capillary pressure model;
3. rate-limited mass transfer to describe NAPL dissolution and volatilization;
4. advective-dispersive transport in both the water and gas phases.
The major purpose of developing such a simulator is that, once compiled, it can
be used to verify the theory describing the physics of the problem and to quantify
parameter sensitivity. In addition, the information derived from verification and
sensitivity analysis can be used to simplify the system (for example, derive and
evaluate simpler constitutive models) which is an important consideration given
the computationally demanding nature of the problem (the solution of coupled
nonlinear equations with sharp-front transport characteristics). Through this type
of development, the resulting simulator, coupled with dedicated pre- and post-
processing software for data input and output visualization, respectively, will be
an efficient and effective engineering tool to be used for field-scale analysis.
1.3. APPROACH
To describe the physical problem mathematically, a set of coupled nonlinear par-
tial differential balance equations (PDE's) which govern the temporal and spatial
variability of the system, and a set of constitutive and thermodynamic equations
which relate physically-based parameters occurring in the PDE's (for example
relative permeability) to the dependent variables (for example fluid pressure and
saturation) must be defined. There are two constitutive models which are of
particular importance with regard to this simulator:
13
-------
1. a model of three-phase relative permeability-saturation-capillary pressure
relationships which includes flow- path-history-dependent functionals (hys-
teresis), fluid entrapment considerations, and functional dependence on fluid
and soil properties;
2. a model of rate-limited interphase mass transfer processes, including NAPL
dissolution and volatilization.
The resulting mathematical interpretation of the physical system is solved
using a numerical solution algorithm which employs the following conceptual tools:
• An implicit-in-time collocation finite element method with Hermite cubic
basis functions is used to generate the systems of algebraic equations.
• A successive substitution iteration scheme is used for nonlinear terms.
• A sequential solution procedure is used to solve the coupled balance equa-
tions to minimize the system matrix order and bandwidth, where the phase
flow equations are solved sequentially using a total velocity formulation, and
given the flow solution, the transport equations are then solved.
1.4. MODEL CAPABILITIES
The mathematical model which is presented in this documentation can be char-
acterized by a list of attributes. This list is intended to provide the reader with a
summary of the capabilities and limitations of the simulator.
1. Problems in one, two and three spatial dimensions (Cartesian coordinate
system) are applicable.
2. The finite element mesh utilizes rectangular elements.
3. The simulator can accommodate as many as three-fluid phases, identified as
a water-phase, a NAPL-phase and a gas-phase, and can model flow of either
one, two or three phases in any combination.
4. Darcy's law is a valid model for quantifying water-, NAPL-, and gas-phase
advection;
14
-------
5. Pick's law is a valid model for quantifying water- and gas-phase diffusion
processes.
6. The water phase is characterized as follows:
• it is incompressible;
• it is made up of two species: water and dissolved NAPL;
• its properties are a function of contaminant concentration only;
• the dissolved NAPL can vaporize into the gas phase and adsorb onto
the solid phase.
7. The NAPL-phase is characterized as follows:
• it is incompressible;
• it is made up of a single chemical species;
• it is able to dissolve into the water phase;
• it is able to volatilize into the gas phase.
8. The gas phase is characterized as follows:
• it is incompressible1;
• it is made up of two species: gas and volatilized NAPL;
• its properties are a function of contaminant concentration only;
9. The porous medium is characterized as follows:
• it is non-deforming;
• it is generally heterogeneous and isotropic;
• it is made up of two species: soil and adsorbed NAPL.
10. Isothermal conditions prevail.
lrrhe model is intended to be applied to physical problems where pressure gradients are
small, and where the time-scale defining gas-phase pressure response to forcing is much shorter
than that which defines gas-phase saturation response. Therefore the gas-phase is considered
effectively incompressible and Darcy's law is assumed to apply.
15
-------
11. Mass transfer relationships are denned by:
• NAPL phase dissolution is governed by a rate-limited, first-order, ki-
netic rule;
• NAPL phase vaporization is governed by a rate-limited, first-order,
kinetic rule;
• dissolved NAPL species vaporization is governed by a rate limited,
first-order, kinetic rule;
• dissolved NAPL species adsorbtion is governed by a linear equilibrium
partitioning rule.
12. The rules and relationships which define the multiphase flow parameters are
as follows:
• the porous medium is isotropic, and it may be saturated with one, two
or three phases;
• when more than one phase is present the relationship between relative
permeability, saturation and capillary pressure (called the k-S-P model)
is based on wettability considerations and two-phase data, where phase
wettability is constrained to follow from most to least: water-NAPL-
gas.
• the k-S-P model is subject to hysteresis due to capillary and fluid en-
trapment effects;
• the van Genuchten (1980) saturation-pressure model is modified to
accommodate hysteresis;
• the Mualem (1976) relative permeability-saturation model is modified
to accommodate hysteresis;
• the capillary pressure between phases can be scaled to accommodate
variable fluid and soil properties.
13. Boundary and external forcing conditions are summarized as follows:
• Dirichlet data for the multiphase flow problem is specified as either
one of the three phase pressures known or all of the primary flow vari-
ables known (i. e., both pressure and saturation). Dirichlet data for
contaminant transport is specified as a known value for the species of
interest.
16
-------
• Non-zero flux conditions for a phase or species within a phase are ac-
commodated by specifying an appropriate point source or sink term (i.
e., a well).
• Diffusive mass flux of the NAPL species in the gas-phase through a
boundary layer at the interface between the ground surface and the
atmosphere is accounted for by the use of a mixed boundary condition.2
14. The numerical model includes the following features:
• The code is written in standard FORTRAN?? with the intent of making
it portable.
• The routines which can run in parallel mode have been coded to do so
using the Silicon Graphics f?7 compiler.
• The code has a re-start facility.
• The code is memory intensive, especially in 3-D mode where its utility
for solving 3-D problems is limited. Therefore, a 2-D version is also
included as part of this simulation package.
• Standard upstream weighting as a means of adding artificial diffusion
to the solution of parabolic equations is not utilized in this simulator.
Instead, a physically-based diffusion term is added to the solution in a
point-wise fashion when necessary, the magnitude of which is based on
a user-defined critical Peclet number.
• Time-step control is provided by two algorithms. One is based on the
number of iterations required for convergence on the nonlinearity, and
one is based on a user-defined maximum local Courant number.
• The simulator has been integrated with a commercially available graph-
ical user interface (GUI), and any interested users can contact the au-
thors for more information.
1.5. ORGANIZATION
This report is organized as follows. The conclusions and recommendations are
provided in Sections 2 and 3, respectively. Section 4 describes the mathematical
2 The two-dimensional version of the simulator includes the mixed-type boundary condition
for gas-phase transport, but the three-dimensional version does not.
17
-------
model which includes development of the governing and constitutive equations.
Section 5 provides a detailed description of the hysteretic k — S — P model used
in the simulator. Section 6 presents the numerical approach used to make the
mathematical model amenable to computational solution. The development in
this Section includes equation discretization, linearization and iterative solution.
Section 7 details the data input and output structure of the simulator called
NAPL. Section 8 details the analysis used to verify the mathematical and numer-
ical models, including convergence and comparison to experimental results. This
Section also details a set of example problems which are tutorial in nature.
18
-------
2. CONCLUSIONS
The movement of NAPL in a porous medium environment can be described by
a coupled set of non-linear partial differential equations that are generated by
combining fundamental balance equations with constitutive and thermodynamic
relationships. Three-phase coupled flow, including hysteresis and mass transfer
effects, is considered. Examination of the literature suggests that the model pro-
vided above is the most physically comprehensive mathematical representation for
the problem of defining the emplacement of NAPL residual. The ability to repro-
duce a physically complex NAPL experiment without calibration demonstrates
that this model is a suitable real-world surrogate.
The nature of the model equations requires the use of a numerical technique
to obtain a solution to physically meaningful problems. The method selected
in this work, collocation finite elements, provides a very accurate solution for the
examples tested, including one field application. The iterative procedure employed
allows for the solution of three-dimensional problems on conventional computer
systems.
19
-------
3. RECOMMENDATIONS
While the work described herein provides an accurate, physically sophisticated
simulator for NAPL , additional work on this software is recommended:
1. The solution algorithm, while adequate for small problems, is inappropriate
for large, three-dimensional field applications. A faster, less computer -
memory intensive solution technique is needed.
2. The new permeability-saturation-pressure relationship should be examined
experimentally.
3. The current model allows for a single species NAPL. A multi-component
NAPL capability would be an important extension.
4. More model verifications using carefully conducted experiments would be
helpful. Comparisons between various currently existing models would be
enlightening.
5. The model should be extended to accommodate various remediation strate-
gies. Consider for example:
(a) the use of surfactants and co-solvents for NAPL mobilization and re-
moval;
(b) the accommodation in the model of biochemical transformations;
(c) the inclusion of the NAPL model into optimal design software;
(d) the effect of thermal forcing (i. e., steam stripping).
20
-------
4. THEORETICAL DEVELOPMENT
4.1. OVERVIEW
In this section a mathematical model based upon the attributes presented in the
previous section will be developed. The mathematical model is made up of two
major components: the mass balance equations which define the distribution in
time and space of the primary variables, and the constitutive equations which
define the inter-relationship between primary and secondary variables1.
The balance equations describe the conservation of mass of each phase (fluid
and solid) and each constituent within a particular phase (species) as they move
and intermingle within a porous medium. One mass balance equation can be
generated for each phase and constituent of interest. Each equation is composed
of terms which define the various components of mass transport at the macroscopic
scale: accumulation, advection, dispersion, external sources and sinks, and mass
transfer between phases. Each term in turn is defined by a set of parameters which
quantifies the physics of the transport process for a particular physical system.
In order to solve the resulting balance equations they must be augmented by a
set of constitutive relationships which relate the primary and secondary variables.
In the discussion which follows, the constitutive relationships are separated into
four categories:
1. those which define the fluid properties density, compressibility and viscosity]
2. those which define fluid flow or advection]
3. those which define non-advective species transport, namely dispersion, and
diffusion;
4. those which define interphase mass exchange.
lrrhe primary variables are those which are advanced in time and space by solving the
mass balance equations. The secondary variables are those variables which are functions of the
primary variables.
21
-------
As will be seen, while the constitutive relationships can be categorized as indi-
cated, they exhibit inter-functional dependence through phase pressure, saturation
and compositional dependence. Therefore, the equation set that is generated by
combining the balance equations and the constitutive relationships is coupled and
nonlinear.
Another outcome of this development is a list of primary variables which are
used to solve the balance equations, and an outline of the simplifying assumptions
employed thereby providing a basis for model applicability. Given proper initial
and boundary conditions, the solution of the equation set will yield the phase
volume and pressure distributions and the constituent concentrations in time and
in space.
In the following discussion, the measurement unit scale representing a par-
ticular measurable quantity is written in brackets. For example, the units for
fluid viscosity are [M/(TL)j, where M represents mass, T represents time and L
represents length.
4.2. MASS BALANCE EQUATIONS
As a point of departure, let us employ the equation development of Finder and
Abriola (1986). First, consider the mass balance law for each fluid constituent, an
ordered pair (L, a) representing a species i in a fluid phase a:
-eSaK"p" = p"Qa+'f\ (4.1)
dt
where the five constituents (L, a) relevant to this simulator are identified as:
(w,W),a water species in the water phase; (n, W), a NAPL species in the wa-
ter phase; (n, TV), a NAPL species in the NAPL phase; (n, G), a NAPL species
in the gas phase; and (g,G), a gas species in the gas phase2. Other symbols
occurring in equation 4.1 are used to represent the following:
e is the porosity of the porous medium.
Sa is the saturation of the a-phase.
The convention used to identify phase and species is that the phase is represented by upper-
case letters [i. e., W (water), N (NAPL), G (gas), and S (solid)], and the components of the
phase (species) are represented by lower-case letters [i. e., w (water), n (NAPL), g (gas)].
22
-------
p" is the mass concentration of species i in the a-phase [M/L3].
va is the mass average velocity of phase a, a vector [L/T].
Da is the dispersion coefficient for the a-phase, a symmetric second-order tensor
Qa is the point source (+) or sink (-) a-phase mass [1/T].
K" is the decay coefficient for species i in the a-phase [1/T].
pf is the source or sink of mass for a species i in the a-phase [M/L3T] due to
interphase mass exchange (i. e., dissolution, volatilization and adsorption).
The exchange of mass for each constituent in equation 4.1 is denned by:
P» = 0 (4.2)
-~-W r~,W pG T?S
Pn = l^n ~ tin/W ~ En/w
^N ( WW , j?G\
Pn = ~ (hn + tin )
^ pG , pG
Pn = tin + &n/w
^G n
Pg = °
where
E^ represents dissolution mass transfer of the NAPL species from the NAPL
phase to the water phase;
^n/w represents volatilization mass transfer of the NAPL species from the water
phase to the gas phase;
E!% represents volatilization mass transfer of the NAPL species from the NAPL
phase to the gas phase;
En,w represents adsorption mass transfer of the NAPL species from the water
phase to the soil.
A sixth mass balance equation is required to describe the NAPL species mass
which is adsorbed onto the soil. This equation is written as:
9 ([1 ~ e] A
-
23
-------
where ps is the density of the soil [M/L3] and tuf is the mass fraction of the
adsorbed NAPL on the solid [dimensionless] . As described in Section 4.8.2, the
balance equation 4.3 is replaced by the following linear equilibrium relationship:
u,f = K^ (4.4)
where K& is a distribution coefficient [L3/M].
To ensure global mass conservation, the following definitions and constraints
on fluid volume, density and mass exchange are employed:
1. The a-phase saturations must sum to one:
Sw + SN + SG = l (4.5)
2. The a-phase mass density, pa [M/L3], is the sum of the species mass con-
centrations in the a-phase:
Pa= . P? ,<* = W,N,G (4.6)
L=w,n,g
3. The sum of mass fluxes of all species L into the a-phase, must equal the total
mass change in the a-phase:
ff,<* = W,N,G (4.7)
4. The total mass change over all phases must be zero:
V V ^Q,
/ p = U
a=W,N,G
5. The sum of the reacting mass must be equal to the sum of the produced
mass:
y" K°p° = 0 , a = W, N, G
/ > i, / i,
A set of fluid phase mass balance equations can be generated by summing
the balance equations 4.1 for each species within the phase, and by incorporating
equations 4.2, 4.6 and 4.7. The three resulting fluid-phase mass balance equations
are:
24
-------
Water-phase:
Q / C
U \t-°WfJ
— r Q W'rWl nWnW , i?W PG nS (A o\
- — -- h V • [ebwp v J = p Q + hn - hn/w - hn/w (4.8)
NAPL-phase:
9 £ P JV * "" (4-9)
Gas-phase:
3
PGQG + EG + EGW (4.10)
With this development, the physical problem can be cast into a mathematical
representation consisting of five mass balance equations. Of the balance equations
written, the following five are used in the simulator:
1 to 3) The three fluid-phase balance equations, equations 4.8, 4.9 and 4.10.
These equations define the temporal and spatial distribution and the flow
properties of the water-, NAPL- and gas-phases throughout the domain.
4 and 5) The two NAPL species balance equations, equations 4.1 with (L, a) =
(n, W) and (n, G). These equations define the temporal and spatial distrib-
ution of the NAPL species as they are transported within and between their
respective phases.
4.3. PRIMARY VARIABLES
Five primary (or dependent) variables are required to solve the balance equations
listed in the previous Section. The five primary variables used in the simulator
are:
{Pw,Sw,STw,p™,p°} (4.11)
where Pw/ is the pressure in the water phase [M/(LT2)], and from equation 4.5
we have defined a new saturation measure,
ST-W = Sw + 5jv = 1 —
-------
called the total liquid phase saturation. These variables in 4.11 are denned such
that they are continuous in time and space. In addition, they are intended to
apply regardless of which phase configuration exists (i. e. one-, two-, or three-
phase flow). The following constraints on the relationship between the pressure
and saturation variables allow for this attribute to be implemented:
1. The property of fluid wettability is defined as, from most to least, water-
NAPL-gas.
2. The property of capillary pressure between immiscible phases is defined as
a function of phase saturation, called herein a saturation-pressure model
(S — P model), where Pco,p = Pa — Pl3 is the capillary pressure between
the nonwetting a-phase and the wetting /3-phase3. In addition when all
three phases are present, the NAPL-phase renders the water-gas interactions
negligible, and three-phase behavior can be gleaned from two, two-phase
S - P models4:
W
PCNW = PCNw(Sw) = PN ~ PW (4-12)
CGN =
PCGN = PCGNSTW = PG - PN
3. The functions PCNW(SW} and PCGN(STW) are defined such that:
(a) they are continuous in time and space regardless of phase configuration;
(b) for the two-phase water-gas case, PCQW is determined from:
PCGW = PCNW + PCGN (4-13)
(c) PcGWi PCNW and PCGN are related through the following scaling rule:
PcNW = ^£E = PcGW (4.14)
&NW &GN &GW
where C^GW> &NW and O"GJV are the interfacial tensions along the inter-
faces between the gas and water phases, the NAPL and water phases,
3 Because the definition of the S-P model is an important component of this simulator, its
development is detailed in Chapter 5.
4 For the case when NAPL is the intermediate wetting fluid, STW is called the total wetting
phase saturation. By convention the S — P relationship is written in terms of the wetting phase
saturation.
26
-------
and the gas and NAPL phases, respectively. In order for equations 4.13
and 4.14 to be compatible, the interfacial tensions are constrained to
be related by:
°GN = 0
An important outcome of applying the capillary pressure scaling rule (equa-
tion 4.14) and the constraint of a neutral spreading coefficient (equation
4.15) is that, of the three capillary pressure relationships required to model
two- and three-phase flow, only one of the three relationships needs to be
measured, with the other two gleaned from the use of equation 4-14-
4. The NAPL and gas-phase pressures are nonlinear functions of Pw , Sw and
STW as indicated from the following definitions:
PN = Pw + PcNw(Sw) (4-16)
pG = Pw + PcNW(sw) + PcGN(STw)
The constitutive models which define the relationships between the primary
and secondary variables are derived with reference to the primary variables in
4.11 and the constraints imposed on the relation between fluid-phase pressure
and saturation.
4.4. FLUID PROPERTIES
The a-phase fluid properties of interest are the density, pa[M/L3], viscosity,
yUa[M/(TL)], and the interfacial tension between immiscible phases a and (3, cra/3
[M/T2]. Under the assumption of isothermal conditions, these parameters are,
in general, a function of the chemical makeup of the phase (i.e. phase composi-
tion) and the applied pressure (to the degree that it effects the phase properties
via phase compressibility). Assuming that natural or induced pressure varia-
tions characteristic of applications involving near-surface soils in contact with the
atmosphere are small (Sleep and Sykes, 1989, Mendoza and Frind, 1990, and
Brusseau, 1991), the dependance of fluid properties on pressure is neglected in this
simulator.
4.4.1. Density
The dependance of fluid density on composition is modeled as follows:
27
-------
pN = pNr = constant
pW =
1-
P
Wr
1-
P
P"
Gr
pNr
(4.17)
(4.18)
where pw/r, pNr and pGr are the mass densities of pure-phase water, NAPL and
gas, respectively.
4.4.2. Viscosity
Fluid viscosities are modeled as a function of composition as follows:
= constant
log/3 =
JVr
P
Nr
(4.19)
(4.20)
where yUW/r, yuArrand n
and gas, respectively.
are the pure liquid phase viscosities of the water, NAPL
4.4.3. Interfacial Tension
The interfacial tensions between fluid phases are assumed to be known constants.
The combination of interfacial tension and contact angle defines fluid wettability.
Regardless of the values of the interfacial tensions, the contact angle is assumed
to be such that the fluid wetting order is constant, and follows, from most to least,
water-NAPL-gas.
The interfacial tensions are used to scale the capillary pressures as indicated in
equation 4.14. In addition, in order to model the general case where the number of
phases saturating the pore space can vary over time between one, two and three,
the simulator employes the algebraic constraints 4.13 and 4.15. Therefore, only
two of the interfacial tensions are independent with the third defined by equation
4-15. The two independent values can be considered fitting parameters.
28
-------
4.5. PHASE ADVECTION
Fluid flow is defined by the parameter va in equations 4.1 and 4.8 through 4.10.
The phase velocity is written in terms of the multiphase extension of Darcy's law:
- TQV^) ,a = W,N,G (4.21)
where Pa is the a— phase pressure [M/(LT2)], 7" = pag is the specific weight of
the phase [M/T2], (7 is the acceleration due to gravity [L/T2], k is the intrinsic
permeability [L2] , considered a scalar herein, and kra is the relative permeability.
The a— phase relative permeability is a scaling factor, 0 < kra < 1, which
accounts for the case where the porous medium is not fully saturated with the
a— phase. This parameter is in general a function of the a— phase saturation.
Given the assumption that the phase wetting-order is constant, and follows, from
most to least, water-NAPL-gas, the following functional dependance is assumed:
krw = krW(Sw) (4.22)
where the relationships listed reduce to their proper two-phase forms when appro-
priate. Since the definition of the relationship between the relative permeability
and saturation (called herein the k — S model) is a major component of this simu-
lator, it is developed from both a conceptual and empirical viewpoint in Chapter
5.
Given that Pw , Sw and STW a-re the primary flow variables, and that the phase
pressures are related to one another through capillary pressure relationships, 4.12
and 4.13, the form of equation 4.21 for each fluid phase of interest is detailed here
as:
(4.23)
P"+PrW ' " '
29
-------
4.6. DISPERSION COEFFICIENT
The dispersive flux of the NAPL species in both the water and gas phases is
defined by the third term in equation 4.1. For the case where phase density is
related to phase composition as defined in equations 4.17 and 4.18 the dispersion
term in equation 4.1 is simplified as:
a S0ar\
pQDQ • V(— ) « I — 1 DQ • Vp" (4.24)
1 pa \paj
This relationship shows that if pa is a constant equal to par, then the standard
definition for dispersive flux applies, but if pa is a function of phase composition,
then dispersive flux becomes a nonlinear function of concentration.
The dispersion coefficient, Da, is a second-rank tensor, and its components in
(x, y, z] Cartesian coordinates are represented as (Scheidegger, 1961):
Daxx = a? K| + K - 4) «)2 / va + TaD*m (4.25)
a -
- a
where a = W and G, a£ is the longitudinal dispersivity [L], af, is the transverse
dispersivity [L], t>", v£ and v" are the components of the interstitial pore water
velocity vector (equation 4.21), va is the mean phase velocity magnitude, D^ is
the a-phase coefficient of molecular diffusion [L2/T], and ra accounts for diffusion
porosity and tortuosity effects [dimensionless] . The approximation for ra used in
the simulator is provided by the Millington and Quirk (1961) model:
= e1/3 (Sa
/3
4.7. NATURAL DEGRADATION
Many NAPL's are subject to biological and chemical degradation in the soil,
and therefore they are not stable in the natural environment. To account for
this natural attenuation phenomenon a constant first-order degradation rate is
30
-------
assumed for account for all degradation processes. With this, the decay coefficients
in equation 4.1, K^ and K^ [1/T], are defined in terms of the half-life, £1/2 [T], of
the NAPL species as:
«f = «£ = In (2) /tl/2 (4.26)
where the assumption is made that £1/2 is the same for the NAPL species in the
gas and water phases.
4.8. MASS TRANSFER
To obtain a perspective on current research applications in mass transfer processes
affecting the fate of NAPLs in near surface granular soils, consider Figures 4.1,
4.2 and 4.3 which summarize relevant publications. Upon review of these figures
the following conclusions can be drawn:
1. Four types of mass transfer processes are important in defining the physics
of the fate of NAPLs in near surface granular soils:
(a) dissolution mass transfer of pure phase NAPL to the water phase;
(b) evaporation mass transfer of pure phase NAPL to the gas phase;
(c) evaporation mass transfer of NAPL species in the water phase to the
gas phase;
(d) adsorption mass transfer of NAPL species in the water phase to the soil
phase.
Adsorption of NAPL species in the gas phase directly to the soil phase is
neglected.
2. Three basic mass transfer models are typically utilized :
(a) Local equilibrium model: (indicated by 'E' in the review figures) Based
on the assumption that equilibrium partitioning dictates NAPL phase
concentrations. This model-type has been shown to have utility in
problems involving low Darcy velocities and homogeneous soil proper-
ties.
31
-------
JD
CL h^
CO
O CO
GO
to
CD
O
CD
Bt
CD'
CO
CO
CD
X
O
O
O
CD
CO
CO
CD
CO
reference
Armstrong et
al. (1994)
Mendoza and
Frind (1990)
Sleep and
Sykes
(1989)
Falta et al.
(1989)
phases and
components
considered
Sw>0, krw=0
SG>0, krG>0
N 0
(MX
Sw>0, krw=0
SG>0, krG>0
SN=0
P° (^
Sw>0, krw>0
SG>0, krG>0
SN>0, krN=0
G W
Dn Dn
Sw>0, krw=0
SG>0, krG>0
SN>0, krN=0
G S
Pn,»n
exchange processes
modeled
AG Henry's law -
Pn/w E and K
AW T-, , v
pn E and K
A s organic carbon
Pn/w partitioning -
E andK
AG Henry's law -
Pn/W E
A s organic carbon
Pn/w partitioning -
AG Henry's law -
Pn/W K
AW „
Pn K
£G Henry's law -
Pn/W E y
AS organic carbon
Pn/w partitioning -
kinetic
process
diffusive
boundary
N/A
diffusive
boundary
layer
N/A
parameter
variability
constant rate coefficients
variable temperature and
soil properties
PG = constant
Sw = constant
p° = f(composition)
isothermal
PG = constant
Sw = constant
pG = f(composition)
|IG = f(composition)
looked at layered soils
constant rate coefficients
uncoupled water-gas flow
PG = constant
S» = constant
pw = constant
p° = f(composition)
get saturated vapor conc-
entration from ideal
gas law
homogeneous soil
S\y= constant
SN = constant
|IG = constant
pG = f(composition)
investigation
type
numerical
model:
advective/
dispersive
gas-phase
transport
numerical
model:
advective/
dispersive
gas-phase
transport
numerical
model:
advective/
dispersive
gas- and
water-phase
transport
numerical
model:
advective/
dispersive
gas-phase
transport
• focus on soil gas venting application
• rate-limited gas-water mass transfer
is the controlling process for the
time-scale chosen
• first-order kinetic is adequate
• cannot calibrate the model with
effluent concentrations
• density-driven gas-phase transport is
important when NAPL source is
present and in high-K soils
• three important transport processes:
diffusion, density-driven advection,
vapor source mass flux
• density-driven gas-phase transport is
important
• neglected the effects of capillarity
• the magnitude of density-driven
gas-phase flow velocity varies
linearly with K and density contrast
between pure and contaminated gas
• need to look at: soil heterogeneity,
multi-component organic liquids,
the effect of gas-phase pressure
forcing
O
O
-------
GO
GO
CD
O
CD
CO
c-t-
CO
O
CO
CD
X
o
CD
o
o
CD
CO
CO
CD
CO
reference
Brusseau
(1991)
Brusseau
Rabideau and
Miller (1994)
Powers et al.
(1991 and
1992)
Miller et al.
(1990)
Imhoff et al.
(1992)
phases and
components
considered
Sw>0, krw=0
SN=0
G W S
Pn , Pn , (On
Sw>0, krw>0
SN>0, krN=0
w s
Pn,»n
Sw=l
w s
Sw>0, krw>0
SG=0
w
Pn
Sw>0, krw>0
SG=0
SN>0, krN=0
w
Pn
exchange processes
modeled
AG Henry's law -
Pn/W E
As organic carbon
Pn/w partitioning -
EandK
pW EandK
Pn/w E and K
JVw EandK
pjf EandK
AW T£
P
kinetic
process
diffusive
boundary
layer
diffusive
boundary
layer
diffusive
boundary
layer
first-order
mass
transfer
first-order
mass
transfer
parameter
variability
constant rate coefficients
variable soil properties
isothermal
PG = constant
Siff = constant
pG = constant
constant rate coefficients
multi-component NAPL
heterogeneous soil due to
residual NAPL
isothermal
SN = constant
constant rate coefficients
layered soils
variable rate coefficients
homogeneous soils
consider variable NAPL
saturation and residual
blob shape
variable rate coefficients
homogeneous soils
consider variable NAPL
saturation and residual
blob shape
investigation
type
numerical
model:
advective
gas-phase
transport
numerical
model:
advective/
dispersive
water-phase
transport
numerical
model:
pump-and-
treat
application
numerical
model:
pump-and-
treat
application
experimental:
determine
rate-
coefficient
parameters
comment
• focus on rate-limited sorption
• differentiate between advective and
non-advective domains
• for gas-phase transport: Darcy's and
Fick's laws apply
• the limiting mass transfer step is solid
to water for water at residual
• focus on rate-limited mass transfer
• heterogeneity due to NAPL residual
greatly affects mass transfer
• rate-limited mass transfer is an
important consideration under
reduced-gradient conditions
• when SN = 0, mass transfer is most
sensitive to spatially variable K and
sorption capacity, and/or sorption
non-equilibrium with heterogeneity
dominating
• rate-limited mass transfer is important
for spills of small areal extent, high
Darcy velocities, large blob sizes,
and low residual SN
• diffusive boundary layer is an
inadequate exchange model
• exchange rate is a strong function of
water velocity and NAPL saturation
-------
GO
CD
O
CD
CO
c-t-
CO
O
CO
CD
X
o
CD
o
o
CD
CO
CO
CD
CO
reference
Sleep and
Sykes
(1993a, b)
Falta et al.
(1992a, b)
Delshad et al.
(1995)
Guarnaccia
etal. (1997)
phases and
components
considered
Sw>0, krw>0
Sr>0 kr>0
^G ' rG
SN>0, krN>0
oG ow cc£
'
Sw>0, krw>0
SG>0, krG>0
SN>0, krN>0
P°,P^
Sw>0, krW>0
SG>0, krG>0
SM>O k M>0
W
Pn
Sw>0, krW>0
SG>0, krG>0
SN>0, krN>0
G W S
Pn , Pn , (On
exchange processes
modeled
pn/w Henry's law - E
AW
Pn k
AG y
Pn
A s organic carbon
Pn/w partitioning -E
pn/w Henry's law - E
Aw E
AG T7
Pn fc
A s organic carbon
Pn/w partitioning -E
AW
Pn K
A G
pn/w Henry's law - K
AW ^
Pn ^
rll
0
AS organic carbon
Pn/Wpartitioning -E
kinetic
process
N/A
N/A
firet
lllSL
order
Kinetic
first-
order
kinetic
parameter
variability
phase densities are a
function of pressure and
composition
phase viscosities are a
function of composition
heterogeneous soil
multi-component NAPL
isothermal
use ideal gas law for gas-
phase properties
gas not a function of P
fluid properties are a
function of composition
and temperature
homogeneous soils
single-component NAPL
use ideal gas law for gas-
phase properties
water properties are a
function of composition
NAPL properties are
constant
heterogeneous soils
single-component NAPL
non-isothermal
water and gas density and
viscosity are a function
of NAPL concentration
fluid properties not a
function of P or T
heterogeneous soils
single-component NAPL
with constant properties
investigation
type
numerical
model:
3-phase, 3D
compositional
simulator
numerical
model:
3-phase
steam
injection
simulations
numerical
model:
3-phase
NAPL
remediation
using
surfactants
numerical
model:
3-phase
NAPL
remediation
comment
• multi-component DNAPL rem-
ediation in heterogeneous soil at
the field-scale
• investigate the effect that infiltrating
water fronts have on DNAPL
vapor transport
• local chemical and thermal
equilibrium
• NAPL saturation allowed to go to 0,
but water and gas saturations forced
to be >0
• model emplacement of residual NAPL
and rate-limited NAPL dissolution
• fluid properties a function of
surfactant concentration
• residual NAPL a function of capillary
and Bond numbers
• compare results with field data
• model emplacement of residual NAPL
and volatilization, dissolution and
adsorption mass transfer
-------
(b) Rate-limited mass exchange model: (indicated by 'K' in the review
figures) Based on the assumption that mass transfer between phases
is limited by diffusive transport across a stagnant boundary layer ac-
cording to Pick's law (diffusive boundary layer). The model has limited
utility because the thickness of the boundary layer cannot be estimated,
and it neglects advective, viscosity- and density-driven processes.
(c) First-order kinetic mass transfer model: (indicated by 'K' in the review
figures) More general than the rate-limited mass exchange model in
that it considers advective- diffusive processes and changing interfacial
contact area. Theoretical considerations and experimental data are
used to define the functional form of the coefficients.
3. The use of rate-limited mass transfer provides predictive flexibility, as local
equilibrium conditions can be simulated by increasing the rate coefficient.
4. Factors which favor rate-limited mass transfer include: heterogeneous soils,
inhomogeneous residual distribution, inhomogeneous blob size distribution,
and high fluid flow rates.
5. Modeling pump-and-treat remediation of NAPL-contaminated soils must
include rate-limited mass transfer to mimic both experimental- and field-
scale data, specifically effluent concentration tailing.
6. With respect to modelling desorption of NAPL from liquid unsaturated
soils, where the transport path of an adsorbed NAPL species is soil-to-
water, water-to-gas, data suggests that one need consider only one of the
mass transfer processes as rate limited since the rate-limited process will
dictate the overall mass entering the gas phase.
7. Henry's law is usually used to define vapor concentrations of a NAPL species
dissolved in the water phase.
8. The amount of NAPL mass which can be adsorbed onto the soil is usually
defined using an organic carbon-based model.
4.8.1. Liquid-Liquid Mass Transfer
When the organic phase is at an immobile residual state, saturation is no longer
considered a function of capillary pressure since capillary pressure becomes unde-
fined. Consider the NAPL phase balance equation 4.9 for the case of an immobile
35
-------
residual with constant phase density, constant porosity, and no external sources
or sinks. For these conditions equation 4.9 reduces to:
£ptfdSN_ = _Ew _ EG (4_27)
This equation states that change in NAPL saturation is due to mass transfer
processes.
The dissolution model defining the mass exchange term, E™, is assumed to
be a first-order kinetic-type reaction of the form:
(4.28)
where C™' \\/T] is the rate coefficient which regulates the rate at which equilib-
rium is reached, and ~p% [M/L3] is the equilibrium concentration of the NAPL
species in the water phase (solubility limit). In the simulator, ~p% is assumed to
be a measurable constant value.
To determine the parametric from of C^ , the work of Imhoff et al. (1992) is
employed. They conducted column experiments designed to study dissolution ki-
netics of residual trichloroethylene (TCE) in a uniform sand by flushing the system
with clean water and tracking the dissolution front as a function of time. Using
a lumped parameter model, they derived the following power-law relationship for
cw =
vw
where/52 ~ 0.5 and /53 ~ 1.0 are dimensionless fitting parameters. The parameter
fii [1/T] is the rate coefficient, and it is fit to available experimental- or field-
scale data.
The volatilization model defining the mass exchange term, E^ , is assumed to
follow a similar model as for dissolution, i. e.:
PG (~iG f-?G _G\ (A qn\
hn = ^n (Pn ~ Pn) (4^U)
where C^ [i/T] is the rate coefficient which regulates the rate at which equilibrium
is reached, and ~p^ [M/L3] is the constant equilibrium vapor concentration of the
NAPL species in the gas phase (vapor solubility limit). The rate coefficient, C%
is assumed to have the form:
(4.31)
36
-------
where/^2 is the same as for the dissolution model, and (3l [1/T] is fit to available
data.
Consider now the volatilization of a dissolved NAPL species in the water phase
to the gas phase. Assuming that the water phase is at residual saturation in the
vadose zone, and that there are no external sources or sinks of mass, then equation
4.8 can be written as:
dpw
£Sw—r^~ = En - En/w - En/w (4.32)
where the exchange term E™' is defined in equation 4.28, and E^/w governs the
volatilization mass transfer of a dissolved NAPL species in the water phase to the
gas phase:
(4-33)
where H is the dimensionless Henry's law coefficient which is defined at equilib-
rium conditions as follows:
H = P^/P^
and C^/w [V-^1 *s the mass transfer rate coefficient which is assumed to be defined
by the power law:
C%/w = (3fv (eSwf* (4.34)
where the fitting parameter (32 'ls assumed to be the same as for the liquid-liquid
mass transfer models, and (3l \\/T\ is fit to the available data.
4.8.2. Liquid-Solid Mass Transfer
Finally, mass exchange due to adsorption, E^,WJ is assumed to be defined by a
linear equilibrium model:
usn = Kdp™ (4.35)
where K& is the distribution coefficient [L3/M] defined as a function of the organic
carbon content of the soil and the relative hydrophobicity of the dissolved NAPL
species:
Kd = focKoc (4.36)
where foc is the mass fraction of organic carbon and Koc is the organic carbon
partition coefficient. Combination of equations 4.3 and 4.35 yields the following
definition for Esn,r:
/W = P
dt
37
-------
where pb = [1 — e\ps is the bulk density of the soil.
4.9. SUMMARY
The physical problem can be cast into a mathematical representation consisting
of five mass balance equations: three fluid-phase balance equations [equations
4.8, 4.9 and 4.10], and two NAPL species balance equations [equations 4.1 with
(i,a) = (n,W] and(n,G)].
The five primary variables used in the simulator to solve the balance equations
are Pw/, Sw-, STW-, P^, and Pn- All other parameters which constitute the balance
equations are assumed to be known physical constants or functions of the primary
variables.
The required physical constants are:
1. The reference fluid properties density and viscosity: pw/r, pNr', pGr, yUW/r,
nNr nGr
H , fj, .
2. The interfacial tensions, CFQWI &NW and &GN- These parameters are used
to scale the capillary pressures as per equation 4.14. Their magnitudes are
constrained to be related by equation 4.15. This constraint is required to
model the general case of two- and three-phase flow. The two independent
values can be considered fitting parameters.
3. The parameters which define water- and gas-phase dispersion: a™ , a^ , D^ ,
af, 4, Dg.
4. The parameters which define the mass transfer include: the rate coefficients:
/3^w, ^N, /3fw and the fitting parameters: /32, /33, the solubility limits:
P« , p^, the Henry's law coefficient: H, and the organic carbon-based par-
titioning parameters: /oc, Koc,
5. Soil properties porosity and permeability, where permeability is assumed to
be a scalar.
The constitutive relationships which relate the primary and secondary vari-
ables are:
1. The saturation-pressure model (S — P model) is predicated on the assump-
tion that phase wetting follows, from most to least, water-NAPL-gas, and
38
-------
when all three phases are present, only two saturation-capillary pressure
functions exist, PCNW(SW} and PCGN(STW)- Requiring that
• PCNW(SW) and PCGN(STW) are continuous in time and space regardless
of phase configuration,
• PCGW = PCNW + PCGN, and
• the capillary pressures are related through interfacial tension scaling,
equation 4.14,
allows the simulator to model the general case of two- and three-phase flow
without changing primary variables. An in depth presentation of the S — P
model is provided in Chapter 5.
2. The saturation-pressure model (k — P model) is given by equation 4.22. A
detailed discussion of these functional dependencies is provided in Chapter
5.
3. The density and viscosity of the water and gas phases are functions of NAPL
species mass concentrations only, equations 4.17, 4.18, 4.19 and 4.20;
4. The phase velocity is defined using Darcy's law, equation 4.21, and it is
considered for the three phases and two species of interest.
5. The phase dispersion coefficient is defined by equation 4.25, and it is used
to define the dispersive mass flux of the NAPL species in the water- and
gas-phases.
6. The following predictive mass transfer models are employed:
• a first order kinetic model for dissolution and volatilization mass transfer,
equations 4.28, 4.30 and 4.33; and
• an equilibrium adsorption/desorption model, equation 4.35.
39
-------
5. HYSTERETIC K-S-P MODEL
With respect to modeling three-phase flow in porous media, a substantial amount
of research defining appropriate constitutive relationships exists in both the petro-
leum reservoir and the water resources literature. From an historical perspective,
the development and application of many of the basic physical models quantifying
the physics of three-phase flow in porous media appeared first in the petroleum
reservoir literature. These models are applied to predictive modeling of oil and
natural gas recovery from petroleum reservoirs. This early work is summarized in
the books of Collins (1961), Aziz and Settari (1979), and Marie (1981).
Relatively recently, the problem of quantifying the spatial and temporal dis-
tribution of NAPL's in the mobile or immobile residual state, due to surface or
near-surface release(s), has become a major concern in the water resources area.
All of the predictive models derived for this application are built upon the ap-
plicable theoretical underpinnings previously derived in the petroleum reservoir
literature. However, because of the substantially different physical problem (ge-
ologic environment, fluid properties and driving forces), and engineering goals
(NAPL recovery techniques) encountered within the two disciplines, there is a
need to augment existing theory for petroleum reservoir applications, and to de-
velop new sub-models to describe physical processes unique to water resources
applications.
Considering this perspective, a review of the state-of-the-art modeling tools
available to describe three-phase flow in near-surface granular soils is provided in
Figures 5.1 and 5.2. These summary figures highlight the major physical model
components, where the S — P model describes the functional relationship between
saturation and capillary pressure, and the k — S model describes the functional
relationship between the relative permeability and saturation. In addition the
following parameters are defined: Sa, a = W (water), N (NAPL), and G (gas), is
the a-phase saturation (the percent of the pore space occupied by the a-phase),
STW = Sw + SN is the total liquid saturation, PCQ/3 = Pa — P'3, is the capillary
pressure between nonwetting and wetting immiscible fluids a and /3, respectively,
40
-------
and kra, a = W, N, G, is the a-phase relative permeability (scales the intrinsic
permeability to account for variable a-phase saturation).
Based on conclusions able to be drawn from review of the models summarized
in Figures 5.1 and 5.2, a relative permeability-saturation-capillary pressure model
(k — S — P model) is derived which:
1. is capable of modeling the simultaneous flow of a water phase (W), NAPL
phase (N), and gas phase (G);
2. is based on wettability considerations and two-phase data;
3. reduces to the appropriate two-phase and single-phase cases;
4. employs capillary pressure scaling to accommodate variable soil and fluid
properties;
5. includes flow-path, history-dependent functionals (hysteresis);
6. includes a mechanism for fluid entrapment as it is displaced from the pore
space.
For completeness and clarity in presentation the concepts which characterize
the k — S — P model that were introduced in Chapter 4 will be reiterated herein.
The presentation of the three-phase hysteretic k — S — P model is organized as
follows. First, a set of physical and model constraints which will allow for the ac-
commodation of attributes 1 through 4 above will be presented. This will provide
a basis for the definition of functional dependency. Second, with this conceptual
model as a foundation, a closed-form empirical k — S — P model which incorpo-
rates attributes 5 and 6 above will be developed as a series of three interrelated
sub-models: entrapment/release, saturation-pressure, and relative permeability-
saturation. In addition, for clarity, the k — S — P model will be developed first
for the two-phase case. The two-phase model will then be augmented to accom-
modate three-phase flow. Third, a section is included which addresses k — S — P
model implementation issues. These issues include details on how certain model
parameters affect the shape of the k — S — P functionals, and details on the def-
inition of a set of computations which will ensure that the functionals generated
are well behaved and amenable to inclusion into a numerical model. Finally, the
use of capillary pressure scaling in the simulator to account for variable fluid and
soil properties is discussed.
41
-------
t± hrj
'
to
g
CO
O
O
O
g
B
o
5T
8-
CO
o
reference
Eckberg and
Sunada (1984)
Abriola(1984)
Pinder and
Abriola(1986)
Faust (1985)
Faust et al.
(1989)
Parker et al.
(1987),
Kuppusamy
et al. (1989),
Lenhard and
Parker (1988),
Lenhard et al.
(1988), Kaluar-
achchi and
Parker (1989)
Parker and
Lenhard (1987)
Lenhard and
Parker (1987),
Lenhard et al.
(1989),
Lenhard(1992)
Ferrand et al.
(1990)
S-P
model
SW(PCNW)
C /T> \
5TwCrcGN)
SW(PCNW)
STW(PCGW)
SW(PCNW)
^TW(PCGN)
^W(PCNW)
STw(PcGN)
SW(PCNW)
Slw(PcGN)
SW(PCNW)
S/D \
TW("CGN)
k-S
model
N/A
krw(Sw)
krN(Sw,SG)
krw(Sw)
krN(Sw,SG)
krw(Sw)
krN(Sw SG)
krG(SG)
krW(Sw)
krG(SG)
N/A
hyster-
esis
no
no
no
no
both
S-P
and
k-S
no
entrap-
ment
no
no
no
no
NAPL
by water,
gas by
water
and
NAPL
no
numerical
model
yes
yes
yes
yes
yes
no
air phase
pressure
constant
constant
constant
constant
constant
constant
PC
scaling
none
none
none
surface
tension
surface
tension
surface
tension
verification
laboratory
experiment
2-phase
data only
field
data
laboratory
experiment
2-phase
S-P model
by lab-
oratory
experiment
laboratory
experiment
comment
Drainage only; Sw(PcNW) was a
good predictor, STw(PcGN) was not.
Assume this model is applicable
when water is mobile.
Field-scale simulations
• analysis of monotonic drainage
• discrepancies in STw(PcGN) funct-
ional due to constant gas presure
assumption
• surface tension scaling is a valuable
tool if real data is not available
• 1 -way transition from W-G to
W-N-G system
• poor mass balance reported
• surface tension scaling provided
qualitatively good results
• need to better determine nonwetting
phase residual saturations
• a wetting phase cannot be trapped
by a nonwetting phase
• numerical experiments indicate that
entrapment processes dominate
the solution
Drainage only; S-P model was
adequate except near residual
-------
comment
§
"ed
O
'S
on
P?l
W g
-S i
, U
03 QH
*«
O
S 15
a 1
i
&
ed -£
II
fe
l-'l
CC ^3
J^ §
^
PH ^
1 O
tD
a
(D
-
(D
"S —
cti | 1 S
t§ oZ £
o
a
^S-
^ C/3 C/3
Is ^ ^6
^? ^S* ^?
1 1
P^ fe
on on
*«
T3 ' — i "S
§ 8^ F3
•^ ^ S
"^ o « ^
£^^c^
a
•s"3 -°
S § B^
• put constraint on interfacial 1
to ensure proper S-P funct:
for 2- and 3-phase flow
• force a constant residual NA
modifying krN functional
£
8 1
2'C
it
8§
ll
u
-g —
ed | 1 5
n O ^ fD
O
a
^s—
C/3 c/3 C/3
^ ^ ^6
i? i? ^?
^ 1
p? fe
x s
^ H
ai ai
— ; ^^
*- "r:-
S ^
£c^
a
'•S o
•g E
Hysteresis, entrapment and va
in hydraulic properties are the
important modeling features
2 =3
2 *
8§fet
ill 1
•g
*3
a
o
o
1
|3
hj C3 >^hj
II- ll
13
1 ^H S
^ 5
^£ ^
i? i?
1 1
PH 9j,
^> 5
^ H
on on
i
ed
"S ^_^
13 m
'S^
S (^
w-
g
•§•3 -^
§ O £u
• put constraint on interfacial t
to ensure proper S-P functi
for 2- and 3-phase flow
• force a constant residual NA
modifying kr» functional
2 "o
Q. '-P
o •"§ a
o > rt
8 c
ll
<0
3
•-H
ed
1
+- „
C$ K-J §
llzl
o
^ fe" 'o
KJ M CC
^ ^ ~0
i? ^? ^
1 1
PH ^S
^> S
;> H
on on
IT)
*^3 -^ ^
^ wf ' — '
&»^ §•
8 ^ o> 8
X x ~~^x
"B 1
•a & ^S 3
•S 2 1 § 1
• surface tension scaling provi
qualitatively good results
• does not account for NAPL t
by gas
• numerical experiments indie;
entrapment processes domi
the solution
• solution sensitive to soil hyd
properties
£
0 |
2'G
II
o o
11
| 1
^Z jzi OD & 2^
o
-3s--
^ C/3 C/3
Is ^£ ^6
.i? .i? .M
1 1
PH fe
^ s
> H
GO 00
_^ w ^
8-3 rf
^ C« c£
>§-
a, t3
00^ S ^
•a o"§ g 2
• agrees well with published ai
• accounta for NAPL trapped 1
• numerical experiments indie;
entrapment processes domi
the solution
• solution sensitive to soil hyd
properties
•8 "c ^
§"0 S"" S
oV^TS.S B
1 1 1 •§ •§ 1
O rt G &^5 T3
8 § <5 ^
2 u S fl
yi -^- TO G
U
3
•£
ed
>
(D
"8 t
g (D Q_!
Cti QH*^ O (D
O 1 *§ 1
O ^ 5^ ^
^ S_
^ cc cc
1 ~S ^2
i? i? i?
1 1
PH° Pj,
x s
^ H
Ol Ol
rt f-^
•g^
C3 ' '
E *~?
S t3
Figure 5.2: Continued: Summary of conceptual and numerical models describing
three-phase flow in granular soils.
43
-------
5.1. CONCEPTUAL MODEL OVERVIEW
If fluid wettability is constrained to follow, from most to least, water-NAPL-gas,
and the intermediate wettability of the NAPL phase is assumed to render the
water-gas phase interactions negligible, then three-phase flow behavior can be
gleaned from readily available two-phase data as follows:
water NAPL gas
(5.1)
S-
k-
P
S
Sw — Si
v (PCNW)
krW = krW (Sw)
SN
krff
1O O
O M/ O (~2
1 ( O O \
— Ky-]\l I iJ p-|/ ^ O (^ 1
^G = ^
^rG =
?G (^cG^v)
^rG (^G)
This representation allows us to combine saturations based on wettability con-
siderations, where the total wetting phase saturation relative to the gas phase
is STW = Sw + SN, implying STW(PCGN) = 1 — SG, and the total nonwet-
ting phase saturation relative to the water phase is STU = SN + SG, implying
SW(PCNW) = 1 ~ STU- The motivation for using equation 5.1 as a representation
for three-phase flow is two fold. First, three-phase k — S — P data is generally not
available. Second, many natural subsurface systems conform to the wettability
and spreading assumption, most notably the cases where the NAPL is a hydro-
carbon derivative and the soil contains a high percentage of quartz and a low
percentage of organic matter. However, as noted in several recent publications
[most notably, Bradford and Leij (1995 and 1996)], there are many natural systems
which do not conform to this assumption, namely those systems where the NAPL
is a mixture of many different chemicals, and where the soil is mineralogically
heterogeneous and contains significant amounts of organic matter. Therefore, the
wettability assumption utilized herein must be recognized as a limiting attribute for
model applicability.
In order to model the general conditions of two- and three-phase flow, a second
set of constraints is imposed. First, for the two-phase water-gas case, PCQW is
determined from the algebraic constraint:
PcGW = PcNW + PcGN (5-2)
Second that the magnitude of the capillary pressure between any two immiscible
fluid phases is proportional to the interfacial tension between those phases (after
Leverett, 1941):
PcNW PcGN PcGW
44
-------
where crow? GNW and &GN are the interfacial tensions along the interfaces between
the gas and water phases, the NAPL and water phases, and the gas and NAPL
phases, respectively. Finally, note that for equations 5.2 and 5.3 to be compatible
the value of the parameters GGWI &NW and O"GJV must be constrained to be related
as:
&GW = &NW + &GN (5.4)
Equation 5.4 defines a neutral spreading coefficient. Therefore, only two of the
interfacial tensions are independent with the third defined by equation 5.4. The
two independent values can be considered fitting parameters. Several recent pa-
pers have questioned the use of constraints 5.3 and 5.4 (Wilson et al, 1990 and
McBride et al., 1992, among others). To date, the sensitivity of the solution to
this type of model constraint has not been formally addressed.
Section 5.5 continues the discussion of the use of capillary pressure scaling in
the simulator.
5.2. TWO-PHASE MODEL
Hysteresis in the k — S — P relationship is assumed to be caused by two main
processes, fluid entrapment effects and capillary and contact angle effects. These
effects are summarized as follows:
1. Fluid entrapment effects - When a fluid is drained from the pore space,
a volumetric fraction is rendered effectively immobile1 either by capillary
isolation in the case of a nonwetting phase, or by capillary and adhesive
forces in the case of a wetting phase. Trapped volumes can be re-mobilized
either by physio-chemical processes (for example, by changing the capillary
or Bond numbers) or by imbibition of the displaced phase. In addition, if
the fluids are slightly miscible, mass-transfer processes can result in a total
reduction of entrapped phase.
2. Capillary and contact angle effects - Assuming that the pore throats control
the wetting-phase drainage process and the pore bodies control the wetting-
phase imbibition process, and that the curvature of the menisci are enhanced
during wetting-phase drainage, the characteristic capillary pressure defining
displacement will be higher for drainage than for imbibition.
lrrhat is, the volume of fluid will no longer respond to a gradient in potential.
45
-------
Figure 5.3 presents a qualitative look at the two-phase hysteretic S — P curve-
type summary for the case when both capillary and entrapment effects are in-
cluded. The nomenclature for curve-type name describes whether the flow-path
is draining (D) or imbibing (I) with respect to the wetting phase, and how the
curve was spawned, primary (P), main (M), or scanning (S). Table 5.1 provides
a descriptive summary of the curve-types considered in this development, where
the curve-type numbering scheme aligns odd numbers with drainage curves and
even numbers with imbibition curves. Note that while the subscripts W and G
indicate water and gas phases, respectively, they can be considered to represent
a generic wetting-nonwetting phase pair. In addition, the parameters Sar are the
residual a-phase saturations, where the residual saturation is defined herein as
the saturation at which the phase effectively no longer responds to a gradient in
potential over the time-scale of interest.
Finally with respect to Figure 5.3, note that the different pressure scales char-
acterizing drainage and imbibition curves are a result of capillary effects, and the
end-points defining curves 1 through 4 are determined by fluid entrapment effects.
Figures 5.4 and 5.5 present a qualitative look at the two-phase hysteretic k — S
curve-type summary for the data presented in Figure 5.3. Hysteresis in krW(Sw)
and AVG^G) is assumed to be caused by fluid entrapment effects only. To see the
connection between hysteresis in the k — S relationships and phase entrapment, let
us define the entrapped phase volume as that volume which will no longer respond
to a gradient in potential. Since the two-phase extension of Darcy's law applies
only to that quantity which is hydraulically connected, the trapped volume must
be regarded as part of the porous matrix. Therefore, kra is not only a function of
Sa but specifically a function of that percentage of Sa which is hydraulically con-
nected. Figures 5.4 and 5.5 show that the magnitude of the relative permeability
at a particular saturation is a function of entrapment considerations. Note that
another definition of residual saturation is the saturation at which the relative
permeability is zero.
5.2.1. Entrapment and Release Sub-Model
Consider the following example of a water-gas displacement process intended to
qualify the salient features of the conceptual model for phase entrapment and
release, where water represents a general wetting phase and gas represents a gen-
eral nonwetting phase. A porous medium is initially saturated with water which
46
-------
0 0.2 0.4 0.6
Swr
water-phase saturation
0.8 1
(l-Sgr)
0 0.2 0.4 0.6
Swr
water-phase saturation
0.8 1
(l-Sgr)
Figure 5.3: Definition plot of the hysteretic relationship between saturation and
capillary pressure employing the empirical model used in the simulator. Curve po-
sition and shape is governed by the mobility status and the magnitude of the phase
saturations when the curve is spawned and whether the displacement process is
drainage (£)) or imbibition (/) with respect to the wetting phase. Primary (P)
and Main (M) curves are spawned when only one phase is mobile [curve numbers
1 and 2]. Scanning (S) curves are spawned when both phases are mobile [curve
numbers 3 to 6]. The curve-type numbering scheme is set such that odd numbers
are aligned with drainage curves [1 = PDC or MDC, 3 and 5 = SDC] and even
numbers are aligned with imbibition curves [2 = PIC or MIC, 4 and 6 = SIC].
47
-------
0.2 0.4
Swr
0.6
0.8 1
(l-Sgr)
1
0.9
.
OJ
>
I 0.5
i-H
OJ
w
f.O-4
f
I 0.3
0.2
0.1
(3
0.2 0.4
Swr
wetting-phase saturation
0.6
wetting-phase saturation
0.8 1
(l-Sgr)
Figure 5.4: The resulting hysteretic k — S functionals for the wetting phase gen-
erated from the empirical model used in the simulator and the data denning the
S — P relationship. Note that the MDC and MIC', shown as dashed lines, are
practically coincident, and that the scanning curves, shown as dashed lines, are
group-lab led because by model definition they are coincident.
48
-------
1
0.9
0.8
120.5
fO.4
0.2
0.1
0
MIC*
(2)*
0 0.2 0.4 0.6 0.8 1
Swr (l-Sgr)
wetting-phase saturation
1
0.9
120.5
8
.a
P-0.4
gf
0.2
0.1
0
0
0.2 0.4 0.6
Swr
0.8 1
(l-Sgr)
wetting-phase saturation
Figure 5.5: The resulting hysteretic k — S functional for the nonwetting phase
using the data denning the S — P relationship and the empirical model used in
the simulator. Note that the curve labeled MIC* is obtained upon reversal from
a PDC where no nonwetting phase was previously trapped. Subsequent reversals
follow the MDC and MIC which are practically coincident. Also note that the
scanning curves are group-labled because by model definition they are coincident.
49
-------
number
1
2
3
4
5
6
name
PDC-Priniary Drainage Curve
MDC- Main Drainage Curve
MDC*
PIC-Priniary Imbibition Curve
MIC- Main Imbibition Curve
MIC*
SDCi- First-order Scanning
Drainage Curve
SIC}- First-order Scanning
Imbibition Curve
SDC2- Second-order Scanning
Drainage Curve
SIC2- Second-order Scanning
Imbibition Curve
qualitative notes
- origin at S\y = l
- origin at S\v = l-SQi- and only the
wetting phase is mobile
- origin at l>S\y>l-SQi- and only
the wetting phase is mobile
- origin at S\y=0
- origin at S\y=S\Vr an
-------
is considered a continuum, and therefore, the initial trapped quantity is denned
as zero. Water is then displaced by the gas until the water effectively ceases to
flow as water films surrounding the soil grains thin to the point that adhesive and
capillary forces dominate pressure gradients. Over the displacement process the
trapped quantity of water rises from zero to a maximum, called the irreducible or
residual. The residual wetting phase saturation, Swn is a function of the porous
medium properties and the chemistry of both the wetting and nonwetting phases.
The concept of a saturation-dependent entrapment mechanism is based on the
interpretation that as the different pore classes are drained (the larger ones first),
some water is rendered effectively immobile as it adheres to the empty pore walls.
This displacement process is called primary drainage (the PDC in Figure 5.3).
With respect to the gas phase, it is assumed that no trapped saturation results
as it displaces the water.
If at this point (that is, at Sw = SVr) the water phase is imbibed, thereby
displacing the gas-phase, two phenomenon occur. First, some of the trapped water
phase is remobilized as the in-flowing water thickens the water film associated
with the smaller pores thus reducing the trapped quantity. Second, as the gas
phase is displaced, a volumetric fraction becomes disconnected from the flowing
volume and isolated in a pore or series of pores by the mechanisms of snap-off
and bypassing [see Wilson et al., 1990, pages 110 to 123]). The resulting trapped
gas saturation rises from zero to a maximum residual value where the gas phase
ceases to flow. The residual nonwetting phase saturation, SGV, is a, function of the
porous medium properties and the chemistry of both the wetting and non-wetting
phases. This process is called main imbibition (the MIC in Figure 5.3). At this
point (that is, Sw = 1 — SGV) the water-phase is the only mobile phase, however,
since the water saturation did not return to unity, some trapped water phase can
still exist (for example, in some of the larger pores predominately filled with gas).
In addition, as the gas dissolves in the water phase over time, the water saturation
approaches unity and, as it does, the trapped quantity approaches zero.
The above displacement experiment involves the situation incurred when full
drainage and imbibition cycles are realized (primary and main curve-types). That
is, displace one fluid until it ceases to flow. Three conceptual attributes are high-
lighted in this description. The first two involve flow-path dependent entrapment
and release mechanisms, where trapped quantities increase as a phase drains from
the pore space and decrease as the phase fills the pore space. The third is that
if the fluids are slightly miscible, mass transfer processes can result in a total
reduction of entrapped phase.
51
-------
The final conceptual attribute has to do with the case of an incomplete dis-
placement process and the spawning of a scanning curve. That is, what happens
when a flow reversal occurs before one of the phases is rendered immobile. For ex-
ample, consider the SICi (curve 4) in Figure 5.3. First, at the point of premature
reversal along the PDC, according to the previous argument the trapped water
phase will have a value somewhere between the extremes, zero and full residual
(5Vi/r), and the gas phase will have zero trapped volume. In general, an imbibing
nonwetting phase tends to fill the larger pore classes first, so from a physical point
of view, it seems reasonable to assume that if incomplete drainage occurred before
reversal, then not all the different pore classes would be filled by the gas-phase,
and upon water imbibition, the full gas-phase residual would not be realized. In-
stead, a value between the extremes would be appropriate, as shown by curve 4
in Figure 5.3.
This concept was first quantified by Land (1968), who developed an empirical
model, based on experimental observation, to estimate the residual nonwetting
phase saturation magnitude. The model utilizes the assumption that the difference
in the reciprocals of the initial and residual nonwetting phase saturations is a
constant for a given sand. That is, for the gas- water system described above, if
one assumes that the maximum possible gas-phase residual, S^y, is obtained for
the case of gas-phase drainage from SG = 1 (i. e., for the PIC, all the pore classes
are initially filled with the gas-phase), then a reduced residual, SQT < SGVJ will
be obtained for the case of gas-phase drainage from SG = SQ^^ < 1 (i. e., for an
MIC or SICi , not all the pore classes are initially filled with gas-phase) according
to the rule (after Land [1968]):
Given this qualitative interpretation of the physics of phase entrapment and
release, let us now derive an analogous quantitative empirical representation. To
facilitate our discussion on fluid entrapment, let us differentiate between the fluid
volume which is 'free' to respond to a gradient in potential (i. e. hydraulically
connected) and the volume which is 'trapped' and cannot respond to a gradient
in potential (i. e. hydraulically disconnected). In terms of a-phase saturation
this differentiation is written as:
Sa = Saf + Sat (5.6)
where the subscripts / and t indicate free and trapped, respectively. It is clear
that free and trapped saturations can not be measured per se. Therefore, equation
52
-------
5.6 must be considered an empirical vehicle upon which to derive a model which
will approximate the qualitative aspects introduced above.
To this end, let us assume that the trapped quantity, Sat, can be described by
a saturation-dependent blending rule of the form:
[Cmax _ c "I
__2 _ 2.
Cmax C*
da ~ ^ar]
where Sat is constrained to lie within predefined, flow-path, history-dependent,
limits, S™tm < Sat < •S'ar- The parameters in equation 5.7 will be described in
turn. S™n is the lower limit of entrapped a-phase. It is intended to quantify
the condition which exists for MDCs and MICs where there is some non-zero
residual a-phase saturation at the time when the a-phase was re-imbibed. Given
the magnitude of the residual at the origin of a main curve, Sar(j), ct = G and j
= 1 for an MDC, and a = W and j = 2 for an MIC, S™n is computed as:
i _ c
u
where it is assumed that a linear relationship exists between S™n and Sa-
The parameter S™3* is the highest a-phase saturation that has occurred since
it was last at immobile residual conditions, i. e., the maximum imbibed a-phase
saturation which is available for displacement. Specifically, Sg3* is the furthest
progression along a PDC or MDC, and 5|jfx is the furthest progression along a
PIC or MIC.
The parameter S^r is the magnitude of the residual a-phase, 0 < S^r < Sar,
at the terminus of an a-phase drainage process. Its value is computed using
the model of Land (1968), modified from equation 5.5 to include the following
additional processes:
• a wetting phase trapped by a nonwetting phase;
• the existence of previously entrapped phase, which acts as a lower limit for
further entrapment.
The resulting empirical relationship defining S^,r is written as:
Cmax _ cmin
C* _
^
at
_i_ Jj> /^Cmax ^Cmin\
^min^ ~1
53
-------
Finally with respect to equation 5.7, e represents a blending parameter which
governs how fast the phase becomes entrapped during drainage flow conditions
or released from entrapment during imbibition flow conditions (e > 0, and for
example, e = I yields a linear relationship between Sat and Sa).
In summary, the trapping model, equations 5.7, 5.8 and 5.9, is designed to
represent any wetting-nonwetting, two-phase system, and it utilizes three fitting
parameters: Swn SQT and e. The remaining parameters, S£r, 5'™ax and S™n, are
functions of flow-path history as discussed above.
5.2.2. Saturation-Pressure Sub-Model
The two-phase hysteretic saturation-capillary pressure sub-model which is de-
scribed below is based on that derived in Luckner et al. (1989). As a point of
departure for this presentation, consider the van Genuchten (1980) saturation-
pressure function for monotonic, non-hysteretic, displacement written here as:
hc(Se) = [(Se)-l/m ~ I]'7" (a)"1 (5.10)
Se = ^—A 0 a
-------
between saturation and capillary pressure for any one of the primary and main
drainage and imbibition curves, / = 1 and 2 (Table 5.1 and Figure 5.3), provided
the following parameter set is known:
S — P M —/• {6(4/r, O(jr, dfi, o,j, r/}
where the notation S — P\M indicates that the S—P model is defined for monotonic
displacement. These parameters are determined either by fitting them to exper-
imental data, approximating them based on soil properties, or a combination of
the two (see for example Nielsen and Luckner, 1992).
Let us now modify the non-hysteretic S — P model defined by equations 5.10
and 5.11 to account for hysteresis caused by the two main processes identified
above: fluid entrapment and capillary and contact angle effects. Fluid entrap-
ment effects are accounted for by making the parameters Sr and Ss functions of
saturation, flow-path history, and the trapping parameters computed in equation
5.7. Hysteresis due to capillary and contact angle effects is accounted for by mak-
ing the scaling parameter, a, a function of saturation and flow-path history. The
parameter 77 is assumed independent of flow-path.
The resulting closed-form empirical relationship between saturation and cap-
illary pressure, assumed to be valid for primary, main and scanning curve-types,
is written as:
Im il1/'? I \-l (r 1Q\
' — 1J (a(f)) (5.13)
0 <&(/) Q Q '
Os(f) - br(f)
where the subscript (/) indicates that the parameter is valid for a specific curve-
type, i. e., / = 1, 2, 3, 4, 5, or 6, as referenced in Table 5.12. This hysteretic
(H) version of equations 5.10 and 5.11 can be computed provided the following
parameter set is known:
S — P\
ff
where the notation S — P\H indicates that the S — P model is defined such that
hysteresis is accommodated, and it is assumed that the parameters Swr, SGV,
a
-------
Table 5.2 provides the details for evaluating the scaling parameter am, where
So(.) is the reversal point saturation at which curve (•) was spawned, and j3 is
a blending parameter. First note that a^ < am < a,^ and for primary and
main curve-types (/ = 1, 2) the limiting values are utilized. Second, note the
use of a blending function to define a(/) for scanning curve-types (/ = 3, 4, 5,
6). The purpose of the blending function is for numerical model implementation
purposes. When a flow reversal results in the generation of a scanning curve, there
is a commensurate jump in the slope of the S — P functional. That is, while the
functional itself is continuous, its slope is discontinuous (for example, see Figure
5.3 curve 4). Since the slope of the curve affects the magnitude of the phase
mass flux, a strong discontinuity can lead to convergence problems in a numerical
scheme as the relationship between the flow parameters and the discrete time- and
space-scales is dramatically altered. For example, a large increase in mass flux may
require a much smaller time step (related to the Courant number) not only to limit
temporal discretization errors, but also to allow for convergence of the linearized
system of equations. Therefore, in order to limit numerical difficulties associated
with functional discontinuities, the blending function has been incorporated to
allow for a smooth transition in slope when a flow reversal occurs.
A blending parameter, (3, has been included to provide modeling flexibility,
where in general J3 > 0. For example, an instantaneous change in am occurs
when J3 = 0, and as J3 gets larger the transition from the reversal point value is
delayed. As an implementation illustration, consider a reversal from a PDC to an
SICi. As shown in Table 5.2, for the new curve 4, «(4)(5Vi/) varies from a^ at the
reversal point, 6*0(4), to a,i at the terminal saturation, Sw = 1 — SG^, according to
the power law defined by the parameter J3.
Given the data quantifying phase entrapment, equations 5.7 and 5.9, and
the value of am., the parameters SV(/) and Ss^ are determined by following the
scaling procedure of Scott et al. (1983) (see also Kool and Parker, 1987, Parker
and Lenhard, 1987, Lenhard and Parker, 1987, Luckner et al., 1989, and Lenhard
et al., 1991). Effectively the model requires that any S — P functional defined by
equation 5.13 pass through two predetermined constraint points chosen to ensure
the following:
constraint (1) That the functional either terminate at the proper residual sat-
uration as defined in equation 5.7 or terminate such that a closed scanning
loop results, whichever is appropriate.
constraint (2) That the functional remain continuous as definitive parameters
56
-------
f
1
2
name
PDC
MDC
MDC*
PIC
MIC
MIC*
bIC2
definition of Ob
\ di=ad
\ a2 = a-i
n n (n n \ \ S°W~Sw ]
UA at (a, ad) [So(3)_Swt\
, / \ T SW—SQ(^\ 1 "
a4 ad | (a, ad) ^^^J
[c c< 1 /5
oocsi— iny
c c
*0(5) -^0(4) J
n n 1 (n n } [ lSw'~5°(6) 1^
«b a0 1 ^4 «oj [s0(5)-So(6)J
computational details
grade change between
span of curve 3
grade change between
span of curve 4
grade change between fl/4\,
the value of CL when the
reversal occured, and d({
over the span of curve 5
grade change between dtfy,
the value of CL when the
reversal occured, and <2-(4)
over the span of curve 6
Table 5.2: DEFINITION OF THE SCALING PARAMETER 'a' FOR THE
TWO-PHASE HYSTERETIC S-P MODEL
57
-------
are altered to reflect changing entrapment and capillary effects.
Constraint f is required first to ensure that residual saturation data is honored.
For example, in Figure 5.3, curve (4) terminates at Ss = (f — «S(jr), the point
where the gas-phase ceases to flow, and curve (3) terminates at Sr = 5^/r, the
point where the water-phase ceases to flow, where SQT and 5|l/r are computed from
equation 5.9. The second part of this constraint, that there be closed scanning
loops, is required so that erroneous pumping effects do not occur during cyclic
flow reversals (Janes, 1984). For example, curve (5), a second order scanning
drainage curve, ties back in with the PDC at the point where the scanning loop
originated.
Constraint 2 is required because a continuous capillary pressure field is required
for mass conservation. It becomes an important consideration when viewed from
the standpoint that the S — P model is to be implemented in a time-discrete
numerical model where saturation is the dependent variable and where parameter
update is lagged in time with respect to saturation. That is, given the saturation
solution at the advanced time step (superscript n + 1), •Sjy'1, the corresponding
capillary head, h"+l is determined from equation 5.13 with the parameter set:
< £>w~l, 57(/), S™(f)i af/)> 'H f; where the superscript n indicates that the saturation-
dependent parameters are dated at the previous time step. To prepare for the next
time step, a,*-, is updated as discussed above, and «SV(/) and Ss(f) are updated such
that not only is constraint 1 satisfied, but that hc does not change value.
Tables 5.3 and 5.4 provide the quantitative details, consistent with the con-
straints listed above, for evaluating the curve-type-specific parameters ST^ and
Ss(f). From a qualitative perspective the following descriptive statements regard-
ing SV(/) and Ss(f) apply:
• For curves PDC, MDC, MDC*, SDCi (/ = 1 and 3): Sr = Swt [from
equation 5.7] so that constraint 1 is satisfied, and Ss becomes a scaling
parameter so that constraint 2 is satisfied, and as such, it has no physical
meaning.
• For curves PIC, MIC, MIC*, SDCi (/ = 2 and 4): S8 = (I - SGt) [from
equation 5.7] so that constraint 1 is satisfied, and Sr becomes a scaling
parameter so that constraint 2 is satisfied, and as such, it has no physical
meaning.
58
-------
f
1
2
3
name
PDC
MDC
MDC*
PIC
MIC
MIC*
SDCi
SIC!
definition of 9
r
1 Sr(i) Swt
I
[
} ^(2) l-SB
(
Sr(3} =Swt
Sw-SeS,w
^(4) l-SB
T5awsw--3.ssw
°'m S^~S'
~SB(.nSw-'SBSs(:r)
^(6) SeU)-Se
computational details
Swt fr°m equation 5.7
computed after s ,_,
has been updated
Swt fr'om equation 5.7
computed after s , .,
has been updated
force closure at the
ning loop originated
3 if imbibition and
spawned from an SDC^
spawned from an SDC2
6 if drainage
Table 5.3: DEFINITION OF THE SCALING PARAMETER 'Sr' FOR THE
TWO-PHASE HYSTERETIC S-P MODEL
• For curves SDC2 and SIC2 (/ = 5 and 6): both Ss and Sr become scal-
ing parameters so that the closed scanning loop part of constraint 1 and
constraint 2 are both satisfied.
• For application in a time-discrete numerical model, the S — P model en-
forces residual saturations, conserves mass, and converges as the increment
in saturation change approaches zero.
Summary and Evaluation Procedure
The empirical hysteretic S — P model is defined by equations 5.7, 5.13 and 5.14.
It requires the input of the following parameter set:
C, Swri SGT-, Q-di ai,r
59
-------
f
1
2
3
4
5
6
name
PDC
MDC
MDC*
PIC
MIC
MIC*
SDCi
SICi
SDC2
SIC2
definition of gs
(
I *«(!) 9 ' ^r(l)
\ *~ e
I
f
1
sw-sr(3} ,
*'s(3) 9 ' *'r(3)
Ss(4) 1— SG*
Sw-Sr(6} ,
•S's(S) 9 ' *'r(5)
SW-S
e
computational details
computed after g , ,
has been updated
S(-'+ from equation 5.7
computed after g ,„,
has been updated
SGf from equation 5.7
computed after ,9 ,
has been updated
computed after g , *
r( )
has been updated
Table 5.4: DEFINITION OF THE SCALING PARAMETER 'Ss' FOR THE
TWO-PHASE HYSTERETIC S-P MODEL
Given these values a sequential, step-wise, procedure is used to compute the
S — P functionals. The computational steps are listed as follows:
1. Given initial data for the parameter set: < 5|J/, 5'™m, S"'(f), a™n [-, compute
the solution for the next time step: Sff1. Given
compute /i™+1 from equation 5.13.
2. Compute Swt and S^t from equation 5.7 with Sw =
3. Determine which curve-type is appropriate given the change in saturation:
co qn+l an
and the current value of /.
/ is currently 1 or 3: If 6Sw > e, where e is a small positive number,
then a reversal is indicated and / is incremented by one, otherwise a
reversal is not indicated and / is not altered.
60
-------
/ is currently 5: A reversal is indicated if 6Sw > e and / is incremented
by one. A closed scanning loop is indicated if 5V < 5*0(4), where 5*0(4)
is the reversal point saturation at which curve 4 was spawned, and /
is reset to equal one. If neither of these things happen then / is not
altered.
/ is currently 2 or 4: If 6Sw < — e, then a reversal is indicated and / is
incremented by one, otherwise a reversal is not indicated and / is not
altered.
/ is currently 6: Since this is the highest-order curve considered in this
simulator, a flow reversal will not produce a new drainage curve. A
closed scanning loop is indicated if either 6^fl > 6*0(5), °r Sw"1 < 6*0(6),
where 6o(.) is the reversal point saturation at which curve (•) was
spawned. For the former, / is reset to equal four, and for the latter, /
is reset to five. If neither of these things happen then / is not altered.
4. Given the updated value of / from part 3, compute the value of ay) from
Table 5.2 with Sw = S$~l.
5. Given /i"+1 from part 1 and a^ from part 4, compute Se using equation
5.13. This step is included for the continuity requirement.
6. Given Se from part 5, compute 6S(/) and SV(/) as described in Tables 5.4
and 5.3 with Sw = 6™/1.
5.2.3. Relative Permeability-Saturation Sub-Model
A closed-form empirical model used to predict relative permeability for a two-
phase system is derived by employing the Mualem (1976) statistical model. The
Mualem model relates the water relative permeability to saturation by a series of
theoretical steps. First, the capillary pressure-saturation functional is assumed to
be analogous to the pore-size distribution function. Second, using the capillary
law, the capillary pressure at which a pore will drain or fill is uniquely related to
the pore radii. Third, the water relative permeability is obtained by integrating
over the contributions of water-filled pores:
krw(Sw) = (6,
eW)
61
eW
9(0,1)
(5.15)
-------
where
-B
g(A,B) = I h-ldSeW
JA
and Sew is an effective water saturation. It is a normalized saturation, analogous
to that denned by equation 5.11, which is a measure of the saturation range over
which water will flow. We identify herein two unique definitions of Sew which
apply to the krw(Sw) functional:
Definition 1. Sew = [reduced saturation] '/[total % pore volume available for
How]
SeW = (Sw - Swt) / (1 - SWt - Set) (5.16)
Definition 2. Sew = [reduced saturation] /[total % pore volume available for
water How]
SeW = (Sw ~ Swt) I (1 - Swt) (5.17)
Note, that the terms Swt and Sot are those which are computed from equation
5.7.
The term in brackets in equation 5.15 represents the ratio of the pores that
are contributing to water flow at a particular effective saturation, g(Q, Sew), to all
pores contributing to water flow, g(0, 1). The term pre-multiplying the integrals
accounts for tortuosity effects and incomplete correlation between pores, and £ is
a pore connectivity parameter.
Substitution of the van Genuchten saturation-pressure relationship (equation
5.10) into equation 5.15 gives (Parker et al., 1987):
g(A, B) = a(l- Al'm} -(I- Bl'm} (5.18)
and equation 5.15 becomes:
l/m~m
where the prescripts on the effective saturations indicate that, for model flexibility,
the user can choose not only which definition of Sew is to be used (i.e., equations
5.16 or 5.17), but can choose a different definition of Sew for each of the two terms
in equation 5.19. Note that this relationship indicates that krw(Sw = 1) = 13
and krw(Sw = Swt) = 0 (where Swt is the residual water saturation in the limit
as Sw approaches Swr, as per equation 5.7).
62
-------
An analogous derivation for the gas phase takes into account that the pores
will fill with the nonwetting phase from the largest to the smallest. As in Luckner
et al. (1989) and Parker et al. (1987), the relationship can be written as:
(5.20)
0(0,1)
where (p is a pore connectivity parameter for the gas-phase, and SeG is the nor-
malized or effective gas saturation which, when applied to the krG(SG) functional,
is assumed to have two admissible definitions:
Definition 3. SeG = [reduced saturation]/[total % pore volume available for How]
SeG = (SG - SGt) I (I - Swt - S^ (5.21)
Definition 4. SeG = [reduced saturation]/[total % pore volume available for gas
How]
SeG = (SG ~ SGt) I (1 - Sot) (5.22)
where the terms Swt and SQt are those which are computed from equation 5.7.
Substituting the van Genuchten saturation-pressure relationship (equation 5.10)
into equation 5.20 gives:
krG(SG) = (aSeGf [l-[l- (bSeG)]1/m] ^ (5.23)
where the substitution, 1 — SeG = Sew? nas been made, and the prescripts on
the effective saturations indicate the user can choose different definitions of SeG
(i.e., equations 5.21 or 5.22) for each of the two terms in equation 5.23. This
relationship indicates that krG(SG = !) = !, and krG(SG = SGt) = 0 (where SQt
is the residual gas saturation in the limit as SG goes to SQT, as per equation 5.7).
In summary, the two-phase relative permeability-saturation sub-model, equa-
tions 5.19 and 5.23, written for a water-gas system, is assumed valid for any
wetting-nonwetting phase system. It includes hysteresis as the result of fluid en-
trapment effects through the definition of the effective saturations. In addition, it
utilizes data from the saturation-pressure sub-model, the parameter TO, and data
from the entrapment/release sub-model, the parameters Swt and SGt. The k — S
sub-model includes two fitting parameters, £ and if>, and an additional model fit-
ting feature is included by allowing the user to choose the definition for effective
saturation (equations 5.16 or 5.17 and 5.21 or 5.22).
63
-------
5.3. THREE-PHASE MODEL
The two-phase k — S — P model derived in Section 5.2 is now adapted for the
three-phase case. To reiterate, the fundamental assumptions which are used to
construct the three-phase model are:
1. Water is the most wetting phase, and it spreads as a film over the soil grains.
NAPL has intermediate wettability, and it spreads as a film over the water.
Gas is the least wetting phase, and it is surrounded by the total wetting
phase, i. e., the water and NAPL. This idealization allows a three-phase
k — S — P model to be generated using two-phase data as shown by the
functional dependencies in equation 5.1.
2. The functionals PCNW(SW} and PCGN(STW) are defined such that PCGW =
PCNW + PCGN regardless of which phases are present in the domain. For
this to be the case then PCNW(SW} and PCGN(STW) must be continuous
in time and space for any phase configuration, and the three fluid pair
capillary pressures must be related through the interfacial tension scaling
relationship 5.3, with the interfacial tension values constrained to yield a
neutral spreading coefficient (equation 5.4).
3. The three-phase k — S functionals are defined such that they reduce to the
appropriate two-phase functionals when appropriate.
As with the two-phase model, the three-phase model can be idealized as three
inter-related sub-models: entrapment-release, saturation-capillary pressure, and
relative permeability-saturation. These sub-models are presented in turn.
5.3.1. Entrapment and Release Sub-Model
The addition of a third phase complicates the entrapment description, but the
following qualitative statements are assumed to hold:
• a-phase drainage results in an increase in trapped a-phase volume;
• a-phase imbibition results in a decrease in trapped a-phase volume;
• an incomplete displacement process will admit residual saturations less than
the maximum measured value.
64
-------
It is this conceptual model we wish to adapt into an empirical model to quantify
phase entrapment and release for three-phase immiscible flow.
Since the water phase is most wetting, we idealize the entrapment mechanism
as reducing to a wetting phase, Sw, being displaced by a nonwetting phase, SV«,
where STU = SN + SG- Therefore, we assume that the magnitude of Swt will be
the same regardless of which phase(s) displaced it. The same argument holds for
the gas phase: since it is most nonwetting, we idealize the entrapment mechanism
as reducing to a nonwetting phase, SG, being displaced by a wetting phase, STW,
where STW = Sw + SN. Therefore, the trapping model, equations 5.7, 5.8 and
5.9, is assumed to describe water and gas phase entrapment (a = W and G)
for both the two-phase water-gas system and the three-phase water-NAPL-gas
system where the NAPL has intermediate wettability..
With respect to the NAPL phase, since the trapping mechanisms are different
for wetting and nonwetting phases, the magnitude of trapped NAPL-phase should
be a function of which phase displaced it. From two-phase displacement experi-
ments, one can measure the maximum NAPL residual in a water-NAPL system,
SjVnr, and in a NAPL-gas system, Sp?wr, where the subscript Nnr indicates NAPL
residual as a nonwetting phase, and Nwr indicates NAPL residual as a wetting
phase. In general, S?jnr ^ Sffwr. Since we are using the same phase entrapment
model for both the water and gas phases, with the fundamental difference being
that the measured residual values are different, consider the following entrapment
model for NAPL. Like the water and gas phases, the NAPL entrapment is defined
using equations 5.7, 5.8 and 5.9, for a = N, and SVr is defined as a linear function
of water and gas saturations (after Payers and Matthews, 1984), i. e.:
i c
— - — — + bNnr
bw + ba/
In summary, the trapping model, equations 5.7, 5.8, 5.9 and 5.24, valid for a
= W, G and N, is designed to represent a three-phase system where the NAPL
has intermediate wettability between water and gas. In addition, it reduces to the
appropriate two-phase model when appropriate. It utilizes five fitting parameters:
Swn $Gr i Spfwn SjVnr and e. The remaining nine parameters, S£r, S™ax and S™11,
a = W, G, N, are functions of flow-path history.
5.3.2. Saturation-Pressure Sub-Model
In Section 5.1 we introduced the fundamental assumption that the three-phase
S — P model can be decomposed into two, two-phase S — P models by virtue of
65
-------
the fact that the NAPL phase is constrained to be the intermediate wetting fluid.
Based on the material presented in Section 5.1, the two relevant S—P relationships
are: PCNW(SW} and PCGN(STW)- As will be seen shortly, these relationships are
interrelated through fluid entrapment considerations.
First let us view each S — P relationship as a unique two-phase system, where
the two-phases are idealized as follows:
functional
PCNW(SW]
PCGN(STW)
wetting phase saturation
Sw
STW = Sw + S N
nonwetting phase saturation
STU = 5jv + SG
SG
It can be seen from this table that in the limit as Sw —» 0, PCGN(SN) represents the
two-phase NAPL-gas system where NAPL is the wetting phase, and that in the
limit as SG —> 0, PCNW(SW} represents the two-phase water-NAPL system where
NAPL is the nonwetting phase. In addition, in the limit as SN —>• 0, STW = Sw
and STH = SG, and PCGW = PCNW(SW} + PCGN(SW} represents the two-phase
water-gas system.
Second, from equation 5.3 we know that the magnitude of the capillary pres-
sure between any two phases is related through interfacial tension scaling. There-
fore, assuming that one of the phase pairs is used to fit the empirical model
parameters a
-------
functional
PCNW(SW)
PCGN(STW)
two-phase model variables used
in equations 5.13 and 5.14 and
in Tables 5.2, 5.3 and 5.4
hc
Sw
Swt
Sat
hc
Sw
Swt
Sat
analogous variable
in the three-
phase model
hcNW
Sw
Swt
Sm + Sat
hcGN
ST-W
Swt + 5jvt
SGt
Table 5.5: THE RELATION BETWEEN TWO- AND THREE-PHASE S-P
MODEL PARAMETERS
is computed by substituting the parameter set: {hc^wj Swr, Swrt, (Sm + ^Gt)}?
and the PCGN(STW) functional is computed by substituting the parameter set:
{hcGN-, STW-, (Swt + Sfft) , Sat}-
In summary, the three-phase hysteretic S — P model requires the following
input data set:
1. the curve-fit parameters: ad, ai; 77, as defined by one of the three fluid pairs
2. the four residual saturations: Swr, Sj\fwr, Sj\fnr, SGV;
3. the blending parameter for entrapment, e;
4. two of the three interfacial tension parameters: CFGW, &NW, °"GJV, with the
third defined such that &GW ~ &NW ~ ^GN = 0.
5.3.3. Relative Permeability-Saturation Sub-Model
The three-phase k — S model is built upon the following fundamental assump-
tion: the most wetting phase fills the smallest pores, the least wetting fills the
largest pores, and the intermediate wetting phase fills the intermediate size pores.
The main outcome of this assumption is that the wetting and nonwetting phases
become spatially segregated to the extent that their relative permeability func-
tions become dependent only upon their respective saturations. In addition it
is assumed that the relative permeability functions apply for two- or three-phase
67
-------
systems. Therefore, for the three-phase system, as long as the wettability assump-
tion holds, the water-phase relative permeability is given by equation 5.19 and the
gas-phase relative permeability is given by equation 5.23 with the definitions for
effective saturations (i.e., equations 5.16, 5.17, 5.21, and 5.22) altered to account
for the presence of trapped NAPL, i. e.,
Definition 5. Se =[reduced saturation] /[total % pore volume available for How]
Sew = (Sw — Swt} I (1 — Swt — Sm — Set) (5.25)
SeG = (SG-SGt)/(l-Swt-SNt-SGt) (5.26)
Definition 6. Se = [reduced saturation] /[total % pore volume available for phase
now]
Sew = (Sw- Swt] I (I -Swt] (5.27)
SeG = (SG-SGt}/(l-SGt) (5.28)
where again the trapped quantities are computed from equation 5.7.
With respect to quantifying the intermediate- wetting NAPL-phase, the appro-
priate Mualem model is given by Parker et al. (1987):
,~ 9Qs
(5.29)
L ylu) -U
where £ is a pore connectivity parameter for the NAPL-phase, and S£N is the
normalized or effective NAPL saturation which, when applied to the krj\f(Sw, SG)
functional, is assumed to have two admissible definitions:
Definition 7. Sefj = [reduced saturation] /[total % pore volume available for How]
SeN = (SN - SNt) I (I - Swt - SNt - SGt) (5.30)
Definition 8. SCN = [reduced saturation] /[total % pore volume available for
NAPL How]
SeN = (SN - SNt) I (1 - SNt) (5.31)
The term in brackets in equation 5.29 represents the ratio of the pores that are
contributing to NAPL flow at a particular combination of effective saturations,
g(Sewi SeTw}i to all pores contributing to water flow, g(0, 1). The parameter S£TW
is the normalized or effective total wetting phase saturation which, when applied
to the krff(Swi SG} functional, is assumed to have three admissible definitions:
68
-------
Definition 9. S£TW = [reduced saturation] /[total % pore volume available for
How]
Sw + SK — Swt — Sm
ew - ~ ~
(i — Owt — ONt — dGt
Definition 10. SeTw = [reduced saturation] /[total % pore volume available for
NAPL and water How]
_ Sw + SN — Swt — Sm , •>
deTw — - 7-j - ^ - ^ - s - (b.66)
(i — Owt ~ ONt)
Definition 11. SeTw = [reduced saturation] /[total % pore volume available for
NAPL How]
Sw + SN — Sm , .x
SeTw = - 7- - ^ - (5.34)
(1 - oNt)
where the trapped quantities are computed from equation 5.7.
Finally, let us substitute (1 — SeTn) for Sew in equation 5.29, where SeTn is the
normalized or effective total nonwetting phase saturation. The parameter SeTn is
assumed to have three admissible definitions:
Definition 12. SeTn = [reduced saturation] /[total % pore volume available for
How]
_ SK + SG — Sfft — Set
OeTn -
(i — Owt — &Nt
Definition 13. SeTn = [reduced saturation] /[total % pore volume available for
NAPL and gas How]
SN + SG — Sm — Sot ,- ocs
SeTn = - 7-j - ^ - ^^ - (5.36)
U — bet — ONt)
Definition 14. SeTn = [reduced saturation] /[total % pore volume available for
NAPL How]
_ SN + SG — Sm , -x
SeTn ~ (5'37)
Substitution of equation 5.18 into equation 5.29 leads to:
f ( r -i/im r i/irra>|2
krN(SW, SG) = (SeN^ { [l ~ (1 - SeTn)^ ~ [l ~ (SeTw)^] } (5.38)
where the effective saturations are as defined above.
In summary the three-phase k — S model is described by the following phase
relative permeability functionals:
69
-------
Water-phase - equation 5.19, with aSew and bSe\v denned either by equation 5.25
or 5.27;
Gas-phase - equation 5.23 with aSec and bSec denned either by equation 5.26 or
5.28;
NAPL-phase - equation 5.38 with Sej\f denned either by equation 5.30 or 5.31,
SeTw defined either by equation 5.32, 5.33 or 5.34, and SeTn defined either
by equation 5.35, 5.36 or 5.37.
The model includes hysteresis as the result of fluid entrapment effects through
the definition of the effective saturations. In addition, it utilizes data from the
saturation-pressure sub-model, the parameter TO, and includes three fitting para-
meters, £, £ and (p. The choice of effective saturation definition provides additional
flexibility. Finally, it reduces to the appropriate two-phase model when appropri-
ate.
5.4. MODEL IMPLEMENTATION
The hysteretic k — S — P model described in this Section is effectively an alteration
of the van Genuchten (1980) and Luckner et al. (1989) empirical models (subse-
quently referred to as the VG model). The alterations include a set of predefined
empirical relationships which describe hysteresis due to capillary and fluid en-
trapment effects. The validity of the resulting model with respect to representing
the hydraulic properties of granular soils must be a consideration. In addition, in
order to be amenable for implementation in a numerical simulator, the resulting
model must generate k — S — P functionals which are well behaved over the full
range of saturations. For example, recall that the formulation requires the S — P
functionals to be continuous in time and space, but that if the model defined by
equation 5.13 were to be used directly, at Se = 0, hc = oo and dhc/dSe = oo, and
at Se = 1, hc = 0 and dhc/dSe = oo. Therefore, additional computational steps
must be taken to ensure that the functionals have finite slope, and that numeri-
cally pathologic scanning curves are not generated. Finally, a necessary attribute
of the time-discrete S — P model is that it generate a unique set of curves as the
time step becomes small, and for mass balance considerations, that the functionals
be continuous. These issues are addressed in the subsections below.
70
-------
5.4.1. Phase Entrapment and Release
One important difference between the VG model and the current model has to
do with the definition of effective saturation. Consider the case of describing
the primary drainage curve (PDC). The PDC is defined by both models using
equations 5.10 and 5.11. The difference between the VG model and the current
model is in the definition of 5V, where, according to the VG model, Sr = Swr
is a constant, and according to the current model, Sr = Swt(Sw) is a variable
function of Sw as defined by equation 5.7. The functional form of Swt(Sw} is
dependent on the blending parameter, e, which governs how fast a phase becomes
entrapped during drainage flow conditions or released from entrapment during
imbibition flow conditions. Note that with respect to the PDC, for e = 0, the
current model reduces to the VG model, and for e > 0, the current model becomes
a modified form of the VG model, where as e gets larger, water entrapment is
delayed with respect to the change in Sw- In Figure 5.6 we consider the effect
that this modification has on the k — S — P functionals. Let us assume that some
soil moisture retention data was fit to the VG model, where the best-fit is given
by the following parameter set: ad = 0.02, a, = 0.04, r\ = 6.5, Swr = 0.17, and
Scr = 0.20. The plots show primary drainage (PDC) and main imbibition (MIC)
for S — P and k — S. The dashed lines represent the VG model (note, also the
current model with e = 0). The solid lines represent the current model with e = I
(a linear entrapment model). The general effect in employing the entrapment
model is that, relative to the VG model, the capillary pressure is under-predicted
during drainage events and over-predicted during imbibition events. In addition,
phase relative permeability is over-predicted during phase drainage. One can
avoid these discrepancies by fitting the parameters defining the present model to
experimental data directly.
5.4.2. S-P Curve Pressure Scale Transition
As discussed in Section 5.2.2, a reversal in flow direction can result in a dramatic
change in slope of the S — P functional as the pressure scale changes due to
capillary and contact angle effects. The change in pressure scale is accounted for
by the curve-type-dependent parameter &(/). Because a discontinuity in slope can
lead to convergence and stability problems in the numerical model, a blending
rule (defined in Table 5.2) is included in the S — P model. This rule governs how
fast a(j) changes as a function of saturation after a reversal occurs. The slower
changes over the span of the new curve, the smaller the discontinuity in slope
71
-------
100
60
MIC
0.5
Sw
0.4
0.2
0.5
Sw
0.8
Imbibition
0.5
Sw
Figure 5.6: An illustration of the effect that the entrapment model blending
parameter e has on the shape of the k — S — P functionals. For the solid curves,
e = 1, and for the dashed curves e = 0.
will be at the flow reversal point.
For example, consider the displacement process depicted in Figure 5.7 where
a flow reversal along a PDC occurs at Sw = 0.40 and an SIC is spawned. The
S — P model is defined by the following parameter set: e = 1, a^ = 0.02, a>i = 0.04,
1] = 6.5, Swr = 0.17, and Scr = 0.20. The SIC described by the dashed curve was
generated with (3 = 0, where the transition is instantaneous and the discontinuity
is a maximum. The SIC described by the solid curve was generated with (3 = 0.2,
where the transition is delayed thereby reducing the discontinuity. Numerical
experiments indicate that (3 = 0.2 yields a sufficiently smooth transition function.
5.4.3. S-P Curve Restriction Parameters
Four auxiliary computational procedures are included in the S—P model to ensure
that the functionals are well behaved over the full range of saturations for which
they are defined.
1. The S — P functionals are approximated by linear extrapolation near their
endpoints to ensure finite values of capillary pressure and slope. Using
equation 5.10 as the base-model, the hc(Se) functional which includes the
extrapolation computation becomes:
72
-------
100
80
60
40
20
PDC
0.5
Sw
Figure 5.7: The effect of using the blending rule described in Table 5.2 to de-
fine (i - s^
where S^r and S^s are critical effective saturations at the residual end and at
the saturated end, respectively, at which point the linearization takes place.
Values for S^r and S^s are of order 0.001.
2. A pathologic S — P functional will be generated if the denominator defining
effective saturation in equation 5.14 is too small. Therefore, a constraint is
included which limits how small the denominator can be, i. e.:
Ss(f) - Sr(f) > span
where span is the specified tolerance. Numerical experiments indicate
span > 0.01 is appropriate. For those S — P functionals where parameter
update would produce a 'span' in violation of this constraint, the current
functional is retained.
(5.39)
73
-------
3. A pathologic scanning curve will be generated if it is spawned too close to
the parent-curve endpoints. Therefore, a constraint is imposed which denies
scanning-curve generation unless the effective saturation of the parent-curve
is between the limits:
e < Se(f) <(l-e)
where e is the specified tolerance. Numerical experiments indicate e > 0.01
and e > span is appropriate.
4. To eliminate premature reversals due to numerical irregularities, a tolerance
in saturation change is defined:
rro+l cmin ^
•-V ~~ ^W > rd
omax rro+1 -- ^
&W °W ^ 'i
where S^,-11 is the minimum wetting phase saturation recorded for the current
drainage curve, S^r3* is the maximum wetting phase saturation recorded for
the current imbibition curve, r& is the tolerance to indicate a reversal from
drainage to imbibition, and r,i is the tolerance to indicate a reversal from
imbibition to drainage. Numerical experiments indicate r values of order
0.001 is appropriate.
5.4.4. Mass Balance and Consistency
By choosing the primary flow variables to be saturation and pressure, the math-
ematical formulation requires that the relationship between saturation and capil-
lary pressure is both continuous and approaches a unique path as the time step
is reduced. To show how the model addresses these issues consider the following
example flow scenario presented in Figure 5.8. A water drainage process has taken
the S — P relationship along the PDC to point 1, at which time a reversal to
imbibition is indicated. An SIC curve is generated (the dashed curve) by forcing
it to pass through the current point 1 and 1 — Sat (the point 1') where Sat is the
entrapped nonwetting phase saturation, currently zero. Over the time step the
S — P model parameters are held fixed and imbibition has progressed to point 2
on the dashed curve. To prepare for the next time step, S — P model parameter
update begins with the definition of phase entrapment, which is computed from
equation 5.7 to be at point 2'. If the remaining parameters were updated by
fitting the curve between points 1 and 2', the dotted line would result, thereby
creating a jump in capillary pressure (a drop in Pc from point 2 to 2*. To avoid
74
-------
Figure 5.8: Graphical representation of how the S—P functional retains continuity
over a time step. Initial drainage is from point 1' to 1. Then a reversal is indicated,
and the SICi is generated (the dashed curve between points 1 and 1'). Over the
time step Sw increases to the point 2 along the dashed curve. Model parameters
are updated to fit the SICi between points 2 (the current point) and 2' (the
updated trapped phase condition, 1 — Set), thereby generating the thin solid
curve.
this, the parameters are computed such that the new S — P curve (the thin solid
curve) passes through points 2 and 2', therefore preserving continuity.
Finally, one can see that as the change in saturation over a time step becomes
small, the discontinuity becomes small and the scanning curve approaches a unique
path.
5.5. CAPILLARY PRESSURE SCALING
Leverett (1941) reasoned that the capillary pressure between two immiscible flu-
ids should depend on the porosity, the interfacial tension between the immiscible
phases, and some sort of mean pore radius. Assuming that the ratio of the perme-
ability to porosity is proportional to the square of the mean pore radius, Leverett
75
-------
defined the dimensionless function of saturation called the J — function. For the
three-phase system considered herein, a J — function can be written for each
capillary pressure:
1/2
J(S) =
-------
was measur
ed, i. e.:
2-phase
system measured
G-W
N-W
G-N
Pc =
PcGW
PcNW
PcGN
a* =
CTGW
&NW
&GN
— (7jYi4/ -|- &GN
=
-------
6. NUMERICAL MODEL DEVELOPMENT
The five mass balance equations (4.1 [(/,, a) = (n, W) and (n,G*)], 4.8, 4.9 and
4.10) coupled with the constitutive and thermodynamic conditions presented in
Chapters 4 and 5 provide a complete description of the mathematical system when
proper initial and boundary conditions are imposed. The system of equations is
represented by the following set of primary variables:
(Pw/(x, t), SW(x, t), SW(x, t), pf (x, t), p£(x, t)}
where (x,t) G il x (0,T), represents the three-dimensional Cartesian spatial co-
ordinates of the domain and temporal scale, respectively.
This section describes the salient features of the numerical model, and it is
organized as follows. First, the governing equations are rewritten in the form
which is used in the numerical model. Then the numerical discretization and
solution strategy used for the set of nonlinear coupled equations is discussed. The
numerical model features a sequential solution of the governing equations and an
implicit-in-time collocation finite-element discretization. Next, in order to close
the system mathematically, the scope of applicable initial and boundary conditions
is presented. In addition, in order to improve the robustness of the algorithm,
a set of numerical tools, which are designed to ensure that the physical problem
and the discrete problem are compatible, is discussed. Finally, the computational
steps associated with the numerical model are summarized.
6.1. FINAL FORM OF THE BALANCE EQUATIONS
By substituting the relevant functional definitions and simplifying assumptions
and using the chain rule of differentiation, the balance equations are rewritten in
the form which will be discretized and solved numerically.
78
-------
6.1.1. Water-Phase Flow and Contaminant Transport
The water-phase flow and transport equations (4.8 and 4.1 [(L, OL) = (n, W)])
are coupled through phase pressure and saturation variables, the composition-
dependence of phase density and viscosity, and the mass exchange terms. Em-
ploying the chain rule of differentiation and rearranging, equation 4.8 is rewritten
as:
__
ti/
P
(6.1)
where qw/ = e5Vi/vw/ is the water flux vector, and pw is the water density as-
sociated with Qw . Substituting equation 4.17 into 6.1 for pw and rearranging
yields:
V • aw = Qw + —- /W /W
p
Wr
p
Wr
where rw = (l — pWr /p
w
Nr
and
(6.2)
is the mass concentration associated with
The dissolved NAPL transport equation (4.1 [(i,a) = (n,W)]) is rewritten
in an analogous way. Using the approximation of the dispersion term defined
in equation 4.24, and making the simplifying assumption that that (pw/r/pw/) is
approximately a constant, the dissolved NAPL transport equation becomes:
dt
P?
dt
(6.3)
Given equation 6.3, the water flow equation 6.2 can be simplified algebraically
by substituting the right hand side of equation 6.3 for the last term in brackets
79
-------
on the right hand side of equation 6.2 to yield upon rearranging:
pW7 _ pG _ f^S
i w/ w/ /W /W
pw
(6.4)
This form of the water-phase flow equation highlights the fact that volumetric
changes in the water-phase will occur as the result of mass exchange processes,
and, when pw/ is a function of p% (i. e., rw/ ^ 0), as the result of NAPL decay
processes and dispersive fluxes.
Algebraic manipulation can also be used to simplify the dissolved NAPL trans-
port equation (6.3) by substituting the right hand side of equation 6.1 for the term
in parenthesis on the left hand side of equation 6.3. This operation yields upon
rearranging:
- V
(Ow _ Ow\
v ' n '
- -
where pn is a known concentration when Qw/ represents a source (QW positive),
and ~pn = pjf when Q^ represents a sink (Qw/ negative) . This equation is further
manipulated by substituting equations 4.28, 4.33 and 4.37 for the mass exchange
terms, and rearranging i. e.,
[cf (1% - P*'} - c*/w
(6.6)
80
-------
This form of the contaminant transport equation highlights the fact that the
relationship between the divergence of the water flux and the mass exchange
processes is nonlinear, and that NAPL species transport in the water and gas
phases is coupled through the exchange term E^/w.
6.1.2. NAPL-Phase Flow
Employing the chain rule of differentiation and rearranging, the NAPL-phase flow
equation (4.8) becomes:
+ V.MV-!M) (6.7)
where pN = pNr is a constant.
6.1.3. Gas-Phase Flow and Contaminant Transport
As with the water-phase, the gas-phase flow and transport equations (4.10 and 4.1
[(L, a) = (n, G*)]) are coupled through phase pressure and saturation variables, the
composition-dependence of phase density and viscosity, and the mass exchange
terms. Following an analogous series of steps as was used for the water-based
equations, the equivalent form of the gas-phase flow equation (4.10) and the NAPL
vapor transport equation (4.1 [(/,, a) = (n, G*)]) are written as follows, respectively:
dt
-r
G J £SGKnp[
pGr
(6.8)
[eSGVG . VpG
pG ' eSGKGPG
(6.9)
81
-------
where rG = (l — pGr/pNr}, Pn is a known concentration when QG represents a
source (QG positive), and J>n = pG when QG represents a sink (QG negative),
and (pGr IpG] is assumed to be approximately a constant. This form of the gas-
based balance equations highlights the nonlinear relationship between the flow
and transport parameters.
6.2. SEQUENTIAL, ITERATIVE SOLUTION
This section details the development of the time-discrete form of the balance
equations, and the iteration method used to linearize the resulting equations.
The phase mass balance equations (6.4, 6.7 and 6.8) define the distribution of
the phases given phase composition. The constituent balance equations (6.5 and
6.9) define the distribution of the NAPL species in the water- and gas-phases,
given the phase distribution and velocity field. In an attempt to minimize the
overall computational effort, an iterative scheme is adopted whereby the balance
equations are solved by sequentially lagging certain dependent variables, such that
the balance equations become uncoupled. This strategy minimizes the size of the
matrix equation which is required to be solved at any one solution step.
The time scale, t, is discretized into a series of finite intervals, with each interval
defining a time step, A t = (tk+l — tfc), where tk is the current time at which the
solution is known and tk+l is the advanced time at which the solution is sought. All
time-dependent parameters are indexed with respect to where they are evaluated
in time, where those terms dated at the current, known, time level are indexed
by the superscript (k) and those terms at the advanced, unknown time level are
indexed by the superscript (k + 1). In addition, let us indicate those variables
which are temporarily fixed in time by the superscript, (k + 1)*, indicating that
the variable is dated at the last available solution step. The time-discrete system
is summarized in the following subsections.
6.2.1. Time-Discrete/Linearized Form of the Flow Equations
The flow equations 6.4, 6.7 and 6.8 are solved using an algorithm based on the
total velocity formulation of Spillette et al. (1973), which requires algebraic ma-
nipulation to generate an elliptic-like pressure equation and two parabolic-like
saturation equations. The key computational feature in applying this method is
that the total velocity is calculated after the solution of the pressure equation,
and this velocity is then used in the saturation equations. This sequential solution
82
-------
strategy has two major attributes (Peaceman, 1977):
• The total velocity obtained after the solution of the pressure equation will
be the same regardless of whether the capillary pressures are evaluated ex-
plicitly or implicitly.
• The saturation equations, written in fractional flow form, are not susceptible
to erroneous changes in a phase's saturation when that phase's mobility is
zero.
The pressure equation is derived by adding the phase mass balance equations
6.4, 6.7 and 6.8 to eliminate the saturation time derivative, yielding:
(6.10)
= (Qw + Q1
ESn X
/ \ (fc+1)*
where qT = qw/ + q^ + qG is the total fluid flux, ( E^,w 1 defines adsorptive
mass transfer (the time-discrete form of equation 4.37):
W(k+l)* _ W(k)
(6.11)
and (rjf) and (FG) define respective changes in p w and p G due to the
presence of NAPL species:
(6.12)
-DG * Vpg]
pGr
(6.13)
Choosing Pw/ as the dependent variable, the phase fluxes used in equation 6.10
are defined from equations 4.23 as:
83
-------
N _ _A(fc+l)* v pW(k+l) ,
q — A
qG = -A* V P+ + P + P - 7+* V*
(6.14)
where AQ = kkra/ [ia (a = W, N, G), is the a— phase mobility scalar.
Given its time-discrete form, the pressure equation 6.10 is a linear elliptic
equation in terms of Pw , and thus it can be solved directly. Let us denote its
solution by pw'(fc+1)*. The total flux is now computed from equations 6.14 with
Pw(k+i) approximated by pw/(fc+1)*. Let us denote this dated variable by qT(fc+1)*.
The saturation equations are chosen to be the water phase balance equation
(6.4) and the gas phase balance equation (6.8). Through algebraic manipulation
of equations 6.14, qw/ and qGcan be rewritten in terms of qT as follows:
+fw {XN [VPcNW+ A fWNVz] + XG [V (PcNW + PcGN) + A fWGVz]}
(6.15)
-/G {XN [VPcGN+ A ~fNGVz] + Xw [V (PcNW + PcGN) + A 7H/C
(6.16)
where
fa=Xa/ (Xw + XN + AG), a = W, G (6.17)
is the a—phase fractional flow function, and A 7Q/3 = 7" — 7^.
With these definitions, consider the following time-discrete forms of the water-
and gas-phase balance equations.
Water phase transport:
A t
prw _ pG _ i?s
-*-Jt-> -*-J~^ lixr -*-J~^ /T.
(6.18)
84
-------
where
qw(fc+i) = jW(fc+i)qT(fc+i)*
G\(fc+1)*
W /1G\\+1* ffilQA
n ~ Pn ) / (6.19)
(fc+1)* r (jf i -i U
-Bf/^ and (r^l/) are denned in equations 6.11 and 6.12,
respectively.
Gas phase transport:
_
Tw + V . qG(fc+l) = gG(fc)
/ \ 6
(fc+1)*
(6.20)
where ^T^ = (1 - SG),
qG(fc+i) = jG(fc+i)qT(fc+i)*
- (fG\ "l(fc+1) [vP(fc+1)+ A
(J AN) [ v -^cGJV ~r A
- (fG\ }(k+1} [V Tp(fc+1) 4-
{J AW) [V \cNW +
-------
• Each equation is linearized by a Picard-type iteration, wherein all functions
of Sw and STW m equations 6.18 and 6.20 are lagged an iteration.
• The uncoupled linearized equations are solved concurrently and in parallel if
possible. That is, both equations are solved before the saturation-dependent
terms are updated.
• Iteration continues until convergence is obtained.
This sequential solution strategy is chosen for the following reasons. First, it
minimizes the computational effort required per solution step. Second, sequential
solution requires that certain primary variables be temporarily lagged in order
to de-couple the system as each balance equation is solved. In particular to our
system, STW is lagged when the water balance is solved for Sw, and Sw is lagged
when the gas balance is solved for STW- This is equivalent to saying that when
solving either balance equation, a change in the dependent variable is due to a
change in SN- Since the NAPL has intermediate wettability between the water
and gas, a symmetry is created in the form of the water and gas balance equations.
This symmetry is shown in the phase velocity equations 6.15 and 6.16, where the
last term in each equation is equivalent, and it represents water-gas interaction
(recall that PCGW = PCNW + PCGN)- It is appropriate then, that this symmetry
be preserved in the evaluation of the water and gas balance equations so that
erroneous NAPL saturations are not created as the result of imperfect iterative
convergence.
To denote iteration level the index 'TO' is used, and terms dated at time level
(k + 1) are double indexed such that (k + l,m) indicates evaluation at the last
known iteration solution, and (k + I,TO + 1) indicates the unknown quantity for
which a solution is desired.
Using this notation, the linearized form of the water-phase transport equation
(6.18) is:
Q(K-f-±^irL-f-±j O\K)
C^W ~ °W , T7 . (fW(k+l,m)QT(k+l)*\
At Vl/ ^ )
+ V • < (/ XM) ' VP-.ATTI/ + A
+V • fw\' V Pm + P'm + A
\ (fc+l)x
«7 _ EG - Es V
i -^n/W ^n/W)
^
86
-------
Note that since Sw is chosen as the dependent variable, VPcjvw m must be ex-
pressed in terms of Sw >m . Consider the following strategy. First, approximate
using a first order Taylor series expansion:
(6.22)
i c n\ K~T±.iii>~f-± i I n\ K~Ti..iii>~Ti. i n\ K~Ti-.ni, I \ • i i • , , • • in i
where oow = ( Sw — Sw I is the iterative increment. Second,
to minimize chaining requirements, simplify the second term on the right hand
side of equation 6.22 as:
(fc+l,m+l)
cNW
(k+l,m)
cNW
'dPf
cNW
dSw
(6.23)
Note that as the solution converges, the iterative increment goes to zero, and
equation 6.23 becomes exact. Following a similar chaining procedure, the gas-
water pressure derivative is written in terms of Sw by using the capillary pressure
scaling rule (4.14), i. e.:
V P
cNW
p(fc+l,m+l)\ _
^cGN J —
V P
D(fc+l,m)
cGN
dSw
(6.24)
This mixed parameter iterative scheme is used in lieu of more traditional pro-
cedures in order to correctly accommodate capillary pressure scaling with respect
to heterogeneous fluid and soil distributions. Finally, substituting equations 6.23
and 6.24 into 6.21, and using SSW >m as the unknown quantity, one obtains
the following linearized version of the water phase transport equation:
'(fc+l,m) T(fc+l)*
A t
W\ \(k+l,m) f^pCfc+l.m)
fW \ XV^TJ-,'"/ I V7plrc-
/ AN) I v/dvw
87
-------
{(/"%)
(fc+l,m)
(fc+1)*
(6.25)
Following an analogous procedure, the linearized form of the gas-phase trans-
r(fc+"
'Tw
port equation (6.20) is written in terms of the iterative increment SS^ 'm
rt(fc+l,m+l)
oTw
A t
V • / ( fG \ \ (fc+1-m) [ v p(fc+1'm)
-V • <^(J AN) [^^cGN
(fc+l,m) p(fc+l,m)
(fc+1)*
(6.26)
After equations 6.25 and 6.26 are solved, the required variables are updated:
Q(fc+l,m+l) _ Q(fc+l,m+l) c
~
_
Tw — Tw — Tw
Iterations continue until SSW 'm and SSj,w 'm are sufficiently small.
6.2.2. Time-Discrete/Linearized Form of the Transport Equations
In general, the fluid phase properties are only slightly affected by the contami-
nant species, thus the flow and contaminant transport equations are only weakly
coupled. Therefore, the solution to the flow problem is used to define the flow
-------
variables used in the transport equations (6.6 and 6.9). Indexing the most recent
flow solution by the superscript (k + 1)*, the following flow variables are defined:
S^+1)*, 4fc+1)*. qw/(fc+1)*, qG(fc+1)*, Dw/(fc+1)* and DG(fc+1)*.
From equations 6.6 and 6.9 it is clear that the transport equations are coupled
and nonlinear as the result of mass exchange processes. In order to solve the
equations, consider the following iterative sequential solution algorithm:
Wfc+l)
Equation 6.6 is solved for pn and equation 6.9 is solved for pn
G(fc+l)
Each equation is linearized by a Picard-type iteration, where all functions
f
of pn
j G(fc+l) , , ., ,.
and pn are lagged an iteration.
• The uncoupled linearized equations are solved concurrently and in parallel
if possible. That is, both equations are solved before the concentration-
dependent terms are updated.
• Iteration continues until convergence is obtained.
Rewriting equations 6.6 and 6.9 in a linearized, time-discrete form consistent
with iterative sequential solution, one obtains:
Pb
Kd
W
Pn
W(k+l,
Pn
pNr
(fc+1
m)\
|
)
m+l) W
~ pn
At
-
pbKd
^
(fc)
Pjf
n rn
rn
(fc+l,m+l) _ G(fc+l,m)*
rn
(6.27)
eS,
G(fc+l,m+l) G(fc)
(k)\ Pn - Pn
G
+q
C(«H-D*
(fc+l,m+l)
89
-------
(6.28)
i Wffc+l.m)* i Gffc+l.rra)* . • i j i r j i i i
where pn and pra represent a weighted average 01 the lagged vari-
ables, i. e.,
pa(fc+l,m)* = 9pa(k+l,m) + (I - 0) pa(k+l,m-l) ^ a = W, G (6.29)
and 9 is a weighting parameter, 0 < 9 < 2.
6.3. SPATIAL APPROXIMATION
The three-dimensional domain is discretized into a finite number of brick-like
elements (i. e., a regular rectangular mesh) with nodes located at element bound-
ary intersections. All spatially varying parameters are represented by a linear
combination of basis functions where nodal values can be interpolated into adja-
cent elements in a continuous manner. These parameters include the dependent
variables and functions of the dependent variables, and soil properties and other
physical constants which are allowed to vary spatially.
Consider first the dependent variables defining equations 6.10, 6.25, 6.26, 6.27
and 6.28: Pw/ , Sw, STW, pv% and pG, respectively. Since these variables require con-
tinuous second-order derivatives, the Cl continuous Hermite cubic interpolation
polynomials are chosen. They are defined in Appendix B.
For three-dimensional discretizations, a general function, tt(x, t) is approxi-
mated over each grid block as:
8
«(x, t] « S(x, t) = {U(t)}T {$(x)}t (6.30)
where u(x., t) is the spatially-discrete approximation to it(x, t). At each of the eight
i nodes in a grid block is defined a vector of eight time-dependent undetermined
90
-------
coefficients:
{U(t)tf = {U(t), Ux(t), UtV(t), Uz(t), Uxy(t), Uxz(t), U,yz(t), U
where the subscript [ , (•)] represents the partial derivative with respect to (•),
and a vector of eight space-dependent Hermite interpolation functions, {$(x)}i,
the form of which is detailed in Appendix B.
Each dependent variable is approximated spatially by an equation of the form
6.30, yielding for example, Pw/(x, t), SV(x, t), SV^x, t), pjf(x, £) and p^(x, t).
These approximations are then substituted into their respective balance equa-
tions yielding five equations in 5[8(/)] nodal unknowns (5 dependent variables, 8
undetermined coefficients per node, and I nodes).
The remaining parameters requiring spatial representation are either functions
of a dependent variable (e.g. capillary pressure and fluid mobility) or spatially
varying physical constants (e. g., soil properties e and k). With respect to capil-
lary pressures, because a mixed parameter iterative scheme (defined in equations
6.25 and 6.26) is employed where, as the iterative increment in saturation ap-
proaches zero, the space operator is defined by capillary pressure gradients, it is
necessary to use a Hermite representation for PCNW(SW] and PCGN(STW}- The
eight nodal coefficients are computed as functions of saturation, and, where ap-
plicable, the chain rule for differentiation is applied (see Appendix B for details).
The remaining coefficients are defined nodally using relevant data, and defined on
each mesh block in terms of the tensor product ordering of the one-dimensional
linear Lagrange polynomials.
The combination of using piecewise continuous soil properties and Hermite in-
terpolated pressure variables has the advantage of yielding a piecewise continuous
velocity field for use in the transport equations.
The nodal unknowns generated from the Hermite discretization are determined
by employing the collocation method, a method of weighted residuals where the
weighting function is the displaced Dirac delta function (see Frind and Finder
[1979] for a detailed discussion). As a result, the residual errors incurred by using
approximated dependent variables in the governing equations are driven to zero
at specified points in the domain, called collocation points. If the Gauss quadra-
ture points are chosen as the collocation points, the method is called orthogonal
collocation.
Using the general function u(x., t) as a surrogate for each approximated depen-
dent variable, a system of linear algebraic equations is generated for each balance
91
-------
equation by imposing the interpolation constraints:
w(xj-,t)=w(xj-,t), j = I,2,...,8(NEL)
where (xj) are the locations of the collocation points, and NEL is the number of
elements. This results in a matrix equation of order [8(NEL)] for each discretized
balance equation. Note that equation generation requires no formal integrations,
and therefore, collocation is computationally analogous to the finite difference
method in that equations are written at points in the domain.
6.4. IMPOSED CONDITIONS
The problem posed in the previous sections is an initial, boundary value problem,
and as such, requires initial and boundary conditions to mathematically close
the system. In addition the method used to impose external flux conditions is
detailed.
6.4.1. Initial Conditions
Initial conditions are specified such that saturation and composition are defined
at t = 0. The following variables require initial data:
Flow variables
SW(x,0) = SWo(x)
SW(x, 0) = SWo(x)
where the specified values SVo(x) and STIUO(X) define the initial k — S — P
functionals according the decision tree provided in Appendix C.
Contaminant transport variables In general the conditions are
In addition, by default, at any node where the initial conditions for satu-
ration indicate that NAPL is present, i. e., (STWQ — 0, the concen-
tration variables are set to their respective solubility limits, i. e., p^0 = ~p^
and p^Q = ~p~n-
92
-------
6.4.2. Boundary Conditions
A sufficient number of boundary conditions must be specified to augment the
collocation equations such that the number of equations equals the number of
unknowns. This is accomplished by specifying as known a group of Hermite coef-
ficients at each boundary node. The grouping is based on the particular boundary
condition-type imposed and the boundary node-type at which the condition is ap-
plied:
1. boundary condition-type: three types of boundary conditions can be imposed
at a given boundary node: Dirichlet (first-type), Neumann (second-type) or
mixed (third-type)1. For Dirichlet data the function value is specified at the
node. For Neuman data the derivative normal to a specific boundary plane
associated with the node is specified. For mixed data, a linear combination
of Dirichlet and Neumann data is specified at the node.
2. boundary node-type: typing is based on the orientation and number of
plane(s) in which the node lies (x-y-plane, y-z plane, and x-z-plane). Along
with the user-defined coefficients identified in (1), the simulator specifies,
by default, those Hermite coefficients which are derivatives of the boundary
data in the direction(s) tangential to the boundary plane(s) associated with
the node, where the node-type is used to identify which coefficients apply.
For example, if Dirichlet data is to be applied to a node which is in an x-y
boundary-plane, then the user specifies the Hermite coefficient U (the Dirichlet
data) and the simulator applies appropriate default conditions on the coefficients
£/a;, U^y and U}Xy (the tangential derivatives of U in the x-y boundary-plane). The
details of which degrees of freedom are associated with which boundary condition-
and node-types are provided in Appendix B. Note that this discussion is specific
to the three-dimensional case. An analogous procedure is applied for the two-
dimensional case where the two-dimensional problem is equivalent to one plane of
nodes in the three-dimensional problem.
This subsection provides a description of the specific boundary conditions used
in the simulator. Imposition of point source and sink terms is discussed in the
following subsection.
JThe current version of the simulator only considers mixed boundary data for the NAPL
species in the gas-phase transport equation.
93
-------
Flow Variables
Two types of flow boundary conditions are considered in this simulator: Dirichlet
and Neumann. A boundary which is considered closed to mass flux is modeled
using a Neumann condition. A boundary which is considered open to mass flux is
modeled using either a Dirichlet condition, or using the combination of a closed
boundary Neumann condition and a point source or sink term representing the
specified flux.
Consider a specification system which identifies a series of boundary-type cases
defined on the basis of the three-phase flow condition. Each case is defined by the
following criteria:
• for each of the phases, specify whether the boundary is open or closed to
flow;
• a boundary which is considered open to flow of a given phase is defined by
specifying either phase pressure, saturation or flux.
Table 6.1 provides a summary of the specification system used in the simulator
for which five cases have been identified. This Table includes information on the
flow status of the boundary, the Dirichlet variable(s) specified where applicable,
and what kind of boundary condition results for each dependent variable. In addi-
tion, Table 6.1 provides information on whether the resulting boundary conditions
are linear or nonlinear with respect to the flow variables. Note that any condition
which involves gas- or NAPL-phase pressure related information is a nonlinear
function of the solution. This is a direct result of the choice in primary variables,
Pw/, Sw and STW (recall equations 4.16).
With respect to this simulator, all the boundary nodes are assumed to be no
flow (Case 1) unless otherwise specified (i. e., no flow is the default condition). In
addition, the nonlinear terms in each boundary condition equation are dated at
the old time level (indexed by the superscript k). As a result, boundary conditions
are in general constant over a time step and time varying (updated after each time
step)
The specific boundary equations corresponding to each of the cases listed in
Table 6.1 are presented below. Note that the superscript (k + I) represents infor-
mation dated at the advanced time level , and that the notation (•) „ indicates
the spatial derivative of the quantity (•) normal to the boundary.
94
-------
case
1
2
3
4
5
boundary
open (o)
closed (x)
to flow
of phase
W
X
X
X
o
o
N
X
X
o
X
o
G
X
o
X
X
o
variables
specified
Pu
PN
Pw
PW sw sTw
conditions for the
primary variables,
1 = Dirichlet,
2 = Ncunian
(L = linear,
NL = nonlinear)
Pw
2(L)
2(L)
2(L)
1(L)
1(L)
S\V
2(NL)
2(NL)
2(NL)
2(NL)
1(L)
STW
2(NL)
2(NL)
2(NL)
2(NL)
1(L)
comment
the default condition
replace with source term
replace with source term
Table 6.1: SUMMARY OF THE SYSTEM USED TO DEFINE BOUNDARY
CONDITIONS FOR THE FLOW VARIABLES
Case 1 The boundary is closed to water, NAPL and gas flow. This is the default
condition used in the simulator, i. e., if no other flow condition is imposed at
a given boundary node, then the simulator assumes it is a Case 1 boundary
node. The boundary equations are derived by setting equations 6.14 equal
to zero and solving for the respective variables, i. e.,
w(k+l)
(Pw}n
i N w\
(P -P )
G pf\
P -P
(6.31)
Case 2 The boundary is made open to gas flow by specifying PG = PGQ. Con-
sider an algorithm which approximates the gas-phase pressure Dirichlet
boundary condition by generating an appropriate source or sink term, the
rate of which results in PG K?> PGO at the boundary node. The algorithm
takes advantage of the fact that a sequential total flow formulation is uti-
lized. First equation 6.10 is solved with the following boundary data which
satisfies the stated condition:
pw =
95
-------
Second, qT(fc+1)* is computed from the solution. Before solving the satu-
ration transport equations, the results from the total flow solution at Case
2 boundary nodes are interpreted as a combination of a no flow boundary
with conditions defined in 6.31 and an active point source or sink with the
rate defined as:
for inflow conditions:
for outflow conditions:
QG =
Qw =
Qa = f
q^* (area}
QN = 0
CK^fC~j~ i j* g~* \ /I firppf] \ (~\, T/T/' ]\T (-1
**-\n ILt/OLtK Ct V V 5 1 V ^ O"^
where q« * is the normal component of the total flux vector at the
boundary node, (area) = l/4(element area about the node in the 2-D plane
through which q^ is defined) is the area associated with the boundary node,
and for outflow conditions the rate is allocated based on a-phase fractional
flow conditions, Ja(fc+1)* (defined in equation 6.17). Note that if inflow con-
ditions prevail at the node, then the boundary is open only to the gas-phase,
but if outflow conditions are indicated, then the rate is allocated as indi-
cated.
Case 3 The boundary is made open to NAPL flow by specifying PN = PN0.
Following the same logic as for Case 2, solve the total flow equation with
the following boundary data which satisfies the stated condition:
r>W _ r>NQ p
f — f — J^cN
Before solving the saturation transport equations the results from the total
flow solution at Case 3 boundary nodes are interpreted as a combination
of a no flow boundary with conditions defined in 6.31 and an active point
source or sink with rate defined as:
for inflow conditions:
for outflow conditions:
QW'=q™Gfc+1)*(area)
g« = f<*(k+i)*(£(>°+V*(area^ a = W,N,G
Case 4 The boundary is made open to only the water-phase by specifying Pw/ =
Plv ° (a linear Dirichlet condition) in combination with the no flow conditions
for the NAPL- and gas-phases as defined in Case 1.
96
-------
Case 5 The boundary is made open to all three phases by specifying the following
linear Dirichlet data:
pW
S\y
OTU
pWQ
Transport Variables
NAPL Species in the Water-Phase The assumption is made herein that,
with respect to water-phase contaminant transport, dispersion is small compared
to advection, and therefore, the applicable conditions reduce to Dirichlet and
homogeneous Neumann.
The default conditions used in the simulator are as follows:
if at a given boundary node, the boundary condition for flow is Case 5, and
SN > 0, then Dirichlet data are specified so that p% = ~p%, otherwise a
homogeneous Neumann condition is applied, i. e., -rp- = 0.
Given the default conditions for p^ at each boundary node, the user needs to
specify only Dirichlet conditions, i. e., pjf(x, t) = p^l0.
NAPL Species in the Gas-Phase Since diffusion in the gas-phase is a signif-
icant transport mechanism a mixed-type boundary condition is included to model
gas-phase mass transport across a boundary which represents the ground surface.
If the atmosphere is considered an infinite sink (i. e., p°, . =0) and the mass
JT \ i i n atm '
flux is defined by diffusion across a stagnant boundary layer of thickness 8 then
the mixed condition is written as:
DG
(pGcf - eSGVG • VpG) • n = -fpG (6.32)
where n is the inwardly directed unit vector normal to the boundary. Note that
as 8 goes to zero the condition becomes homogeneous Dirichlet, and that as 8 gets
large the condition becomes homogeneous Neumann.
The default conditions used in the simulator are as follows:
if at a given boundary node, the boundary condition for flow is Case 5, and
5V > 0, then Dirichlet data are specified so that pG = p^, otherwise a
homogeneous Neumann condition is applied, i. e., -j^ = 0.
97
-------
Given the default conditions for p^ at each boundary node, the user has the choice
of specifying either a Dirichlet condition by denning p^(x, t) = p^0, or a mixed
condition by denning the value of 6.
6.4.3. External Flux Conditions
External flux conditions are point source or sink terms represented in equations
6.4, 6.7, 6.8, 6.3 and 6.28 as Qa [1/T]. The user specifies Qa at a node in units of
[L3/T] and the simulator normalizes that flow rate by the volume associated with
the node (1/8 of the volume of the elements surrounding the node). By convention
positive Qa represents a source and negative Qa represents a sink. Numerically,
Qa is defined at a node and interpolated using linear Lagrange basis functions.
Source forcing terms are applied with known composition, and are thus incor-
porated directly into the mass balance equations. For sink forcing terms, however,
fluid composition is not known, and the distribution of the phases must be deter-
mined in order to define the phase flow rates. The procedure used in the simulator
to allocate flow rates to each of the phases as the solution evolves through time is
based on fractional flow defined by phase mobilities. During any given time step,
(k) to (k + 1), the outflow phase fluxes will be allocated using information from
the last time step (explicit formulation). Given that total flux is known at the
current time level, QT = Qvv + QN + QG = (QT) , the phase fluxes are defined
as:
(QQ)(fc) = p« (QT) (fc) , a = w, N, G (6.33)
6.5. DIAGNOSTIC TOOLS
Solving for the temporal and spatial distribution of NAPL contaminants results
in a problem with moving fronts. In addition, because the model considers the
emplacement of residual NAPL in the vadose- and saturated-zones, it needs to
consider problems which evolve to quasi-static flow conditions, with the distrib-
ution of phase saturations governed by system forcing, soil and fluid properties,
and the nonlinear flow parameters. In a heterogeneous distribution of physical
properties and forcing terms, the static solutions involve different length-scale
boundary layers where phases appear and disappear. Therefore, any numerical
model derived to solve this type of physical problem must be able to resolve a
complex distribution of fluid saturations. As such, it is imperative that the spatial
and temporal discretizations utilized be compatible with the flow parameters so
98
-------
that discretization errors do not adversely affect the quality of the solution, and
that the case of phase appearance and disappearance can be accommodated.
In light of this discussion a series of diagnostic numerical tools have been devel-
oped which are designed to ensure this compatibility. These tools can be separated
into two categories: those which assess the spatial discretization and those which
assess the temporal discretization. The spatial discretization is assessed through
the use of a Peclet constraint. The temporal discretization is assessed through the
use of a Courant constraint and a dynamic time-stepping scheme based on the
number of iterations required for convergence on the nonlinearity.
Effectively, this section describes the diagnostic numerical tools designed to
make the simulator more robust.
6.5.1. Peclet Constraint
It is well known that application finite element numerical methods to immiscible
flow problems with sharp fronts will generate solutions which exhibit oscillations
when the spatial truncation errors are unable to correctly propagate the short
wave-length parts of the solution (Allen, 1983 and Mercer and Faust, 1976). Os-
cillations can lead to violation of the maximum principle, which results in non-
physical values of the primary variables and potential instabilities in evaluating
the nonlinear coefficients. In addition the combination of the inability of the ap-
proximation to propagate short wave-length oscillations with the sensitivity of the
frontal velocity with respect to those short wave-length components, can result in
convergence to the wrong solution (see Allen, 1983 for details).
Assuming that some physical capillarity exists for the problems under con-
sideration herein, the relationship between the flow parameters which define the
shape of the front and the grid spacing which defines the resolution scale, is quan-
tified by the dimensionless group of parameters called the Peclet number, Pe,
defined as:
p (aduection) (characteristic length scale) /r> o A\
^e ~ (diffusion) ! (V.04)
where in this case the advection parameters are dominated by gravity forces, the dif-
fusion parameters are dominated by capillary forces, and the characteristic length
scale is related to the local grid spacing. Large values of Pe indicate a sharp front
relative to the grid spacing, and small values indicate a diffused front relative
to the grid spacing. Consider the definition of a critical working Pe, called here
Pecrrt, for a particular finite element method. Pecr", determined from numer-
ical experiment by successively refining the mesh for a given problem, is that
99
-------
Pe where the numerical solution exhibits oscillations at or below some tolerance.
Therefore, given the flow parameter values, nonphysical numerical oscillations in
the finite element solution can be minimized by using a spatial discretization for
which Pe < Pecrit.
However, because mesh refinement is computationally expensive, it is in gen-
eral an impractical option. As a result, a second option is to alter the flow para-
meters given the grid spacing so that Pe < Pecrit. Specifically one adds artificial
diffusion sufficient to smooth the solution, and thus, eliminate the undesirable
truncation errors.
Adding artificial diffusion alters the physical problem, resulting in a smooth
but overly smeared front. Therefore, the term should maintain the following
attributes:
1. be applied only where needed, and in minimum magnitude;
2. be a function of the local grid spacing, flow parameters, and solution evolu-
tion;
3. not be a function of coordinate rotation; and
4. provide a convergent solution as the grid is refined.
Upstream weighting is a common class of methods used to add diffusion by
introducing a low-order spatial truncation error, the size of which is meaningful
only near sharp fronts. Because the size of the non-physical error term introduced
is in general dependent on the orientation of the front with respect to the grid,
upstream weighting methods do not meet requirement (3.) for problems involving
large mobility-ratio displacements (the viscosity of the invading fluid is much less
than that of the displaced fluid, for example gas displacing NAPL).
The goal of this analysis is to define a physically-based diffusion term to be
added to the governing equations which meets the three attributes identified
above, and which is applicable for the general multiphase flow problem where
a non-zero physical capillarity is inherent in the system..
The differential equations which describe phase transport can be classified as
nonlinear parabolic advection-difiusion equations, with the advection term dom-
inated by gravity forces and the diffusion term defined by capillary forces. The
water- and gas-phase saturation equations can be written in advection-difiusion
form by using the fractional flow formulation described in Section 6.2.1. For clarity
in presentation and without loss in generality, consider the physically simplified
100
-------
case of constant fluid density, isotropic porous medium, and no sources or sinks
of mass. For this case the water and gas equations are written as:
9SW
dt
fV . (fV + [(fw\N) A TWW + (/W/AG) A 7H/G] V*)
fV • [AWN (fwXN) VPcNW + AH/G (/W/AG) VPcGH/]
= 0 (6.35)
dt
+V • (/GqT - [(fGXN] A 7jVG + (/GAH/) A 7H/G] Vz
- [\NG (fGXN) VPcGN + AH/G (f°Xw) VPcow]
= 0 (6.36)
On the left-hand-side of equations 6.35 and 6.36 the second term represents the
advective component of mass transport which consists of two sub-components:
advection due to total fluid velocity, and advection due to gravitational forces,
and the third term represents the diffusive component of mass transport which
arises from capillary forces where PCGW = PCNW + PCGN, and the terms A;vw,
AJVG and AGW are new parameters which govern how much artificial physical
diffusion is added to the system. Specifically, Aa/3 > 1 (a(3 = WN, NG, WG),
and the case when AQ/3 = 1 represents no additional diffusion. The focus of this
development is to derive an algorithm to calculate the terms A;vw, AJVG and AGW
based on an appropriate definition of grid Peclet number for two- and three-phase
flow.
Section 6.2.1 provides the details of the sequential solution procedure, in-
cluding the time-discrete nonlinear iteration algorithm. For the purposes of this
discussion, note that the total flux vector (qT) is computed from the total flow
equation and is then considered a known variable for the solution of the saturation
equations. In addition, recall from Section 6.3 that all spatially-dependent para-
meters are defined nodally and interpolated into the elements using appropriate
interpolation functions. Therefore, a Peclet number can be defined at each node
i in the domain, i. e.,
_ max{|q
-------
where P&i is the Peclet number denned at node i. The numerator is a measure of
the advective component of transport at node i, and it is denned by the largest
component of the product of the advection vector and the grid spacing in each of
the s-directions (s = x, y, z). The denominator is a measure of the scalar diffusion
coefficient at node i.
To derive the expressions for grid Peclet number in the form of equation 6.37,
equations 6.35 and 6.36 must be rewritten in terms of the dependent variables
Sw and STW by using the chain rule for differentiation. After performing the
necessary chaining operations (see Appendix D for details) the Peclet equations
for the water- and gas-phases are denned at each node i in the domain as:
where the notation
and
) Tw represents the partial derivative of the func-
tion (•) with respect to Sw and STW, respectively, and 9S indicates the angle
that the s— direction takes to the direction of gravity. In addition, when chaining
the terms involving gas- water interactions, the following approximations based on
capillary pressure scaling (equation 4.14) are used:
\PcGW\W —
(P \ ~
\*cGW),Tw =
— -
\Pc
NWW
cGN),Tw
Note that in equation 6.38 the Peclet number is denned by the largest spatial
component. Finally, the derivatives with respect to saturation are in general
computed numerically (see Appendix D for details).
Because the terms denning the discrete form of the balance equations used in
the simulator are not the same as those which define equation 6.38, the Peclet
numbers computed from equation 6.38 provide only an estimate of the actual
discrete values. Realizing the approximate nature of the analysis, a second Peclet-
like parameter is defined using a different set of discrete parameters in an effort
102
-------
to provide an alternative measure (after El-Kadi and Ling, 1993):
Pe, = max
= x,y,z
= W,N,G
(6.39)
where the subscript (, s) represents the partial derivative in the s—dimension, and
a Pej is denned for each phase. This expression represents a measure of the shape
of the a—phase saturation front relative to the mesh based on the fractional flow
function. Equation 6.39 is made amenable to computation in the discrete model
by chaining the numerator in terms of the saturation variables, i. e.:
Pe, = max
fa
= x,y,z
where the three Pe values are defined by the largest spatial component. Note
that the derivatives (Sw},s and (STW},S are available from the nodal Hermite coeffi-
cient vector, and that the fractional flow derivatives are available from computing
equation 6.38, therefore no additional computations are required to generate this
information.
The parameters A.NW? A-NG and A.QW are computed at each node i in the do-
main by requiring that Pecnt is not violated. The procedure is not so straight
forward because the parameter AGW occurs in both balance equations, thus cou-
pling the system. Therefore, consider the following algorithm which is based on
numerical experiment and is designed to simplify the analysis:
AWN. = max jl, Pew/Pecrit, PeW/Pecrit\
ANGi = max jl, PeG/Pec"*, P^G/Pec"*,TeN / Pecrit\
I J i
r i-C1
for fr >
AWN. = max {1, Pew/Pecrit, Pe /Pec"*, Pe /Pecrit}
> r - J{
ANG. = max I 1, PeG/Pecrit, Pe /Pecrit
for /f < f?
(6.40)
where the particular combination is chosen based on the relative magnitudes of
the gas- and water-phase fractional flow functions. Note that the equations in
6.40 reduce to the proper two-phase equations when applicable.
103
-------
In summary, before each solution step for the saturation equations, equation
6.40 is computed at each node in the domain using the current nodal saturation
data. This equation is based on two Peclet-like constraints and a user-defined
critical value, Pecrit. The result is that capillary diffusion is augmented where
necessary, thereby ensuring that the flow variables are compatible with the spatial
discretization. Finally, three important points should be highlighted:
1. based on numerical experiments conducted for a broad range of problem-
types, Pecnt = 2 appears to be the optimal value for the finite-element
formulation chosen.
2. if Pecnt is made large enough, then no additional diffusion will be added to
the system;
3. because diffusion is added in a point-wise manner and in a known magnitude,
one can use this information in a diagnostic sense to either better qualify
the solution or to design the mesh in a more appropriate manner.
6.5.2. Time Step Control
Stability, accuracy and computational efficiency considerations require that a dy-
namic time stepping control be included in the simulator. In general, the value of
a given time step is based on one of three criteria:
1. a Courant constraint;
2. the number of iterations required to solve any one equation over a time step;
3. how much the soultion has changed over the last time step.
Based on a Courant constraint
One way to ensure that time truncation errors are not adversely affecting the flow
solution as it evolves in time is to define a maximum time step size based on
a Courant constraint. Advection-diffusion partial differential equations, such as
those which describe phase mass transport, can be characterized by the dimen-
sionless group of parameters called the Courant number, Co,
^-i (advection)(characteristic time scale) /^ A~, \
(characteristic space scale) \ ' '
104
-------
As with the Peclet number, a discrete point- wise measure of the Courant number
can be denned at a node i as:
Coi = max {\advection\s / A s , s = x,y, z}. A t (6.42)
where At is the time step and the advective component of transport at node i is
denned by the largest component of the product of the advection vector and the
grid spacing in each of the s-directions (s = x, y, z).
Given a constant spatial discretization and a user-defined upper limit for Co,
called C*oc™*, expression 6.42 can be used to provide a measure of the maximum
time step as a function of flow variables and space step. Specifically, using an
analogous derivation as was used for the Peclet constraint, given the current solu-
tion to the problem, the maximum time step allowed for the next solution interval
is approximated as:
Co
crit
/As,
= x,y,z
.=1
(6.43)
where the denominator represents the maximum value of the quantity at all the
/ nodes in the mesh.
Based on numerical experiments conducted for a broad range of problem-types,
Cocnt < 2 appears to provide an effective constraint on time step for the finite-
element formulation chosen. In general, Cocnt should be smaller for three-phase
flow problems (Cocrit < 1).
Based on iterative convergence
While equation 6.43 provides a measure of the maximum time-step allowable over
a time step, an adaptive time-stepping scheme is employed to determine an appro-
priate time-step size based on the number of iterations required for convergence
of the nonlinear problem. This provides an effective tool which minimizes the
cumulative number of iterations required over the course of the simulation.
Atfc+i = e/ A tk for m < mmax, e/ > 1
= £R A tk for m > mmax, 0 < eR < I
105
-------
Based on change in solution
As discussed in subsection 5.4.3, when considering hysteresis, the denominator of
equation 5.14 can become very small when scanning curves are generated. There-
fore, small changes in saturation can yield large changes in effective saturation
and capillary pressure which in turn can cause instabilities. As a result, a con-
straint is put on the maximum change in either Sew °r SCTW over a time step, i.
e., if during iterative convergence,
where 6SfSK is the maximum allowable change, then the time step is restarted
with a new At equal to the old At times a factor e < 1. Based on numerical
experiments, the following values provide useful limits: if considering hysteresis
then 8Sfax ~ 1.0, if not considering hysteresis then 8Sfax ~ 0.1.
6.5.3. Phase Discontinuities
With respect to the solution of the phase flow equations, the appearance and
disappearance of phases is accommodated by using the total flow formulation
in conjunction with the requirement that the pressure variables are continuous in
time and space. The result is that equations 6.10, 6.25 and 6.26 are representative
of the system regardless of phase configuration.
With respect to the solution of the contaminant transport equations, note that
equations 6.27 and 6.28 are undefined when Sw = 0 and STW = 0, respectively.
Therefore, a special algorithm is required to accommodate these conditions as the
solution evolves through time. First assume that the variables p" , a = W, G*,
are continuous in time and space. Second, given the phase saturation distribution
from the solution to the flow problem identify those nodes where Sw < e and
SG < e, where e is a small positive number of order 0.001. Third, a moving
boundary condition problem is set up by specifying Dirichlet data for pn
at any node where Sw < e, and for pn at any node where SG < e.
Therefore, the matrix equation for each of the NAPL species contains only those
algebraic equations written at the collocation points associated with non-Dirichlet
condition nodes.
106
-------
6.6. SOLVING THE SYSTEM OF LINEAR EQUATIONS
The system matrix generated by the Hermite collocation finite-element method
has the attributes of not being symmetric, positive definite, or diagonally domi-
nant (Dyksen and Rice, 1986). In addition, the fact that the degrees of freedom
and the collocation equations are defined at different locations in the domain
(at the nodes and at the Gauss quadrature points, respectively), the collocation
system matrix has a structure which depends strongly on the numbering scheme
chosen. For regular meshes (as are used in this simulator) two desirable num-
bering schemes have been identified. First is termed finite-element ordering by
Lai et al. (1994). In finite element ordering the collocation points are numbered
consecutively within each element and each element is numbered in the shortest
directions first, while the degrees of freedom are numbered consecutively by node
and each node is numbered in the shortest directions first (see Frind and Finder
[1979] for an example in 2-D). This numbering scheme has the attribute of creat-
ing the most compact block-banded matrix possible, where the matrix bandwidth
is given by:
2(4nelsm(nelmed + 1) + 11) + 1
where nelsm is the number of elements in the smallest dimension and nelmed is
the number of elements in the medium or next smallest dimension. It also has the
attribute of creating a matrix with mostly zeros on the diagonal. As a result, the
finite-element ordering scheme is not generally amenable to iterative solution. It
is however amenable to direct solution, and this version of the simulator provides
an option to use a banded LU decomposition with partial pivoting (LAPACK
driver routine DGBSV, 1993).
The second numbering scheme termed tensor-product ordering by Lai et al.
(1994), numbers degrees of freedom and equations in the following way (see Lai
et al. [1994] for an example in 2-D):
1. associate each unknown with a nearest-neighbor collocation point.
2. sweep along lines in the mesh marching in the shortest directions first.
3. sweep along each line four times, each time numbering two unknowns and
associated collocation points, where the pairing of unknowns is based on
their association with Dirichlet and Neumann data (see Appendix C).
107
-------
This numbering scheme creates a matrix with a block tri-diagonal structure,
and non-zeros on the diagonal, which makes solution by an iterative method possi-
ble. This version of the simulator provides an option to use a preconditioned GM-
RES (generalized minimum residual) method (SLATEC driver routine DSLUGM,
1989 [incomplete LU GMRES]).
6.7. SUMMARY OF SEQUENTIAL ITERATION
The sequential solution iteration procedure can be summarized by the following
series of computational steps:
1. Given Sw, STW, P^ and pG and boundary data at time-level k:
compute Pw/ using the total flow equation (6.10);
compute qT, qw/, qG from equations 6.14;
for boundary nodes with Dirichlet conditions 2, 3 or 4, convert the com-
ponent of qT normal to the boundary to the appropriate source/sink
terms, Qw , QN, QG, as described in Section 6.4.2.
2. Given Sw and STW, qw/, qG, Qw and QG from step 1:
compute Dw/, DG;
compute pjf and pG via the concurrent iteration algorithm described in
Section 6.2.2.
3. Given the current values of qT, Qw , QN and Q from step 1, and the current
values of p^ and pG from step 2:
compute Sw and STW via the concurrent iteration algorithm described in
Section 6.2.1.
4. Repeat steps 2 and 3 as required for convergence.
108
-------
7. SIMULATOR DOCUMENTATION
The simulator has been written in standard FORTRAN 77 and does not require
linkage to any external libraries. Input and output (I/O) files consist of a series of
standard ASCII files. In general this section details I/O file content and format.
The input files which provide the necessary data to run the simulator can be
defined either manually or by using suitable pre-processing software designed to
write to the necessary files in the required format. The output files can be read
by appropriate post-processing software. This section is written so that the user
can adapt appropriate pre- and post-processors in an efficient manner.
The documentation in this section consists of the following:
• Data Input: Detailed description of data input files, data requirements, and
data format, including references to sections detailing important algorithmic
aspects.
• Input and output file description: A listing of input and output files associ-
ated with the simulator and a brief explanation of file content and format.
Note, Chapter 8 provides a detailed description of several illustration problems
intended to provide the user with examples of model capability and problem set-
up and implementation strategies. Reference is made to the disk location of the
appropriate input files used to set up and run each problem, and the output files
used to analyze the results. As such, the description of simulator I/O contained
in this Chapter can be augmented using the examples found in Chapter 8.
Finally, note that Appendix E contains a listing of the names of the source
code FORTRAN files which are compiled and linked to create the executable
application and a brief explanation of each file's contents and computational task.
7.1. DATA INPUT
With respect to data input, consider the following important notes:
109
-------
parameter
permeability (k)
porosity (e)
soil dry bulk density (pb)
fraction of organic carbon (foc)
residual saturations (Swr, Sj\fnr, Sj\fwr, S(jr)
S — P model fitting parameters (a
-------
3-D
2-D
Figure 7.1: Default axis orientation for the 3-dimensional and the 2-dimensional
versions of the simulator. Note that the x-axis is aligned with gravity in both
versions. Also illustrated are the definitions of the grid rotation option and the
input parameter nface which is used in setting boundary conditions
9. The default axis orientation for both the three-dimensional and two-dimensional
versions of the simulator is illustrated in Figure 7.1, where the x—direction
is vertical and increasing downward, the y—direction is horizontal and in-
creasing to the right, and the z—direction is horizontal and increasing into
the plane of the paper. The two-dimensional version uses the same (x, y)
orientation and ignores the z—dimension. As described in more detail in
subsection 7.1.2, an option for grid rotation is available where the user de-
fines the angels to rotate the grid in the counter-clockwise direction about
the z— and y—axes (9Z and dy, respectively in Figure 7.1).
10. The rectangular mesh is defined by the union of rectangular elements where
the total number of elements is defined by the product of the number of el-
ements in the x—, y—, and z—directions (nex, ney, and nez, respectively).
Note, the user must make sure that the dimension statements in the IN-
CLUDE file, include.f, reflect these values (see subsection 7.1.1 for details).
11. The global node numbering scheme which is required for all input requiring
a node number follows the rule: increment first in the y—direction, then in
x—direction, and then in the z—direction, i. e.,
Ill
-------
nodel: x(I), y(I), z(I) [e.g. (0,0,
node 2: rr(l), y(2), 2(1)
node nny: x(l), y(nny\ z(I)
node nra/ + 1: x(2)
node nra/ nnrr: x(nnx), y(nny\ z(l)
node nny nnrr + 1: x(l), |/(1), z(2)
..., etc.
where nnrr = (nex + 1), nra/ = (ney + 1) and nnz = (nez + 1) are the
number of nodes in the x—, y—, and z— directions respectively.
12. As discussed in subsection 6.4.2, the simulator assumes that all boundary
nodes are no flow with respect to the phases and homogeneous Neumann
with respect to the NAPL species, unless otherwise specified. When apply-
ing a Dirichlet or mixed condition at a boundary node it is necessary to
provide the simulator with information qualifying the axis (in 2-D) or plane
(in 3-D) that is normal to the applied condition. This is because the sim-
ulator uses Hermite cubic polynomials to interpolate the primary variables,
and as discussed in subsection 6.4.2 and Appendix B, the group of Hermite
coefficients which are associated with a given boundary condition is a func-
tion of the axis or plane which is normal to the direction from which the
condition is applied. Using Figure 7.1 as a definition sketch, the simulator
uses the following code to identify the axis or plane type, where for the 3-D
simulator the code is: nface = 1 if the face is an x — y plane, nface = 2 if
the face is a y — z plane, and nface = 3 if the face is an x — z plane, and for
the 2-D simulator the code is: nface = 1 if it is the x — axis and nface = 2
if it is the y — axis. Finally, note that this information is technically only
required at corner nodes in 2-D and edge and corner nodes in 3-D. For ex-
ample, in 2-D, one can apply two different conditions at a corner node, such
as a constant head condition on the horizontal side (nface = 2) and a no
flow condition on the vertical side (nface = 1).
112
-------
In the discussion which follows for reference purposes, all variable names writ-
ten in italic font are those used in the FORTRAN source code, and all file names
are written in bold font.
7.1.1. Set up for the include.f file
The include.f file defines the array dimensioning and common block definition,
and for the 2-D version houses a flag which determines which linear algebraic
equation solver to use.1 There are three groups of PARAMETERS which must
be set to reflect the particular problem one wishes to run.
1. Maximum mesh definition.
For the two-dimensional version, set the PARAMETERS mnndJ, the max-
imum number of nodes in the long dimension, and mnnd-S, the maximum
number of nodes in the short dimension, where
mnndJ > max[(nea; + 1), (ney + 1)]
mnnd_s > min[(nea; + 1), (ney + 1)]
and nex and ney are the number of elements in each spatial dimension of
the model mesh.
For the three-dimensional version, set the PARAMETERS mgnd_x, mgndjy,
and mgnd-z such that
mgnd-X > (nex + 1)
mgnd.y > (ney + 1)
mgnd-Z > (nez + 1)
NOTE, setting these parameters to the minimum values minimizes run-time
memory requirements.
2. Flags to choose which linear algebraic solver to use. For the 2-D version
only. Two options exist, use a direct solver (LU with partial pivoting), or
use an iterative solver (GMRES with incomplete LU preconditioning). See
subsection 6.6 for details on solver description. The flags are itsolJ, and
itsoLf where:
JThe 2-D version allows one to choose between a direct solver and an iterative solver. The
3-D version employes only the iterative solver.
113
-------
it sol.t: if = 1, then use the iterative solver for the transport-like equations
(the equations for Sw, STW, P% and /)„), if = 0, then use the direct
solver.
it sol. f: if = 1, then use the iterative solver for the total flow equation (the
equation for Pw), if = 0, then use the direct solver.
In general both these flags should be set equal to one. Note, the 3-D version
employes only the GMRES solver.
3. For the GMRES iterative solver, the maximum number of search-direction
vectors to be saved and orthogonalized against before restart is set:
nsave > I
In general nsave should be as large as the maximum number of iterations
allowed, however, this maximizes memory requirements. As nsave increases,
the number of iterations required for convergence decreases, but the memory
requirement increases. Results from numerical experimentation suggests
nsave = 10 to 50. NOTE on the use of the GMRES solver. The combination
of tensor-product ordering of the equations and unknowns (described in
subsection 6.6) and the use of the GMRES solver, requires the following
constraint-rules for problem setup: for the 3-D code, the first 2-D plane
of nodes numbered (the smallest dimension boundary plane) must have at
least one node with a Dirichlet condition, and for the 2-D code the first 1-D
line of nodes numbered (the smallest dimension of x and y) must have at
least one node with a Dirichlet condition applied for the method to work.
This does not apply to the direct solver.
No other changes to the file include.f are required to be made.
7.1.2. Data Input Driver
The file sni.in is the 'driver' for data input. It either contains the required data
which is read from the FORTRAN file niain.f, or it defines the 'call' to other data
input files. The other data input files define the mesh, the initial and boundary
conditions, and all spatially varying parameters (Table 7.1). A description of
these input files can be found in subsection 7.2.1.
The first 74 lines of sni.in are fixed, and each line requires a fixed set of one
or more parameters as specified below. All data must be input whether the data
114
-------
is used or not (e. g., if the problem is set up so that contaminant transport is
switched off, then some 'dummy' numbers must be specified for all transport-
specific parameters).
This section describes the input requirements for the file sm.in by line number.
For each line number, the parameter list is specified in italics, where the names of
these parameters are those which are used in the FORTRAN source code (provided
for reference purposes). Following the parameter list is a brief description of the
meaning of the parameter(s), where appropriate documentation subsections are
referenced for additional detail.
I. title
provide a character string (a maximum of 68 characters) for the title of the
output to appear in file echo.out which echoes the data input (see section
7.2.5 for details on what echo.out contains).
2. iphase
Determines which phases are to be modeled. The following code is used.
Let water = 1, NAPL = 2, gas = 3, then you have three options:
iphase = 12, water-NAPL flow problem (solve two flow equations and
SG = 0)
iphase = 13, water-gas flow problem (solve two flow equations and SN = 0)
iphase = 123, water-NAPL-gas flow problem (solve three flow equations)
Admissable values: 12, 13, or 123. If iphase does not equal one of these
three values then an error message appears on the screen and the simulator
stops. If one of the admissible values are input then a warning message is
output to the screen to tell that you which phases are to be considered in
the simulation.
3. iscr
A switch to turn screen output on/off. If iscr=l then iteration, solution
and time information is output to the screen during the simulation. If iscr
is not equal to one then screen output is suppressed. Iteration and solution
information is output to the screen in the order that the computations are
performed. This output is diagnostic in nature and provides the user with
run-time information on simulator performance. See subsection 7.2.4 for
details on the contents of the screen output.
115
-------
4. iherm
Option to write Hermite data to file at print-intervals (see input line 8 for
definition of print interval). If iherm=l then output the data, if iherm is
not equal to one then do not. If enabled, then at each print interval in a
separate file for each Hermite parameter, a time header is printed followed
by a node list of the Hermite coefficients for that parameter. The files thus
generated are: pa.out (Pw data), sw.out (Sw data), st.out (STW data),
oa.out (p% data), og.out (p^ data). The data for each print interval over
the course of the simulation is concatenated. See subsection 7.2.5 for details
on the contents of these files.
5. mass
Option to perform and output mass balance information. If mass=l then
compute and output the data, if mass is not equal to one then do not.
If enabled, then after each time step a material mass balance is performed.
The files thus generated are: mass.out (summary), cmass.out (cumulative
mass balance as a function of time), massw.out (water mass balance over
each time step), massn.out (NAPL mass balance over each time step),
massg.out (gas mass balance over each time step), masst.out (total phase
mass balance over each time step). Poor mass balance is an indication of
incorrect time and/or space discretization, or incorrect specification of the
problem forcing terms. See subsection 7.2.5 for details on the mass balance
computations and output files.
6. nex, ney, nez
The number of elements in the x—, y— and z—dimensions (rectangular
mesh). Note, make sure that the PARAMETER statements in the IN-
CLUDE file include.f reflect these values, i. e., for the 2-D simulator:
mnndJ > max[(nea; + 1), (ney + 1)]
mnnd-S > min[(nea; + 1), (ney + 1)]
and for the 3-D simulator:
mgndjx > (nex + 1)
mgnd.y > (ney + 1)
mgnd-Z > (nez + 1)
116
-------
7. time
The time associated with the initial conditions (e. g., the reference time).
For example, if one is restarting the simulator from a previous run, this
input variable can be used to date the data associated with the restart files
(see subsection 7.2.2 for details on restart file content and print frequency).
8. tlpr, tmprnt, tmax
Define print interval and maximum simulation time, where:
tlpr - first print when elapsed simulation time equals time + tlpr
tmprnt - the time interval between prints after the first print.
tmax - the maximum simulation time.
Note, the time step is adjusted so that print times and tmax are honored.
The following files are written to at each print interval and the data for each
print interval, identified by its time, is concatenated:
Hermite data (if output is enabled, see input line 4): pa.out (Pw data),
sw.out (Sw data), st.out (STW data), oa.out (p% data), og.out (p^
data)
General solution: sat.out (SV, 5jv, SG, PCNW, PCGN, PVn and Pn)
Mass balance (if output is enabled, see input line 5): mass.out (summary)
Flux vector data: velg.out (gas), veln.out (NAPL), velw.out (water)
Restart files: all files with the extension rs (data used for simulator restart,
see input line 68)
Reference subsection 7.2.5 for details on output file contents.
9. itincs, itincc, tmult
Time step control based on the iterative solution of the nonlinear saturation
and contaminant transport equations. In general larger time steps require
more iterations to converge on the nonlinearity. After each solution inter-
val, if the number of iterations required for convergence of each nonlinear
equation is less than the values specified, where
itincs - saturation iterations (suggestion: 3 < itincs < 15)
117
-------
itincc - contaminant transport iterations (suggestion: 3 < itincc < 15)
then increase the time step for the next solution interval, i. e., A tnew =
tmult A t0id, where
tmult - factor to increment (tmult > 1.0, and tmult PS 1.15 is an appropri-
ate value)
See subsection 6.5.2 for more detail on dynamic time-step control.
10. itreds, itredc, tdiv
Time step control based on the iterative solution of the nonlinear saturation
and contaminant transport equations. After each solution interval, if the
number of iterations required for convergence of any one nonlinear equation
is more than the value specified, where
itreds - saturation iterations (itreds > itincs and suggestion: 8 < itreds <
20)
itredc - contaminant transport iterations (itredc > itincc and suggestion:
8 < itredc < 20)
then reduce the time step for the next solution interval, i. e., A tnew =A
iv, where
tdiv - factor to reduce (tdiv > 1.0, and tdiv PS 1.25 is an appropriate value)
See subsection 6.5.2 for more detail on dynamic time-step control.
11. ithangs, ithangc, tdivh
Time step control based on the iterative solution of the nonlinear saturation
and contaminant transport equations. If during a solution interval, the
number of iterations required for convergence of any one nonlinear equation
exceeds the specified value, where
ithangs - saturation iterations (ithangs > itreds and suggestion 10 <
ithang < 25)
ithangc - contaminant transport iterations (ithangc > itredc and sugges-
tion 10 < ithanc < 25)
118
-------
then the current time step is restarted with a reduced time step, A tnew =A
where
tdivh - factor to reduce (tdivh > 1.0, and tdivh ~ 1.5 is an appropriate
value)
This flag usually trips when the iterative scheme will not converge as the
iteration solutions repeat with a period > 2 (a ringing phenomenon).
12. itermx
When the GMRES solver is used, this parameter represents an upper limit
on the number of iterations allowed to solve any one matrix equation. Typi-
cally, this variable applies to the elliptic pressure equation which is the most
difficult to solve iteratively. This number may be based on the number of
iterations required to solve the pressure equation for the first time step. Af-
ter setting a reasonable value, if for a given time step, the pressure equation
will not converge within this maximum number of iterations the simulator
will stop without saving the current solution and an error message will be
output to the screen. To restart after a 'crash' caused by a maximum iter-
ation violation, one must go back to the time of the last output interval as
specified by the print interval chosen (see input line 8). This type of 'crash'
is indicative of an excessively coarse convergence criterion for the saturation
equations. A number < 200 is appropriate.
13. dslim
An upper limit for the change in Sew and SBTW over a time step. As discussed
in subsection 5.4.3, when considering hysteresis, the denominator of the
effective saturation (equation 5.14) can become very small as the S— P curve
is fit between the predefined constraint points. Therefore, small changes in
saturation can yield large changes in effective saturation. If during iterative
convergence this threshold is reached, then restart the time step with a new
At equal to the old At divided by the factor tdivh (defined in input line
11). Suggested working values: if considering hysteresis then dslim ~ 1.0,
if not considering hysteresis then dslim ~ 0.1. See also the discussion in
subsection 6.5.2.
14. co
119
-------
The critical Courant number, Cocrit as denned in subsection 6.5.2. Used
to define an upper limit on the time step size based on time truncation
error. Appropriate values are case-specific, however, in general co < 4 and
smaller values of co are appropriate for three-phase flow problems (e. g.,
co < 1). NOTE: make co a large number to avoid this upper limit time step
constraint.
15. dtO
The time step used to start the simulation. Typically the initial time step is
relatively small so that the simulator can propagate boundary and forcing
terms which are initially shock-like. Let the dynamic time stepping tools
(defined by the parameters listed in input lines 9 through 14) bring the time
step up to its 'optimal' value.
16. tsmax, tsmin
The upper and lower limit for time step, respectively. Time steps which
are too large or too small can lead to instabilities as numerical round-off
errors dominate the solution. The value tsmax overrides all other time step
incrementing computations. The maximum time step is problem specific
and can be estimated or obviated by employing the Courant constraint. If
the time step is cut to below tsmin due to lack of iterative convergence, then
the simulation stops and the last solution is output as per a typical print
interval. This type of 'crash' is an indication that the discrete problem and
the physical problem are not compatible.
17. grf_ on, grinc, ngrch, fgrch, gmax
Animated graphics output control. For the 2-D version of the simulator these
files are formatted specifically for the Jacquard graphics software package2.
For the 3-D version of the simulator these files are formatted specifically for
the SciAn graphics software package3 according to the STF format.
The simulator generates files with the extension stf and jin to be read
by the post-processing software for animated display of saturation (sw.stf,
st.stf), concentration (oa.stf, og.stf) and velocity (velw.stf, veln.stf and
An EPA sponsored pre- and post-processing graphics software package written specifically
for Silicon Graphics computer hardware.
3 SciAn is a public domain post processing software package. For information on obtaining
SciAn, send email to: scian-info@scri.fsu.edu.
120
-------
velg.stf). The stf files contain time-dated concatenated output where the
output interval is defined as follows:
grf_ on - if = 1 then enable graphics output, =0 then turn off.
grinc - the initial time increment to print
ngrch - the number of prints at current increment before changing the
increment
fgrch - factor to change the current increment
gmax - maximum increment
18. nloop
The number of sequential iteration loops between the pressure equation and
the saturation equation (s) during sequential solution. Setting this value
to zero results in sequential solution, values > 0 results in an iterative se-
quential solution. Use a value > 0 only for problems where separate phase
velocities must be well defined. See section 6.2 for details on the sequential
solution algorithm.
19. erip
GMRES convergence criterion for the pressure equation. Experience sug-
gests that values of order 1.0* 10~6 provide accurate solutions with relatively
low iteration requirements. Use mass balance as your ultimate guide.
20. eris, eras
Convergence criteria for the saturation equations. The parameter eris is
for the GMRES solver (inner iterations), and the parameter eros is for
convergence on the nonlinearity (outer iteration). Experience suggests that
eris xi eros pa 1.0 * 10^3 provides accurate solutions with relatively low
iteration requirements. Use mass balance as your ultimate guide.
21. erit, erot
Convergence criteria for the contaminant transport equations. The parame-
ter erit is for the GMRES solver (inner iterations), and the parameter erot
is for convergence on the nonlinearity (outer iteration). Experience suggests
that erit ?a 1.0* 10^5 and erot ?a 1.0* 10^3 provides accurate solutions with
relatively low iteration requirements. Use mass balance as your ultimate
guide.
121
-------
22. idxdy
Mesh definition flag. The main program now reads information from the
file space.in and the flag idxdy defines what must be in the file. If idxdy =
0 then the simulator will generate a uniform mesh based on the number
of elements specified in line 6 and the maximum dimensions listed in file
space.in, i. e.,
xmax, ymax, zmax
If idxdy = I then the simulator will generate a mesh based on a catenated
x—, y— and z—dimension node spacing list in file space.in, i. e.,
x—node spacing list, x(i), i = 1, nex + 1
y—node spacing list, y(i),i = 1, ney + I
z—node spacing list, z(i\ i = 1, nez + 1
23. grav
The magnitude of the gravity vector [L/T2].
24. th_z, th_y
Orient the grid at a specified angle to the horizontal, where th_ z is the angle
to rotate the grid in the counter clockwise direction about the z — axis,
and th_y is the angle to rotate the grid in the counter clockwise direction
about the y — axis. For example, th-z = 0 and thjy = 0 yields the default
orientation with +x = depth (see Figure 7.1), and th-z = 90 and thjy = 0
yields — y = depth. Note, the two-dimensional version of the simulator only
considers the angle th_z.
25. permb
The value of the intrinsic permeability scalar [L2] applied by default at every
node in the domain. Spatial variability is defined in the next input line.
26. ndev
The number of nod-specific deviations from the global value permb. If
ndev > 0 then the simulator opens the file perm.in and reads ndev lines of
information, where each line lists the node number and the value of perme-
ability (nod, perm)
122
-------
27. porb
The value of the soil porosity scalar [dimensionless] applied by default at
every node in the domain. Spatial variability is denned in the next input
line.
28. ndev
The number of node-specific deviations from the global value porb. If
ndev > 0 then the simulator opens the file por.in and reads ndev lines of
information, where each line lists the node number and the value of porosity
(nod, por).
29. bulkb
The value of the dry bulk density [M/L3] applied by default at every node
in the domain. This parameter is used to model adsorption mass transfer
(see subsection 4.8.2). Spatial variability is defined in the next input line.
30. ndev
The number of node-specific deviations from the global value bulkb. If
ndev > 0 then the simulator opens the file bulk.in and reads ndev lines
of information, where each line lists the node number and the value of the
dry bulk density (nod, bulk).
31. vw_r, vn_r, vg_r
The viscosity (pure phase) [M/(TL)] of the water-, NAPL- and gas-phases,
respectively.
32. rw_r, rn_r, rg_r
The density (pure phase) [M/L3] of the water-, NAPL- and gas-phases,
respectively.
33. siggw, signw, siggn
The fluid-fluid interfacial tension [M/T2] between the water and gas, the
NAPL and water, and the gas and NAPL phases, respectively. These para-
meters are used to scale the capillary pressures for three-phase flow. Recall
that the simulator assumes that the S — P model parameters were fit to
a two-phase displacement experiment (the type of which is defined in the
next input line), and that it uses the interfacial surface tension data to scale
123
-------
the S — P model to represent the other fluid pairs. See sections 4.3, 5.1
and 5.5 for the details regarding the use of capillary pressure scaling in this
simulator. Note, any set of self-consistent units can be used for these values
since they are used to generate dimensionless ratios.
34. n_phase
This entry tells the simulator which phase pair was used to obtain the S — P
model fitting parameters. This information is used to define the capillary
pressure scaling relationship. Given the following form of equation 4.14:
the J3 parameters are defined as a function of the value of n-phase, i. e.,
njphase = 1
njphase = 2
njphase = 3
Paw = 1
PGN = (°GN + &NW) I°GN
PNW =
PGW =
PGN =
PNW =
Admissable values: 1, 2, or 3. If n_phase does not equal one of these three
values then an error message appears on the screen and the simulator stops.
See sections 4.3, 5.1 and 5.5 for the details regarding the use of capillary
pressure scaling in this simulator.
35. swr, snnr, snwr, sgr
Global definition of residual saturations applied by default at every node in
the domain, where:
swr = residual water phase (Swr)
snnr = residual NAPL as a nonwetting phase (Sffnr)
snwr = residual NAPL as a wetting phase (S?jwr)
sgr = residual gas phase (Scv)
124
-------
Spatial variability of these parameters is defined in the next input line.
36. ndev
The number of node specific deviations from the global values swr, snnr,
snwr, sgr. If ndev > 0 then the simulator opens the file residual, in and
reads ndev lines of information, where each line lists the node number and
the four residual values swr, snnr, snwr, sgr. Use this option if the domain
contains different soil types which warrant the specification of unique S — P
model fitting parameters. This file is usually set up in concert with the files
shape.in (input line 38) and base.in (input line 40).
37. asd, asi, eta
The S — P model curve shape parameters (equation 5.13) applied by default
at every node in the domain, where:
asd = the pressure scale for drainage (ad) [l/L]
asi = the pressure scale for imbibition (a 0 then the simulator opens the file shape.in and reads ndev lines
of information, where each line lists the node number and the three fitting
parameters asd, asi, eta. Use this option if the domain contains different
soil types which warrant the specification of unique S — P model fitting
parameters. This file is usually set up in concert with the files residual.in
(input line 36) and base.in (input line 40).
39. perm_ b, por_ b
The soil-type parameters upon which the residual saturations (input line
35) and the curve shape parameters (input line 37) are based. Applied by
default at every node in the domain, where:
perm-b = intrinsic permeability scalar [L2]
porJ) = porosity
125
-------
These parameters are equivalent to k* and e* denned in section 5.5, equation
5.41. Note that the soil parameters permb and porb denned in input lines
24 and 26 are equivalent to k and e in equation 5.41. The combination of
these soil property definitions defines capillary pressure scaling relationship
with respect to soil properties.
40. ndev
The number of node specific deviations from the global values perm_ b and
por_ b. If ndev > 0 then the simulator opens the file base.in and reads
ndev lines of information, where each line lists the node number and the
two parameters k* and e*. Use this option if the domain contains different
soil types which warrant the specification of unique S — P model fitting
parameters. This file is usually set up in concert with the files residual.in
(input line 36) and shape.in (input line 38).
41. alfw, nsewl, nsew2
These parameters define the krw(Sw} functional (equation 5.19), where,
alfw is the connectivity parameter £
nsewl is a flag which determines the definition of aSew, i.e.,
nsewl
nsewl
= I
= 2
a c
a c
Oew
= (Sw - S
= (Sw - S
W*) / (1
W*) / (1
— Swt
- Swt
-sNt
)
-Sat)
nsew2 - a flag which determines the definition of bSew, i-e-5
nsew2
nsew2
= I
= 2
b Q
bo
OeW
= (Sw - S
= (Sw - S
Wt) 1 (1
wo / (i
- SWt - S,
— Swt)
Nt ~ S(jt)
If these flags do not equal one of the numbers indicated then an error message
is printed and the simulator stops.
42. alfn, nsenl, nsen2, nsenS
These parameters define the krN(Sw, SQ) functional (equation 5.38), where,
alfn is the connectivity parameter £
126
-------
nsenl is a flag which determines the definition of S£N, i.e.
nsenl
nsenl
= I
= 2
SeN —
SeN =
(SN — SNt]
(SN — SNt]
/ (1 - Swt
/ (i - sm)
-sNt
-SGt)
nsen2 - a flag which determines the definition of SeTn, i.e.
nsenl = I
nsen2 = 2
nsenl = 3
SeTn —
q t>N+t>G — t>Nt—t>Gt
^eTn — —/-, o
nsenS - a flag which determines the definition of SeTw, i.e.
nsen3 = I
nsen3 = 2
nsen3 = 3
If these flags do not equal one of the numbers indicated then an error message
is printed and the simulator stops.
43. alfg, nsegl, nseg2
These parameters define the
G) functional (equation 5.23), where,
alfg is the connectivity parameter £
nsegl is a flag which determines the definition of aSeG, i.e.,
nsegl
nsegl
= I
= 2
a o
iJeG —
a o
JeG —
(SG
(SG
- SGt) /
- Sot) 1
(1 — Swt
(1 - Sot)
-sNt
-SGt)
nseg2 - a flag which determines the definition of bSea, i.e.,
nsegl = 1
nsegl = 2
6 C
b C
G = (SG
G = (SG
~ S^ 1 (1
- sGt) 1 (i
- Swt - S,
-Sat)
vt ~ Sat)
If these flags do not equal one of the numbers indicated then an error message
is printed and the simulator stops.
127
-------
44. se_sl, se_rl
Linearly extrapolate the S — P functionals at effective saturations near one
and zero to make them amenable to numerical model implementation as
discussed in subsection 5.4.3, Part 1, where the meaning of these parameters
is as follows:
if Se > (I — se_sl) then linearize the S — P functional (se_sl = S^,s in
equation 5.39)
if Se < sejrl then linearize the S — P functional (sejrl = S^r in equation
5.39)
Numerical experiments suggest that values between 0.001 and 0.01 have no
noticeable effect on the solution.
45. sfact_ kr
Force kra(Se} = 0 for Se < sfact-kr. This parameter eliminates the effect
that small oscillations in the saturation and concentration solutions have on
residual saturations. Numerical experiments suggest that sfactjkr pa 0.01
has no noticeable effect on the solution.
46. nhyst
Switch to turn on/off the hysteresis option. If nhyst = 1, then use the
hysteretic k — S — P model definition as defined in Chapter 5, equations
5.13 and 5.14. If nhyst = 0, then use the model for monotonic displacement,
equations 5.10 and 5.11 with Ss = 1 and Sr = SWT.
47. e_ r
Phase entrapment/release definition - the blending parameter (e) defined
in equation 5.7 which governs how fast the phase becomes entrapped dur-
ing drainage flow conditions or released from entrapment during imbibition
flow conditions (see also subsection 5.4.1). ejr = I yields a linear entrap-
ment/release model.
48. b_a
Blend the S — P model scale parameter (am in equation 5.13) between
reversal points. This is the blending parameter (3 as used in Table 5.2 which
defines how fast the van Genuchten scale parameter (am.) changes after a
128
-------
reversal occurs. In general 6_a > 0.2. See subsection 5.4.2 for additional
detail.
49. sp_ min, sr_ min
Restriction parameters for S — P curve generation and update, where:
sp_min is the minimum span (the denominator denning effective satura-
tion), for which a new set of scaling parameters will be calculated. See
subsection 5.4.3, Part 2, where sp_min = span. Numerical experiments
suggest a range: 0.025 > sp_min > 0.1
sr_ min is a minimum tolerance away from the S—P curve end points above
which a scanning curve will be generated after a flow reversal, i. e., if
sr_min < Se < (I — sr_min) then a scanning curve will be generated
after a flow reversal otherwise the computation will be suppressed. See
subsection 5.4.3 Part 3, where sr_min = e. Numerical experiments
suggest a range: 0.025 > sr_min > 0.1 and sr_min > sp_min.
50. factd, facti
The tolerance in saturation change to indicate a reversal from drainage to
imbibition (factd) or from imbibition to drainage (facti). See subsection
5.4.3 Part 4, where factd = r^ and facti = r-i. In general factd = facti
> 0.001.
51. pe_w, pe_g
The critical Peclet numbers for the water and gas phases, respectively. The
parameter Pecrit in equation 6.40, where pe_ w is used for the water equation
and pe_g is used for the gas equation. In general pe_w = pe_g > 2. See
section 6.5.1 for details.
52. pg_ref
The reference gas phase pressure [M/(LT2)]. Usually set to atmospheric
pressure (e. g., 1 * 106dynes /cm2).
53. ntr_ ow, ntr_ og
Switches to turn on/off constituent mass transport and mass transfer com-
putations, where:
129
-------
ntr_ ow - if =1 then model dissolved NAPL contaminant transport (solve
for pjf); if =0, then bypass
ntr_ og - if =1 then model NAPL vapor contaminant transport (solve for
p^); if =0, then bypass
As indicated below, when either or both of the contaminant transport com-
putation switches is (are) turned off, then the input parameters required
for contaminant transport are set internally by the simulator to appropri-
ate values, thereby overriding user denned input. If either ntr_ ow = 0 or
ntr_ og = 0 then a warning message is output to the screen to tell that you
have turned off appropriate transport mechanisms.
54. theta
The projection parameter used to model the coupling between the water-
and gas-phase contaminant transport equations due to mass transfer of the
NAPL species between phases. The parameter 9 in equation 6.29. This
parameter is used only if ntr-ow = ntr_og = 1.
55. along, atran, diffw, diffg
Dispersion tensor definition (see subsection 4.6), where:
along - longitudinal dispersivity - the parameter a^ [L] in equation 4.25
where a^ is the same for both phases.
atran - transverse dispersivity - the parameter a-p [L] in equation 4.25 where
-------
focb is the fraction of organic carbon in the soil, the parameter /oc, applied
by default to every node in the domain.
Note, these parameters are used only if ntr_ ow = 1.
57. ndev
The number of node specific deviations from the global value focb. If
ndev > 0 then the simulator opens the file o_c.in and reads ndev lines
of information, where each line lists the node number and the value of the
fraction of organic carbon in the soil (nod, foe). Note, if ntr_ow = 0 then
ndev is set to zero internally.
58. d_ layer
The thickness of the boundary layer which separates the boundary between
the domain and the atmosphere. Used for the mixed boundary condition on
the p^ transport equation, where d_ layer defines a constant parameter 8
(equation 6.32). See subsection 6.4.2 for more detail. Note, this parameter
is used only if ntr_ og = 1 and a third-type boundary condition is specified.
59. bow_ 1, bow_ 2, bow_ 3
Definition for the first-order mass transfer kinetic for the NAPL-phase dis-
solving into the water-phase (see subsection 4.8.1), where
bow_l: the rate coefficient (3l \\/T\ in equation 4.29
bow_2: the coefficient (32 in equation 4.29
bow_3: the coefficient j3-A in equation 4.29
Note, if ntr_ ow = 0 then bow_l is set to zero internally.
60. parow
The equilibrium concentration of the NAPL species in the water phase (sol-
ubility limit), the parameter pjf [M/L3] in equation 4.28. Note, if ntr_ow
= 0 then parow is set to zero internally.
61. bog_l, bog_2
Definition for the first-order mass transfer kinetic for the NAPL-phase va-
porizing into the gas-phase (see subsection 4.8.1), where
131
-------
bog_ 1: the rate coefficient (3l [i/T] in equation 4.31
bog_2: the coefficient (32 m equation 4.31
Note, if ntr_ og = 0 then bog_I is set to zero internally.
62. parog
The equilibrium concentration of the NAPL species in the gas-phase (solu-
bility limit), the parameter ~pP [M/L3] in equation 4.30. Note, if ntr_og =
0 then parog is set to zero internally.
63. bowg_l, bowg_2
Definition for the first-order mass transfer kinetic for the dissolved NAPL
species (p%) vaporizing into the gas-phase (see subsection 4.8.1), where
bowg_ 1: the rate coefficient /5fw/ [1/T] in equation 4.34
bowg_2: the coefficient /32 in equation 4.34
Note, if either ntr_ ow = 0 or ntr_ og = 0 then bowg.I is set to zero inter-
nally.
64. e_ henry
The dimensionless Henry's law constant. The parameter H in equation 4.33.
Note, used only if both ntr_ ow = 1 and ntr_ og = 1.
65. t_half
The effective half-life of the NAPL species in the water and gas phases [T]
(the parameter t-\_m in equation 4.26). By default the decay coefficients, K^
and K^, are set to zero if the input parameter tJialf < 0. If tJialf > 0
then equation 4.26 is used to define the decay coefficients.
66. swinit
The global initial condition for water saturation (the parameter SWQ in
subsection 6.4.1).
67. stinit
The global initial condition for total wetting phase saturation (STWQ = 5
-------
iphase = 12 (see input line 2) and stinit ^ 1
iphase = 13 (see input line 2) and stinit ^ swinit
swinit > stinit
See Appendix C for the decision tree which defines initial k — S — P model
given initial saturations.
68. roainit
The global initial condition for dissolved contaminant concentration (the
parameter p^0 in subsection 6.4.1). Note, that if ntrjow = 0 then roainit
is set to zero internally, and if ntr_ow = 1, (stinit — swinit) > 0, and
swinit > 0, then roainit is set to the solubility limit, parow, internally.
These operations override the user-defined input.
69. roginit
The global initial condition for vapor contaminant concentration (p^0). Note,
that if ntrjog = 0 then roginit is set to zero internally, and if ntr_og = 1,
(stinit — swinit) > 0 and (I — stinit) > 0, then roginit is set to the solubility
limit, parog, internally. These operations override the user-defined input.
70. ncont
Restart option switch. If ncont = 1 (yes) = 0 (no). If yes, then the de-
fault initial values defined in lines (64 through 67) are overwritten and the
initial conditions are read in from the set of files with the rs extension (for
restart). These files are in binary format, and they contain output from the
last output interval of the previous simulation. See input line 8 for output
interval definition and subsection 7.2.2 for description of restart files.
71. ndev_s
Number of deviations from initial global saturation definition. If ndev_s > 0
then the simulator opens the file sat.in and reads ndev lines of informa-
tion, where each line lists the node number and the values of Sw and STW
(nod, sw, st). This data over-writes all previous data. Note an error message
will appear on the screen and the simulator will stop if any of the following
input combinations occur:
iphase = 12 (see input line 2) and st
133
-------
iphase = 13 (see input line 2) and st ^ sw
sw > st
Also note, at all nodes where (st — sw) > 0, if ntr_ow = I and sw > 0,
P%Q is set internally to the solubility limit, parow, and if ntrjog = I and
1 — st > 0, p^Q is set internally to the solubility limit, parog. Finally, see
Appendix C for the decision tree which defines initial k — S — P model given
initial saturations.
72. ndev_roa
Number of deviations from initial global definition. If ntr_ow = 0 and you
input ndevjroa > 0 then the code prints a WARNING to the file echo.out
alerting you that transport is off before setting ndevjroa to zero internally.
Otherwise, if ndevjroa > 0 then the simulator opens the file roa.in and
reads ndevjroa lines of information, where each line lists the node number
and the value of p^0 (nod, roa). This data overwrites all previous data.
73. ndev_rog
Number of deviations from initial global definition. If ntrjog = 0 and you
input ndevjrog > 0 then the code prints a WARNING to the file echo.out
alerting you that transport is off before setting ndevjrog to zero internally.
Otherwise if ndevjrog > 0 then the simulator opens the file rog.in and
reads ndevjrog lines of information, where each line lists the node number
and the value of p^0 (nod, rog). This data overwrites all previous data.
74. nbcroa
The number of Dirichlet boundary conditions for the parameter pjf (recall
from subsection 6.4.2 that since homogeneous Neumann conditions are ap-
plied as the default at all boundary nodes, the user needs to specify only
Dirichlet conditions). If ntrjow = 0 and you input nbcroa > 0 then the code
prints a WARNING to the file echo.out alerting you that transport is off
before setting nbcroa to zero internally. Otherwise, if nbcroa > 0 then the
simulator opens the file bc_roa.in and reads nbcroa lines of information,
where each line lists the following:
nface, nod, roa
where,
134
-------
nface is the mesh face which is normal to the direction from which the
condition is applied. See figure 7.1 for the definition of nface for the
2- and 3-D simulators.
nod is the node number for specified boundary condition
roa is the value of the known concentration [M/L3]
75. nbcrog
The number of Dirichlet boundary conditions for the parameter p^ (recall
from subsection 6.4.2 that since homogeneous Neumann conditions are ap-
plied as the default at all boundary nodes, the user needs to specify only
Dirichlet and mixed conditions)4. If ntrjog = 0 and you input nbcrog >
0 then the code prints a WARNING to the file echo.out alerting you
that transport is off before setting nbcrog to zero internally. Otherwise,
if nbcrog > 0 then the simulator opens the file be rog.in and reads nbcrog
lines of information, where each line lists the following:
nface, nod, rog
where,
nface is the mesh face which is normal to the direction from which the
condition is applied. See figure 7.1 for the definition of nface for the
2- and 3-D simulators.
nod is the node number for specified boundary condition
rog represents one of two things. If rog > 0, then rog is the value of the
known concentration [M/L3] at the node. If rog < 0, then the 2-D
simulator assumes that the condition is of the mixed type and applies
equation 6.32 at the node. The 3-D simulator will not accept rog < 0
as a valid entry, and it will print an error and then stop.
76. ncondf
The number of flow-variable boundary conditions different from the de-
fault no flow condition. If ncondf > 0 then the simulator opens the file
bc_flow.in and reads the data for ncondf boundary nodes, where each
node requires the following:
4The 2-D version allows both Dirichlet and mixed conditions. The 3-D version allows only
Dirichlet conditions.
135
-------
n/ace, nod, ncode
where
n/ace is the mesh face which is normal to the direction from which the
condition is applied. See Figure 7.1 for the definition of nface for the 2-
and 3-D simulators, and the for the significance of this parameter refer
to the discussion in Section 7.1, Item 12. Note, while this information
is only important at corner and edge boundary nodes, a value must be
specified even at internal boundary nodes. If at any boundary node
there is a discrepancy between the user-defined value of n/ace and the
admissible values, then an error message will appear on the screen and
the simulator will stop.
nod is the node number for the specified boundary condition.
ncode defines the condition-type, where the code follows from the number-
ing scheme used in Table 6.1. In summary:
if ncode = 2 then prescribe gas-phase head (hG = PG/7G)
if ncode = 3 then prescribe the NAPL-phase head (hN = PN/JN)
if ncode = 4 then prescribe the water-phase head (/iw/ = Pw/ /7W/)
if ncode = 5 then prescribe h , Sw and St
'G
Refer to subsection 6.4.2 for additional details on admissable boundary
conditions. Note if ncode does not equal one of these values, then an
error message will appear on the screen and the simulator will stop.
Depending on the value of ncode the simulator expects the following
information on the following input line:
If ncode = 2, list the gas head value [L], h-g
If ncode = 3 list the NAPL head value [L], /i_n
If ncode = 4 list the water head value [L], h-w
If ncode = 5 list the parameter n_ opt, and then the values for water
head [L], water saturation and gas saturation values,h-w, sjw, s_g,
where
n_ opt tells the simulator whether to use the Dirichlet data listed
(n_opt = 0), or to define the Dirichlet data using the initial con-
ditions from a previous run or stress period (n_opt = 1). Note,
by setting n_ opt = 1 the code assumes that the current solution
at the node is the Dirichlet data, therefore, n_ opt = 1 can only
136
-------
be used in conjunction with either ncont = 1 (see input line 68)
or after a stress period change (see input line 68). If neither form
of initial data is available then an error message will be printed to
the screen and the simulator will stop. Also note that if any of the
following input parameter combinations are used an error message
will be printed to the screen and the code will stop:
ncode = 2 and iphase = 12
ncode = 3 and iphase = 13
77. nwella
Specification of point sources and sinks (wells), where nwella is the number
of nodes at which a point source or sink condition is to be applied. If
nwella > 0 then the simulator opens the file well.in and reads nwella lines
of information, where each line lists the following:
nod, q_ tot, ff_ u>, ff_ g, roa, rog
where
nod is the node number where the well is idealized
q_tot is the total well flow rate [L3/T] (QT in subsection 6.4.3), where
injection is positive (+), extraction is negative (-).
ff_ w is the fractional flow of water of the injected fluid, where the simulator
internally computes the water flow rate as5: Qw = ff-w qjtot.
ff_ g is the fractional flow of gas of the injected fluid, where the simulator
internally computes the gas flow rate as2 : QG = ff~g q-tot.6
roa is the concentration of the NAPL species in the water-phase in the
injected fluid when ntr_ ow = 1, and no NAPL is present.2 If NAPL
is present and ntr_ ow = 1, then roa is set internally to the solubility
limit parow. If ntr_ ow = 0 then roa is set to zero internally.
rog is the concentration of the NAPL species in the gas-phase in the injected
fluid when ntr_ og = 1, and no NAPL is present.2 If NAPL is present
°Note that this only applies for injection wells (i. e., qJt.ot > 0) where the makeup of the
source fluid is known. At extraction wells (i. e., qJ.ot < 0) the the makeup of the sink fluid is a
function of the solution. See equation 6.33.
6At injection wells, given qJLot, ff~w and ff-g, the NAPL flow rate is known since QN =
qjkot-Qw-QG.
137
-------
and ntr_ og = 1, then rog is set internally to the solubility limit parog.
If ntr_ og = 0 then rog is set to zero internally.
78. timec, iph_new, dtnew, co_new, dtmax
Define a new stress period: change the boundary and external flux conditions
at specified times during the simulation without stopping, where
timec - time for new stress period.
iph_ new - define number and types of phases modelled (redefines the pa-
rameter iphase initially defined on input line 2)
dtnew - new initial time step.
co_new - new Courant constraint (redefines the parameter co initially
defined on input line 14)
dtmax - new maximum time step (redefines the parameter tsmax initially
defined on input line 16).
If the simulation time is greater than or equal to timec then what follows is
a re-specification of boundary conditions starting at (72.). Therefore, repeat
input lines 72 through 76 for each stress period change.
This concludes the description on data input. Examples of the actual input
files can be found in Chapter 8.
7.2. INPUT- AND OUTPUT-FILE DESCRIPTION
7.2.1. Input Files
All the files are written in standard ASCII format. Reference to specific input
lines refers to those described in subsection 7.1.2. Input files are listed by category.
1. Main data driver
sm.in - defines all data input for the simulator. Defines calls to other data
input files.
2. Grid spacing
138
-------
space.in - grid spacing information, list first nodal x-coordinates (nnx
entries) and then nodal y-coordinates (nny entries)
3. Initial conditions
roa.in - list of nodal exceptions to global specification of p^ as defined in
input line 70 of file sm.in. Lists the node number and the specified
value.
rog.in - list of nodal exceptions to global specification of /)„ as defined in
input line 71 of file sm.in. Lists the node number and the specified
value.
sat.in - list of nodal exceptions to global specification of Sw and STW as
defined in input line 69 of file sm.in. Lists the node number and the
specified Sw and STW values.
4. Spatially-varying physical constants (Table 7.1)
bulk.in - list of nodal exceptions to global specification of the soil dry bulk
density as defined in input line 30 of file sm.in. Lists the node number
and the value of the soil bulk density.
o_c.in - list of nodal exceptions to global specification of the fraction of
organic carbon in the soil as defined in input line 55 of file sm.in. Lists
the node number and the value of organic carbon content.
perm.in - list of nodal exceptions to global specification of permeability
scalar as defined in input line 26 of file sm.in. Lists the node and the
permeability value.
por.in - list of nodal exceptions to global specification of porosity as defined
in input line 28 of file sm.in. Lists the node and the porosity value.
shape.in - list of nodal exceptions to global specification of S — P curve
fitting parameters ad, a,L and r\ as defined in input line 38 of file sm.in.
Lists the node number and the values of ad, a,!f and r\.
5. Boundary and external forcing conditions
be flow.in - list of nodal Dirichlet data for flow for each stress period.
Data for each stress period is concatenated.
139
-------
be roa.in - list of nodal Dirichlet data for dissolved NAPL concentration.
Data for each stress period is concatenated.
be rog.in - list of nodal Dirichlet or mixed data for NAPL vapor concen-
tration. Data for each stress period is concatenated.
well, in - list of nodal point source and sink data for each stress period.
Data for each stress period is concatenated.
7.2.2. Restart Files
Restart files contain un-formatted double precision output written at the last
specified output interval (as determined by the print command), and they consti-
tute all the necessary nodal information required to restart the simulator from the
time of the last specified output interval. Used for either a planned continuation,
or after an unexpected program stoppage. The files contain nodal solution data,
and they are categorized into two groups:
1. those files associated with the hysteretic k — S — P model:
a.rs -the scaling parameters a^ used in equation 5.13 and defined in Table
5.2.
nhc.rs - the integer curve-type indicator /.
s_max.rs - the highest a-phase saturation that has occurred since it was
last at immobile residual conditions, S™3* (equation 5.7).
sb.rs - the saturation at which curve / was spawned, «Sb(/) used in Tables
5.2, 5.3 and 5.4.
sr.rs - the minimum wetting phase saturation for curve-type /, Sr(f) used
in equation 5.13 and defined in Table 5.3 5V(/)-
ss.rs - the maximum wetting phase saturation for curve-type /, Ss(f) used
in equation 5.13 and defined in Table 5.4.
trap c.rs - the current trapped quantity, Sat, defined in equation 5.7.
trap mn.rs - the lower limit of entrapped a-phase S™11 used in equation
5.7 and defined in equation 5.8.
trap mx.rs - the magnitude of the residual a-phase at the terminus of
an a-phase drainage process, S^r, used in equation 5.7 and defined in
equation 5.9.
140
-------
2. those files associated with the solution variables:
oa.rs - the nodal Hermite pjf data
og.rs - the nodal Hermite p1^ data
pa.rs - the nodal Hermite Pw/ data
sw.rs - the nodal Hermite Sw data
st.rs - the nodal Hermite STW data
7.2.3. Compilation Files
The simulator uses the FORTRAN 77 INCLUDE statement. The file include.f
is used to set the dimensions for the arrays used in the simulator. See subsection
7.1.1 for proper dimensioning requirements.
After changing the include.f file one needs to compile the code. A file called
makefile contains the UNIX system commands which create and link the object
files to yield the executable file napl.
7.2.4. Screen output
If the output flag iscr = 1 (set in line 3 of file sm.in ) then iteration and solu-
tion information is output to the screen in the order that the computations are
performed (see subsection 6.7 for a summary of the computational steps). This
output is diagnostic in nature and provides the user with run-time information
on simulator performance.
• P - GMRES (number of iterations) or P - DIRECT. The total flow
equation was just solved for Plv. Note, this equation is linear, there-
fore iteration is necessary only when the GMRES iterative solver is
employed.
• FLUX CONVERSION FOR BCs 2, 3 and 4. At each boundary node
where a Case 2, 3 or 4 flow boundary condition is applied, the bound-
ary fluxes were just converted into point source/sink conditions. The
output lists the phase flux data for each node:
node number Qlv QN QG
where Qa [L3/T] is the equivalent volumetric flux of the a—phase at
the node.
141
-------
• Rog - GMRES (number of iterations) or Rog - DIRECT. .
Roa - GMRES (number of iterations) or Roa - DIRECT.
The contaminant transport equations were just solved for p^ and p^,
respectively. Note, this information is displayed only if contaminant
transport is turned on (i. e., ntr_ow = 1 and/or ntr_og = 1 in input
line 48 of file sm.in), otherwise it is suppressed. Also, if contaminant
transport is on, then since the concentration equations are nonlinear,
the iteration information for nonlinear convergence is displayed in the
next line, where the respective entries from left to right are:
Equation-type : NL_C,
iteration number,
node where the largest change occurred,
the values of p% and p^ at that node,
the LI norm of the iterative increment,
the infinity norm of the iterative increment.
• STw - GMRES (number of iterations) or STw - DIRECT.
Sw - GMRES (number of iterations) or Sw - DIRECT.
The saturation equations were just solved for STW and Sw, respectively.
Note, if iphase = 12, then only Sw is solved, and if iphase = 13,
then only STW is solved (see input line 2 in file sm.in). Also, since
the saturation equations are nonlinear, the iteration information for
nonlinear convergence is displayed in the next line, where the respective
entries from left to right are:
Equation-type : NL_S,
iteration number,
node where the largest change occurred,
the values of Sw, SN and SG at that node,
in parentheses, the maximum change in Sew or S£TW over the time
step,
the L2 norm of the iterative increment,
the infinity norm of the iterative increment.
After a convergent solution is computed for the time step, the time information
is printed, including from left to right,
142
-------
the elapsed time (in the units specified) since the beginning of the simula-
tion,
the time step value for the last time step,
the number of times the time step had to be restarted because of lack of
iterative convergence (see input line 12),
the number of time steps taken since beginning of the simulation.
7.2.5. Output files
Output files are grouped according to the information they provide. We identify
three groups:
1. Output files to be used by post-processing graphics software. These files are
generated only if the graphics output is enabled (see input line 17 in file
sm.in):
solution vectors (see input line 16 for output interval definition)
sw.stf and st.stf (saturations Sw and STW)
oa.stf and og.stf (concentrations, p^ and p^ )
velw.stf, veln.stf and velg.stf (fluid velocities vw/ , v^ and VG )
soil properties
soil.stf - column 1 is nodal values of permeability, column 2 is nodal
values of porosity.
mesh information
mfile.jin - contains a list of nodal coordinates and element-to-node
connectivity
2. Output files containing mass balance and simulation performance informa-
tion, this information is output only if the mass balance computations are
turned on. There are six files that contain mass balance information.
cmass.out - cumulative mass balance calculation after each time step.
Mass balance for each phase at time T > 0 is defined as:
(the cumulative phase mass that has crossed the boundary from t = 0 to
T)
143
-------
/ (the total phase mass in the domain at time T ).
The file lists the following nine entries for each time step:
elapsed time,
water phase mass error,
NAPL phase mass error,
gas phase mass error,
A t,
number of Pw iterations,
number of Sw iterations,
number of SG iterations,
number of iteration hang-up/restarts.
mass.out - for each specified solution output (as defined by the print
interval), summarizes mass, time step, and iteration performance.
massg.out - gas mass balance error per time step. Mass balance over the
time step is defined as:
(the change in the amount of phase mass in the computational domain over
At)
/ (the phase mass that has crossed the boundary over At)
Perfect mass balance is indicated by a ratio of one. The following six
entries are listed for each time step:
elapsed time
gas mass in over boundary over At
gas mass out over boundary over At
change in gas mass over the boundary over At
change in gas mass in the domain over At
mass balance ratio defined above.
masso.out - same as massg.out except for NAPL phase.
massw.out - same as massg.out except for water phase.
masst.out - total fluid mass balance error per time step. For each time
step, lists the elapsed time and the ratio:
(the change in the amount of total phase mass in the computational domain
over At)
144
-------
/ (the total phase mass that has crossed the boundary over At)
3. Output files containing information on the solution. The frequency of the
output is defined by the print interval (see input line 8 in file sm.in), and
the output is concatenated.
If Hermite data output is enabled (see input line 4 in file sm.in), then the
following Hermite coefficient information at each node in the domain
is provided:
sw.out, st.out (saturations Sw and STW)
pa. out (water pressure)
oa.out and og.out (concentrations p% and p^}
Referenced by time information, each print interval lists the solution
data for each node:
node number, values of the 4 (8) degrees of freedom.7
echo.out - echoes parameter input, and lists WARNING messages.
sat.out - Solution summary. Referenced by time information, each print
interval lists the following solution data for each node:
Sw, SN, SG, PCNW, PCGN, Pln anci p%
velg.out - the computed gas flux components at each node.
veln.out - the computed NAPL flux components at each node.
velw.out - the computed water flux components at each node.
7 Four degrees of freedom for the 2-D simulator and eight degrees of freedom for the 3-D
simulator.
145
-------
8. MODEL TESTING AND EXAMPLE
PROBLEMS
The focus of this Section is on testing the model for self-consistency and its ability
to simulate experimental procedures and results. Self-consistency is established
by investigating spatial and temporal convergence attributes and mass balance
performance. A series of pseudo-one-dimensional example problems are presented
in order to evaluate convergence and mass balance, and to give the user an in-
dication of appropriate discretization for a given set of input data. In addition
to addressing self-consistency, four example problems were designed to simulate
specific physical experiments:
1. a three-phase LNAPL spill and redistribution experiment (Van Geel and
Sykes, 1995);
2. a three-phase DNAPL spill and redistribution experiment conducted at the
EPA's Subsurface Protection and Remediation Division of the National Risk
Management Research Laboratory in Ada, OK;
3. an experimental investigation of the dissolution of residual DNAPL in a
saturated sand (Imhoff et al, 1992);
4. an experimental investigation of DNAPL vapor transport in an unsaturated
sand (Lenhard et al., 1995).
As discussed in subsection 8.3, the data sets for these problems 1, 2, and 4 are
included with the software in the appropriate dedicated directory.
8.1. CONVERGENCE AND MASS BALANCE
8.1.1. Compatibility of the grid and the flow model
Because we are modeling the emplacement and dissolution of NAPL in the vadose
zone, we need to consider problems which evolve to quasi-static flow conditions.
146
-------
As a result, it is imperative that the spatial discretization be compatible with the
flow parameters, especially those defining the S — P model, so that oscillations
in the saturation solutions are minimized. This point is illustrated in Figure 8.1
which compares the capillary rise for two different soils.
The problem is to simulate water drainage in a one-dimensional soil column,
1.2 meters long and initially saturated with water. The relevant soil properties,
S — P model parameters, and resulting van Genuchten PDC are shown in Fig-
ures 8.la [for a relatively fine sand] and 8.1c [for a relatively coarse sand]. The
boundary conditions are: at the top, open to the atmosphere, and at the bottom,
specified water head (fine sand = 10 cm, coarse sand = 60 cm). The columns are
then allowed to drain under the influence of gravity until quasi-static conditions
prevail. Figure 8.1b shows that for the fine sand, an appropriate grid spacing is
approximately 10 cm, and Figure 8. Id shows that for the coarse sand an appropri-
ate grid spacing is approximately 2 cm. It must be noted here that the simulator
includes no implicit mechanism to add artificial diffusion, or capillarity, to the sys-
tem (for example, an upstream weighting algorithm is not employed). Artificial
capillarity must be included explicitly either by altering the S — P curve fitting
parameters or by utilizing the Peclet criterion (see Section 6.5.1 for details).
With respect to time step, we note the following. Because the Dirichlet pres-
sure data is a nonlinear function of saturation, and because no saturation Dirichlet
data has been specified, the time step represents the only explicit mechanism to
damp oscillations in the pressure solution as the system approaches steady-state.
The appropriate size of the time step is problem-dependent, and in general, the
more nonlinear the S — P functional the smaller the quasi-steady-state time step
must be. For example, the quasi-steady-state At for the problem in Figure 8.1b
is of order 1000 seconds, while that for the problem in Figure 8. Id is of order
25 seconds. It is suggested that as the system approaches steady state, the gas
pressure boundary condition be changed to Dirichlet conditions on Pw, Sw and
SQ (a case 5 boundary condition, see Section 6.4.2).
8.1.2. Analysis of the three-phase hysteretic k-S-P model
Here we consider convergence and mass balance attributes of the simulator when
the three-phase hysteretic k — S — P model is employed. Consider the model
problem presented in Figure 8.Id, and the result therein to be the initial conditions
for a DNAPL spill simulation. Relevant model parameters are presented in Table
8.1. These parameters mimic those to be used in an artificial aquifer experiment
147
-------
120
100
80
60
40
20
0
120
110
,100
90
80
70
60
50
Experimental S-P relationship
n=6.49
a = 0.0203/cm
Swr = 0.17
K = 8.64m/d
0.2 0.4 0.6 0.8
water saturation
Experimental S-P relationship
n= 11.434
a = 0.0849/cm
Swr = 0.08
K=170m/d
120
100
80
60
40
20
0
Computed steady-state moisture profile
0 0.2 0.4 0.6 0.8 1
water saturation
0.2 0.4 0.6 0.8
water saturation
120
110
100
90
80
70
60
50
Computed steady-state moisture profile
(d)
dx = 5 cm
dx = 2cm
0 0.2 0.4 0.6 0.8 1
water saturation
Figure 8.1: Analysis of appropriate grid spacing to compute capillary rise for
different soil-types. Parts (a) and (b) are for a relatively fine sand, and parts (c)
and (d) are for a relatively coarse sand.
148
-------
Fluid properties
pw = 0.9982 g/cm3
yUW/ = 0.01 poise
S-P model definition
ad = 0.0849/cm
SWr = 0.08
GGW = 72.75 dynes /cm
k-S model definition
aSew from eq 5.25
aSeG from eq 5.26
Se]\[ from eq 5.30
C = V? = £ = 0.5
Field properties
e = 0.3115
pN = 1.626 g/cm3
^N = 0.0093 poise
a>i = 0.12 /cm
SGr = 0.16
&NW = 39.5 dynes/cm
bSew from eq 5.25
"S^c from eq 5.26
S£TW from eq 5.32
k = 170 m/d
p(> = 0.00129 g/cm3
/jG = 0.0002 poise
TJ = 11.434
SNnr = 16, SW = 0.08
&GN = 31.74 dynes/cm
SeTn from eq 5.35
Table 8.1: PARAMETERS USED TO MODEL THREE-PHASE PCE MIGRA-
TION IN OTTAWA SAND
being conducted at RSKERL involving a controlled release of a DNAPL. The
boundary conditions for the current problem are given as: /iw = 60 cm at the
bottom boundary for all time, and, at the top boundary, the following time-varying
conditions are applied:
1. For time = 0 to 100 s, DNAPL is injected at a constant volume rate of 0.03
cm3 / s.
2. For time = 100 s to the end of the simulation, P = atmospheric.
This problem is considered a severe test of the numerical model for two reasons:
1. The S — P functionals are extremely nonlinear.
2. The incompressible gas-phase assumption (while not appropriate for this
simulation) leads to a result in which, as the DNAPL is injected, the water
table is depressed. Then, after the source is removed and the top re-opened
to the atmosphere, the water table rebounds, forcing the gas-phase to per-
colate up through the descending DNAPL.
149
-------
Simulation results for two different discretizations are presented in Figure 8.2.
A time step of order 2 seconds was required to obtain a solution for this problem, as
larger time steps caused convergence problems. Figure 8.2a shows the total liquid
saturation solution at initial conditions (T = 0) and at T = 200 s (100 s after the
DNAPL source was removed). One can see that while both discretizations capture
the sharp DNAPL front, the x = 2.5 cm solution exhibits oscillations behind the
front. Figure 8.2b presents saturation results at time = 5000 s, after the DNAPL
has migrated to near static, residual state. It is apparent that for these model
parameters, a grid spacing of approximately 1 cm is required. Figures 8.2c and
8.2d present mass balance results for this simulation. The definition of the mass
balance ratio used in the figure is:
(the change in the amount of total phase mass in the computational domain over /At)
/(the total phase mass that has crossed the boundary over /At)
Perfect mass balance over a time step is indicated by a ratio of one. One can see
that in general the model performs well with respect to mass balance except when
boundary forcing is changed, and several time steps are needed to accommodate
the discontinuity imposed.
8.1.3. Analysis of the mass transfer model
Here we investigate the convergence attributes of the kinetic mass transfer model.
Consider the following one-dimensional water flow and contaminant transport
problem where the domain is the same as that defined in Figure 8.1c, and the
model parameters are those given in Table 8.1. The initial conditions are set such
that the domain is saturated, and there is a zone of residual DNAPL, SN = 0.15,
uniformly distributed from x = 25 cm to the bottom. The boundary conditions are
set such that there is a constant influx of clean water at the top at a rate of 0.008
cm/s, and an equivalent efflux of contaminated water at the bottom. Relevant
mass transfer and transport parameters are: (32 = 0-5 and /33 = 1.0, avl = I
cm, and pjf = 0.001 gm/cm3. The results for different values of the exchange
rate coefficient, /3fN\ are presented in Figures 8.3a and 8.3b. As shown in the
Figure, a distinct dissolution front is created, the shape of which is a function of
the size of {3™/N. High values effectively approximate the equilibrium partitioning
approximation and produce a sharp front, while low values produce a broad front.
From a numerics standpoint, the dissolution front should be resolved over several
elements to minimize oscillations in the solution which can cause erroneous NAPL
saturations upstream of the source area.
150
-------
120
110
100
Jl 90
f 8°
70
60
50
T = 0, dx=2.5cm
T = 0, dx=lcm
T = 200, dx=2.5cm
T = 200, dx=lcm
0.2 0.4 0.6 0.8 1
total liquid saturation
120
100
0 0.2
(b)
Sw, dx = 2.5 cm
Sw, dx = 1 cm
Sn, dx=2.5 cm
Sn, dx= 1cm
0.4 0.6 0.8
saturation
z.
1.5
1
1
0.5
A
0
1C
z
EARLY TIME
dx= 1cm
dt=2s 0 1.5
Is
I
80.5
T=100
i Q
0 200 300 400
time
(d)
LATE TIME
dx= 1cm
dt=2s
1000 2000 3000 4000 50C
time
Figure 8.2: Results of a one-dimensional, three-phase, DNAPL injection and re-
distribution simulation, highlighting spatial convergence and mass balance.
151
-------
0.2 0.4 0.6 0.8
NAPL saturation/O.I 5
0 0.2 0.4 0.6 0.8 1
Dissolved concentration / 0.001
0 0.2 0.4 0.6 0.8 1
NAPL saturation/O.I 5
0 0.2 0.4 0.6 0.8 1
Dissolved concentration / 0.001
Figure 8.3: Computational analysis of the dissolution model. Parts (a) and (b)
illustrate the effect that the rate constant (ex in the figure, expressed in units
of I/day) has on the solution. As the dissolution front sharpens, oscillations
appear indicating that a finer grid spacing is required. Parts (c) and (d) illustrate
spatial convergence for ex = 24/d. For the parameters chosen a grid spacing of
approximately 5 cm is appropriate.
152
-------
Fluid properties
pw = 0.9982 g/cm3
yUW/ = 0.01 poise
S-P model definition
ad = 0.0203/cm
SWr = 0.17
Paw = 1
k-S model definition
aSew from eq 5.25
aSeG from eq 5.26
Se]\f from eq 5.30
C = v? = £ = 0.5
Field properties
e = 0.374
pN = 0.6858 g/cm3
^N = 0.00409 poise
cii = 0.0271/cm
SGr = 0.20
/3NW = 0.5128
^eH/ from eq 5.25
bSeG from eq 5.26
S£TW from eq 5.34
k = 1.02a;10-7 cm2
p« = 0.00129 g/cm3
/jG = 0.0002 poise
7] = 6.49
5jVnr = «SWiur = 0.18
PGN = 0.27397
S'eTn from eq 5.37
Table 8.2: PARAMETERS USED IN THE LNAPL SPILL PROBLEM
Spatial convergence is illustrated in Figures 8.3c and 8.3d. For a constant
rate coefficient, /3^'N = 24/rf, the model exhibits convergence as the mesh is
refined. As with the previous example, oscillations in the saturation solution
appear behind the front when the grid scale is too large. As with diffusion in
contaminant transport, the size of the exchange coefficient is limited by the spatial
discretization.
8.2. COMPARISONS TO EXPERIMENTAL RESULTS
8.2.1. LNAPL Spill
Here we compare current model results to the experimental and simulation results
of Van Geel and Sykes (1995 a and b). The model parameters are presented in
Table 8.2 (note, these are the same parameters as those used in the problem
defined in Figures 8.la and 8.1b). A plot of the S — P model primary and main
curves for the simulation is shown in Figure 8.4. The two-dimensional problem
domain and boundary conditions are defined in Figure 8.5.
With reference to Figure 8.5, the forcing conditions for the problem can be
separated into three stress periods:
1. From time = 0 to 63000 s (part a), the initially saturated sand is allowed to
153
-------
100
I 60
OH
s
S3 40
1-1 20
d< zu
I
0
0
0.2
0.4 0.6
water saturation
0.8
100
60
40
20
0.2
0.4 0.6
total liquid sat.
0.8
Figure 8.4: Plot of the primary and main S-P functionals denned by the current
model for the LNAPL spill simulation, where the drainage curves are represented
by the thick lines and the imbibition curves are represented by the thin lines.
Here the fitting parameters, assumed to be valid for a water-gas system, have
been scaled to represent the water-NAPL and NAPL-gas systems.
(a) Stress periods 1 and 3
Time = 0 to 63000 s and
64120 to 66000s
Initial condition: S\y = 1
gas pressure = atmospheric
Ocm —
WT
114.5 cm-
water head = 8.7 cm
0 cm 73 cm
(b) Stress period 2
Time = 63000 to 64120 s
NAPLheadl
= 3 cm \ no
V flow
0 cm —
gas pressure =
atmospheric
cm
cnj
instrumented
section
WT
water head = 8.7 cm
Figure 8.5: Definition sketch for the LNAPL spill simulation, showing spatial scale
and boundary conditions.
154
-------
— — — current
VGS model
-60 -40 -20 0
water head (cm)
~ 40
60
80
100
experiment
current
VGS model
0.2 0.4 0.6 0.8
water saturation
Figure 8.6: Comparison of results from the physical experiment, the current
model, and the model used by Van Geel and Sykes (VGS). Part (a) shows the
vertical distribution of water pressure head and Part (b) shows the vertical dis-
tribution of water saturation.
drain, thereby creating a quasi-steady-state moisture profile.
2. From time = 63000 to 64120 s (part b), an LNAPL source is applied to the
top left corner. Specifically, a constant LNAPL head of 3 cm is applied for
1120 s to yield a cumulative infiltrated volume of 2 L.
3. From time = 64120 to 66000 s (part a), the LNAPL source is removed, and
the resulting infiltrated volume is allowed to redistribute.
The results of this comparison study are presented in Figures 8.6 and 8.7.
Figures 8.6a and 8.6b show the water pressure head and saturation solutions,
respectively, at the end of the initial water drainage regime. The current model
results match the experimental data.
Figures 8.7 illustrate the results during and after the LNAPL spill, where
the times are given in elapsed time after the LNAPL was first applied (i. e.,
T = 0 at the start of stress period two). The plots on the left show the LNAPL
distribution computed from the current model, and the plots on the right compare
the three simulation results along the instrumented section shown in Figure 8.5.
The current model captures the experimental data quite well.
The following additional details of this study are provided to highlight specific
aspects of the current model:
155
-------
16
• 32
' 48
- 65
100
114.5
16
' 32
' 48
- 65
100
114.5
6 16 26 36 46 56 73
distance from center-line
6 16 26 36 46 56 73
distance from center-line
6 16 26 36 46 56 73
distance from center-line
0
20
40
60
80
100
0
20
40
60
80
100
experiment
— — — - current
VGS model
0.4 0.6 0.8
NAPL saturation
0.4 0.6 0.8
NAPL saturation
experiment
— — — - current
VGS model
0.4 0.6 0.8
NAPL saturation
20
40
60
80
100
:__-—-£'"
• ™"
time = 600 s
_ _ _ .
experiment
current
VGS model
Figure 8.7: Comparison of results for the LNAPL spill problem. The plots on the
left show the NAPL saturation contours as computed by the current model at the
times indicated from the initiation of the LNAPL spill. The plots on the right
compare results taken along the instrumented vertical section (the vertical dotted
line in the plots on the left).
156
-------
2500
1,2000
£
| 1500
1
.> 1000
—
"3
| 500
u
0
0
model (boundary)
model (interior)
500
1000
1500
elapsed time (seconds)
2000
2500
3000
Figure 8.8: A comparisson of the cumulative LNAPL mass which has entered the
domain as a function of time. At time = 1120 s the LNAPL source was removed.
The solid line is the computed cumulative mass which has crossed the boundary.
The dashed line is the change in mass in the domain. The dash-dot line is the
experimental data.
1. Mesh definition (number of elements [spacing]): horizontal - 3 (2 cm), 2 (2.5
cm), 11 (5 cm), and 1 (7 cm); vertical- 15 (4 cm), 10 (5 cm), and 1 (4.5
cm), into page - 1 (6 cm).
2. Memory requirements (using the GMRES solver): 16 Mb RAM.
3. Time stepping information: stress period (1) - A tmax = 2000s, stress period
(2) - A tmax = 10 s, , stress period (3) - A tmax = 100 s
4. The LNAPL mass balance computation is illustrated in Figure 8.8 which
plots the cumulative LNAPL mass entering the domain as a function of
time. The solid line is the computed cumulative mass which has crossed
the boundary, while the dashed line is the computed change in mass in the
domain. A zero mass balance error is indicated when the lines are coincident.
The dash-dot line represents the experimental data for which the LNAPL
mass equals 2000 ml at time > 1120 seconds.
5. The close match with the experimental data effectively verifies the algorithm
used to impose gas and NAPL pressure boundary conditions.
6. The Van Geel and Sykes model (VGS model) utilizes the hysteretic k — S — P
model of Parker and Lenhard (1987) and Lenhard and Parker (1987). This
157
-------
Fluid properties
pw = 0.9982^/cm3
yUW/ = 0.01 poise
S-P model definition
ad = 0.04/cm
SWr = 0.12
GGW = 72.75 dynes /cm
k-S model definition
aSew from eq 5.25
aSeG from eq 5.26
Sgfj from eq 5.30
C = v? = £ = 0.5
Field properties
e = 0.37
pN = l.Gg/cm3
^N = 0.009 poise
a>i = 0.06/cm
SGr = 0.02
&NW = 31.74 dynes /cm
bSew from eq 5.25
bSeG from eq 5.26
S£TW from eq 5.32
k = 3.5xlO"7cm2
pG = 0.00129^/cm3
/jG = 0.0002 poise
ri = 10
5jVrar = 5jVujr = 0.16
°"GJV = 47.5 dynes 1 cm
SeTn from eq 5.35
Table 8.3: PARAMETERS USED IN THE DNAPL SPILL PROBLEM
particular k — S — P model does not account for NAPL entrapment as it is
being displaced by the gas-phase (a wetting phase cannot be trapped by a
nonwetting phase). This explains the under-prediction of the NAPL-phase
saturation behind the front. The comparison results support the empirical
hysteretic k — S — P model described herein.
7. Some of the k — S — P model parameters are different for the two empirical
models. Specifically, the definition of residual saturation and the capillary
scaling term (3GW.
8.2.2. DNAPL Spill
An artificial aquifer experiment was conducted by Mikhail Fishman at the EPA's
Subsurface Protection and Remediation Division of the National Risk Manage-
ment Research Laboratory in Ada, OK. The model parameters are presented in
Table 8.3, and the details of the experimental setup can be found in Appendix F.
The purpose of the experiment was to gather quantitative and qualitative data
on DNAPL migration through a variably saturated homogeneous sand. The
DNAPL used in the experiment is called tetrachloroethylene (PCE, a common
chlorinated hydrocarbon used in the dry cleaning industry). As is detailed in
158
-------
Appendix F, the data from the experiment consists of several types:
• soil and fluid properties;
• moisture retention data for the sand used in the experiment;
• experimental initial and boundary conditions;
• DNAPL influx data;
• video images of the box at various points in time showing the areal extent
of the DANPL which is dyed to maximize contrast.
There are three main types of physical experimental results which are available
for model validation:
1. The steady-state moisture profile (Table F.I). A comparison plot of the
experimental and model moisture profile data is provided in Figure 8.9.
This data represents the initial condition for the DNAPL flood.
2. The volume of PCE infiltrated as a function of time (Table F.2). The
DNAPL source was applied until 200 err? infiltrated into the aquifer. Figure
8.10 provides a comparison plot of the cumulative volume of PCE infiltrated
as a function of time. In addition, with respect to the model results, we pro-
vide a volume balance check by superimposing the plots for the cumulative
PCE volume which has crossed the boundary and the change in PCE vol-
ume in the domain as a function of time (for a perfect volume balance the
curves would be coincident).
3. A series of video frames at specific times (e. g., Figures F.7, F.8, F.9 and
F.10). Figures 8.11, 8.12, 8.13 and 8.14 provide a comparison between the
experimental and model results for the times indicated. It appears that the
model had the most difficulty in simulating the behavior in the vicinity of
the capillary fringe where all three phases have a meaningful mobility. This
indicates that the three-phase relative permeability model is missing some
important information.
159
-------
70
60
50
40
.230
20
10
dashed (computed): a= 0.04/cm, n= 10, Swr=0.12
solid (exp. data normalized by porosity = 0.37)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
water saturation
Figure 8.9: A comparison plot of the experimental and model moisture profile
used as the initial condition for the DNAPL flood.
160
-------
300
-§250
£
| 200
3
1 150
v
•a 100
—
1 50
u
0
early time
250
---- model
experiment
model (boundary)
(interior)
0 20 40
late time
60
80 100 120
elapsed time (seconds)
140
160
180
200
• experiment
• model (boundary)
model (interior)
500
1000
1500 2000 2500
elapsed time (seconds)
3000
3500
4000
Figure 8.10: A comparison plot of the experimental and model data quantifying
the cumulative volume of PCE infiltrated as a function of time. Specific to the
experiment, it took 143 seconds for 200 cm3 of PCE to infiltrate.
161
-------
Time = 143 second
Experiment Model
DNAPL Saturatioi
(shading)
Water Saturatioi
(contours)
Figure 8.11: Comparison between experiment and model results at time = 143
seconds (the time when the DNAPL source was removed).
162
-------
Time = 283 second
Experiment Model
Figure 8.12: Comparison between experiment and model results at time = 283
seconds after the DNAPL source was first applied.
163
-------
Time = 1195 second
Experiment Model
Figure 8.13: Comparison between experiment and model results at time = 1195
seconds after the DNAPL source was first applied.
164
-------
Time = 3595 second
Experiment Model
Figure 8.14: Comparison between experiment and model results at time = 3595
seconds after the DNAPL source was first applied.
165
-------
8.2.3. DNAPL Dissolution
The mass exchange/transport portion of the simulator was verified by numeri-
cally simulating a laboratory experiment which was designed to study dissolution
kinetics of residual trichloroethylene (TCE) in a uniform sand column by flushing
the system with clean water and tracking the dissolution front as a function of
time. Problem definition and model results can be found in Guarnaccia et al.
(1992).
8.2.4. DNAPL Vapor Transport
Here we consider the simulation of a two-phase, water-gas, flow and contaminant
transport experiment conducted by Lenhard et al. (1995). The experiment is
described as follows. A one meter deep by 2 meter long by 7.5 cm wide experi-
mental box is filled with a relatively coarse, homogeneous, sand. A water table
is maintained near the bottom of the sand column with a small head differential
imposed such that water flows from right to left. The experiment is set up such
that the water saturation profile is initially in static equilibrium. A zone near
the top center of the column is excavated, and a container filled with sand and
residual TCE is placed in the void to act as the contaminant vapor source. The
spatial and temporal distribution of TCE vapor concentration is measured using
a regularly spaced assemblage of sampling points. Please refer to Lenhard et al.
(1995) for additional detail.
The idealization of the experimental setup for simulation purposes, including
the model domain dimensions and initial and boundary conditions is presented
in Figure 8.15. Table 8.4 provides the relevant physical parameter data. With
respect to Table 8.4, note that several simulation parameters (identified by an
asterisk) had to be estimated either because the data was not reported or because
the model problem setup required an augmented data set. For example, the
porosity was estimated from the reported value of the soil bulk density (pb = 1.4
g/cm3), i.e., e w 1 - pb/ps = 1- 1.4/2.65 = 0.47.
With respect to the DNAPL vapor source, the authors state that the vapor
concentration as a function of time is known (i.e., measured experimental data),
but they did not report these data. As a result, we mimic the DNAPL vapor
source by placing a volume of residual NAPL saturation in the vicinity of the
experimental source (see Figure 8.15). Volatilization and dissolution of this resid-
ual saturation occurs as the result of kinetic mass transfer processes. The rate
coefficients are defined in Table 8.4.
166
-------
top boundary : PG= Patm
water and gas water and gas
concentration = 0 concentration = 0
o
o
a
pW patm
i
i low^^Hi|
kand £ ^^L
zones ^initial residual SN = 0.20
initial Sw = static saturation profile
100 cm
iwT V
water flow
no flow
pW= patm
— PW=l5.5cm
Ocm
200cm
Figure 8.15: An illustration of the DNAPL vapor transport model problem domain
including initial and boundary conditions.
We mimic the presence of a container with the sides impermeable and the top
and bottom permeable by placing low permeability and porosity zones on either
side of the source as shown in Figure 8.15.
Figure 8.16 provides a comparison plot of the experimental and model gas-
phase concentration distribution after 12 hours of source loading. Because the
source forcing conditions for the two simulation experiments are different, it is
not surprising that only the relatively low concentration values could be matched.
8.3. SOFTWARE
This section describes the material available on disk. The material can be cate-
gorized into four file-types:
1. README files (documentation-like files)
2. FORTRAN source code files
3. Executable files
4. Data input files
167
-------
Fluid properties
pw = 0.9982 g/cm3
p^' = 0.01 poise
S-P model definition
a^ = 0.156/cm
SWr = 0
&GW = 72.75 dynes /cm
k-S model definition
aSew from eq 5.25
aSea from eq 5.26
C = V? = £ = 0.5
Transport parameters
D% = 0.00001cm2/s *
fi" = I/day *
H = 0.236 *
Field properties
e = 0.41 *
p^ = 0.00553 g/cm3
^N = 0.0002 poise
Oi = 0.156/cm
SGr = 0
&NW = 31.74 dynes /cm
bSew from eq 5.25
bSeQ from eq 5.26
D% = 0. 009cm2 /s
fiw = I/day *
-ft = 0.00052 g/cm3
k = 2.lxlO-Gcm2
p° = 0.00117 g/cm3
/^G = 0.0002 poise
r] = 4.26
5jVrar = iJNwr = 0.20
UON = 47.5 dynes /cm
P?w =Q
-pW = 0.0011 g/cm3
Table 8.4: PARAMETERS USED IN THE DNAPL VAPOR TRANSPORT EX-
PERIMENT (note, parameters with an asterisk are estimated)
10
2,0
Time = 12 ti<
Figure 8.16: A comparison plot of the experimental (solid line) and model (dashed
line) results at time = 12 hours.
168
-------
The main directory (the highest level) houses an information README file
and two sub-directories:
2D - the two-dimensional version of the simulator
3D - the three-dimensional version of the simulator
Each of these sub-directories house a relevant set of FORTRAN, executable, and
data input files. Specifically, sub-directories 2D and 3D each house five sub-
directories:
For - the FORTRAN source code files and makefile:
For/UNIX - the UNIX version
For/Dos - the DOS version (the makefile is set up to use the Lahey F77
compiler)
Exe - the executable application:
Exe/UNIX - the UNIX version (compiled on the Silicon Graphics plat-
form)
Exe/Dos - the DOS version (Lahey F77 compiler)
Ex_l - the complete set of data input files required to run the problem described
in sub-section 8.2.1.
Ex_2 - the complete set of data input files required to run the problem described
in sub-section 8.2.2:
Ex_2/IC - generate the static moisture profile to be used as an initial
condition for the DNAPL flood experiment.
Ex 2/Flood - simulate the DNAPL flood and redistribution experiment.
Ex 3 - the complete set of data input files required to run the problem described
in sub-section 8.2.4:
Ex 3/IC - generate the static moisture profile to be used as an initial
condition for the DNAPL vapor transport experiment.
Ex_3/Force - simulate the DNAPL vapor transport experiment.
169
-------
BIBLIOGRAPHY
[1] Abriola, L. M., Multiphase migration of organic compounds in a porous
medium, Lecture notes in engineering, number 8, Springer-Verlag, New York,
1984.
[2] Allen, M. B., Collocation simulation of multiphase porous-medium flow,
SPEJ of AIME, pages 135-142, February 1983.
[3] Armstrong, J. E., E. O. Frind, and R. D. McClellan, Non-equilibrium mass
transfer between the vapor, aqueous, and solid phases in unsaturated soils
during vapor extraction, Water Resour. Res., 30(2), 355-368, 1994.
[4] Bear, J., Dynamics of Fluids in Porous Media, Dover Publications, Inc., New
York, 764 pages, 1972.
[5] Bradford, S. A., and F. J. Leij, Fractional wettability effects on two- and
three-fluid capillary pressure-saturation relations, J. Contaminant Hydrology,
20, 89-109, 1995.
[6] Bradford, S. A., and F. J. Leij, Predicting two- and three-fluid capillary
pressure-saturation relationships of porous media with fractional wettability,
Water Resour. Res., 32(2), 251-259, 1996.
[7] Brusseau, M. L., Transport of organic chemicals by gas advection in struc-
tured or heterogeneous porous media: Development of a model and applica-
tion to column experiments, Water Resour. Res., 27(12), 3189-3199, 1991.
[8] Brusseau, M. L., Rate-limited mass transfer and transport of organic solutes
in porous media that contain immobile immiscible organic liquid, Water Re-
sour. Res., 28(1), 33-45, 1992.
[9] Demond, A. H., and P. V. Roberts, Effect of interfacial forces on two-phase
capillary pressure-saturation relationships, Water Resour. Res., 27(3), 423-
437, 1991.
170
-------
[10] Delshad, M., G. A. Pope, and K. Sepehrnoori, A compositional simulator for
modeling surfactant enhanced aquifer remediation, J. Contaminant Hydrol-
ogy, accepted for publication, 1995.
[11] Dyksen, W. R., and J. R. Rice, The importance of scaling for the Hermite
bicubic collocation equations, SIAM J. Sci. Stat. Comput., 7(3), 707-719,
1986.
[12] Eckberg, D. K., and D. K. Sunada, Nonsteady three-phase immiscible fluid
distribution in porous media, Water Resour. Res., 20(12), 1891-1897, 1984.
[13] El-Kadi, A. I. and G. Ling, The Courant and Peclet number criteria for the
numerical solution of the Richards equation, Water Resour. Res., 29(10),
3485-3494, 1993.
[14] Essaid, H. I., W. N. Herkelrath, and K. M. Hess, Simulation of fluid distri-
butions observed at a crude oil site incorporating hysteresis, oil entrapment,
and spatial variability of hydraulic properties, Water Resour. Res., 29(6),
1753-1770, 1993.
[15] Falta, R. W., I. Javandel, K. Pruess, and P. A. Witherspoon, Density-driven
flow of gas in the unsaturated zone due to the evaporation of volatile organic
compounds, Water Resour. Res., 25(10), 2159-2169, 1989.
[16] Falta, R. W., K. Pruess, I. Javandel, and P. A. Witherspoon, Numerical
modeling of steam injection for the removal of nonaqueous phase liquids
from the subsurface: 1. Numerical formulation, Water Resour. Res., 28(2),
433-449, 1992a.
[17] Falta, R. W., K. Pruess, I. Javandel, and P. A. Witherspoon, Numerical
modeling of steam injection for the removal of nonaqueous phase liquids
from the subsurface: 2. Code validation and application, Water Resour. Res.,
28(2), 451-465, 1992b.
[18] Faust, C. F., Transport of immiscible fluids within and below the unsaturated
zone: A numerical model, Water Resour. Res., 21(4), 587-596, 1985.
[19] Faust, C. F., J. H. Guswa, and J. W. Mercer, Simulation of three-dimensional
flow of immiscible fluids within and below the unsaturated zone, Water Re-
sour Res., 25(12), 2449-2464, 1989.
171
-------
[20] Payers, F. J. and J. D. Matthews, Evaluation of Normalized Stone's methods
for estimating three-phase relative permeabilities, SPEJ, pages 224-232, April
1984.
[21] Ferrand, L. A., P. C. D. Milly, G. F. Finder, and R. P. Turrin, A comparison
of capillary pressure-saturation relations for drainage in two- and three-fluid
porous media, Adv. Water Res., 13(2), 54-62, 1990.
[22] Forsyth, P. A., and B. Y. Shao, Numerical simulation of gas venting for
NAPL site remediation, Adv. Water Res., 14, 354-367, 1991.
[23] Frind, E.O., and G.F. Pinder, A collocation finite element method for poten-
tial problems in irregular domains, Int. J. Num. Meth. Engrg., 11, 681-701,
1979.
[24] Guarnaccia, J. F., P. T. Imhoff, B. C. Missildine, M. Oostrom, M. A. Celia, J.
H. Dane, and G. F. Pinder, Multiphase chemical transport in porous media,
EPA Environmental Research Brief, EPA/600/S-92/002, 19 pages, 1992.
[25] Guarnaccia, J. F., G. F. Pinder, and M. Fishman, NAPL, Simulator Docu-
mentation, Final Report, EPA Cooperative Agreement No. CR-820499, 1997.
[26] Huyakorn, P. S., S. Panday, and Y. S. Wu, A three-dimensional multiphase
flow model for assessing NAPL contamination in porous and fractured media,
1. Formulation, J. Contaminant Hydrology, 16, 109-130, 1994.
[27] Imhoff, P. T., P. R. Jafie, and G. F. Pinder, An experimental study of the dis-
solution of trichloroethylene in saturated porous media, Princeton University
Water Resources Program Report: WR-92-1, 1992.
[28] Janes, D. B., Comparison of soil-water hysteresis models, J. Hydrol., 75,
287-299, 1984.
[29] Kool, J. B., and J. C. Parker, Development and evaluation of closed form
expressions for hysteretic soil hydraulic properties, Water Resour. Res., 23
(1), 105-114, 1987.
[30] Kueper, B. H., and E. O. Frind, Two-phase flow in heterogeneous porous
media: 1. Model development,Water Resour. Res., 27 (6), 1049-1057, 1991.
172
-------
[31] Kaluarachchi, J. J., and J. C. Parker , An efficient finite element method for
modeling multiphase flow, Water Resour. Res., 25(1), 43-54, 1989.
[32] Kuppusamy, T., J. Sheng, J. C. Parker, and R. J. Lenhard, Finite-element
analysis of multiphase immiscible flow through soils, Water Resour. Res.,
23(4), 625-631, 1987.
[33] Lai, Y., A. Hadjidimos, E. N. Houstis, and J. R. Rice, General interior collo-
cation methods for second-order elliptic partial differential equations, Applied
Numerical Mathematics, 16, 83-200, 1994.
[34] LAPACK driver routine DGBSV, version 1.1, Univ. of Tennessee, March
1993.
[35] Land, C. S., Calculation of imbibition relative permeability for two- and
three-phase flow from rock properties, Trans. Am. Inst. Min. Metall. Pet.
Eng., 243, 149-156, 1968.
[36] Lenhard, R. J., and J. C. Parker, A model for hysteretic constitutive rela-
tions governing multiphase flow: 2 Permeability-saturation relations, Water
Resour. Res., 23(12), 2197-2206, 1987.
[37] Lenhard, R. J., and J. C. Parker, Experimental validation of the theory of ex-
tending two-phase saturation-pressure relations to three-fluid phase systems
for monotonic drainage paths, Water Resour. Res., 24(3), 373-380, 1988.
[38] Lenhard, R. J., J. H. Dane, J. C. Parker, and J. J. Kaluarachchi, Mea-
surement and simulation of one-dimensional transient three-phase flow for
monotonic liquid drainage, Water Resour. Res., 24(6), 853-863, 1988.
[39] Lenhard, R. J., J. C. Parker, and J. J. Kaluarachchi, A model for Hysteretic
constitutive relations governing multiphase flow: 3. Refinements and numer-
ical simulations, Water Resour. Res., 25(7), 1727-1736, 1989.
[40] Lenhard, R. J., J. C. Parker, and J. J. Kaluarachchi, Comparing simulated
and experimental hysteretic two-phase transient fluid flow phenomena, Water
Resour. Res., 27(8), 2113-2124, 1991.
[41] Lenhard, R. J., Measurement and modeling of three-phase saturation-
pressure hysteresis, J. Contaminant Hydrology, 9, 243-269, 1992.
173
-------
[42] Lenhard, R. J., M. Oostrom, C. S. Simmons, and M. D. White, Investigation
of density-dependent gas advection of trichloroethylene: Experiment and a
model validation exercise, J. Contaminant Hydrology, 19, 47-67, 1995.
[43] Leverett, M. C., Capillary behavior in porous solids, Trans. AIME, 142, 152-
169, 1941.
[44] Luckner, L., M. T. van Genuchten, and D. R. Nielsen, A consistent set of
parametric models for the two-phase flow of immiscible fluids in the subsur-
face, Water Resour. Res., 25 (10), 2187-2193, 1989.
[45] McBride, J. F., C. S. Simmons, and J. W. Gary, Interfacial spreading effects
on one-dimensional organic liquid imbibition in water-wetted porous media,
J. Contaminant Hydrology, 11, 1-25, 1992.
[46] Mendoza, C. A., and E. O. Frind, Advective-dispersive transport of dense or-
ganic vapors in the unsaturated zone: 1. Model development, Water Resour.
Res., 26(3), 379-387, 1990.
[47] Mercer, J. W., and C. R. Faust, The application of finite-element techniques
to immiscible flow in porous media, Proc., Intl, Conf. in Water Resources,
Penteuch Press, London, pages 1.21-1.57, 1976.
[48] Miller, C. T., M. M. Poirier-McNeill, and A. S. Mayer, Dissolution of trapped
nonaqueous phase liquids: mass transfer characteristics, Water Resour. Res.,
26(11), 2783-2796, 1990.
[49] Mualem, Y., A new model for predicting the hydraulic conductivity of un-
saturated porous media, Water Resources Res., 14(2), 325-334, 1978.
[50] Nielsen, D. R. and L. Luckner, Theoretical aspects to estimate reasonable
initial parameters and range limits in identification procedures for soil hy-
draulic properties, in Indirect Methods for Estimating the Hydraulic Proper-
ties of Unsaturated Soils, edited by van Genuchten et al., published by the
University of California, Riverside, 147-160, 1992.
[51] 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(4), 618-624, 1987.
174
-------
[52] Parker, J. C. and R. J. Lenhard, A model for hysteretic constitutive relations
governing multiphase flow: 1 Saturation-pressure relations, Water Resour.
Res., 23(12), 2187-2196, 1987.
[53] Peaceman, D. W., Fundamentals of Numerical Reservoir Simulation, Devel-
opments in Petroleum Science, 6, Elsevier, New York, 176 pages, 1977.
[54] Pinder, G. F., and L. M. Abriola, On the simulation of nonaqueous phase
organic compounds in the subsurface, Water Resour. Res., 22 (9), 109s-119s,
1986.
[55] Powers, S. E., C. O. Loureiro, L. M. Abriola, and W. J. Weber, Theoretical
study of the significance of nonequilibrium dissolution of nonaqueous phase
liquids in subsurface systems, Water Resour. Res., 27(4), 463-477, 1991.
[56] Powers, S. E., C. O. L. M. Abriola, and W. J. Weber, An experimental
investigation of NAPL dissolution in saturated subsurface systems: Steady-
state mass transfer rates, Water Resour. Res., 28, 2691-2706, 1992.
[57] Rabideau, A. J., and C. T. Miller, Two-dimensional modelling of aquifer re-
mediation influenced by sorption nonequilibrium and hydraulic conductivity
heterogeneity, Water Resour. Res., 30(5), 1457-1470, 1994.
[58] SLATEC driver routine DSLUGM, Lawrence Livermore National Laboratory,
Livermore, Ca., 1989.
[59] Scott, P. S., G. J. Farquhar, and N. Kouwen, Hysteretic effects on net infiltra-
tion, in Advances in Infiltration, American Society of Agricultural Engineers,
163-170, 1983. Sleep, B. E. and J. F. Sykes, Modeling the transport of volatile
organics in variably saturated media, Water Resour. Res., 25(1), 81-92, 1989.
[60] Sleep, B. E., and J. F. Sykes, Modeling the transport of volatile organics in
variably saturated media, Water Resour. Res., 25(1), 81-92, 1989.
[61] Sleep, B. E., and J. F. Sykes, Compositional simulation of groundwater con-
tamination by organic compounds: 1. Model development and verification,
Water Resour. Res., 29(6), 1697-1708, 1993a.
[62] Sleep, B. E. and J. F. Sykes, Compositional simulation of groundwater conta-
mination by organic compounds: 2. Model applications, Water Resour. Res.,
29(6), 1709-1718, 1993b.
175
-------
[63] Sleep, B. E., A method of characteristics model for equation of state com-
positional simulation of organic compounds in groundwater, J. Contaminant
Hydrology, 17, 189-212, 1995.
[64] Spillette, A. G., J. G. Hillestad, and H. L. Stone, A high-stability sequen-
tial solution approach to reservoir simulation, Soc. Pet. Eng. 48th Annual
Meeting, Las Vegas, Nev., S.P.E. Paper No. 4542, 14 pages, 1973.
[65] van Genuchten, M. Th., A closed form equation for predicting the hydraulic
conductivity of unsaturated soils, Soil Sci. Soc. of Am. J., 44, 892-898, 1980.
[66] Van Geel, P. J., and J. F. Sykes, Laboratory and model simulations of a
LNAPL spill in a variably-saturated sand medium: 1. Laboratory experiment
and image analysis techniques, Journal of Contaminant Hydrology, 17(1), 1-
26, 1995a.
[67] Van Geel, P. J., and J. F. Sykes, Laboratory and model simulations of a
LNAPL spill in a variably-saturated sand medium: 2. Comparison of labo-
ratory and model results, Journal of Contaminant Hydrology, 17(1), 27-53,
1995b.
[68] Wilson, J. L., S. H. Conrad, W. R. Mason, W. Peplinski, and E. Hagan,
Laboratory investigation of residual liquid organics from spills, leaks, and
the disposal of hazardous wastes in groundwater, EPA/600/6-90/004, April
1990.
176
-------
A. PARAMETER LIST
- the curve shape fitting parameter for curve-type /
- the curve shape fitting parameter for water drainage curves
- the curve shape fitting parameter for water imbibition curves
% - the NAPL- water mass transfer rate coefficient [1/T]
™ - the NAPL-water mass transfer rate coefficient [1/T]
' ^ne NAPL species water-gas mass transfer rate coefficient [1/T]
D - the gas-phase dispersion tensor, [L2/T] (equation 4.25)
Dw/ - the water-phase dispersion tensor, [L2/T] (equation 4.25)
e - the blending parameter which governs how fast the phase becomes entrapped
during drainage flow conditions or released from entrapment during imbibi-
tion flow conditions (e > 0) [see equation 5.7]
E™' represents dissolution mass transfer of the NAPL species from the NAPL
phase to the water phase;
E^IW represents volatilization mass transfer of the NAPL species from the water
phase to the gas phase;
E^ represents volatilization mass transfer of the NAPL species from the NAPL
phase to the gas phase;
E^,w represents adsorption mass transfer of the NAPL species from the water
phase to the soil.
/ - the curve-type index, takes on a value from 1 to 6 (see Table 5.1)
177
-------
F - the maximum number of curve-types considered (= 6, see Table 5.1)
foe - the mass fraction of organic carbon [dimensionless]
fa - a—phase fractional flow function (equation 6.17)
g - gravity [L/T2]
H - Henry's law coefficient
hc - the capillary pressure head: hc = Pc/(p^ g)
k the intrinsic permeability scalar [L2]
k (as a superscript) - indicates the time level of parameter evaluation.
k* - the intrinsic permeability magnitude [L2] of the soil used to measure the S-P
model curve fit (see equation 5.41)
Kd - the distribution coefficient [L3 /M] (equation 4.36)
Koc - the organic carbon partition coefficient [dimensionless] (equation 4.36)
kra - the relative permeability of the a—phase, a = W, N, G [dimensionless]
krQ - the relative permeability of the gas-phase [dimensionless]
krff - the relative permeability of the NAPL-phase [dimensionless]
krw - the relative permeability of the water-phase [dimensionless]
m = 1 — 1/77 - a fitting parameter used in the k — S — P model
m (as a superscript) - indicates the iteration level of parameter evaluation.
MDC - the main drainage curve
MIC - the main imbibition curve
PDC - the primary drainage curve
PIC - the primary imbibition curve
Pa - pressure of the a-phase, a = W, N, G [F/L2]
178
-------
PCNW - capillary pressure for the NAPL-water system [F/L2]
PCGN capillary pressure for the gas-NAPL system [F/L2]
Qa - point sources (+) or sinks (-) of the a— phase, a = W, N, G [1/T]
t> ( Q GmiiA ~1 /i omiiA
Ka = (dar - bed ) ~ I1 ~ ^crf J
Se - the effective wetting-phase saturation for a given S — P curve (equation 5.10)
Se(f) - the effective water saturation for a specific S — P curve type / (equation
5.13)
S'g - the effective saturation as computed from equation 5.13 given hc and a^
(used during parameter update)
Sec - the effective gas saturation used to define the kr(j functional (equations
5.26 and 5.28
°'Sea - the effective gas saturation used to define the connectivity term of the kr(j
functional (see equation 5.23)
bSea - the effective gas saturation used to define the integral of the kra functional
(see equation 5.23)
Sepf - the effective NAPL saturation used to define the kr^ functional (equations
5.30 and 5.31)
- the effective total wetting-phase saturation used to define the kr^ func-
tional (equations 5.32, 5.33 and 5.34)
- the effective total nonwetting-phase saturation used to define the kr^
functional (equations 5.35, 5.36 and 5.37)
Sew - the effective water saturation used to define the krw functional (equations
5.25 and 5.27)
aSew - the effective water saturation used to define the connectivity term of the
krw functional (see equation 5.19)
bSew - the effective water saturation used to define the integral of the krw func-
tional (see equation 5.19)
179
-------
Sr - the minimum saturation for a given S — P curve
SV(/) - the minimum saturation for a specific S — P curve type /
Ss - the maximum for a given S — P curve
Ss(f) - the maximum saturation for a specific S — P curve type /
Scfr - a fitting parameter representing the maximum residual gas-phase saturation
Sffr - the residual NAPL-phase saturation computed from equation 5.24
Sffnrj a fitting parameter representing the maximum residual NAPL-phase sat-
uration in a two-phase NAPL-water system (NAPL as a nonwetting phase)
SNWT, a fitting parameter representing the maximum residual NAPL-phase sat-
uration in a two-phase NAPL-gas system (NAPL as a wetting phase)
SWT - a fitting parameter representing the maximum residual water-phase satu-
ration
Sa - the a—phase saturation, a, = W, N, G
Saf - the a—phase saturation which is free to flow
Sat - the a—phase saturation which is trapped and unable to flow
5"™11 - the lower limit of entrapped a-phase (equation 5.8)
S'ar - the magnitude of the residual a-phase, 0 < 5^r < Sar, at the terminus of
an a-phase drainage process (equation 5.9)
^™ax ~ the highest a-phase saturation that has occurred since it was last at im-
mobile residual conditions, i. e., the maximum imbibed a-phase saturation
which is available for displacement
Sar(j) - the magnitude of the residual at the origin of a main curve
SQT - the maximum Sot value for the system
STW - total wetting phase saturation = Sw + «SW
va - mass average a—phase velocity vector [L/T], a = W, N, G (equation 4.21)
180
-------
VT = vw/ + v^ + VG - the total velocity
co an+1 on
Oow — &W ~ &W
f32 - dimensionless fitting parameter for mass transfer
/33 - dimensionless fitting parameter for mass transfer
0l - the rate coefficient for NAPL-water mass transfer, [1/T]
/3GN - the rate coefficient for NAPL-gas mass transfer, [1/T]
0l ' - the rate coefficient for gas-water mass transfer, [1/T]
C - the pore connectivity parameter for the krW functional [equation 5.19]
(p - the pore connectivity parameter for the kra functional [equation 5.23]
£ - the pore connectivity parameter for the krfj functional [equation 5.38]
77 - a curve fitting parameter for the S — P model [see equations 5.10 and 5.13]
e - porosity of the porous medium [dimensionless]
e* - the porosity of the soil used to measure the S-P model curve fit (see equation
5.41)
Y* _ potg jg ^g Specinc weight of the phase [M/T2], a = W, N, G
Fjf - defines the change in pw due to the presence of NAPL species (equation
6.12)
FG - defines the change in pG due to the presence of NAPL species (equation
6.13)
AQ = kkra/n™ (a = W, N, G), is the a—phase mobility scalar [L^/FT]
Ha -the phase viscosity [FT/L2]
tuf is the mass fraction of the adsorbed NAPL on the solid [dimensionless]
pa - a-phase mass density [M/L3] a = W, N, G
ps - the density of the soil [M/L3]
181
-------
p^ - the mass concentration of NAPL in the water-phase, [M/L3]
Pn - the mass concentration of NAPL in the gas-phase, [M/L3]
p" - the mass concentration of a species i in the a—phase, [M/L3]
~p™ [M/L3] - the equilibrium concentration of the NAPL species in the water
phase (solubility limit)
~p~n [M/L3] - the constant equilibrium vapor concentration of the NAPL species
in the gas-phase (vapor solubility limit)
pf - the source or sink of mass for a species i in the a-phase [M/L3T] due to
interphase mass exchange (i. e., dissolution, volatilization and adsorption).
A o/ ., — <-va — V3
^ Ja/3 — I I
GGW ' the interfacial tension between the gas and water phases
GGN - the interfacial tension between the gas and the NAPL phases
&NW - the interfacial tension between the NAPL and water phases
182
-------
B. PARTICULARS OF HERMITE
COLLOCATION
This appendix contains details regarding the three-dimensional, Hermite colloca-
tion formulation.
B.I. Nodal Degrees of Freedom
Each node has eight degrees of freedom (DF), and the vector representing the
values, U, is numbered as:
12345678
u ux uy uz uxy uxz uyz uxyz (0-L)
where the subscript [ , (•)] represents the partial derivative with respect to (•).
B.2. Basis Function Definition
The Cl continuous Hermite cubic interpolation polynomials are denned on a gen-
eral one-dimensional interval (£o)£i) as:
defined at £ = £0
1 - 2-r3- ) defined at f = ^
defined at £ = £0
defined at £ = ^
(B.2)
where there are four functions defined on each one-dimensional interval (element).
183
-------
The three dimensional version of the Hermite cubic can be derived from the
tensor product ordering of the one-dimensional Hermite basis functions. At each
of the i nodes in a three-dimensional grid block, i = 1,2, ...8, the following eight
Hermite polynomials are denned:
where the ordering is consistent with the degree of freedom ordering in B.I.
B.3. Hermite Interpolation of Capillary Pressure
Given capillary pressure as a function of saturation, and saturation interpolated
by Hermite cubics, then the Hermite cubic interpolation coefficients of PC(S) are
generated using the chain rule, i. e.,
4
P'S.,
5
p'o
-1 ^
6
p'o
-<££
7
P'S.
yz
p'
1
where P'C = dPc/dS, and terms 5 through 8 have been simplified by assuming
the chained terms which include the derivatives d2Pc/dS2 and d3Pc/dSs are small
compared to the other terms and can thus be neglected.
B.4. Boundary Condition Specification
The set of DF's which are specified at a boundary node is defined as a function of:
node-type and boundary condition-type. There are three boundary-node-types to
be considered:
node- typo
rnidsido (on a plane)
edge (intersection of two planes = line)
corner (intersection of three planes = point)
number of
DP specified
4
6
7
number of
DP free
4
2
1
184
-------
The following lists the DF's (referenced by the nodal degree of freedom num-
bering system shown in B.I) that are specified for each node-type and boundary
condition-type (Dirichlet or Neumann).
Midside Nodes (lie in a plane) - The DF's specified depend on the plane and the
boundary condition imposed:
x-y plane y-z plane x-z plane
Dirichlet
Neumann
1, 2, 3, 5
4, 6, 7, 8
1, 3, 4, 7
2, 5, 6, 8
1, 2, 4, 6
3, 5, 7, 8
Edge Nodes (lie along a line) - The DF's specified depend on what boundary
condition is imposed on each intersecting plane:
• When the edge is parallel to the x-axis (intersection of the x-y and x-z
planes):
Dirichlet on x-y side Neumann on x-y side
Dirichlet on x-z side
Neumann on x-z side
1, 2, 3,4, 5, 6
1, 2, 3, 5, 7, 8
1, 2, 4, 6, 7, 8
3, 4, 5, 6, 7, 8
When the edge is parallel to the y-axis (intersection of the x-y and y-z
planes):
Dirichlet on x-y side Neumann on x-y side
Dirichlet on y-z side
Neumann on y-z side
1, 2, 3, 4, 5, 7
1, 2, 3, 5, 6, 8
1, 3, 4, 6, 7, 8
2, 4, 5, 6, 7, 8
When the edge is parallel to the z-axis (intersection of the x-z and y-z
planes):
Dirichlet on x-z side Neumann on x-z side
Dirichlet on y-z side
Neumann on y-z side
1, 2, 3, 4, 6, 7
1, 2, 4, 5, 6, 8
1, 3, 4, 5, 7, 8
2, 3, 5, 6, 7, 8
Corner Nodes - The DF's specified depend on what boundary condition is im-
posed on each intersecting plane. There are 8 combinations:
Dirichlet on x-y, Dirichlet on x-z, Dirichlet on y-z: 1, 2, 3, 4, 5, 6, 7
185
-------
Neumann on x-y, Dirichlet on x-z, Dirichlet on y-z: 1, 2, 3, 4, 6, 7, 8
Dirichlet on x-y, Neumann on x-z, Dirichlet on y-z: 1, 2, 3, 4, 5, 7, 8
Dirichlet on x-y, Dirichlet on x-z, Neumann on y-z: 1, 2, 3, 4, 5, 6, 8
Neumann on x-y, Neumann on x-z, Dirichlet on y-z: 1, 3, 4, 5, 6, 7, 8
Neumann on x-y, Dirichlet on x-z, Neumann on y-z: 1, 2, 4, 5, 6, 7, 8
Neumann on x-y, Dirichlet on x-z, Neumann on y-z: 1, 2, 3, 5, 6, 7, 8
Neumann on x-y, Neumann on x-z, Neumann on y-z: 2, 3, 4, 5, 6, 7, 8
186
-------
C. INITIALIZE TRAPPING PARAMETERS
Given the initial phase saturation distribution, SWQ, SNQ and SGQ, this is the
decision tree which determines the initial phase trapping parameters.
• if SNO = 0
if (SGQ < SGT) (gas is at residual)
Sat = SGQ, Swt = 0, and / = 1 for both S — P functionals
else if (SWQ < Swr} (water is at residual)
Swt = SWQ, Set = 0, and / = 2 for both S — P functionals
else (both phases are mobile)
Swt from equation 5.7, Set = 0, and / = 1 for both S — P functionals
endif
• else if SGO = 0
Set = 0, and / = 1 for STW(PCGN) functional
if (Sm < SNnr) (NAPL is at residual)
Sm = SNO, Swt = 0, and / = 1 for SW(PCNW] functional
else if (Swo < Swr} (water is at residual)
Swt = SWQ, Sfft = 0, and / = 2 for SW(PCNW] functional
187
-------
else (both phases are mobile)
Swt from equation 5.7, S^t = 0, and / = 1 for SW(PCNW) functional
endif
• else if Sjvo > 0 and SGQ > 0
if (SGQ < SGV) (gas is at residual)
Set = SGQ, and / = 1 for STW(PCGN) functional
if (Sjvo < SWnr) (NAPL is at residual)
Sm = SNQ, Swt = 0, and / = 1 for SW(PCNW} functional
else if (SWQ < «SW) (water is at residual)
Swt = SWQ, Spft = 0, and / = 2 for SW(PCNW] functional
else (both water and NAPL phases are mobile)
Swt from equation 5.7, Sm = 0, and / = 1 for SW(PCNW} functional
endif
else if (SWQ ^ Swr} (water is at residual)
Swt = SWQ and / = 2 for SW(PCNW} functional
if (SNQ < Sffwr) (NAPL is at residual)
Sm = SNO, Sa = 0, and / = 2 for STW(PCGN) functional
else (both gas and NAPL phases are mobile)
Set from equation 5.7, Sm = 0, and / = 2 for STW(PCGN} functional
endif
else (mobile water and gas after full drainage and imbibition cycles)
Swt — Swr and JGt — '--'Gr
if (SVo ^ Sffnr + Sffwr) (NAPL is at residual)
SVt from equation 5.24
188
-------
else (NAPL is mobile)
endif
endif
• endif
189
-------
D. PECLET CONSTRAINT
To derive the expressions for grid Peclet number, equations 6.35 and 6.36 must
be written in terms of the dependent variables Sw and STW by using the chain
rule for differentiation:
dt
*N\w A TWW + (fWXG),w
V • { [AWN (fw\N) (PcNW\w + AWG (fw\G) (PcGW\w
V • { [AWG (fwXG] (PCGW),TW\ VSTw}
-/H/V-qT (D.I)
dt
+ (
-------
and
( f W \ \ _ fW(\ \ I fN(\ \ fWfN
(J *N),W — J (AN),W + J (*w),w ~ J J
1 f W \ \ _ fW(\ \ fWfN(\ \
(J *N),TW — J (^N),TW — J J (AT),TW
(fGXN\TW = f° (\tf\Tu, + fN(^G\TW ~ fGfN(
(fW\G\W = (fG\V\W = f°(\W\W ~ f°fW
f fW \\ _ ( fG\ \ _ fW f \ \ f G fW
(I *G),TW — (I AW),TW — J (^G),TW - I I
(AT))H/ = (AH/))H/ + (\N),w
(XT),w = (^G},TW + (^N),TW
The derivative definitions (D.3) are chosen such that:
(fW),a + (fN),a + (fW),a = 0 , « = W,
which is necessary in order that equations D.I and D.2 sum to the NAPL balance
equation. In addition the derivatives are evaluated using numerical differentiation,
where e is a small positive number of order 10^6.
Finally, the terms (Pccw},w and (PCGW),TW arc computed from equation 4.14
as:
cGW),Tw = - ,-T cGN),Tw
\(?GN
191
-------
E. SOURCE FILE DESCRIPTION
The following provides a list of the FORTRAN program routines which make up
the simulator and a brief description of what each routine does. Programs are
listed in alphabetical order
BASIS.F: Set up arrays for Hermite and Lagrange basis functions evaluated at
local orthogonal collocation points.
Subroutine BASIS_3D: 3-D Hermite and Lagrange basis functions.
Subroutine BASIS_2D: 2-D Hermite and Lagrange basis functions.
BCSET.F - Set up boundary conditions.
Subroutine BCSET: Set up default Neumann conditions for all dependent
variables.
Subroutine BC_FLOW: Set up Dirichlet conditions for phase flow equa-
tions.
Subroutine BC_OA: Set up Dirichlet conditions for dissolved NAPL species
contaminant transport equations.
Subroutine BC_OG: Set up Dirichlet conditions for NAPL vapor species
contaminant transport equations.
Subroutine BC_UP: Update nonlinear boundary conditions after each time
step.
Subroutine PR_BC: Set flux terms for gas and NAPL pressure conditions
after total flow solution.
Subroutine PR_SAVE: for flow boundary conditions 2 and 3 save the PG
and PN values at the beginning of each time step.
Subroutine NO_FLOW: For no flow boundary nodes, set the S — P curve
slope to a high value to mimic a linear no flow condition.
192
-------
DIRECT_SOLVE.F The LAPACK library routine dgbsv.f (LU with partial piv-
oting)
DSLUGM.F The SLATEC library routine of the same name (incomplete LU
preconditioned GMRES)
EQ_NUMBER.F Number the non-boundary data degrees of freedom for each
dependent variable (column of the system matrix)
Subroutine DF_NUM_P: Degree of freedom numbering for Pw/
Subroutine DF_NUM_S: Degree of freedom numbering for Sw and STW
Subroutine DF_NUM_T: Degree of freedom numbering for the NAPL
species.
Subroutine NUMBER: Set the number for the degrees of freedom.
EXCHANGE.F Compute the mass exchange terms for the water and gas phases.
Subroutine EXCH_W: exchange terms for the water-phase
Subroutine EXCH_G: exchange terms for the gas-phase
FUNCTIONS.F Set nodal values of all functions of the dependent variables
Subroutine FLOW_FUN Set nodal values of all functions of saturation
Subroutine WATER_PROP Compute water-phase properties based on
composition.
Subroutine GAS _ PROP Compute gas-phase properties based on compo-
sition.
Subroutine FUN_AT_N Set nodal values of all functions dated at the
current time level.
Subroutine VEL_W Compute the nodal values of water velocity for the
transport equation
Subroutine VEL_G Compute the nodal values of water velocity for the
transport equation
Subroutine VEL_N Compute the nodal values of NAPL velocity for mass
balance and output
193
-------
Subroutine V_TOT Compute the nodal values of total velocity for the
saturation equations
Subroutine DISP_W Compute the dispersion tensor for the water trans-
port equation
Subroutine DISP_G Compute the dispersion tensor for the gas transport
equation
Subroutine DFDS3 Compute the terms which make up the Peclet and
Courant constraints.
HYST.F Set nodal values of all the hysteresis variables
Subroutine HYST_IC Initialize hysteresis variables
Subroutine TRAP_UP Update nodal values of trapping variables
Subroutine SW_PC Update nodal values denning SW(PCNW} functionals
Subroutine ST_PC Update nodal values denning STW(PCGN) functionals
ICSET.F: Set up initial conditions for the simulation.
Subroutine ICSET Set the global initial conditions for all the variables.
Subroutine IC_SAT Set the node-specific initial conditions for the satura-
tion variables.
Subroutine IC_ROA Set the node-specific initial conditions for the dis-
solved NAPL species variables.
Subroutine IC_ROG Set the node-specific initial conditions for the NAPL
vapor species variables.
INCLUDE.F: Account for all parameter definition, array dimensioning, and com-
mon block definition.
MAIN.F: Driver program taking care of I/O for the simulator.
MBAL.F: Mass balance computations.
Subroutine MBAL Determine mass in the system based on Gauss quadra-
ture integration of the solution.
194
-------
Subroutine MBALB Determine the mass in and out of the system along
the boundaries of the domain.
Subroutine PORE_VOL Compute the pore-volume of the domain.
PE_CO.F Compute the Peclet and Courant constraints.
Subroutine MP_DIFF Compute the Peclet and Courant constraints.
Subroutine DT_CNTRL Time step control.
POINT_SOURCE.F: Definition of source and sink terms.
Subroutine PT_SRCE: Read from input file and set initial nodal point
source/sink data
Subroutine QOUT: For outflow wells, distribute the rate between phases
based on the saturation solution at the node where the well is idealized.
Subroutine QOUT_ADJ: For outflow wells defined from case 2 and 3 flow
boundary conditions. Distribute the rate between phases based on the
saturation solution at the node where the well is idealized.
Subroutine WELL_MASS: Compute contribution to mass balance from
point source/sink data
POINTER.F Generates pointer vectors to map between default numbering and
shortest-direction- first numbering
Subroutine POINT Set the pointer arrays
Subroutine RELEM Generate element index list
PW_SOL.F: Driver for the total flow equation solution algorithm.
RESTART.F: Setup to restart the time step because of non-convergence
SAT_SOL.F: Driver for the saturation equation solution algorithm.
Subroutine SAT_SOL: Takes care of I/O for saturation solution and com-
putes convergence
Subroutine SW_SOL: Driver for the Sw equation solution algorithm
195
-------
Subroutine ST_SOL: Driver for the STW equation solution algorithm
SATPAR.F: Compute the functions of saturation: capillary pressure and relative
permeability
Subroutine PC_GN: Compute the Sw(Pcww} functional.
Subroutine PC_NW: Compute the STW(PCGN) functional.
Subroutine SFUNKW: Compute the krw(Sw} functional
Subroutine SFUNKN: Compute the krN(SwiSc) functional
Subroutine SFUNKG: Compute the krQ(Sa) functional
Subroutine LEV_SW_P: Interfacial tension scaling.
Subroutine LEV_ST_P: Interfacial tension scaling.
SYS_OA.F Set up and solve system of equations for the dissolved NAPL species
contaminant transport equation.
Subroutine SYS_OA_I: Version for iterative solver.
Subroutine SYS_OA_D: Version for direct solver.
SYS_OG.F Set up and solve system of equations for the NAPL vapor species
contaminant transport equation.
Subroutine SYS_OG_I: Version for the iterative solver.
Subroutine SYS_OG_D: Version for the direct solver.
SYSTEM_G.F Set up and solve the system of equations for the linearized STW
equation.
Subroutine SYSTEM_GI: Version for the iterative solver.
Subroutine SYSTEM_GD: Version for the direct solver.
SYSTEM_P.F Set up and solve the system of equations for the total flow equa-
tion.
Subroutine SYSTEM PI: Version for the iterative solver.
196
-------
Subroutine SYSTEM_PD: Version for the direct solver.
SYSTEM_W.F Set up and solve the system of equations for the linearized Sw
equation.
Subroutine SYSTEM_WI: Version for the iterative solver.
Subroutine SYSTEM_WD: Version for the direct solver.
TRAN_SOL.F: Driver for the contaminant transport equation solution algo-
rithm.
Subroutine TRAN_SOL: Takes care of I/O for the transport solution and
computes convergence
Subroutine ROA_SOL: Driver for the p% equation solution algorithm
Subroutine ROG_SOL: Driver for the /)„ equation solution algorithm
WRITE_OUT.F Takes care of most file I/O
Subroutine READ_RS: Read restart files
Subroutine WRITE_RS: Write restart files
Subroutine PRINT: Write to files *.out
Subroutine GRAPH: Write to files for animated graphical display
Subroutine JAQ_SET: Write grid and soil definition for animated graph-
ical display
Subroutine ECHO: Echo data input to file echo.out
197
-------
F. NAPL PROJECT
F.I. Purpose
The purpose of this project is to:
1. build and calibrate a simulation model which approximates the transient
behavior of a two-dimensional, three-phase DNAPL flood experiment;
2. obtain a sense of how certain physical parameters affect the solution;
3. obtain a sense of how discretization errors affect the solution.
In addition, the analysis thus conducted can be used to verify whether the
mathematical representation of three-phase flow through porous media that is
used in the simulator represents reality. In other words, is the simulator a sur-
rogate for reality or are there aspects of the mathematical interpretation of the
physics which are either erroneous or simply incomplete?
F.2. Scope
An artificial aquifer experiment was conducted by Mikhail Fishman at the EPA's
Subsurface Protection and Remediation Division of the National Risk Manage-
ment Research Laboratory in Ada, OK (formerly RSKERL) to gather quantitative
and qualitative data on DNAPL migration through a variably saturated homo-
geneous sand. The DNAPL used in the experiment is called tetrachloroethylene
(PCE, a common chlorinated hydrocarbon used in the dry cleaning industry). As
will be detailed in subsection F.3, the data from the experiment consists of several
types:
• soil and fluid properties;
• moisture retention data for the sand used in the experiment;
198
-------
• experimental initial and boundary conditions;
• DNAPL influx data;
• video images of the box at various points in time showing the areal extent
of the DANPL which is dyed to maximize contrast.
At this point it may be clear that this information represents only a subset of
that required to model the experiment. Thus, part of this project is to identify
the data requirements for the model, to generate a cause and effect relationship
between specific parameters and simulation results, and to quantify the physical
limits of specific parameters.
The overall project has been separated into a series of Tasks, with each sub-
sequent Task using information generated from the previous one. Each Task
highlights an important aspect of model development, and leads the user in a
sequential manner to developing the full three-phase model. The Tasks are sum-
marized as follows:
Task 1 Define the S — P model fitting parameters ad, r\ and Swi-
Task 2 Define appropriate model discretizations (i. e., node spacing and time
step). Identify and quantify those input parameters which are used to ensure
compatibility between the physical and the discrete problems.
Task 3 Set up the two-dimensional model for the DNAPL flood. Define appro-
priate boundary conditions. Identify a subset of input parameters which
can be used to fit the model to the experimental data. Estimate a physi-
cal range of acceptable values for these data, and qualify how the solution
should respond to a change in a particular parameter or set of parameters.
Identify how the experimental data can be used to calibrate the model.
Task 4 Model the experiment and calibrate the model.
F.3. Experimental Setup and Data Base
A photograph of the actual experimental apparatus is provided in Figure F.I,
and an idealization showing dimensions is provided in Figure F.2. The inside
dimensions of the box (defining the volume of the sand) are 67 cm deep, 49 cm
wide, and 2 cm thick. The top boundary is open to the atmosphere and the bottom
199
-------
boundary is impermeable. The vertical sides are constructed such that water can
flow across the boundary but not air. A constant phreatic surface is defined by
specifying appropriate water source/sink ports along the vertical sides. The box is
filled with a uniform medium-grained sand in a manner which is assumed to result
in a homogeneous, isotropic porous medium. The soil properties are reported to
be:
permeability
porosity
3.5:rl(r7cm2
0.37
The sand-filled box is imbibed with water to the top and allowed to equilibrate to
create an initial condition where the system is in static equilibrium and Sw = 1
throughout. The following fluid properties are provided:
density (g/cm3)
viscosity (p)
water
1.0
0.01
NAPL
1.626
0.0093
gas
0.00129
0.0002
interfacial tension (dynes / cm)
(jNW = 39.5 (7GN = 31.74 OGW = 72.75
From this initial condition, three sequential displacement experiments are run:
1. The phreatic surface is lowered to elevation 35.5 cm from the top of the box
and, the system is allowed to return to equilibrium conditions. The data
reported which defines the moisture content as a function of depth is given
in Table F.I.
2. Given the initial condition from Part 1, the PCE source is applied as shown
in Figure F.2. That is, a 0.5 cm head of PCE is applied uniformly over
a 10 cm2 surface at the center/top of the box until 200 cm3 enters the
domain. Note that for the experiment this took 143 seconds. For this
forcing condition, Table F.2 provides the data which relates the cumulative
volume of PCE entering the domain to the elapsed time. Other information
qualifying this phase of the experiment is a series of video images of the front
face of the box similar to the one shown in Figure F.I. Video images are
available for the following times (elapsed time in seconds since the DNAPL
was applied: 5, 14, 62, 143.
3. Given the initial condition from Part 2, that is, the data at time = 143s, the
PCE source is removed and the system is allowed to return to equilibrium for
a period of 3452 seconds (total elapsed time since the DNAPL was applied is
200
-------
depth from top (cm)
0
2.5
5
7.5
10
12.5
15
17.5
20
22.5
25
27.5
30
32.5
35
67
moisture content
.0514
.063
.0788
.1273
.1822
.2543
.2927
.3088
.3171
.3237
.3229
.3292
.3476
.3619
.37
.37
Table F.I: EXPERIMENTAL DATA - MOISTURE CONTENT AS A FUNC-
TION OF DEPTH
ELAPSED TIME (s)
0
1.88
10.84
23.12
31.41
60.29
83.47
102.72
113.41
143.0
> 143.0
CUMULATIVE VOLUME IN
(cm3)
0
22.0
67.5
87.5
101.0
145.0
165.5
179.5
187.0
200.0
200.0
Table F.2: EXPERIMENTAL DATA - PCE VOLUME ENTERING THE DO-
MAIN AS A FUNCTION OF TIME
201
-------
Figure F.I: Photo of the experimental apparatus just after the PCE was removed,
showing dimensions and vertical constant head boundaries. The PCE (dark grey)
has been dyed red to maximize contrast.
3595 seconds) . Information qualifying this phase of the experiment consists
of a series of video images for the following times (elapsed time in seconds
since the DNAPL was applied: 285, 185, 1195, 1795, 3595.
F.4. TASK 1
F.4.1. Purpose
Given the initial static moisture profile as defined in Table F.I, determine the
appropriate S — P model parameters: a^,?7, Swi- That is, fit the van Genuchten
S — P model (equation 5.10) to the experimental data.
202
-------
open
boundary
center
line
2 cm
DNAPL Source (PCE)
0.5 cm head
for 143 seconds
(total volume = 200 ml)
22cm
31.5cm
Figure F.2: Idealization of the experimental setup superimposed on a video image
of the box (5 seconds after the PCE source was applied).
203
-------
= 30
top. of im?g?.
dashed (computed): a= .I/cm, n = 2, Swr = .22
solid (exp. data normalized by porosity = 0.37)
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
water saturation
Figure F.3: Plot of the initial static saturation profile. The computed curve is fit
to the experimental data by altering the S-P model curve fitting parameters.
F.4.2. Procedure
Use a trial-and-error sequential approach wherein you systematically vary the
parameters 0^,77, and Swi until a qualitatively "good" fit is obtained. Equation
5.10 is used because the experiment represents a uni-directional displacement
process.
For each parameter set chosen, superimpose the plots of equation 5.10 and the
data from Table F.I. An example of such a plot is provided in Figure F.3.
This is a heuristic analysis, that is, you should try to fit the model to the
data by iteratively choosing different parameter sets and "eye-balling" the fit.
Given some knowledge of how the model is affected by the three parameters,
your analysis should converge quickly. In addition, please note that there is no
'correct' answer, and that this data represents only a piece of the calibration
puzzle. Therefore, it behooves you to generate a suite of parameter sets which
give a good fit so that you can make more meaningful decisions when trying to
fit the model to the three-phase flow problem.
204
-------
F.4.3. Results
In addition to providing a plot of the fitted data, please answer the following
questions:
1. What type of S — P curve is defined by this data?
2. From a qualitative perspective, how is the S — P curve affected by the
parameters a^ and 77?
3. What are the ranges of the values of a
-------
pG = patm
Figure F.4: Illustration of the one-dimensional water-gas displacement experiment
showing boundary and initial conditions and dimensions. The mesh has one 2 cm
element in the horizontal direction, and the vertical direction is discretized in an
appropritate manner.
appropriate discretization to resolve the vertical direction. An illustration of the
computational domain including initial and boundary conditions is provided in
Figure F.4. Run the simulation for 5000 seconds (this should be near steady
state) .
Focus your attention on the following aspects of model input:
1. GRID DEFINITION - The grid must be able to resolve not only the steady-
state moisture profile, but also the transient gas-phase front as the water
drains from the soil. A general rule of thumb is that the front should be
resolved in by about four elements (an increase in resolution will have little
effect on the solution, and a decrease in resolution will result in oscillations).
Associated with grid definition is the Peclet number (file sm.in input line
51). By using a small enough value of Pecnt (e. g., Pecnt < 2), one alters
the problem by increasing the amount of the capillary force, thereby fitting
the physical problem to the fixed grid.
2. TIME STEP CONTROL - The size of the time step is defined in file sm.in,
input lines 8 through 15. Time step definition is based on:
206
-------
the initial time step - must be relatively small in order to resolve the upper
boundary condition which is a shock to the system.
the number of iterations required for convergence over any given time step - there
is a direct relationship between the time step size and the number of
iterations it takes to converge on the nonlinear problem.
the Courant constraint - a check on the time truncation error.
the change in effective saturation over the time step - the effective saturation is
a normalized saturation, and for 'tight' scanning curves small changes
in saturation can lead to large changes in effective saturation.
the maximum and minimum specified values - the time step cannot be too
large or too small because at these extremes round-off errors can dom-
inate the solution and cause instabilities.
3. K-S-P CURVE DEFINITION - You already have the parameters ad, 77, and
Swi from Task 1. The remaining parameters which constitute the k — S — P
model are denned in file sm.in, input lines 33 through 50.
F.5.3. Results
Please respond to the following:
1. From the screen output, note the initial and final total flux across the hori-
zontal boundaries.
2. Run the model with different uniform grids using no additional diffusion
(i.e., set Pecrit > 100), use Az = 10cm, Az = 5cm and Az = 2.5cm.
Note how the model resolves the moisture profile (if it does). Is the solution
converging to something? How do these profiles compare the one defined in
Task 1?
3. Try running the model with added diffusion, Pecnt = 2, and Az = 5cm.
See how the solution is altered.
4. Write Darcy's law for the water and gas phases for the case of static equi-
librium.
5. What role does the relative permeability play in defining the static moisture
profile?
207
-------
6. What role do the following parameters play in this experiment: S?fnr, S?fwr, S(jr
and
-------
1. set up the two-dimensional model for the DNAPL flood experiment. Specif-
ically, you need to define the following modeling parameters:
• appropriate boundary conditions,
• discretization in the horizontal dimension,
• the remaining parameters which define the hysteretic k — S — P model1,
i. e.,
«i, SNnr, SNwr, SGr, <;, £, and (p
2. identify which experimental data can be used to calibrate the model and
how it can be used; and
3. identify a subset of physical parameters for which the solution is most sen-
sitive and with which you can calibrate the model.
F.6.2. Procedure
As introduced in Section F.3, think of the DNAPL experiment model as a series
of three sequential sub-models, i. e.,
1. Starting with the initial condition of full water phase saturation, drop the
water table to match the experimental condition, and allow the system to
approach steady-state conditions. This is the 2-D version of the model you
set up for Task 2.
2. Given the initial conditions from sub-model 1, apply the DNAPL source for
the specified time period (i. e., 143 seconds).
3. Given the initial conditions from sub-model 2, remove the DNAPL source
and allow the system to re-equilibrate for 3452 seconds.
Let us itemize the procedure into three steps: (1) define the horizontal dis-
cretization, (2) set up the appropriate boundary conditions for the three sub-
experiments, and (3) get the feel of the physics.
1 One should estimate a physical range of acceptable values for these data, and qualify how
the solution should respond to a change in a particular parameter or set of parameters.
209
-------
-r
? =
C>i '
ci—
P
o
no flow
P-
C/5"
(Jl
§
no flow
no flow
9
O
23
O
Case 5 Dirichlet data
use hw, S\v and SG from solution of Experiment 1
no flow
Case 5 Dirichlet data
use hw, S\y and SG from solution of Experiment 2
O
II
»
3
;.'
bd ^>
•8-
3 ^-K
-sr
•
•
¥
1 1
1
1
no flow
O
o
o
W
TS
2
3'
rs
^^>
s>
n
s-
n
x
5"
o'
9
i— '
T-,
O
II
3
»«
1'
rt-
h^
I
H
TS
fB
ft
rt-
u
I
H
II
in
I
Lft
I
H
00
(Jl
Figure F.5: Illustration of the 2-D model setup for the three sequential simu-
lations, showing boundary and initial conditions and the time frame for each
experiment.
210
-------
Horizontal Discretization
There are two important issues to consider with respect to spatial discretization:
computational efficiency and numerical accuracy. As the number of grid nodes
increases linearly, the computational effort increases exponentially. Therefore,
when designing the mesh consider the following:
1. The experimental setup is symmetrical about the center line, and you should
take advantage of this fact from a computational efficiency point of view.
2. The mesh should be refined in the vicinity of the NAPL source so that the
model can accurately represent it. The mesh can coarsen away from the
source.
3. In general the mesh scale in the horizontal dimension can be coarser than
that in the vertical dimension because gravitational forces are absent and
capillary forces dominate (i. e., the saturation front is more diffused in the
horizontal that in the vertical because of gravitational effects).
Figure F.5 provides a diagrammatic representation of the first two considera-
tions, where the three experiments shown refer to the three sequential sub-models
referenced above.
The only physical constraint with respect to defining the horizontal discretiza-
tion, aside from defining the appropriate width of the domain, is that relating to
the definition of the NAPL source area. If we realize that the 2-D model has unit
length in the omitted dimension, then the area of the surface associated with a
node normal to the direction from which the condition is applied is equal to 1/2
the length of the elements surrounding the node times unit length in the omitted
dimension. To mimic the experimental conditions then, the grid spacing at the
source must be set such that the NAPL is applied over an area equal to 2.5cm2.
An example of an appropriate horizontal mesh is provided in Figure F.6.
Boundary Conditions
The model boundary conditions should be set up to mimic the experimental
boundary conditions to the extent possible. This task may not be so straight-
forward especially when the physical boundary conditions do not match the sim-
ulator capabilities.
211
-------
§
a
- hG/
'
1 cm
1 cm
1 cm
2 cm
2 cm
Figure F.6: Illustration of the uniform distribution of forcing conditions about a
node. For this example, NAPL is ponded with head HN over an area equal to 2.5
cm
For the sequence of experimental conditions described above, consider the
example of boundary condition specification shown in Figure F.5. Let us look at
each experiment in turn. In Experiment 1 the vertical sides are set as no flow
boundaries2, and for the conditions specified on the top and bottom sides, the
problem becomes essentially 1-D (in fact, the same one you solved in Task 2).
With regard to the setup for Experiment 2, in order to make the model condi-
tions mimic the experimental conditions, first the bottom side is set as a no flow
boundary. Second, the right-vertical side is opened to flow by assuming that the
data for the dependent variables (Pw, SV, •S'c) from Experiment 1 at each of these
boundary nodes will not change over the course of the second simulation. That is,
the data from Experiment 1 at each of these boundary nodes becomes Dirichlet
data. Note, by specifying Dirichlet data at a node for a given variable we are
effectively solving for the gradient normal to the boundary of that variable. Using
Darcy's law to compute flow, the boundary is therefore open to flow of the phase
if the phase is mobile (i. e., the relative permeability is nonzero). By making this
side a Dirichlet boundary, we are assuming that it is far enough away from the
action so that changes in S^v, ^G and 5jv remain effectively zero and that changes
in vertical gradients are negligible..
Finally, with respect to Experiment 2, let us consider the top boundary. As
2The fact that symmetry in the solution is utilized here requires, by definition, that the
left-vertical boundary be set as no flow.
212
-------
introduced in the discussion on the horizontal discretization, in order to mimic the
NAPL source condition3, we need to apply a 0.5 cm NAPL head uniformly over
2.5 cm2. Figure F.6 provides an example of an appropriate boundary condition.
Note, that there are several nodes between the two different head conditions which
are set to no flow. This separation is employed because the NAPL- and gas-phase
head conditions are nonlinear4, and numerical experiments suggest that, for the
updating algorithm used in the simulator, putting a no flow buffer between them
is more computationally robust.
Given the setup for Experiment 2, the setup for Experiment 3 is straight
forward. As shown in Figure F.6, the NAPL head conditions are replaced by no
flow conditions.
It must be emphasized that the boundary condition setup used in Figure F.6
is only one of several admissable strategies.
Get to know the physics
Everything you need to know regarding the cause-and-effect relationship between
the model parameters and the model solution can be found in the simulator doc-
umentation. So by all means, read relevant sections.
F.6.3. Results
Please respond to the following:
1. Define an appropriate horizontal discretization and set up and run Exper-
iment 1 in Figure F.5. This result is your initial condition for the next
simulation. The files with the extension rs contain all the data necessary
for the code to restart the simulation (see see file sni.in, input line 70).
Since these files are over-written at each print interval, copy them into an-
other directory so that, if you have to, you can go back to this simulation
time (e. g., 5000s).
2. What role do the interfacial tension parameters play in modeling the three-
phase flow problem? How can you use them as calibration parameters?
3 The experimental source area is 5 cm by 2 cm. The model is assuming symmetry, and the
2-D domain has unit length in the omitted dimension. Therefore, the model source area is 2.5
cm by 1 cm.
4Actually the code converts the conditions to Pw, L e., Pw = PCNW(SW)+PCGN(STW)+P°
for the gas condition, and Fw = PCNw(Sw) + PN f°r the NAPL condition.
213
-------
3. Assuming that your discretization is denned such that numerical errors are
not significantly affecting the solution, focus your attention on the following
'short list' of calibration parameters:
• permeability, k
• residual saturations:Swr,
-------
F.7. TASK 4
F.7.1. Purpose
The purpose of this Task is to simulate the two-dimensional DNAPL spill exper-
iment using the numerical model, and after comparing results, to consider alter-
native data sets which may provide a qualitatively better match to the available
experimental data.
F.7.2. Procedure
From Task 3, you should be set up to run the DNAPL spill experiment. Choose
an appropriate parameter set for the three-phase flow problem. Run the model
as a series of sub-models, as described in Task 3 and as illustrated in Figure F.5,
where you use the simulator's restart facility to set the initial conditions as the
solution of the previous problem. It is suggested that after each successful model
experiment you save a copy of the restart files in another directory so that you
can restart from that point if necessary.. For example, much of your calibration
efforts will be focused on Experiments 2 and 3 in Figure F.5, therefore, being able
to use the restart data from Experiment 1 several times is appropriate.
Given an appropriate set of boundary conditions, consider the following pro-
cedural summary for calibrating the model:
1. Choose an appropriate set of values for the following 'short-list' parameters:
• permeability, k
• residual saturations :SVi/r, Sfjnr, SjVmr? $Gr
• S — P model fitting parameters: a^, a-i-, r\
• k — S model fitting parameters:
-------
compare the computed solution with the experimental result (i. e., Figure
F.7). Compare the volumetric rate of PCE infiltration predicted by the
model with the experimental results (i. e., Table F.2).
4. If the match is not qualitatively good, especially with respect to the infil-
trated PCE volume, then adjust one or more of the 'short-list' parameters
and re-run step 3 until you think the solution is acceptable. Note, depending
on which parameters you change you may have to re-run Experiment 1 to
generate appropriate initial conditions.
5. After capturing the PCE infiltration phase of the experiment (Experiment
2), run the PCE re-distribution phase of the experiment (Experiment 3)
using as initial conditions the solution to Experiment 2. Visually compare
the computed solution with the experimental results (i. e., Figures F.8, F.9
and F.10).
6. If the match is not qualitatively good, especially with respect to the vertical
distribution of PCE, then adjust one or more of the 'short-list' parameters
and re-run step 5 until you think the solution is acceptable. Note, depending
on which parameters you change you may have to re-run Experiments 1 and
2 to generate appropriate initial conditions for Experiment 3.
F.7.3. Results
Please respond to the following:
1. What are the values of the 'short-list' parameters which gave the best re-
sults?
2. Plot the model solutions for the following times (where T = 0 is assumed to
represent the initial conditions for Experiment 2: 143, 283, 683, 1195, 1795
and 3595 seconds.
3. For which 'short-list' parameters is the model most sensitive with regard to
the cumulative volume of PCE infiltrated? Briefly discuss why.
4. For which 'short-list' parameters is the model most sensitive with regard to
the areal distribution of the PCE over time? Briefly discuss why.
216
-------
5. After the DNAPL penetrates the capillary fringe, the rate of vertical migra-
tion of the model solution should equal that of the experimental results. For
example, even if the model DNAPL front is ahead of the experimental front
in the liquid saturated zone, the difference should not increase remarkably
as time progresses. What parameters define the DNAPL front speed in the
liquid saturated zone?
6. When trying to calibrate the model to Experiment 2, which ^ short-list^ pa-
rameters can you change without necessarily having to re-run Experiment
1.
7. Given your understanding of the mathematical model (i. e., the physical
problem is translated into a series of interrelated mass balance equations
and constitutive relationships), and the results from the model calibration
exercise, does it seem to you that the model is ''missing'1 some of the physics?
For example, no matter how you adjust it, the model can not capture the
shape of the PCE front as it moves through the capillary fringe, therefore,
something must be missing. Recall that, when building the mathematical
model, a series of major simplifying assumptions were incorporated either to
make the physical problem tractable from a computational point of view, or
to effectively 'fill in the blanks' with respect to data which can not currently
be measured (e. g., three-phase relative permeability). Briefly discuss which
major simplifying assumptions may contribute to the to the fact that the
model is physically flawed?
217
-------
Time = 143 seconds
67 cm
35.5 cm
5.5cm
Figure F.7: Video image at time = 143 seconds just after the DNAPL source
was removed (i. e., at the end of Experiment 2), showing dimensions, where the
superimposed grid is for reference purposes (elements are 5.5cm by 5.5 cm)
218
-------
Time = 283 seconds
Time = 683 seconds
Figure F.8: Video images at time = 283 and 683 seconds after the DNAPL spill
began (i. e., take Time = 0 as the initial condition for Experiment 2).
219
-------
Time = 1195 seconds
Time = 1795 seconds
Figure F.9: Video images at Time = 1195 and 1795 seconds after the DNAPL
spill began.
220
-------
Time = 3595 seconds
Figure F.10: Video image at Time = 3595 seconds after the DNAPL spill began.
221
------- |