Citation

## Material Information

Title:
Adding microscale effects to a macroscale SDE model of NAPL flow in heterogeneous porous media
Creator:
Barnhart, Kevin S
Place of Publication:
Denver, Colo.
Publisher:
Publication Date:
Language:
English
Physical Description:
xii, 97 leaves : ; 28 cm

## Thesis/Dissertation Information

Degree:
Master's ( Master of Science)
Degree Grantor:
Degree Divisions:
Department of Mathematical and Statistical Sciences, CU Denver
Degree Disciplines:
Applied Mathematics
Committee Chair:
Lodwick, Weldon
Committee Members:
Bennethum, Lynn
Dean, Dave

## Subjects

Subjects / Keywords:
Nonaqueous phase liquids ( lcsh )
Porous materials -- Fluid dynamics ( lcsh )
Differential equations ( lcsh )
Granular materials ( lcsh )
Differential equations ( fast )
Granular materials ( fast )
Nonaqueous phase liquids ( fast )
Porous materials -- Fluid dynamics ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

## Notes

Bibliography:
Includes bibliographical references (leaves 93-97).
General Note:
Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Kevin S. Barnhart.

## Record Information

Source Institution:
Holding Location:
Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
66462116 ( OCLC )
ocm66462116
Classification:
LD1193.L622 2005m B37 ( lcc )

Full Text
Wi
ADDING MICROSCALE EFFECTS TO A MACROSCALE SDE MODEL
OF NAPL FLOW IN HETEROGENEOUS POROUS MEDIA
by
Kevin S. Barnhart
Bachelor of Science, Montana Tech, 2002
A thesis submitted to the
University of Colorado at Denver and Health Sciences Center
in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Mathematics
2005

This thesis for the Master of Science
degree by
Kevin S. Barnhart
has been approved
by
/
Dr. Weldon Lodwick
fj/)V. 2\$jMOS
Date"7

Barnhart, Kevin S. (M.S., Applied Mathematics)
Adding Microscale Effects to a Macroscale SDE Model of NAPL Flow in Het-
erogeneous Porous Media
Thesis directed by Professor Dr. Weldon Lodwick
ABSTRACT
A three-dimensional stochastic differential equation to model the flow of a
non-aqueous phase liquid in a saturated heterogeneous porous media developed
by D. Dean and T. Russell1 is explained. This model aims to closely describe the
physics of the fluid flow at the macroscale while taking into account phenomena
at the meso and micro scales. An additional microscale addition is presented
which is intended to simulate irregularities in the pore structure. From empirical
data, lognormal distributions of grain size are used to generate a probabilistic
percolation at the macroscale. This methodology has a high degree of agree-
ment with new experimental water retention curves suggesting a powerful link
between the distribution of grain sizes and how the two-phases move through
a porous media. A brief discussion is given on how the computational model
is parallelizeable along with numerical results showing the efficacy of efforts in
implementing the parallelized framework.
1D. Dean is a senior researcher at the UCDHSC Center for Computation Mathematics
where T. Russell served as a tenured faculty (he is currently a director at the NSF). The
model is still being refined by D. Dean in collaboration with T. Illangasekare, M. Mathews,
T. Sakaki and others at the Colorado School of Mines Environmental Science and Engineering
department.
ill

This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
Dr. Weldon Lodwick
IV

DEDICATION
To my parents. This thesis is a production of their values and support.

ACKNOWLEDGMENT
I cannot give enough thanks to (Dr.) Dave Dean for his support. He has
unselfishly given his time and expertise to help me to learn the construction
of his model, overcome my mathematical deficiencies, guide my research, and
explain many of the subtleties of doing applied mathematics. Working with him
has been the most edifying experience that I have had as a graduate student at
UCDHSC.
I would like to thank Dr. Lodwick for his support and more than anything
his extreme patience. It has been nice to work with someone who wishes me
success in however I may define it-even if it doesnt coincide with his own desires.
Many thanks to Dr. Bennethum whose enthusiasm and understanding of the
I would like to acknowledge Dr. Cozzens who is at least somewhat respon-
sible for my involvement in the Applied Mathematics graduate program.
Dr. Illangasekare, Dr. Sakaki, and Dr. Mathews of the Colorado School of
Mines have been helpful by providing feedback and supplying me with valuable
data.
Finally, my family and friends have made earning this degree a much easier
(or even possible) process. I cannot express enough gratitude to them.

CONTENTS
Figures .............................................................. x
Tables............................................................ xii
Chapter
1. Introduction...................................................... 1
1.1 Overview ......................................................... 1
1.2 Summary of Results ............................................... 3
1.3 Nomenclature...................................................... 3
2. Preliminaries .................................................... 6
2.1 Fundamentals of Groundwater Fluid Mechanics.................... 6
2.1.1 Darcys Law.................................................... 6
2.1.2 Representative Elementary Volumes for Sands and Fluids .... 8
2.1.3 Saturation and Volumetric Fraction................ 11
2.1.4 Continuity Equation........................................... 13
2.1.5 Young-Laplace Equation........................................ 14
2.1.6 Water Retention Curves........................................ 17
2.1.6.1 Brooks-Corey Theory ....................................... 17
2.1.6.2 Other Water Retention Curve Theories........................ 18
2.2 Measurability................................................ 19
2.2.1 Measurable Sets .............................................. 19
2.2.2 Measurable Functions ......................................... 22

2.3 Probability Theory.................................................. 23
2.3.1 Random Variables and Distributions.................................. 27
2.3.2 Expectation and Conditional Expectation............................. 30
2.4 Stochastic Processes................................................ 33
2.4.1 Definitions ........................................................ 33
2.4.2 Martingales......................................................... 35
2.4.3 The Weiner Process ................................................. 36
2.5 Stochastic Differential Equations .................................. 40
2.5.1 The Ito Integral.................................................... 40
2.5.2 Ito Stochastic Differential Equations (SDEs)........................ 45
2.5.2.1 Itos Lemma....................................................... 46
2.5.3 Fokker-Planck Equation.............................................. 48
2.5.4 Fokker-Planck PDE to Ito SDE........................................ 50
2.6 Heaviside Step Functions and the Delta Function.................. 51
3. The SDE Model........................................................... 55
3.1 Motivation......................................................... 55
3.1.1 The Saturation-Probability Relationship............................. 58
3.2 Fokker-Planck Equation for Fluid Flow............................... 59
3.3 PDE to SDE ......................................................... 61
3.3.1 Specifics........................................................... 61
3.3.1.1 qt Velocity Field................................................. 63
3.3.1.2 Interface Control with the Jump Term......................... 64
3.3.2 Computational Model................................................ 69
4. Macro Percolation....................................................... 71
viii

4.1 Particle Size Distributions...................................... 71
4.2 Pore Size Distributions.......................................... 75
4.3 Lognormal Water Retention Curves................................. 76
4.4 Water-NAPL Water Retention Curves .............................. 78
4.4.1 Macro Percolation in the SDE Model............................. 79
4.5 Beowulf.......................................................... 85
4.6 Future Work ..................................................... 90
5. Conclusion......................................................... 92
References............................................................ 93

FIGURES
Figure
1.1 An Application: Pollutant Spill................................. 2
2.1 Fluid Volume and Density Characteristics....................... 10
2.2 Contact Angle of Wetting and Non-Wetting Fluids................ 15
2.3 Young-Laplace Equation Example Setup........................... 16
2.4 Water Retention Curve Experimental Data and Brooks-Corey Model 18
2.5 Example Path of a Stochastic Process .......................... 33
3.1 Experimental Tank Setup........................................ 55
3.2 LNAPL Plume in Experimental Tank............................... 56
3.3 Premature Breakthrough of Plume from Numerical Diffusion .... 57
3.4 Microscale Effects that are Stochastic at the Macroscale....... 58
3.5 Dc Df versus Non-Wetting Phase Saturation ......... 66
3.6 Entry Pressure and Corresponding Saturation.................... 67
3.7 Dc Df with Brooks-Corey Theory .............................. 68
3.8 Example Run of Computational SDE Model......................... 70
4.1 Sands Used in Experiments...................................... 71
4.2 Lognormal CDF Curves Fit to Grain Size Data.................... 74
4.3 Lognormal PDF Curves Fit to Grain Size Data.................... 75
4.4 Lognormal WRC Compared with Experimental and Brooks-Corey
WRCs........................................................... 79

4.5 Example Run of SDE Model with a Well-Defined Interfaces .... 80
4.6 Interfacial Layer Between Fine and Coarse Sands............... 81
4.7 Variation in Grain Size Causes Variation in Pore Size......... 81
4.8 Brooks-Corey and Lognormal Curves for #16 and #30 Sands ... 83
4.9 Dc Df with Brooks-Corey and Macro Percolation............... 84
4.10 Before and After Macro Percolation is Applied to the SDE Model . 85
4.11 UCDHSC CCM Beowulf Cluster ....................................... 86
4.12 Computation Time versus Number of Particles for Movement and
Pressure Calculations........................................ 87
4.13 Grid for Finite Element Method............................... 88
4.14 Decrease in Computation Time from Parallelizing Pressure Field . 89
4.15 Fingering in LNAPL Plume..................................... 91
xi

TABLES
Table
4.1 Grain Size Data from Sand Suppliers................................ 72
4.2 Lognormal Distribution Parameters, SSE Fit, and Distribution Statis-
tics ................................................................... 73

1. Introduction
1.1 Overview
High standards of living in industrialized nations have a major negative
impact on the environment. Consumerism is destroying habitat through mining,
deforestation, and pollution at an alarming rate. The demand is unarguably
unsustainable and as time goes on more resources will need to be spent on
mitigating the most immediate of these issues.
Beginning in the mid-to-late eighties there has been increased research
in dense groundwater pollutants. This classification includes polychlorinated
biphenyl oils, creosotes, and chlorinated solvents. These fluids have a wide va-
riety of uses and are fairly inexpensive to produce, making them common in
consumer products and in the production of consumer products.
Because these fluids are so hazardous to biological life, expensive precautions
are required for even controlled laboratory experiments. Field experiments are
even more difficult to perform because of environmental concerns. Unfortunately
such care is not always taken in the industry which sometimes leads to the
contamination of groundwater systems (see Figure 1.1).
1

Figure 1.1: An Application: Pollutant Spill
For the above reasons, there is an increasing need for hydrological researchers
to produce models (both mathematical and computational) that describe how
these contaminants will flow in (water) saturated porous soils. Despite twenty
or more years of research on this particular problem and the similarity of this
problem to the extraction of petroleum from depressurized oil reserves, there are
very few sureties. Uncertainties in the geological structures from the size and
shape of sand grains to the rock types in a large aquifer are certainly stumbling
blocks, but, also, the very nature of the flow of two fluids is being re-thought.
This may be somewhat surprising since hydrologists have been using some of
these laws of fluid flow for over a century.
2

This thesis describes a three-dimensional probabilistic model by Dean and
Russell [12] of two-phase flow in porous medium. While other models may be
computationally faster, this one attempts to capture the physics of the situation
as accurately as possible; resulting in a complex model because the scenario at
hand is also quite complex. To this end, the author has made a contribution
to the model motivated from microscale (i.e. pore scale) phenomena in order to
better the multi-scale representation.
1.2 Summary of Results
The proofs of Theorem 2.13, Lemmas 2.32 and 2.33, and Proposition 2.35
are created by the author, but are not novel constructions.
Improvements to the model are given at the end of Section 3.3.1.2 and
throughout Chapter 4.
1.3 Nomenclature
The following is a list of commonly used symbols and what they represent.
a.s. almost surely
Bd the set of the Borel sets on Rd
B(I) the Borel sets generated by the interval /
Ac the complement of the set A
6 delta function
0 the empty set
T event space
Tt a filtration
T* the natural filtration generated by the stochastic process X
3

9 gravitational constant [Length/Time2]
1 specific weight [Force/Length3]
H Heaviside step function
i
k (intrinsic) permeability [Length2]
K hydraulic conductivity [Length/Time]
A fitting parameter for Brooks-Corey model
9 viscosity [Force x Time/Length2]
n porosity
Q sample space
IX 3.14159...
Pd entry pressure [Force/Length2]
p pressure [Force/Length2]
Pc capillary pressure [Force/Lengtli2]
P probability measure/distribution function
V(S) the power set (set of all subsets of S)
Q specific discharge [Length3/Time]
Q Darcy flux [Length/Time]
Rd the d-dimensional space of real numbers
P density [Mass/Volume]
4

r.v. random variable
S saturation
Se effective saturation
Srw residual saturation
Ac B A is a strict subset of B
ACB A C B or A = B
t time
9 volumetric fraction
0 contact angle
* V velocity [Length/Time]
V volume [Length3]
w.r.t with respect to
wt Weiner process at time t
i
5

2. Preliminaries
2.1 Fundamentals of Groundwater Fluid Mechanics
Within this section, many of the basics of fluid mechanics in porous media
are reviewed. For a more comprehensive explanation, the reader is referred to
the classic Dynamics of Fluid in Porous Media by Jacob Bear [3].
This study is only concerned with the flow of immiscible fluids or phases.
These fluids are not soluble in each other and, consequently, a tension is created
at their interface. This tension creates a pressure difference known as the cap-
illary pressure between the two phases. Fundamental laws that govern the flow
of these fluids in porous media are summarized below.
2.1.1 Darcys Law
In 1855 and 1856, Henry Darcy formulated what is now called Darcys Law
by experimenting with water flow in sand. The law is a constitutive equation.1
which describes the flow of fluid through a porous medium, and is still generally
accepted when used to describe the viscous (slow) flow of fluids.
Darcys Law is an expression of the conservation of momentum which makes
it analogous to Ohms law in the theory of electrical networks. It is given by
the equation:
Q = ~KAha~kb (2.1)
Â¥
where Q is known as the specific discharge (the volume of fluid flowing per unit
time), K is the hydraulic conductivity (a value which describes the ease with
1A constitutive equation is specific to a material or substance.
6

which the fluid can move through the porous media), A is the cross-sectional
area to the flow, and Vh = (ha hb)jL is the hydraulic head drop or hydraulic
gradient where ha and h,b are the heights at points a and b, respectively, and L is
the distance between the two points. See [3] for a derivation. Define h = ha hb
as the hydraulic head which has the corresponding pressure
P = hpg, (2.2)
where p is the density of the fluid and g is the acceleration due to gravity.
If (2.1) is divided by A then an equation for the Darcy flux is obtained:
q = -KVh (2.3)
which describes the flow per unit area.
Since the fluid is only moving through the pores in the medium, the Darcy
flux should be divided by the porosity of the medium in order to find a rule for
the velocity that a tracer would actually experience. This relationship is known
as the Dupuit-Forchheimer Equation, where n is the porosity:
v = -. (2.4)
n
Another relationship is that between the hydraulic gradient and the piezometric
S7h = -V0*.
The piezometric head is defined as the elevation above some datum plus the
pressure head, so it is expressed as:
(j)* z-\-h = z + (2.5)
99
7

Vfi* = (Vp pgz) where z =
P9
v-v
As mentioned above, the hydraulic conductivity measures the ease with
which a fluid can move through a porous medium. To keep with the electrical
circuit analogy, hydraulic conductivity is akin to 1/resistivity. It is related to
the intrinsic permeability or just permeability, k, through the equation:
K = (2.6)
P
where
7 = P9
(2.7)
is the specific weight of the fluid and p is the viscosity. Thus, hydraulic con-
ductivity depends on the properties of both the fluid and the porous medium,
whereas permeability is a property of the porous media alone. Put simply, per-
meability is just a measure of the ability of a material to transmit fluids through
it.
From the equivalences given above Darcys Law is restated as
Q = (Vp 72). (2.8)
ft
2.1.2 Representative Elementary Volumes for Sands and Fluids
This section describes the concept of a Representative Elementary Volume
(REV) used in Continuum Mechanics. The explanation follows that given by
Bear in [3].
8

t
Fluids are composed of molecules which bounce off each other and the solids
around them. While one may be able to deterministically solve the direction
and momentum of a few molecules, we still lack the computing power to reliably
predict the path of even a few hundred molecules. So, instead of trying to
predict fluid flow using the properties at a molecular level, a statistical approach
Bear describes a fluid particle as
...an ensemble of many molecules contained in a small volume.
Its size is much larger than the mean free path of a single molecule. It
should, however, be sufficiently small as compared to the considered
fluid domain that by averaging fluid and flow properties over the
molecules included in it, meaningful values, i.e., values relevant to
the description of bulk fluid properties, will be obtained.
The fluid property most often used to determine the size of the fluid particle
is density the ratio of the mass of fluid, m, and the volume, V, in which it is
contained. Let x be a point in the fluid and consider a volume (e.g. a sphere),
V, for which x is the centroid. The relationship between the size of the volume
and the density can empirically be shown to be similar to Figure 2.1
9

Figure 2.1: Fluid Volume and Density Characteristics
The point Vo in Figure 2.1 indicates a volume where the density is fairly
stable. Let Vo to be the volume of the fluid particle. All points in the fluid
are then associated with a fluid particle with this density, and, in this way, the
fluid particle will have more definite dynamic and kinemetic properties which
can greatly simplify fluid flow models.
The same concept of the REV of a fluid (the fluid particle) can be applied
to porous media. The analogous property to density in porous materials is the
porosity the mass of solid per unit volume. The relationship between volume
and porosity is quite similar to that of volume and density and so a REV is
chosen in the same manner.2
2The notion of a REV is somewhat limited in application by the medium. For some
materials it is very difficult to obtain a reliable REV for which the desired property is stable
due to unreliable factors such as the packing and shape of the sand grains.
10

2.1.3 Saturation and Volumetric Fraction
Consider a REV of a homogeneous porous material with total volume V.
Let Vv represent the volume of the void in V, and Va represent the volume of
the a-phase3 in V.
Saturation is the percent of void space that is filled by the fluid. Accordingly,
is the saturation of the ct-phase in V.
When a 11011-wetting fluid displaces a wetting fluid (such as water) there will
always exist some residual water content (or wetting fluid) in the soil because
water becomes trapped in the crevices of the porous medium. The residual
saturation, Srw, is related to the effective saturation, Se:

S Sj'
1 Srw
(2.9)
where S is the measured saturation. The effective saturation is scaled to range
from zero to one by treating the volume of the residual water as a in-displaceable
substance (such as the porous media). The following equality is usually assumed
to hold in a two-phase system:
Smv + Sw 1-
(2.10)
The volumetric fraction is the fraction of total volume that is filled by a
fluid:
Va
V '
3The term a-phase is used in multi-phase discussions to indicate any one of the phases.
11

It is possible for the volume of the liquid to fill all the void space, so Sa is a
value between 0 and 1. However, the fluid can never fill the entire volume since
part of it is filled by the porous medium, and so 9a must vary between 0 and n
(the porosity)
v; = nV.
Using this and the definitions given above, the saturation and volumetric fraction
are related by
Particular to this study, is how the concept of fluid particles move in porous
media, so saturation and volumetric fraction is described with these in mind.
Each REV of the porous media has a finite volume and can therefore only contain
a finite number of fluid particles.
Let VfP be the volume of a single fluid particle of the non-wetting phase.
Given that p is constant, the number, Np, of a-phase particles needed to com-
pletely saturate an REV is
N
P V
So, if the REV contains x o-phase fluid particles its saturation is
9a = nSa.
(2.11)
and the corresponding volumetric fraction must be
(2.12)
12

2.1.4 Continuity Equation
The theory of continuum mechanics postulates the global balance law given
by
4 f at Jv(t) JS(t)
which states that the time rate of change of the property

of the fluxes, r, across the surface S(t) plus the sum of the forces, /, acting
on the body of interest [4], This law is said to be localized when the derivative
is moved inside the integral a multi-dimensional form of Leibniz Rule, called
the Reynolds Transport Theorem.4 Further, the divergence theorem allows the
second integral to be taken over the boundary of V.5 Assumptions must be
made to make this localization and about the smoothness of the integrand (see
[4]) to provide the local balance law as a differential equation,
^ + V ( As stated above, Darcys Law is specific form of conservation of momentum.
The other governing law is a statement of conservation of mass known as the
continuity equation:
^ + V (pv) 0
4Leibniz rule shows how to take the derivative of a definite integral. It is stated as
~Â§t Sa /Db= f Â§t^x + f(b> Â£)ft f(a, Off where a b are functions of t.
5The divergence theorem maintains that a property within a volume, V, must remain
constant unless there is flow across the boundary, dV, of that volume, f (V-tp)dV = fdv ip-da.
13

where p is again the density of the fluid, t is time, and v is the velocity of the
fluid.6 This equation conies directly from the local balance law, equation (2.13),
if the property to be conserved is density (ip = p, t = 0, and / = 0).
Assuming that the solid is incompressible, the equivalent form of the con-
tinuity equation inside a porous medium uses volumetric fraction and Darcy
flux:
DO
Â¥ + V.? = 0. (2.14)
This form of the continuity equation is prevalent in two-phase flow applications.
2.1.5 Young-Lap lace Equation
As described in detail by Bear in [3], the interfacial tension is the amount of
work required to separate a unit area of substance a from substance (3, which we
denote as (f>ap. Within a particular fluid, molecules have an (inward) attraction
to each other. This creates tension within the fluid. A similar tension exists at
the interface between two substances thus work would be required to separate
them.
If two phases (say a liquid and a gas) come in contact with a solid, they
will form a contact angle, 0, between the solid and the interface of the two
phases. Conventionally, the angle is measured through the denser fluid. A
wetting fluid will create a smaller contact angle than a non-wetting fluid because
of its tendency to wet the surface (see Figure 2.2). The wetting fluid will be
denoted by subscript W and the non-wetting by subscript ATT.
6The notation V x where x is some vector is the divergence of x and is shorthand notation
c)5f 3x dx dx
for ----h -----h . In the literature pertaining to this project this is also written as
OX\ OX 2 OX 3 OXi
where x is understood to be a vector.
14

Wetting Fluid
Non-Wetting Fluid
Figure 2.2: Contact Angle of Wetting and Non-Wetting Fluids
Consider an experimental setup of a capillary tube placed in a bucket of
water (see Figure 2.3). Water is a wetting liquid; consequently, the water level
will rise inside the tube until an equilibrium is attained. Suppose the tube
has a radius, r, so that its circumference is 2irr. In accord with the discussion
above, the water will form a contact angle with the tube which is dependent on
the interfacial tension between the water and air, 7 The weight of the water
will create a downward force in the tube (dependent on the height of the water
column, h) given by
F -jvr2hpwg-
The surface tension at the top of the water column has an upward component
of (f> cos 0 which creates the following force:
Fs = 27rr0cos 0. 7
7Actually, three interfacial tensions are at work to form the equilibrium equation cos 0 =
(sg 4>si)/ 15

atmospheric
pressure
<|> Cos
atmospheric
pressure
0
e
I
Pgh
,4>
2r>-
Figure 2.3: Young-Laplace Equation Example Setup
The atmospheric pressure between the water level in the bucket and the top of
the water column is negligible. So, in a static state F ~ Fs in order to maintain
equilibrium. Setting these equal, the capillary pressure head is given by the
Young-Laplace Equation:
, 2(f> cos 0
hc =-----------,
pwgr
(2.15)
or by (2.2),
Pc =
2(f> cos 0
r
(2.16)
16

In 1966, Brutsaert [9] found that for an air-water-soil system the capillary pres-
_ A 0.149 cm2
c r r
2.1.6 Water Retention Curves
An important characterization of soil moisture content is the water reten-
tion curve (WRC). The WRC describes the relationship between the amount of
water in the soil (either in terms of saturation or volumetric fraction) to the soil
capillary pressure. Many models have been proposed to describe this relation-
ship since it is often both expensive and time consuming to generate a curve
experimentally. Perhaps the most famous have been developed by Brooks and
Corey [8], van Genuchten [51], and Russo [43].
Here, drainage curves are described; these are water retention curves that
show how the non-wetting phase displaces the wetting phase. Imbibition curves
are also studied, which describe the reverse process.
2.1.6.1 Brooks-Corey Theory
In 1964, Brooks and Corey [8] proposed the following relationship between
the effective water saturation, Se, and the capillary pressure:
Pc = pdS7lx. (2.17)
The parameter pd is called the entry pressure. Along with A, the two parameters
are usually obtained by fitting (2.17) to experimental data.
The Brooks-Corey (BC) model is intended to describe phenomena at the
macroscale. Note that when Se = 1, pc = p,i- However, experimental data
indicates that this in not the case. (Refer to Figure 2.4.)
17

Figure 2.4: Water Retention Curve Experimental Data and Brooks-Corey
Model
The primary assumption of the BC model is that experimental data near
full saturation is an experimental artifact and should be ignored. Thus, at the
macroscale the fluid must be at or above the entry, 'pd, before it may enter into
a fully saturated porous medium.
2.1.6.2 Other Water Retention Curve Theories
There are many other WRC models. Many disagree with Brooks and Coreys
assumption at high saturations. For instance, the following S-shaped model was
proposed by van Genuchten in [51]:
se = i/{i + (obcirr
18

where n and a are the empirical parameter and m = 1 1 /n.
Also commonly employed is a two-parameter model by Russo [43]:
Se(Pc) = {(1 + P\Pc\)exp(-/3\pc\)}2/{2+k\
Fractal models have been extensively explored to model the WRC (e.g. [50]
and [11]) because at the microscale the fluid behaves like a fractal as it moves
through the network of pores (cite [22] and [23]).
More recently (e.g. [30]) these and other WRCs have been criticized because
their parameters must be re-fit for any simple change in the setup. Also, the
parameters do not appear to have real physical significance. Ideally the WRC
parameters should be determined using accessible physical characteristics of the
porous medium and the phases. This thesis will emphasize this approach.
2.2 Measurability
2.2.1 Measurable Sets
Modern analysis is built on concepts from set theory. On the real line, one
may desire to work with intervals (closed, open, and half-open), canonical sets
(e.g. integers, rationals) or some other set. An important construction is an
algebra of sets, A, which is a collection of sets where: if A is in A then Ac is
also in the collection, and if A and B are in A then their union (and therefore
the intersection by De Morgans laws) is in A. as well. This may be extended to
an infinite, but countable, number of intersections and unions in which case
the collection of sets is a a-algebra. Formally,
19

Definition 2.1 A a-algebra, A, is a collection of sets such that
1. if A Â£ A then Ac Â£ A; and
OC
2. if (Ai) is a countable collection of sets then [jAi is also in A
i= 1
What sort of sets are produced under countable intersections and unions? It
can be shown that the union of any collection (countable or otherwise) of open
sets is open, and that the intersection of a collection of closed sets is again a
closed set[42]. However, the opposite is not necessarily true. Thus, a er-algebra
of open sets may contain sets that are not open. This necessitates the definition:
Definition 2.2 The collection of Borel sets, B, is the smallest (i.e. it contains
the fewest sets) a-algebra that contains all of the open sets.
It turns out that such a collection exists over R and that it must also contain
all closed sets, open intervals, and half-open intervals[42].8
It is natural to seek the size (or measure) of these sets. In K, the size of an
interval is intuitively its length (the difference of its end points). Based on this,
the outer measure, m*, of a set is defined as
where (In) is a collection of intervals whose union will cover (contain) A and /(/)
is the length of the interval In. Such a covering exists for every set (trivially,
(oo, oo) covers every set on M). Unfortunately, the outer measure turns out to
be only sub-additive on any countable collection of sets of real numbers. This
8In fact, the collection of Borel sets can be constructed from any of these types of sets.
20

restriction is removed and additivity will be attained if the collection of sets is
measurable:
Definition 2.3 The collection of measurable sets, A4, contains all sets, E,
where
m*(A) = m*(A n E) + m*(A n Ec)
for any conceivable set A. If E is such a set then it is called measurable,9
In words, the definition states that a set is measurable if it breaks up nicely.
This is not a surprising construction given the foreknowledge that integration
and differentiation is all about breaking sets up into infinitesimal pieces. The
next two theorems give more examples of the advantages of measurability:
Theorem 2.4 The collection of measurable sets, M, is a a-algebra.
Proof: See an introductory book on measure theory, e.g. [42],
Theorem 2.5 Every Borel set is measurable.
Proof: Again, see [42],
When it is understood that only measurable sets are being used, then the
outer measure is referred to as the Lebesgue measure, m.
9This is known as Caratheodorys Condition.
21

2.2.2 Measurable Functions
Functions that are defined over measurable sets also have nice properties.
Definition 2.6 An extended real-valued function, f, is said to be Lebesgue
measurable or just measurable if its domain is measurable and one of the
following is true:
1. For every real number a the set {x : f(x) < a} is a measurable set.
2. For every real number a the set {x : f(x) < a} is a measurable set.
3. For every real number a the set {x : f(x) > a} is a measurable set.
4. For every real number a the set {x : f(x) > a} is a measurable set.
Of course, the concept of measurability can be generalized to almost any measure
(such as the probability measure given below). Because the set of measurable
sets is a cr-algebra, the above definition is given in general as
Definition 2.7 A function, f, is said to be measurable from (Q, T) to (ib, Q)
(with respect to some measure) if
f-\B) = {u : f{u>) Â£ B} T
for each B Q. FI ere, T is a measurable a-algebra over il and Q is a measurable
cr-algebra over 'h.
22

Example 1. The function / : R > R is Borel measurable if = {x :
f{x) eB)eB for every B Â£ where B and C are Borel cr-algebras in the range
and domain, respectively.
The function of interest may be defined on a product space or product
topology in which case the function is measurable over the a-algebra of the
product space and with respect to a product measure (see, e.g., [42]).
2.3 Probability Theory
This introduction to the theory of probability will follow [27] and [28].
Probability theory is based on the notion that knowing how frequently some-
thing has occurred in the past will give some indication on how frequently it
will occur in the future. For instance, consider a six-sided die that is not evenly
weighted. In order to determine which sides of the die are weighted more heavily
than others, one would roll the die a sufficient number of times and monitor
the relative frequency of each side (e.g. perhaps the sixth side came up 202
times in 1000 rolls or about 20% of the time). Each roll is an experiment that
produces an outcome. Probability theory states that if we roll the die an infinite
number of times then the relative frequency of each side will approach some
asymptotic value: the probability. Even though it is impossible to know in
outcome based on past experience.
There have been many different constructions of Probability Theory. A.N.
Kolmogorovs approach, given in 1933, is still the most common, and will there-
fore be presented here. It is mainly just an extension of Measure Theory.
A random experiment is simply any action which has uncertain outcomes.
23

Definition 2.8 The set of all possible outcomes of a random experiment is
known as the sample space, Q.
Perhaps we desire to know if the larger numbers, {4, 5, 6}, on the die occur
more commonly then the lower numbers, {1, 2, 3}. In this case, we are interested
in the occurrence of subsets of Q (equivalently, elements of V(Q)).
Definition 2.9 The event space, T, is a set of subsets of LI.
An element of E T is called an event. If Q T (as is usually required) it
is called the certain event. The empty set, 0, is known as the impossible event.
If E G Q as well, then it is an elementary event. In general we are interested
in event spaces that are a er-algebras because we wish to define measurable
functions on them.
Consider the relative frequency, fn(E), which returns the number of times
event E occurred in n experiments divided by n. The probability of E can be
roughly defined as lim fn(E) ~ P(E).10 The following definition attempts to
noc
capture this intuition.
Definition 2.10 Let T be a a-algebm on fl. A probability measure is a
function P : T > [0,1] where:
10In some cases, the relative frequency function can be used in the definition of the proba-
bility measure. The limit must exist for each E in order to set lim /(Â£?) = P(E).
1. P(0) = 0 and P(Q) = 1.
2. Ei,E2, ... G T and Vi,j, Ei n Ej = 0 then P

71 1
OO
24

Definition 2.11 A probability space is the triple (0.,F,P).
The probability of rolling the die and none of the sides showing is zero,
while the probability that any of the sides does show is one. It also makes sense
that the probability of obtaining a 4, 5, or 6 is the sum of the probabilities of
obtaining a 4, 5, or 6 individually.
Definition 2.12 A probability space is complete if T contains all subsets of
il that have P-outer measure zero. That is, all sets A c Cl where P*(A) :=
inf{P(F);Fe F,Ac F} = 0.
We still desire the event space to be a cr-algebra, so in order to complete a
probability space, sets with a P-outer measure of larger than zero must also be
included. Actually, any incomplete probability space can easily be completed as
is shown in the next theorem.
Theorem 2.13 If (Q,F,P) is a probability space that is not complete, then we
can construct a complete probability space, (Q,F0,Po), such that:
1. Pc F0;
2. E T => P(E) Po(E); and
3. E e F & E = AuB where B e F, A C C, C e F, and P(C) = 0.
Proof: The last item in the theorem basically states how to construct Fq.
Let
F0 = {E : E = A U P, B G T, A C C, C E F, P(C) = 0}.
25

This construction makes sense from the definition if we recall that 0 Â£ T by
definition, and so A U 0 will produce completeness since A C C with P(C) = 0.
Where B ^ 0 we are filling out the rest of JF0 to make it a cr-algebra.
For the rest of the thesis take T to be a completed cr-algebra. From Defi-
nition 2.10, the probability measure is monotonic. That is, if A,B Â£ T where
AC B, then P(A) < P(B).
It is desirable to find the conditional probability, P(A\B), of event A knowing
that event B has already occurred. If knowing that B has occurred does not
have any affect on the probability of A occurring, then P(A ft B) = P(A)P(B),
and A and B are said to be independent. Formally,
Definition 2.14 A (possibly infinite) collection of events (Et),iej is an inde-
pendent collection if VJ C /, J < oo there is multiplicity under intersections:
P
This definition is different from pairwise independence where multiplicity
exists between pairs and not necessarily between entire subsets.
probability that A will occur, known as dependence. In this case, it makes no
sense to consider where A and Bc both occurred since we are certain B has hap-
pened and because both B and Bc cannot occur at the same time (P(BPi Bc) =
P(0)). This motivates the informal definition: P(A\B) ~ lim ^ . Now,
fn{B)
formally:
Definition 2.15 Let A and B be two events in T where P(B) / 0. The con-
P{AHB)
ditional probability of A given B is P(A\B) =
P(B)
26

2.3.1 Random Variables and Distributions
The previous sections have captured the idea of a random experiment by
specifying its possible outcomes (Cl), specifying subsets of outcomes (T), and
giving a measure of the chance of the outcome (P). Next, a numerical value
(other than the probability) is associated with random events. For instance, let
the sample space be made up of all possible pairs of latitude and longitude in
the continental U.S. If a location is chosen at random, the height above sea level
may be associated with that event. A map from the sample space to a set of
numerical values is needed: a random variable. Formally,
Definition 2.16 Given a probability space (Cl, T, P), a function X : Q is
a random variable (r.v.) if it is a measurable function from (Cl,J-) to (K, B):
Example 2. A simple example of a random variable is the indicator function,
A random variable can be described by its probability distribution which is
defined in the abstract sense:
(uefl: X(u) for all AeB.
27

Definition 2.17 Let (Ll,P, P) be a probability space, and let X be a r.v. on
that space, then the distribution or law of X is:
PX(B) = P(X-\B)) = {Po X~1)(B) = P(u Â£ Ll : X(u>) Â£ B) = P(X Â£ B)
for all B Â£ B.
Working with a set function may not be ideal, so a distribution function is
defined that works with points in M.d:
Definition 2.18 Let (Ll,P,P) be a probability space, and let X be a r.v. on
that space, then the (cumulative) distribution function of X is:
Fx(x) = Px{{-oo,x)) = P(ujÂ£tt: X(u) < x)
Example 3. The distribution function for the indicator function can be given
as:
0
Fia(x) =
< l-P(A)
1
x < 0
0 < x < 1
1 < x
A random variable is called a continuous random variable if its distribution
function is continuous. It necessarily follows that P({lo Â£ Q : X(u>) = x}) = 0
since P must be a countably additive measure and P{Ll) 1. Next, if the
distribution function is differentiable, then we can also use a density function to
describe a r.v.:
28

Definition 2.19 Let (Ll,X, P) be a probability space, and let X be a r.v. on
that space with (differentiable) distribution function Fx, then the probability
density function (pdf) or just density function of X is given by:
f(x) = Fx(x)-
Two of the most common continuous random variables are given below.
Example 4. The distribution function:
/
0 x < a
< xa a < x < b
ba
l b < x
and the associated density function:
Fx(x) = f(x) =
0 x [a, b]
x e K b1
describe the uniform random variable, X, that takes on values in the interval
[a, 6] with equal likelihood.
Example 5. Consider the density function:
/M =
y Zita
whose corresponding distribution function is differentiable everywhere (theoret-
ically), but must be evaluated using numerical techniques. A r.v. with such a
distribution is said to follow a normal or Gaussian distribution which is indi-
cated by the notation: X ~ N(p,,a). The parameter // is the mean or average
29

value of the distribution, and a is the spread of the distribution. The distri-
bution function will look bell-shaped for any choice of p and a. The normal
distribution is often used to model measurement error around the true value
(/a). The cumulative distribution function (right-tailed) is
2.3.2 Expectation and Conditional Expectation
The notion of expectation is most intuitive when considering the discrete
random, variable X (on (fi,.F, P)) that can take on values Xi,X2,x3,.... The
expectation of X is the weighted sum:
where At = {u> G ft : X(uj) = ay} are events. X is expected to take on the value
ay with probability P(Af), so by summing these values, the expected value for
X is obtained.
Expectation is defined for the continuous case as:
Definition 2.20 Let X be a continuous random variable on the probability space
P) that takes on values in R. The expectation or expected value of
X is defined to be:
OC
(2.18)
E(X)=^2xiP[Ai)
i> 1
n
R
30

The integrals in the definition are Lebesgue integrals (see [42], for example).
This definition is analogous to the discrete case given above if one thinks of the
density function Px{x) as estimating the probability P(X = x) even though
this probability is always zero for the continuous case.
Due to the properties of the Lebesgue integral (again, see [42]), expectation
is both linear and monotonic: for r.v.s X and Y and constants a and (5
It is possible to find the expected value of a function of X as well.
Theorem 2.21 Let X be a r.v. on (Q.,F,P) and let g : (Q, T) > (R, B) be a
measurable function then
Definition 2.20 provides the expected value over the range of all possible
values of the r.v.. If, however, it is known that the r.v. will only take on a set of
values that correspond to a sub-er-algebra of the event space, of a conditional expectation can be defined:
Definition 2.22 Let (Q, X, P) be a probability space and let X : Q Rd be a
random variable with finite absolute expectation (T?[|X|] < oo). Let S C T be a
a-algebra. The conditional expectation of X given S, E[X\S], is defined
as any r.v. (from Q to that satisfies:
E[aX + /3Y] = aE[X\ + pE[Y],
and when X < Y
E[X) < E[Y].
Proof: Theorem 9.5 of [27].
31

1. E'fXItS] is S-measurable; and
2. I E[X\S]dP = f XdP for all S eS.
s
s
The existence of Â£?[X|S] is guaranteed by the Radon-Nikodym theorem (see
[42] or [27]). Briefly, define a measure v as:
By the Radon-Nikodym theorem, there exists a unique function Epf | measurable w.r.t. Some of the properties of the conditional expectation are given in the next
theorem.
Theorem 2.23 Let X and Y be random variables on the probability space
(Ll,X, P), where both X and Y have finite absolute expectation. Let a,b 6l be
constants, then the following are true:
s
u(S) = / E[X\S]dP.
s
1. E[aX + bY\S] = aE[X\S} + bE[Y\S};
2. E[E[X\S}} = E[X];
3. Â£J[X|S] = X if X is S-measurable; and
4. 2?[X| Proof: See [27] or [35] for details.
32

2.4 Stochastic Processes
2.4.1 Definitions
This section was written by following the work given in [17], [35], [28], and
[38]-
A stochastic process can be thought of as a collection of random variables,
each of which apply to a position of a particle at some moment in time.
Definition 2.24 Let (fl,T,P) be a probability space. Define the function
X(t,u) : / x Q > Ed as a Ed-valued stochastic process if for each t G I,
X(t) G Md (where I is an interval in Ej.
A particular X(-, u>) is called a path or realization of X (see Figure 2.5). The
term path is particularly descriptive in this application since it might describe
the possible paths that a fluid particle (u) could travel through time (t).
Figure 2.5: Example Path of a Stochastic Process
33

Not surprisingly, stochastic processes that are measurable are convenient to
work with.
Definition 2.25 A stochastic process X is measurable or T-measurable
(from (/ x ft,Â£(/) x X) to (Rd,Bd)) if X~\B) = {(t,ur) : X(t,uj) B} e
B(I) x T for all B G Bd.
Continuity is another valuable property to have in a stochastic process:
Definition 2.26 Let X be a stochastic process. If for every w SI the func-
tion X(t, lo) is right/left continuous or continuous, then the process X is also
right/left continuous or continuous.
As in measure theory, continuity and other conditions are sometimes relaxed
so that the property holds except on a set of measure zero. If this is the case
the property is said to exist almost everywhere (a.e.), but in probability theory
it is more common to say that the property holds almost surely (a.s.). That
is, there may be a set (or sets) N C Q where P{N) = 0 where the property
does not hold (e.g. the process is not continuous on the set N but is continuous
everywhere else). Thus, the condition will hold almost surely since P(Q \ N) =
P(Q) P(N) = 1 0=1 which follows directly from the additivity axiom of
P.
Often, a process takes on events in a specific subset or collection of subsets
of the event space. It is useful to consider a subset that can occur up to a finite
time, t.
34

Definition 2.27 A filtration on (Â£2, T) is an increasing collection of a-
algebras, {Tt)t>o, where
Xs C Tt C T, 0 < s Given a stochastic process, Xt, the natural filtration, T* w.r.t. the process
is the a-algehra of
U <*(*.)
0 where a(Xs) is the a-algebra on the range of Xs.
We then consider stochastic processes that are measurable with respect to
a filtration:
Definition 2.28 A stochastic process X is said to be Tt-adapted if X(t) E Tt
for all t > 0U
2.4.2 Martingales
Consider a game where \$1 is wagered each round. Let the stochastic process
Xt represent the total winnings after t games. If the game is started with Xs
winnings there will be Xt Xs net winnings at some time t > s. If all possible
events that could occur (a filtration) up to time s are known, it would be nice
to make some sort of a prediction about how much money could be expected at
time t. If the expectation at time t is that there is no net gain (or loss) compared
to the amount at time s (Xt Xs = 0) then the game is fair. If, instead, more 11
11 Being measurable with respect to a filtration is not sufficient for the process to be measur-
able with respect to T. Also, there are processes that are measurable, but whose sub-domain
is not included in the filtration (and therefore are not adapted).
35

money is expected, then the game is super-fair. Having less would be an unfair
game. With this in mind, a martingale is defined as follows:
Definition 2.29 Let X be a real-valued, stochastic process and let Tt be a filtra-
tion on its probability space. X is a martingale with respect to Tt if:
1. Â£[|X(t)|] < oo, Vf > 0,
3. E[X(t)\Ts] = X(s) a.s., Vs < t.
To define a sub-martingale and a super-martingale replace the equality in
requirement (3) with a > and a <, respectively. In the context of the game,
a martingale is a random variable that describes a fair game. In other words,
E[Xt = 0. Sub-martingales and super-martingales corresponds to
the super-fair (E[Xt XS\TS\ > 0) and unfair (E[Xt XS\TS\ < 0) games,
respectively.
2.4.3 The Weiner Process
The most common stochastic process used in practice is a Brownian motion
or Weiner process. It is named after Robert Brown who studied the irregular
motion of pollen grains in a liquid.12 The Weiner process was used to model the
position of a pollen grain, u>, over time.
Definition 2.30 A stochastic process X that has the following properties is
called a (one-dimensional) Brownian motion or Weiner process:
12Technically, the process is called a Weiner process after the mathematician; while, the
term Brownian motion refers to the movement of the particles. Weiner is responsible for the
mathematical characterization of this motion.
36

1. X starts at zero i.e. A(0) = 0.
2. X is a.s. continuous (that is, its paths are a.s. continuous).
3. X has independent increments. That is, Xtl,Xt.2 Xtl, ...,Xtk Xtk_1 are
independent for all 0 < t\ < t2 < < tk.
4. X has independent Gaussian increments, that is, the random variable Xt
Xs ~ N{0, t s), t > s.
Actually, the above definition is a standard Weiner process since it starts at
zero. The general definition allows the process to begin anywhere, as long as it
does so with a probability of one (P[A(0) = x] = 1 where x is any valid starting
point). For the purpose of this thesis, the standard Weiner process is sufficient,
so it will be referred to as simply the Weiner process.
An important property of the Weiner process is that its paths can be chosen
to be continuous a.s. (Condition (2) above is not needed but is usually given
anyway). This follows from the following theorem:
Theorem 2.31 (Kolmogorovs Continuity Theorem) Let X be a stochas-
tic process, and suppose that for all T > 0 there exist positive constants a, ft, D
such that
E[\Xt-Xa\a] < D \t s|1+/?
where 0 < s and t Proof: See Kolmogorovs Lemma in [38].
We show that the Weiner process satisfies the hypothesis of this theorem.
37

Lemma 2.32 Let Wt be a Weiner process on K, Wq = 0. Then, Vm e I
E[etuWt] = e"5u2t.
Proof: Recall that Ws-Wt~ N{0, s t), so Wt = Wt W0 ~ N(0, t).
E [juWt] = J(2irt)-hiuxe^dx
R
f , ,1 2tiuac-l2-t2u2-H2u2
= / (2irt)~2e------sdx
R
/ I -(tixi-x)2 -u2t2
(27rt) 2e 2t e 2* dx
R
-n2t /* . . 1 -(tiu-x)2
= e 2 / (27rt) 2e a* dx
R
_u2f
= e 2 -1 since we integrated over the CDF for N(iut,t)
-u2t
= e 2 .
Lemma 2.33 Let lFt fee a Weiner Process on M, W0 = 0. Then,
E[W?]=-^i?;keN.
Proof: We begin with the result from the Lemma 2.32:
E [eiuWt] =
Then, by expanding e into the Maclaurin series and by linearity of E[-}:
n!
E
n=0
(mt-Fy
n!
38

E
71=0
inunE[Wll]
(-1)nu2n^
n\
= E
t J \ ^ 2K
n=0
n!
In particular, consider those terms in the expansion with the same power of
n = 2A;:
{2k)\ = k\
After isolating E[Wfk},13 we obtain the desired result:
E [wl = Â¥%tt-
Theorem 2.34 There exists a continuous version of the Weiner process.
Proof: From the Lemma 2.33, if a = 4, = 1, and D = 3 is chosen, then by
Kolmogorovs Continuity Theorem we have the result.
Whats interesting is that even though we can work with the Weiner process
as a continuous process, it turns out that its paths are nowhere differentiable.
This is true because the sample paths have unbounded variation (see Section
2.4 in [28]).
Many solutions to differential equations that contain stochastic processes
require that the process is a martingale.14 Fortunately, the Weiner process is
such a process.
Proposition 2.35 Let Wt be a Weiner process, then Wt is a martingale with
respect to the natural filtration J-fi generated by {Ws : s < t}.
13Note that i~k = (l)fc which is how we rid ourselves of the signs.
14 Some require only a semi-martingale.
39

Proof: We must show the three properties given in Definition 2.29. Prom the
Cauchy-Schwarz inequality, we can show that the expected value of the absolute
value is finite. For any t,
E[\Wt\] < yjE[W?]
= yjE[(Wt Wo)2}
= v^
< oo.
Consider the natural filtration T]v. Then, for s > t
E\Wa\jf] = E[WS -Wt + Wt\F}
= E[Wa Wt\T} + E[Wt\T)
= E[WS Wt] + E[Wt\^v]
= 0 + Wt
- Wt.
We get that E[WS Wt |JC] = E[WS Wt\ because Ws Wt is independent
of anything in T. The identity E[Wt\E^ ] = Wt is established since Wt is
measurable. Both of these follow from Theorem 2.23. Trivially, Wt is Tf -
adapted since Tf; was built from Wt.
2.5 Stochastic Differential Equations
2.5.1 The Ito Integral
The impetus behind developing the theory for stochastic differential equa-
tions is to take a differential equation that serves as a model and add a noise
40

term in order to make the model more robust. In the simplest case, consider a
first-order ordinary differential equation (ODE):
which can be rewritten in symbolic differentiable form (due to the flexibility of
Leibnitz notation) as:
Alternatively (and arguably more correctly) this can also be expressed in the
integral form:
A unique solution exists provided that the function a is Lipschitz continuous (a
sufficient condition see, for example, [13]). Einstein added a second (Gaussian)
term to the equation which was later refined to be a Weiner process by Weiner
leaving the integral equation:
The problem encountered here is that the second integral cannot be a Lebesgue
or Riemann integral because Ws is not differentiable and does not have bounded
variation. So, a new integral needed to be created, and many were. This appli-
cation uses Itos integral.
Consider the case where b(t,x) = b is a constant. In order to appeal to the
existing notion of calculus we woidd like to have
dx = a(t,x)dt.
t t
to
to
41

Thus, as is done in elementary measure theory, the integral on [a, b] for a step
function b(t,u>) = b3 on tj < t < tJ+\ where [a, b] is divided into subintervals
a = to < ti < < tn+1 = b can be defined to be:
0
f b(t,u>)dW,(u) = Y,b, - Wt,(u)).
i j=>
This works well. Note that in this case b(t,u>) is nonrandom since it does not
depend on u. When a; is a factor, i.e. b(t,ui) = bj(ui) problems arise. Consider
the two cases:
1
1
/n
b(t,u)dWs(u) (Wtj+1(u;) wtj(u))i and
n i=i
i
2. f b(t,u)dWs{u) E bj+i(u) (Wtj+l(uj) Wti(o;)).
n .7 = 1
o i-
Both approximations seem reasonable as long as the discretization is fine enough.
However, if bj(u>) Wt{u) it is shown in [35] (Example 3.1.1) that by properties
of the expectation and Weiner process
T
Y.E lWi M M wh M)] = 0; and
7=1
1. E
1
J b(t,uj)dws(u)
2. E
1
J b(t,L,)dW,(u>)
Y E fttiM - W>M))] = T
7=1
This indicates that the paths of the Weiner process have such large variations
that a Riemann-type integral will not suffice.15
15Actually, the paths of the Weiner process have a.s. infinite variation[35].
42

It turns out that one possible (Itos) solution to the above problem is to
restrict the class of functions, b(t,uj), to a specific set and to only approximate
functions by their left end-points in each discrete interval.
Definition 2.36 Let V = V(S,T) be the class of functions
f(t, uj) : [0, oo) xQ->R
such that
1. f(t,uj) is CxE-measurable, where C denotes the Borela-algebra on [0,oo).
The second condition guarantees that the function is nonanticipative that
its value at time t only depends on events at or before that time. The third
condition states the / must have finite second moments, which is a necessary
requirement in the extension of the integral from simple functions to those func-
tions in V. To define the Ito integral, begin with all simple functions, d, in
V:
where each bj must be Ttj-measurable because d E V. For these functions,
define the integral as
T
3. E / f(t,uj)2dt < oo.
,s

j
T
j> 0
S
43

Theorem 2.37 (Ito isometry) Let {>{t,uj) be abounded simple function. Then,
I Tr V T /*
E = E / tf{t,uj)2dWt{uj)
\s J J _s
Proof: Lemma 3.1.5 in [35].
Approximate a function / G V by simple functions to define the Ito inte-
gral.16
Definition 2.38 Let f G V(S,T) and let (dn) be a sequence of simple functions
that converge to f in the following manner
T
(f{t,uj) dn(t,ui))2dt 0 as n > oo.
Define the ltd Integral by
T T
I f(t,u)dWt(uj)= lim f dnft, u)dWt(u>).
s s
where the limit is taken in the space of L2(P).17
An important property of the Ito integral is that it is a martingale.
16This is quite similar to how the Lebesgue integral is constructed.
17
It
1
turns out that Jdn(t,u>)dWt(aj) forms a Cauchy sequence in Lr{P). As review, define
the L2 -norm of a r.v; X : Q
as
||X||2 = I j \X(uj)\2dP(uj)
1/2
. Then the space is
L2(P) = L2(Cl) = {X; ||AT||2 < oo}. More details on the convergence can be found in Chapter
S of [35].
44

Theorem 2.39 Let f(t,u) 6 V(0, T) for all T. Then,
t
o
is a martingale w.r.t. the natural filtration ofWt.
Proof: Oksendal proves this as a corollary on page 33 of [35].
This means that its expected value is zero between any two times s and t:
Corollary 2.40 Let f(t, u) G V(0, T) for all T. Then,
for 0 < S < T.
Proof: This follows from the definition of the martingale (Definition 2.29).
2.5.2 Ito Stochastic Differential Equations (SDEs)
Now that the Ito integral has been presented, we are able to define a dif-
ferential equation that describes the motion of a small particle suspended in a
moving liquid. The uncertain bombardments on this particle at the molecular
level invite us to describe the diffusion of the particle as a stochastic process. In
this case, we use a Weiner process as was originally selected by Einstein.18
Definition 2.41 A (vector) ltd stochastic differential equation or ltd
diffusion has the form
18A different process (e.g. Poisson) might actually be a better description of diffusion. Very
little research has been done into the distribution of diffusion, so almost by default the Weiner
process is chosen.
dX(t) = a(t, Xt) dt + B(f, Xt) dW(t)
45

which is to be interpreted as
X(t)-X{t0)= f a(T,x)(r)dT+ f B(r, X(r)) dW{r)
Jta V ' Jto
where
Xt is the position of a particle at time t;
a describes the convection or drift of the particle;
B describes the diffusion of the particle; and
W is a Weiner Process.
The reader is referred to [35], [28], or [17] for more properties of the Ito
diffusion.
2.5.2.1 Itos Lemma
In ordinary calculus problems the definition of the Riemann integral is not
usually directly useful. Instead, the fundamental theorem of calculus and the
chain rule are relied on heavily to find solutions. The same situation exists with
Ito integrals, but there is no differentiation theory available (due to the nature
of the Weiner process) to provide something as powerful as the fundamental
theorem of calculus. However, ltd was able to construct an integral version of
the chain rule now known as ltds formula.
It turns out that the ltd integral can be defined over a larger class of functions
than V by relaxing the second and third conditions in Definition 2.36.
Definition 2.42 Let W = W(S,T) be the class of functions
f(t, u) : [0, oo) x Q > R.
46

such that
1. f(t,cj) is CxT-measurable, where C denotes the Borel o-algebra on [0, oo);
2. f(t,u) is Jit-adapted where 7it is an increasing family of a-algebras for
which Wt is a martingale;19 and
From this we define a general stochastic differential equation known as the
Ito process:
Definition 2.43 Define the ltd process as
where v EW and u is Jit-adapted (as in Definition 2.42).
Finally, the one-dimensional version of Itos formula (called ltds Lemma)
can be given.
Theorem 2.44 (Itos Lemma in one-dimension) Let Xt be an ltd process,
and U(t,x) be a function in C2([0, oo) x R). Then Yt = U(t,Xt) is again an ltd
process given by the formula
3. P / f(t,u>)2dt < oo = 1.
_s
0
0
Proof: Theorem 4.1.2 in [35].
The vector version of the lemma is used in this thesis:
19This is a relaxation since C Tit-
47

Theorem 2.45 (Itos Lemma) For a sufficiently smooth function
U C1i2([0,T] x Rd)
U : [0, T\ x
ik and Yt = U{t, Xt)
If Xt is an n-dimensional ltd process, Yt is given componentwise by,
dYf> =
(dUq(t, Xt) , ~sdUq(t,Xt)
^ dt + ^ ( Xt) dXi
d m
z i,j=1 Z=1
m d
8'XidXj

, dxi
i=i i=i
Proof: See, for example, Theorem 4.2.1 in [35].
Q= 1,2,- ,k
2.5.3 Fokker-Planck Equation
The Fokker-Planck Equation describes the evolution of the probability den-
sity function of the position and velocity of a particle through time. Named after
Dutch physicist Adriaan Fokker (1887-1972) and German physicist Max Planck
(1858-1947), it was first used to describe the Brownian motion of a particle in
a fluid. The equation was developed by Einstein in a monograph on Brownian
motion.
Consider the following heat equation for functions f(x,t) G C^q (IT):
/(*t)
9f(x,t) 1 d2f{x, t)
dt 2 dx2
(2.19)
Think of a fixed point in three-dimensions where heat is being generated. As
time progresses the heat will dissipate into the surrounding area. This rate
48

of diffusion is described by the heat equation where f(x,t) is the probability
density function of the percent of heat at point x at time t. It turns out that
the density is Gaussian.
The heat equation is purely a diffusive action. Adding a convective term
yields the Fokker-Planck Equation (in divergence notation):
Definition 2.46 Define the Fokker-Planck Equation as
+ V [(a(t, x) VTDT(t, x)) P(x, t)} V D(f, x)VP(x, t) = 0
where
a : [0, T] xl^ is vector function called the drift coefficient;
D : [0, T] x Md > Kdxm is the d x rn-matrix function called the diffusion
coefficient; and
P(x, t) is the probability density function.
The Fokker-Planck Equation is related to the Ito SDE through Itos Lemma.
Some sources refer to the Fokker-Planck Equation as the generator of Xt. The
formal definition is
Definition 2.47 The (infinitesimal) generator, G, of the ltd diffusion Xt is
defined as
Gf(x) lim
qo
for the set of functions f : Kn
R
t
> R where the limit exists.
49

Example 6. Consider the Ito diffusion
dXt = dWt
then the generator of Xt is the heat equation (2.19). This can be verified using
Itos Lemma.
2.5.4 Fokker-Planck PDE to Ito SDE
Using Itos Lemma the Particle Trajectory SDE can be transformed into a
Conditional Probability Density PDE (a Fokker-Planck Equation). This section
outlines the process of the transformation.
Begin with a multi-dimensional ltd SDE with unknown convection function a
and diffusion function B (Definition 2.41). For the sake of notation, let C(t, x) =
|B(t,T)BT(f,T). By Itos Lemma we are allowed to choose a function U from
which we can create a new process Yt = U(t,Xt) so long as U comes from
the class of functions C'1,2([0,T] x M.d). For this case, let U{t,x) = g(x) where
g : R1 and g ^ 0. Because g is not a function of t the first term in (2.45)
will be zero and the rest simplifies to
a(t, Xt) Vg(Xt) + Y, Xt)dxidxjg(Xt) J dt + Vg(Xt) B(t, X)t)dWt.
i /
(2.20)
The expected value of the Ito integr al is zero (Corollary 2.40), so by taking
expectations of both sides the term Xg(Xt) R(t, X)t)dWt will drop out. The
expectation of the left hand side is by definition
E
' d_
dt
9(Xt)
/ li, ((AS)P(x,t)dx) = J g(x)^Â£(x,t,)dx
50

where P(x,t) = P(:r, t] Xo, to) is the conditional probability density of the posi-
tion of a fluid particle, Xp, at time t which started its journey at (x0,to). Also
from definition of expected value (and linearity of the integral) the right hand
size of (2.20) can be rewritten to give (2.20) the form
j g(x)^-(x,t)dx = j a(t,Xt)-Xg(Xt)P(x,t)dx
+ f Cij(t,Xt)dxidxjg(Xt)P(x,t)dx.
i,j=1 J
Apply integration by parts to the first term to the right of the equal sign, letting
u = a(t,x)P(x,t) and Vu = Xg(x). By assuming that the surface terms are
zero and moving everything to one side yields
/
dP
(x,t) + V (a(t,x)P(x,t)) ^2 Cij(t,Xt)dxidxjP(x,t)
i,j=1
g(x)dx = 0.
According to concepts in measure theory (see [42]), this integral can only be
zero if g(x) is zero (which cant be since g was chosen to be nonzero) or if
dP d -
(f,t) + V (a(t,,x)P(x,t)) - Xt)dxidxjP(x,t) = 0.
i,j=1
This is the Fokker-Planck PDE and is identical to Definition 2.46 if written in
divergence notation:
dP(2 ^ + V [(a(t, x) VTCT{t, x)) P(f, t.)] V C(t, x)VP(Â£, t) = 0.
2.6 Heaviside Step Functions and the Delta Function
In the control of fluid particle movement across an interface where one
porous media meets another, a Heaviside (or unit) step function is encountered.
51

Definition 2.48 A Heaviside step function is the real-valued function:
H(t)
1 t > 0
0 t < 0
t R
which is optionally defined at zero.
The function can be scaled:
Vo + oqH (to t)
where yo is the y-value before time t0 when a jump of magnitude ao will take
place.
It is sometimes the case that the derivative of the Heaviside step function
is needed. Clearly, the derivative is zero when x 0, but it is not directly
differentiable when x = 0. The following shows how to express the derivative of
a function with jump discontinuities (such as the Heaviside step function) as a
distributional derivative.
The theory is developed using two types of functions: test functions, and
generalized functions. Test functions are infinitely differentiable functions that
have compact support (they are only nonzero on a compact set20), denoted here
as u(t) G Cq^R1). Generalized functions are a generalized notion of regular
functions, and are often used to make discontinuous functions appear as contin-
uous functions. Technically, generalized functions can be defined as continuous
linear functionals over the test functions such that all continuous functions have
20Since this project works over the real numbers, by the Heine-Borel theorem a compact set
is one that is both closed and bounded.
52

derivatives which are themselves generalized functions.21
Given a test function, i/(t), and generalized function, u(t):
I u'(t)u{t)dt = f
J M1 J R1
u(t)v'{t,)dt + u(oo)is(oo) u(oo)u(oo)
by Greens Formula. Notice that the last terms equal zero because v(oo) =
u{oo) = 0 since v has compact support. What is left is
(2.21)
which is often given as the definition of a generalized function.
Definition 2.49 Define the 5-function as the function satisfying the equality
j 5(t)v(t)dt = zz(0).
Jr1
It can be shown that H is a generalized function, so by (2.21)
/ H\t)u{t) = - / H(t)u'(t)dt,
J R1 Jr1
then by (2.48) and since u{oo) = 0
/ H\t)u{t) = - ( H(t)v'(t)dt
J K1 Jr1 OC /*
= / v'(t)dt
o
= WO).
By Definition 2.49, uniquely
H\t) = 5(t) (2.22)
21 That is, a generalized function maps a test function into M. It must be continuous, and it
must be linear. Finally, for every continuous function there exists a derivative of that function
which is also a generalized function
53

is the distributional derivative of H.
Informally the 5-function can be thought of having the value of oo when
t = 0 and is zero elsewhere. Of course, this not actually a function and so it can
only be meaningfully discussed when it is multiplied with a test function and
integrated (as above). The delta function is itself useful to model an impulse
(such as the force transfer from a bat to a ball) or point charge.
54

3. The SDE Model
3.1 Motivation
The experimental simplification of an aquifer is a tank filled with different
kinds of sand (heterogeneous) and then saturated with water (see Figure 3.1).
Experimental Tank
Figure 3.1: Experimental Tank Setup
The flow of pollutants in the aquifer can then be modeled by injecting either
a LNAPL (Light Noil-Aqueous Phase Liquid lighter than water) or a DNAPL
(Dense Noil-Aqueous Phase Liquid heavier than water) into the tank (see
Figure 3.2).
55

Figure 3.2: LNAPL Plume in Experimental Tank
Partial differential equations are used to model the saturation of two-phase
flow in porous media. This originated from work by Bernard Kueper and Emil
Frind in [32] with the development of the following equations:
0
DS\y
~dT
d_
dxi
k-ij kj*w
PW
dPw
dxj
PwQ
dz
9xj)_
= 0,
0
ds
NW
d,
d_
dxi
kij krNW
pNW
dP,
NW
dz
dxj
Pnw9z
dxj 7 j
= 0.
Here, the saturation of the wetting phase (Sw) and the non-wetting phase
(Smw) are found by solving the equations simultaneously often, using a Fi-
nite Element Method (FEM). A direct FEM is known to (typically) have the
following problems:
Mass loss. Due to numerical roundoff the total mass in the system may
56

not remain constant.
Numerical diffusion. Buildup of numerical error can cause innappropriate
movement. For example, the NAPL is allowed to move into a fine sand
before it should (see Figure 3.3).
Figure 3.3: Premature Breakthrough of Plume from Numerical Diffusion
Dean and Russell (refer to [12]) proposed using a stochastic differential equa-
tion to model the movement of individual fluid particles. The idea is to combine
deterministic elements Darcys law, the continuity equation, the Brooks-Corey
model, and a PDE of the velocity of the fluid into a unified stochastic descrip-
tion which better describes the underlying physics. Their approach mitigated
the above issues:
No mass loss. A finite number of fluid particles are tracked each of which
has a constant mass.
Interface problems can be avoided. Fluid particle movement can be micro-
managed to avoid unrealistic movement.
57

Additionally, a stochastic approach is a better conceptualization of macroscale
movement since diffusion is inherently uncertain at this scale. At the microscale
a fluid particle is subject to erratic behavior because it collides with other ob-
jects solids, or other fluid particles (see Figure 3.4). At the macroscale there
is no way to deterministically predict these diffusive actions, thus suggesting a
probabilistic solution. Naturally, a SDE model will also have drawbacks, but
the above arguments were provocative enough to try this approach.
Figure 3.4: Microscale Effects that are Stochastic at the Macroscale
3.1.1 The Saturation-Probability Relationship
Let V be a Representative Elementary Volume (REV) of the porous media,
Vv be the void space in V, and Np be the number of fluid particles needed to fill
the REV (see Section 2.1.3). The volume of the o-phase in V is then
Va = xVfp = 0 where Vfp is the volume of a fluid particle. If Nt is the number of total particles
in the system then the saturation of the n-phase is shown to be
58

Therefore, the probability of the ct-phase being in the REV containing x at time
t is:
Pa(x,t) = Sa(x,t). (3.1)
The underlying assumption in this model is that the positions of the fluid par-
ticles are stochastic, and Pa describes the law.
3.2 Fokker-Planck Equation for Fluid Flow
In this section, a Fokker-Planck Equation is derived that describes the flow
of a fluid particle of the a-phase. The formulation is based on Darcys Law
(2.8) and the continuity equation (2.14). Building such an equation is desir-
able because of the intrinsic link between the Fokker-Planck Equation and Ito
Stochastic Differential Equation.
Recall from equation (2.8) that the Darcy flux of the n-phase can be de-
scribed as
Qa = (Vpa 7az).
Pa
This can be written in terms of the hydraulic conductivity tensor K using the
relationship between itself and the intrinsic permeability given by (2.6):
qa = (Vp 7az).
la
By assuming that the density, pa, of the fluid is constant1, the specific weight
will also be constant. Multiplying through by the specific weight gives
qa = I 1 Assuming that the density of the fluid is constant is equivalent to calling the fluid incom-
pressible. This is an idealization used to simplify analysis. In reality all fluids are compressible
to some extent.
59

where the term ^ is denoted using the symbol ipa. pja is dependent on the
volumetric fraction, 0a (as the volumetric fraction increases the saturation of
the phase must increase which causes an increase in pressure of that phase):
Pa i.0a^)
^Pa{0 a)
7a
Take the gradient to obtain a formula for \/ipa (which exists in Darcys Law):
dtpa(0a)dda dipa d0a
dxi dxi
d0a dxi
=
where is longhand for V#Q.* 2
OXi
The hydraulic conductivity Ka is also physically dependent on the volu-
metric fraction, so it is now written as Ka(0a). With substitutions, (3.2) now
appears as
Qct
Here, the term
'dip'' d0a j ) VK + Ki)z (3.3)
) = - (3.4)
is referred to as the capillary diffusivity. It will be used (in part) to determine
dt/t
d\$a
the magnitude of diffusion of the fluid particle. In this application, the term
is the slope of the Brooks-Corey curve.
Insert Darcys Law as given in (3.3) into the continuity equation, (2.14),
and use the relationship between saturation and volumetric fraction, (2.11), to
find
dsa
dt
l
+ V -Ka(0a)zSa ) V D(dQ)VSQ = 0
t/n
(3.5)
dtp
t dxd
2The notation is shorthand notation for the gradient and is identical to I 8x2 1 ip.
dXi
60

From Section 3.1.1 this is related to the saturation in the REV and is used to
obtain a Fokker-Planck PDE for the flow of the a-phase
+ V (^Ka(0a)zPa(x,t)^j V D(0a)X(Pa(x,t)) = 0. (3.6)
3.3 PDE to SDE
By comparing the PDE for the a-phase, (3.6), to the general form of the
Fokker-Planck Equation, (2.46), the formula for the vector function a and the
d x ra-inatrix function B can be obtained.
a(t,x) = -j Ka(0a(t,x))z + VTBr(0a(t,x)) (3.7)
0a{t,x)
C(t,x) = ^B(t,x)BT(t,x) = D(0a(t,x)) (3.8)
where D is the capillary diffusivity, and K is the hydraulic conductivity tensor.
By the relationship given in Section 2.5.4 there is an associated SDE for the
flow of a fluid particle of the a-phase:
dXt= ( ...-1 ^ Ka(0a(t,Xt))z + VTDT(6a(t,Xt))) dt + B(t,Xt)dWt. (3.9)
\0a{t,Xt) J
3.3.1 Specifics
Further manipulation of the SDE is needed to actually compute particle
trajectories. Here, the convective term of the SDE,
1 Ka(0a(t, Xt))z + VrDT(0a(t, Xt)),
0a(t, Xt)
is expanded into pieces that more simply describes the physics.
61

Observe the following two equalities:
<7t qw + <1nw, and (3.10)
Pc = Pnw Pw (3.11)
Equation (3.10) states that the total flux of the system, qt is the sum of the
wetting and non-wetting phase flux. Equation (3.11) asserts that the capillary
pressure, pc, is the difference of the pressure of the non-wetting phase, p^w, and
the pressure of the wetting phase, pw-3 These equations may be combined with
Darcys Law. Letting \a = ,
Qn = fn where fn = \w+xN an(l A = Xw+\N are fracdtonal flow coefficients.
Substitute into the continuity equation, (2.14), which creates a convection-
diffusion PDE for the saturation of the non-wetting phase:
dS.
N
dt
+ V-
1
T
N
(.fnQt Ayu/ki* + \~fNkz)SN
-V-
Ak dpc
~^dS~
N
VSN = 0. (3.12)
Recall that the saturation is related to the probability distribution function
for the fluid particles through (3.1). By using this relationship, (3.12) can be
converted into a Fokker-Planck PDE which (by Itos Lemma) is related to a new
non-wetting Phase SDE:
dXt = (J- (Ayvkz- A7wkÂ£+ fnqt) + VTDT^) dt + B dWt. (3.13)
3This equivalence is only believed to be true at an equilibrium state. Currently, researchers
are trying to determine how to model capillary pressure in a system that is not stable.
62

The volumetric, fraction can be found easily since the model keeps track of
the placement of fluid particles, and the size of the system is known. Also, the
fractional flow coefficients, specific weights, and permeability tensor can all be
calculated from input data. What remains is the value for the total velocity
field, qt, and the jump term VTDT.
3.3.1.1 qt Velocity Field
The total flux or velocity field, qt, of the system is obtained by solving a
PDE for the piezometric head of the wetting phase. Take the gradient of the
piezometric head, (2.5), and multiply through by the wetting-phase density to
define
vv = lw<\>w = Iwz + Pw
Vvc = Vpw ~ Iwz.
This exists in Darcys law, (2.8), so by substitution
_ K w, T.
w
Iw
This, along with the equality, Sw + Snw = 1 can produce a new form of the
continuity equation:
-V-
. 7 w
V(j>w
- V
K
NW
In
V(j>NW
= 0.
By (3.11) we are left with the PDE:
-V
K w Ky
Iw In
Vd'.u, = V -Vpc V KNz + V
7 N
(3.14)
IN J
63

Once this is solved for vv the pressure of the wetting phase, pw, can be found.
This in turn substituted into (3.11) to find the non-wetting pressure, and then
Qn and qw can be calculated from Darcys Law. Finally, qt = Qn + qw delivers
the the velocity field.
3.3.1.2 Interface Control with the Jump Term
The control of the flow of particles across an interface between two sands of
different permeabilities is accomplished using the jump term
' V^ dOn n dSn
i
J
0
In a homogeneous domain, D is assumed constant for a specific time and so
VrDT = 0. In a heterogeneous domain, D will have a different constant value
for each type of sand. At the interface between two sands a jump occurs which
looks like a Heaviside step function in one-dimension. Consider the 2-direction
and assume that D is diagonal4, then
D33(s, Xa) = K(s) + a0(s)H(Xs zx)
where K(s) is the value of the diffusivity of the sand that the material is in at
time s, ao(s) is the difference in diffusivities (do = Dc D/ in which c is for
coarse and / stands for fine), z\ is the coordinate at the interface, and Xs
is the position of the particle.
The model requires that derivative of D33 is found at the interface. Using
the theory introduced in Section 2.6, this is only possible using distributional
4 Assuming that D is diagonal is the same as stating that the permeability of the material,
k is diagonal. This is convenient, but not always true since there may be directional coupling
caused by e.g. grains which do not have coordinate symmetry.
64

derivatives (if it is multiplied by a test function and then integrated). Let
v E Co(f2) be a test function, d(il) be the boundary of the interface, and let
nt be a vector normal to the interface. By Greens formula5
f DzzdzudVt + f v dzDzzdtl = (f v Â£>33n+d7. (3.15)
Jn Ju JdQ
This can written in distributional notation:
(D33,dzu) + (dzD3 3,1/) = (l[D33}}nt8(dn),v) (3.16)
where [[Z)33]] = D33,c Dzzj. According to the usual definition (Definition
2.21) of a generalized function, the first term is equal to ~(dzD:iZ), u) where dz
represents the distributional derivative so that (3.16) is also
-{dzDii,dzv) + (dzDzz,u) = ([[Â£>33]]n+<5(dQ), v)-
In this way,
dzDzz = dzDzz [[D33]]n^(cK2).
On either side of the interface D is constant so dzDzz = 0. Finally,
(BzD33) = -([[Dzz]}n+S(dn)). (3.17)
If this is completed for the other directions
f VTDn
Jo
da =
L[[Du(t)R^
-[[D22{t)}}n+
^-[[Aj3(*)]]w+ }
d{dm-
5Greens formula: / q VvcKl + / v V qdQ = vq- n+dy.
Jo. Jq JdQ.
65

The effect is that VTDr will be in the opposite direction of (A7jvkz A7wkz + fnqt).
Thus, the particle jumps back into the coarser sand. The normal components
jump the particle back at an angle consistent with a physical collision (for ex-
ample, think the collision of a pool ball and the edge of the pool table).
Recall that D is dependent on the saturation (for constant porosity, other-
wise volumetric fraction) of the non-wetting phase. The relationship is given in
Figure 3.5.
Non-wetting Phase Saturation
(in course sand, fine saturation is zero)
Figure 3.5: Dc Df versus Non-Wetting Phase Saturation
As seen in the figure, the diffusivity coefficient allows unrestricted flow once
full saturation is reached. According to Brooks-Corey theory unrestricted flow
will commence once the non-wetting fluid reaches the entry pressure. This corre-
sponds to a saturation less than one (see Figure 3.6). So, once the fluid particle
reaches the entry pressure of the fine sand, the jump term is taken to be zero
66

and the particle can move freely into the fine sand.
Figure 3.6: Entry Pressure and Corresponding Saturation
A recommended addition to the model (other than those proposed in Chap-
ter 4) is to multiply the jump term by another step function to accurately reflect
the WRC theory and bring it explicitly into the model. To be more exacting, let
Pd,c and Pdj be the entry pressure of the coarse sand and fine sand, respectively.
Following the visual in Figure 3.6 and using (2.17) and (2.9) define the entry
saturation, S^w as
c*
oNw
= 1 -
Pd,c
Pd.f
(1 Sr\y ) T Sr]\?
(3.18)
67

The one minus is taken so that this is a non-wetting phase saturation since the
Brooks-Corey formula gives the wetting saturation by (2.10). A step function
can then be formed as
K-(Snw)
(3.19)
0 Smw > S
NW
and then multiplied into the jump term
/
f
V D da fC(SNW)
[{D22(t, 5Vw)]]fhj
+ \
S(d(Q)).
^-[[-DasCT'S'ivw)]]^ }
The jump term is effectively truncated at the entry saturation as is shown
in Figure 3.7.
Non-wetting Phase Saturation
(in course sand, fine saturation is zero)
Figure 3.7: Dc Df with Brooks-Corey Theory
68

3.3.2 Computational Model
The above model was realized in C code by D. Dean. The computational
model is still under construction (as is the model) so it does not boast exceptional
speed or accuracy. Currently, a standard Galerkin method is being used to solve
the PDE for the velocity field which is known to be inefficient, especially in
three dimensions. The stochastic differential equation is being solved by finding
a weak solution6 using a truncated Taylor approximation method:
m
U+i = P + Yn)A + Y.
j = 1
The details for this method can be found in [28] and [29]. Figure 3.8 shows an
example run of the model.
6A strong solution to a SDE is one where the Weiner process has been chosen ahead of
time. A weak solution is found if any Weiner process can be chosen to solve the problem.
Therefore, all strong solutions are weak solutions, but there exist examples where a strong
solution does not exist where a weak: solution does. See [35] or [28] for further discussion. A
weak solution is sufficient in this case since there is no motivation to use one Weiner process
over another. Actually (see page 269 of [17]), no strong solutions exists for the Fokker-Planck
Equation.
69

Figure 3.8: Example Run of Computational SDE Model
70

4. Macro Percolation
4.1 Particle Size Distributions
The experiments at the Colorado School of Mines, from which this study
is based, use sand-blasting sands. These are much more uniform than what is
typically found in nature. The sands range from a #8 size to a #140 size which
indicates the (average) size of the grains (the lower the number the larger the
grains as seen in Figure 4.1).
Suppliers of sands typically conduct a sieve analysis of their product. For
instance, they will obtain a sample of their #8 sand and then sift it through
a number of different sieves (each with a smaller hole size than the previous)
71

I
Table 4.1: Grain Size Data from Sand Suppliers
Purported # and Percent that Passed Each Sieve
# Sieve Size #8 #16 #30 #50 #70 #110 #302
8 2.360mm 96.5 100.0 100.0 100.0 100.0 100.0 100.0
16 1.180mm 26.0 87.1 99.9 100.0 100.0 100.0 99.9
20 0.850mm 3.2 27.2 99.6 100.0 100.0 100.0 99.7
30 0.600mm 1.1 8.8 86.3 100.0 100.0 100.0 90.9
40 0.425mm 0.3 1.7 15.9 97.0 99.8 100.0 42.9
50 0.300nnn 0.1 0.4 1.8 57.0 97.9 99.5 20.2
70 0.212mm 0.0 0.1 0.6 18.0 59.2 95.5 6.4
100 0.150mm 0.0 0.0 0.3 4.0 22.9 77.5 1.5
140 0.106mm 0.0 0.0 0.1 1.0 3.7 33.5 0.4
200 0.075mm 0.0 0.0 0.0 0.0 0.6 8.5 0.0
270 0.053mm 0.0 0.0 0.0 0.0 0.0 0.5 0.0
Pan >0.053mm 0.0 0.0 0.0 0.0 0.0 0.0 0.0
and then compute what percent passed through each sieve (see Table 4.1). This
provides a crude CDF of the grain distribution. Usually, experimentalists just
want to know that most of the sand is of the diameter they wish to use. In this
case, we wish to estimate the actual distribution.
Under the assumption that the particle sizes follow a lognormal distribu-
tion1, lognormal CDFs were fit to the data from Table 4.1 (data obtained from
[44]). This was done using the Nelder-Mead optimizer from the R analysis pack-
age with a SSE objective function.2 Table 4.2 shows the fitted parameters and
1 Other right-tailed distributions have been used for the particle size distribution, but the
lognormal is the most common. See footnote below on pore-size distributions
2SSE stands for Sum of Squared Errors which is given by: YltiVt Vif where y is the
actual data value and y is the predicted value from the distribution function. The Nelder-
Mead optimizer is a fairly robust, but not exceptionally fast, optimization method.
72

Table 4.2: Lognormal Distribution Parameters, SSE Fit, and Distribution
Statistics
111 dyyi Hd SSE Median Mode Mean S.D.
#8 -1.9568 0.28068 0.00012 0.14144 0.13072 0.14712 0.04212
#16 -2.3522 0.19828 0.00664 0.09515 0.09148 0.09704 0.01943
#30 -2.9937 0.16510 0.00035 0.05010 0.04875 0.05079 0.00844
#50 -3.5742 0.27793 0.00267 0.02804 0.02595 0.02914 0.00826
#70 -3.9528 0.30512 0.00437 0.01920 0.01749 0.02011 0.00628
#110 -4.4248 0.31472 0.00073 0.01198 0.01085 0.01258 0.00406
#30&50 -3.1599 0.33794 0.01399 0.04243 0.03785 0.04492 0.01563
statistics for the lognormal distributions of the various sand types.3
Figure 4.2 displays the curves and data points for the lognormal CDFs. The
fit seems reasonable despite a few outliers near the lower end of the distributions.
Also, noticing that the medians and means in Table 4.2 are relatively close
indicates that the data may be more normal than lognormal (when the mean
and the median are the same then there is no skew in the distribution). Knowing
the parameters gives us the PDF of the grain size data (see Figure 4.3).
3The optimization was done using values in centimeters instead of millimeters since the
model uses centimeters for other calculations.
73

Cumulative Density
Particle Size (cm)
Figure 4.2: Lognormal CDF Curves Fit to Grain Size Data
74

o
o
#8
#16
#30
#50
#70
#110
1 I I T
0.00 0.05 0.10 0.15 0.20 0.25
0.30
Particle Size (cm)
Figure 4.3: Lognormal PDF Curves Fit to Grain Size Data
4.2 Pore Size Distributions
Pore sizes (the size of the gaps in between the grains of sand) are also
commonly assumed to follow a lognormal distribution.4 This sort of information
is difficult to verify even in the most controlled conditions.5
4Sometimes a gamma distribution (e.g [36]), a beta distribution (e.g. [39]), a Wei bull
distribution, or some other distribution is used. What is generally recognized is that the true
distribution is skewed slightly to the right.
Experimentalists who work at this level often work with glass beads which are more
uniform than any soil or sand sample.
75

Rouault and Assouline (1998) in [40] proposed that a nonlinear relationship
exists between the grain size and the radius of the pore size:
d = urv (4.1)
where d is the diameter of the sand gr ain, r is the radius of the pore, and u and
v are constants that are related to the packing of the particles, the shape of the
particles, and the shape of the pores. If reliable values of u and v can be found
for different porous media, this model would be an extremely useful relationship
for two reasons: one, radius distributions are useful for flow predictions of fluids
through the media; and two, grain size distributions are relatively easy to obtain.
Assuming that the pore size distribution is lognormal then implies the par-
ticle size distribution is also lognormal by (4.1). This connection was originally
hypothesized by Shirazi and Boersma in [47] and was later experimentally vali-
dated by Shizawa and Campbell (see [46]).
4.3 Lognormal Water Retention Curves
In [30] and [31], Kosugi assumes that pore sizes are distributed lognormally
and then suggests that the Water Retention Curve is also given by a lognormal
distribution. His derivation is based on the Young-Laplace Equation, (2.15).
Given a capillary pressure head the effective saturation can be computed from
the following:
S.Qnh) = F.(^l^=) (4.2)
where
- Se is the effective saturation (defined as Se = (S Sr)/( 1 Sr));
76

- Fn is the normal CDF, (2.18); and
- hm and are the mean and standard deviation of In h.
This is a two-parameter WRC, whereas his first model in [30] used a bubbling
pressure as a third parameter.6
In 2003, Hwang and Powers (see [25]) combined the relationship between
grain size and pore size with Kosugis model. The end result is a WRC that
is based on the distribution of grain sizes. They used soil data for sandy soils
from the UNSODA soil database (see [34]) which is a collection of hydraulic
properties of undisturbed soil samples. The parameters u and v were used as
fitting parameters to optimize a RMSE (Root Mean Squared Error) objective
function.
The parameters hm and 07, are derived from the lognormal grain distribution.
The mean of the pore size distribution is found through the Young-Laplace
relation
. In dm In u
lnrm =-------------.
v
By definition of the standard deviation:
or = yjE [(lnr lnrm)2]
\( In d In u In dm In u VI
[l V V )\
In d In d
v
6In [31] the bubbling pressure is assumed to be zero since it has little affect on the accuracy.
77

Furthermore, by the Young-Laplace equation
In hm = In _ lnrm = A hi r
V Pwg )
which is the first parameter and then
crh = E [(In h In hm)2]
= E
{(A lnr) {A- hi/im))2
E [(lnr lnrm)2]
v
(4.3)
(4.4)
is the second parameter.
4.4 Water-NAPL Water Retention Curves
The work by Kosugi, Hwang and Powers focused on unsaturated soils. That
is, they developed water retention curves for an air-water system. Water reten-
tion curve experimental data is considerably easier to obtain for air and water
compared to water and NAPL. Because of this, researchers hope that WRCs
can be generated for any two phase flows by scaling an air-water WRC. Cur-
rently, work is being done at the Colorado School of Mines (CSM) to test this
hypothesis. For various sands, they are generating WRCs for air-water and
water-soltrol7 systems. The scaling being used is the ratio of the interfacial ten-
sions. For a given saturation, the pressure of the soltrol-water system, psw, is
obtained by multiplying the interfacial tension ratio times the pressure of the
air-water system, Paw'-
4>sw
7Soltrol is a chlorinated solvent that is sometimes polluted by industry. It functions as a
LNAPL.
78

This seems to be scaling well for a #30 sand (see Figure 4.4).
Using the lognormal parameters found in Table 4.2, we fit a WRC to the
new data from CSM (also shown in Figure 4.4).
Figure 4.4: Lognormal WRC Compared with Experimental and Brooks-Corey
WRCs
4.4.1 Macro Percolation in the SDE Model
As mentioned in Section 2.1.6 the behavior of the fluids is not well known
at low and high saturations. According to the Brooks-Corey model, the non-
wetting phase cannot enter into the fine sand until the entry pressure is reached.
This works well in the SDE model to create distinct interfaces. Figure 4.5
shows the tank setup and a resulting run of the model. In this example, the
79

contaminant plume pools along an interface between two sands until it reaches
a high enough saturation to invade the finer sand.
Experiment 9 Tank Design
O
ii
=&=
T5
G
fO
CO
#30
s and # 3
Sand
Sand #16
Sand #30
Sand #16
Sand TF7fi
O
i
=#=
T3
G
CO
Saturation Scale
Figure 4.5: Example Run of SDE Model with a Well-Defined Interfaces
80

At the microscale this sort of clean interface is artificial for two reasons.
First, based on how two sand types are packed together, the finer sand is likely
to fill in voids in the coarser sand creating an interfacial layer (see Figure
4.6). Secondly, if the sand is not uniform then there wall occasionally be larger
pores for fluid particles to move through before the entry pressure is reached
(see Figure 4.7).
Figure 4.6: Interfacial Layer Between Fine and Coarse Sands
t
t
O
O
o
o
Figure 4.7: Variation in Grain Size Causes Variation in Pore Size
81

There are mainly three types of models that describe what happens at the
microscale: Diffusion-Limited Aggregation (DLA) (see [37], [1], [53], [20], and
[24]), percolation (see [45], [52], [5], [16], [14], [49], and [41]), and fractal (see
[33], [22], and [23]). Of these, the percolation models are best motivated by
the physical structure of the soil. Typically, different pore geometries are con-
structed which are connected in a lattice. The fluid can then percolate through
the network entering each pore once it attains enough pressure. In sample scale
models the pore sizes can be randomly distributed ([39] uses a Beta distribution
and a Gamma distribution is used in [36]) to estimate the uncertainty of the
pore structure.
It would be computationally unreasonable to create such a three-dimensional
pore lattice on a scale of decimeters (or larger) since this would involve millions
of pores. As a microscale addition to the SDE model classic percolation modeling
is not realistic. Instead we propose that every time a fluid particle tries to cross
an interface, a random pressure is generated from the grain size distribution of
the fine sand using the theory given in Section 4.3. If that pressure is lower than
the current pressure at the fluid particle, then the model allows the fluid particle
to cross the interface. Thus, a few fluid particles are allowed to percolate into
the fine sand before the entry pressure is reached on the Brooks-Corey curve.
Let Vt be a stochastic process comprised of lognormal random variables
which takes on values from the pressure distribution. From (4.3) and (4.4), we
have that for a fixed t
Vt ~ N(hi hm, ah)pwg
82

where In hm and ah are parameters of the pressure distribution for the fine sand.
The parameters u and v are fit so that the lognormal WRC agrees with the
Brooks-Corey curve as shown in Figure 4.8.
o
o
o ~
o
0.0
0.2
0.4
0.6
0.8
1.0
Saturation of Wetting Phase
Figure 4.8: Brooks-Corey and Lognormal Curves for #16 and #30 Sands
Define the process
similar to what was done in (3.18). The idea is again to use a step function and
apply it to the jump term. When a particle reaches an interface, if a random
saturation is generated that is less than the current saturation of the non-wetting
phase, then the new step term will retard the jump term to allow percolation.
Specifically, from (3.19) define
0 Snw > S*
1 0 < S'jVH' 5: <5t*
83

and then let
/
VtDt da = 1C(Snw)k(Snw)
^ -[[Â£>ii (t, SW)]]rcÂ£ ^
[[D22(t, 5jvwO]]Wy
y-[p33^,S'ath/)]]^ y
wm
Visually, Figure 3.5 has been transformed into Figure 3.7 by the Brooks-
Corey theory. Now Figure 3.7 becomes Figure 4.9 (see Figure 4.9). In this
figure the dotted line represents the fact that a particle will not be forced back
by the jump term if it encounters a large grain in the fine sand which will
correspond to a large pore size into which the particle is allowed to move.
Non-wetting Phase Saturation
(in course sand, fine saturation is zero)
Figure 4.9: Dc Df with Brooks-Corey and Macro Percolation
The question is whether or not the microscale grain size variations will have
an impact on the macroscale model. From sample runs of the computational
model, Figure 4.10 shows what can happen with macro percolation. In the first
84

picture, the plume has broken into the finer sand by reaching a high concentra-
tion. The second picture indicates where additional movement as taken place
through percolation.
Figure 4.10: Before and After Macro Percolation is Applied to the SDE Model
4.5 Beowulf
A major goal is to parallelize the SDE model. In general, parallelization
is becoming a standard request for large computational models. Techniques
used are often very specialized to a specific algorithm, so custom work is almost
always done for any new model.
We initially ported the computational SDE model to run on a single node
of UCDHSC CCMs Beowulf cluster. The cluster consists of 36 separate nodes
(computers) connected together using a Dolphin Interconnect backbone.8 Each
8The SCI Interconnect is manufactured by Dolphin and installed as a 6x6 double torus. It
yields 220MB/s actual bandwidth with a 5/xs latency. A Scab MPI implementation is used.
85

node is running the Fedora Core 3 linux distribution on two Pentium-Ill pro-
cessors operating at 933MHz with 2GB of memory (see Figure 4.11).
Figure 4.11: UCDHSC CCM Beowulf Cluster
An algorithm is parallelized by first finding what portions of the code are
most costly. The cost is a function of how often the piece of code is called times
how much time it takes to run. The most expensive model components are then
analyzed to determine if they can be parallelized.
In initial runs of the model9, approximately 98% of all computation time
was spent computing pressures for the velocity field. The time to make this
computation appeared to be stable even when more particles were inputed into
the tank. The other 2% of the time was spent moving the particles which
9Initial runs were performed for a period of 30 virtual minutes where the velocity field was
recalculated every 10 virtual seconds. For the setup used (i.e. the type of sands, size of tank,
input rate, etc.) approximately 1500 fluid particles were injected into the tank.
86

seemed to increase linearly as more particles were added. So, in the long run
(i.e. hundreds of thousands of particles are in the system), the time to move the
particles would surpass the time to compute the velocity field (see Figure 4.12).
We first concentrated on parallelizing the velocity field code since, even with a
large number of particles, it could still have a significant impact on the overall
performance.
Figure 4.12: Computation Time versus Number of Particles for Movement
and Pressure Calculations
The pressure field (from which the velocity field is found) is calculated by
solving the PDE given by (3.14). As explained in Section 3.3.2 this is currently
done using a standard Galerkin FEM method. Essentially, the tank is broken
into a grid of computational cells. In each cell a system of equations is solved;
the results of which are inputed into a global system of equations (see Figure
87

4.13).10 The pressure is found for each corner of each cell by solving the global
system. The pressure inside the cell is then interpolated by the pressures at
the corners. This is a very rough description and the reader is referred to any
elementary book on the subject (e.g. [7]).
Global System
We assigned an even number of computational cells to each processor. When
the pressure field needs to be calculated the task of solving the local system of
equations for the individual cells can be done simultaneously. Overhead cost of
sending pertinent information to each of the processors and receiving the solu-
tions back to the main (root) processor will cut into the increased performance.
Despite this, we were able to decrease the overall computation time of the pres-
sure field by over 50%. This is shown in Figure 4.14 where the percent decrease
^Computational cells are not necessarily (or usually) cubic. Various shapes are employed
to give more accurate solutions.
88