Project Summary

United States

Environmental Protection

Agency

National Risk Management

Research Laboratory

Ada, OK 74820

Research and Development EPA/600/SR-97/102

NAPL: Simulator Documentation

Joseph Guarnaccia, George Pinder, and Mikhail Fishman

October 1997

A mathematical and numerical model

is developed to simulate the transport

and fate of NAPLs (Non-Aqueous Phase

Liquids) in near-surface granular soils.

The resulting three-dimensional, three

phase simulator is called NAPL. The

simulator accommodates three mobile

phases: water, NAPL and gas, as well as

water- and gas-phase transport of NAPL

contaminants. The numerical solution

algorithm is based on a Hermite

collocation finite element discretization.

Particular attention has been paid to the

development of a sub-model that

describes three-phase hysteretic

permeability-saturation-pressure (k-S-P)

relationships, and that considers the

potential entrapment of any fluid when it

is displaced. In addition rate-limited

dissolution and volatilization mass

transfer models have been included. The

overall model has been tested for selfconsistency using mass balance and

temporal and spatial convergence

analysis. The hysteretic k-S-P and mass

exchange models have been tested

against experimental results. Several

example data sets are provided,

including a setup of the artificial aquifer

experiments being conducted at the

EPA’s Subsurface Protection and

Remediation Division of the National Risk

Management Research Laboratory in

Ada, OK.

This Project Summary was developed

by EPA's National Risk Management

Research Laboratory's Subsurface

Protection and Remediation Division,

Ada, OK, to announce key findings of

the research project that is fully

documented in a separate report of the

same title (see Project Report ordering

information at back).

INTRODUCTION

The physical problem which is addressed

is the contamination of a pristine porous

medium as the result of releases of organic

liquids, commonly referred to as NonAqueous Phase Liquids (NAPLs), in nearsurface heterogeneous granular soils. The

organic liquids can be either lighter than

water or heavier than water. The soil domain

is idealized as consisting of three

interrelated zones: a vadose zone which is

in contact with the atmosphere, a capillary

fringe zone, and a water-table aquifer zone.

A particular problem of interest may include

all three zones or a subset thereof. Granular

soils are stable (non-deforming) and

relatively chemically inert (the soil particles

do not interact with the soil fluids).

Therefore, the soil is idealized as containing

a high percentage of quartz particles and

only a minor percentage of clay particles

and organic matter.

A conceptual illustration of surfacerelease-generated NAPL migration in the

vadose, capillary fringe and aquifer zones

is provided in Figure 1. There are three

fundamental mechanisms for NAPL

migration. First, the NAPL infiltrates into

the soil and migrates both vertically and

laterally under the influence of gravitational

and capillary forces. The distribution of the

NAPL liquid is a function of fluid properties

(density, viscosity, interfacial tension,

wetting potential and variable chemical

composition), soil properties (grain size

distribution, mineral content, moisture

content, porosity, hydraulic conductivity and

spatial heterogeneity), and system forcing

history. If the source is periodic in nature,

then during drying periods, not all the NAPL

will drain from the pore space, leaving
2

behind an immobile residual, held in place

by capillary forces. If the NAPL is more

dense than water, it will migrate through the

capillary fringe and continue its vertical

migration until either the mobility becomes

zero (all the NAPL liquid is at the immobile

residual state) or the NAPL front encounters

an impenetrable geologic horizon.

The second contaminant transport

mechanism is dissolution and consequent

advection in the downward-flowing waterphase, with precipitation providing the water

source in the vadose zone. In the case of

a DNAPL, flowing ground water picks up

dissolved NAPL constituents.

The third transport mechanism is

transport as a vapor NAPL constituent in

the soil gas, where the increased gasphase density induces downward

movement. Partitioning between the gasand water-phase contaminants further

enhances the migratory potential of the

NAPL constituents.

DESCRIPTION OF MODEL

Mathematical Statement of the

Model NAPL

The model NAPL is designed to solve the

following system of governing equations

along with initial and boundary conditions.

Mass Balance Equations

Following the development of Pinder and

Abriola (1986), the mass balance law for

each fluid constituent, an ordered pair (i,a)

representing a species i in a fluid phase is:

∂ε ρ

∂

ε ρυ

ε ρυ

ρ

ρ

ε κρ ρ ρ

α ι

α

α ι

α α

α ι

αα α ι

α

α α

α

ι

α

ι

α α

ι

α

S

t

S

SD S Q

a f +∇⋅ −∇⋅

⋅∇ + = + F

H IK

L

N

M O

Q

P $

(1)

where the five constituents (i,a) relevant to

this simulator are identified as: (w,W ), a

water species in the water phase;(n,W), a

NAPL species in the water phase; (n,N), a

NAPL species in the NAPL phase; (n,G), a

NAPL species in the gas phase; and (g,G),

a gas species in the gas phase.

Other symbols occurring in equation 1 are

used to represent the following:

e is the porosity of the porous medium.

Sa is the saturation of the a-phase.

ri

a is the mass concentration of species in

the a-phase [M/L3].

na is the mass average velocity of a-phase,

a vector [L/T].

Da is the dispersion coefficient for the aphase, a symmetric second-order tensor

[L2/T].

Qa is the point source(+) or (-) a-phase

mass [1/T].

ki

a is the decay coefficient for species in the

a-phase [1/T].

ρ$ι

α is the source or sink of mass for a

species in the a-phase [M/L3T] due to

interphase mass exchange (i.e.,

dissolution, volatilization and

adsorption).

The exchange of mass for each constituent

in equation 1 is defined by:

$

$

$

$

$

ρ

ρ

ρ

ρ

ρ

w

W

n

W

n

WG S

n

N

n

W

n

G

n

G

n

G G

g

G

EEE

E E

E E

n

W nW

n

W

=

=−−

=− +

= +

=

0

0

c h

(2)

where

En

W represents dissolution mass transfer

of NAPL species from the NAPL phase

to the water phase;

11111111111111111111111111

11111111111111111111111111 precipitation

NAPL source area

residual NAPL 

dissolved NAPL

vaporized NAPL

mobile NAPL

groundwater flow

water table capillary fringe zone

soil

moisture

profile

−water head

vadose zone

aquifer zone water content

Figure 1. Definition illustration of NAPL contamination in near-surface soils due to an intermittent release.
3

En

W

G represents volatilization mass transfer

of the NAPL species from the water

phase to the gas phase;

En

G represents volatilization mass transfer

of the NAPL species from the NAPL

phase to the gas phase;

En

W

S represents adsorption mass transfer

of the NAPL species from the water

phase to the soil.

A sixth mass balance equation is required

to describe the NAPL species mass which

is adsorbed onto the soil. This equation is

written as:

∂ ερω

∂ ερκω

1

1 − +− =

S

n

S

S

n

S

n

S S

t

En

W

c h

(3)

where

ρS is the density of the soil [M/L3

] and

ωn

S is the mass fraction of the adsorbed

NAPL on the solid [dimensionless].

The balance of equation 3 is replaced by

the following linear equilibrium relationship:

ω ρ n

S

d n

W = K (4)

where

Kd is a distribution coefficient [L3

/M].

To ensure global mass conservation, the

following definitions and constraints on fluid

volume, density and mass exchange are

employed:

1. The a-phase saturations must sum to

one:

SSS WNG ++=1 (5)

2. The a-phase mass density [M/L3] is the

sum of the species mass

concentrations in the a-phase:

ρ ρα α

ι

α

ι

= =

=

∑ $ , ,, wng , ,

WNG (6)

3. The sum of mass fluxes of all species

into the a-phase must equal the total

mass change in the a-phase:

 

$ $ , ,, , ,

ρ ρα α

ι

α

ι

= =

=

∑wng

WNG (7)

4. The total mass change over all phases

must be zero:

$

, ,

ρα

α=

∑ = WNG

0 (8)

5. The sum of the reacting mass must be

equal to the sum of the produced mass:

κρ α ι

α

ι

ι

α

=

∑ = = wng

WNG , ,

0, ,, (9)

A set of fluid phase mass balance equations

can be generated by summing the balance

equation 1 for each species within the

phase, and by incorporating equations 2,6

and 7. The three resulting fluid-phase

mass balance equations are:

Water-phase:

∂ε ρ

∂ ε ρ

ρ

S

t

S v

QEEE

W

W

W

W W

W W

n

WG S n

W nW

c h +∇⋅

= +−−

(10)

NAPL-phase:

∂ε ρ

∂ ε ρ

ρ

S

t

S v

QEE

N

N

N

N N

N N

n

W

n

G

c h +∇⋅

= −−

(11)

Gas-phase:

∂ε ρ

∂ ε ρ

ρ

S

t

S v

QEE

G

G

G

G G

G G

n

G Gn

W

c h +∇⋅

= ++

(12)

With this development, the physical problem

can be cast into a mathematical

representation consisting of five mass

balance equations. Of the balance

equations written, the following five are

used in the simulator:

1 to 3) The three fluid-phase balance

equations, equations 10,11 and 12. These

equations define the temporal and spatial

distribution and the flow properties of the

water-, NAPL- and gas-phases throughout

the domain.

4 and 5) The two NAPL species balance

equations, equation 1 with (i,a) = (n,W) and

(n,G). These equations define the temporal

and spatial distribution of the NAPL species

as they are transported within and between

their respective phases.

PHASE ADVECTION

Fluid flow is defined by the parameter va

in equations 1 and 10 through 12. The

phase velocity is written in terms of the

multi phase extension of Darcy’s law.

v

k k

S

P z

W, N,G

α α

α

α

α α

ε μ

γ

α

=− ⋅ ∇ − ∇

=

r c h,

 (13)

where

Pα is the a-phase pressure [M/(LT2)],

γ ρ α α = g is the specific weight of the

phase [M/T2

],

g is the acceleration due to gravity [L/T2],

 k is the intrinsic permeability [L2],

considered a scalar herein, and

kra is the relative permeability.

The a-phase relative permeability is a

scaling factor, 0 1 ≤ ≤ krα , which accounts

for the case where the porous medium is

not fully saturated with the a-phase. This

parameter is in general a function of the

aphase saturation. Given the assumption

that the phase wetting-order is constant,

and follows, from most to least, water-NAPLgas, the following functional dependance is

assumed:

k kS

k k SS

k kS

rW rW W

rN rN W Tw

rG rG Tw

=

=

=

b g

b g

b g

,

(14)

where the relationships listed reduce to

their proper two-phase forms when

appropriate. The form of equation 13 for

each fluid phase of interest is detailed here

as:

v kk

S

P z

v kk

S

PP z

v kk

S

PP P z

rW

W

W

W W

N rN

N

N

W

cNW

N

G rG

G

G

W

cNW cGN

G

W

=− ⋅ ∇ − ∇

=− ⋅ ∇ + − ∇

=− ⋅ ∇ + + − ∇

( )

( ) ( )

( ) ( )

ε μ

γ

ε μ

γ

ε μ

γ

(15)
4

MASS TRANSFER

Four types of mass transfer processes

are important in defining the physics of the

fate of NAPLs in near surface granular

soils:

(a) dissolution mass transfer of pure phase

NAPL to the water phase;

(b) evaporation mass transfer of pure

phase NAPL to the gas phase;

(c) evaporation mass transfer of NAPL

species in the water phase to the gas

phase;

(d) adsorption mass transfer of NAPL

species in the water phase to the soil

phase.

Adsorption of NAPL species in the gas

phase directly to the soil phase is neglected.

Liquid-Liquid Mass Transfer

When the organic phase is at an immobile

residual state, saturation is no longer

considered a function of capillary pressure

since capillary pressure becomes

undefined. Consider the NAPL phase

balance equation 11 for the case of an

immobile residual with constant phase

density, constant porosity, and no external

sources or sinks. For these conditions

equation 11 reduces to:

ε ρ

∂

∂

N N

n

W

n

S G

t =− − E E (16)

This equation states that change in NAPL

saturation is due to mass transfer processes.

The dissolution model defining the mass

exchange term, En

W , is assumed to be a

first-order kinetic-type reaction of the form:

E C n

W

n

W

n

W

n

W = − c h ρ ρ (17)

where

Cn

W [1/T] is the rate coefficient which

regulates the rate at which equilibrium is

reached, and

ρn

W [M/L3] is the equilibrium concentration

of the NAPL species in the water phase

(solubility limit).

In the simulator ρn

W is assumed to be a

measurable constant value.

To determine the parametric form of Cn

W ,

the work of Imhoff et al. (1992) is employed.

They conducted column experiments

designed to study dissolution kinetics of

residual trichloroethylene (TCE) in a uniform

sand by flushing the system with clean

water and tracking the dissolution front as a

function of time. Using a lumped parameter

model, they derived the following powerlaw relationship for Cn

W :

C Sv n

W WN

N

W =β ε β β

1

2 3

b g (18)

where

β2 ≈0 5.

 and

β3 ≈10.

 are dimensionless fitting parameters.

The parameter β1

WN [1/T] is the rate

coefficient, and it is fit to available

experimental- or field-scale data.

The volatilization model defining the mass

exchange term En

G is assumed to follow a

similar model as for dissolution, i. e.:

E C n

G

n

G G

n

G = − c h ρ ρ (19)

where

Cn

G [1/T] is the rate coefficient which

regulates the rate at which equilibrium is

reached, and

ρn

G [M/L3

] is the constant equilibrium vapor

concentration of the NAPL species in

the gas phase (vapor solubility limit).

The rate coefficient Cn

G is assumed to have

the form:

C S n

G GN =β ε N

β

1

2

b g (20)

where

β2 is the same as for the dissolution model,

and

β1

GN [1/T] is fit to available data.

Consider now the volatilization of a

dissolved NAPL species in the water phase

to the gas phase. Assuming that the water

phase is at residual saturation in the vadose

zone, and that there are no external sources

or sinks of mass, then equation 10 can be

written as:

ε ∂ρ

∂

S

t W EE E

W

n

WG S n

W nW =− − (21)

where

the exchange term En

W is defined in

equation 15 and En W

G governs the

volatilization mass transfer of a dissolved

NAPL species in the water phase to the gas

phase:

E CH n

W n W

G G

n

W

n

G = − c h ρ ρ (22)

where

H is the dimensionless Henry’s law

coefficient which is defined at equilibrium

conditions as follows:

H n

G

n

W =ρ ρ (23)

and

Cn

W

G [1/T] is the mass transfer rate

coefficient which is assumed to be

defined by the power law:

C S n W

G GW =β ε w

β

1

2

b g (24)

where

the fitting parameter β2 is assumed to be

the same as for the liquid-liquid mass

transfer models, and β1

GW [1/T] is fit to the

available data.

Liquid-Solid Mass Transfer

Finally, mass exchange due to

adsorption, En W

S is assumed to be defined

by a linear equilibrium model:

ω ρ n

S

d n

W =K (25)

where

Kd is the distribution coefficient [L3/M]

defined as a function of the organic

carbon content of the soil and the relative

hydrophobicity of the dissolved NAPL

species:

K fK d oc oc = (26)

where

foc is the mass fraction of organic carbon

and

Koc is the organic carbon partition

coefficient. Combining equations 3 and

25 yields the following definition for En

W

S

:

E K n

W

S b

d

n

W

n

S

n

W = + F

H

G I

K ρ J ∂ρ

∂ κ ρ t (27)

where

ρ ερ b S = −1 is the bulk density of the

soil.
5

MODEL TESTING AND

EXAMPLE PROBLEMS

A series of pseudo-one-dimensional

example problems were presented in order

to evaluate convergence and mass balance,

and to give the user an indication of

appropriate discretization for a given set of

input data. In addition to addressing selfconsistency, four example problems were

designed to simulate specific physical

experiments:

1. a three-phase LNAPL spill and

redistribution experiment (Van Geel

and Sykes, 1995a, 1995b);

2. a three-phase DNAPL spill and

redistribution experiment conducted at

the EPA’s Subsurface Protection and

Remediation Division of the National

Risk Management Research

Laboratory in Ada, Oklahoma;

3. an experimental investigation of the

dissolution of residual DNAPL in a

saturated sand (Imhoff et al, 1992);

4. an experimental investigation of DNAPL

vapor transport in an unsaturated sand

(Lenhard et al., 1995).

The first example simulates water

drainage in a one-dimensional soil column,

1.2 meters long and initially saturated with

water. The boundary conditions are: at the

top, open to the atmosphere, and at the

bottom, specified water head (fine sand =

10 cm, coarse sand = 60 cm). The columns

are then allowed to drain under the influence

of gravity until quasi-static conditions

prevail. Figure 2 presents the results for

the simulation in the first example.

The second example is summarized as

follows:

1. For time = 0 to 100 s, DNAPL is rejected

at a constant volume rate of 0.03 cm3/s

2. for time = 100 s to the end of the

simulation, PG = atmospheric.

Simulation results for two different

discretizations are presented in Figure 3. A

time step of order 2 seconds was required

to obtain a solution for this problem, as

larger time steps caused convergence

problems. Figure 3 shows the total liquid

saturation solution at initial conditions (T =

0) and at T = 200 s (100 s after the DNAPL

source was removed). One can see that

while both discretizations capture the sharp

DNAPL front, the x = 2.5 cm solution exhibits

oscillations behind the front. Figure 3

presents saturation results at time = 5000

s, after the DNAPL has migrated to near

static, residual state. It is apparent that for

these model parameters, a grid spacing of

approximately 1 cm is required. Figures

3.2c and 3.2d present mass balance results

for this simulation. In general the model

performs well with respect to mass balance

except when boundary forcing is changed,

and several time steps are needed to

accommodate the discontinuity imposed.

Mass Transfer Model

Here we investigate the convergence

attributes of the kinetic mass transfer model.

Consider the following one-dimensional

water flow and contaminant transport

problem. The initial conditions are set such

that the domain is saturated, and there is a

zone of residual DNAPL, SN= 0.15, uniformly

distributed from x=25 cm to the bottom.

The boundary conditions are set such that

there is a constant influx of clean water at

the top at a rate of 0.008 cm/s, and an

equivalent efflux of contaminated water at

the bottom. Relevant mass transfer and

transport parameters are:

β2 =0 5. ,

β3 =10. ,

and

aL

W =1cm,

 and

ρn

W =0 001 3 . gm cm .

The results for different values of the

exchange rate coefficient, β1

WN , are

presented in Figures 4a and 4b. As shown

in the Figure, a distinct dissolution front is

created, the shape of which is a function of

the size of β1

WN . High values effectively

approximate the equilibrium partitioning

approximation and produce a sharp front,

while low values produce a broad front.

From a numerics standpoint, the dissolution

front should be resolved over several

elements to minimize oscillations in the

solution which can cause erroneous NAPL

saturations upstream of the source area.

Spatial convergence is illustrated in

Figures 4c and 4d. For a constant rate

coefficient, β1 24 WN = / d , the model

exhibits convergence as the mesh is refined.

BIBLIOGRAPHY

Imhoff, P.T., P.R. Jaffe, and G.F. Pinder,

An experimental study of the dissolution

of trichloroethylene in saturated porous

media, Princeton University Water

Resources Program Report: WR-92-1,

1992.

Lenhard, R.J., M. Oostrom, C.S. Simmons,

and M.D. White, Investigation of densitydependent gas advection of

trichloroethylene: Experiment and a

model validation exercise, J.

Contaminant Hydrology, 19, 47-67,

1995.

Pinder, G.F., and L.M. Abriola, On the

simulation of nonaqueous phase organic

compounds in the subsurface, Water

Resour. Res., 22 (9), 109s-ll9s, 1986.

Van Geel, P.J., and J.F. Sykes, Laboratory

and model simulations of a LNAPL spill

in a variable-saturated sand medium: 1.

Laboratory experiment and image

analysis techniques, Journal of

Contaminant Hydrology, 17(1), 1-26,

1995a.

Van Geel, P.J. and J.F. Sykes, Laboratory

and model simulations of a LNAPL spill

in a variable-saturated sand medium: 2.

Comparison of laboratory and model

results, Journal of Contaminant

Hydrology, 17(1), 27-53, 1995b.
6

0 0.2 0.4 0.6 0.8 1

0

20

40

60

80

100

120

elevation (cm)

water saturation

Experimental S−P relationship

n = 6.49

a = 0.0203/cm

Swr = 0.17

K = 8.64m/d

(a)

0 0.2 0.4 0.6 0.8 1

0

20

40

60

80

100

120

dx = 5 cm

dx = 10 cm

dx = 20 cm

(b)

WT

water saturation

Computed steady−state moisture profile

0 0.2 0.4 0.6 0.8 1

50

60

70

80

90

100

110

120

elevation (cm)

water saturation

Experimental S−P relationship

n = 11.434

a = 0.0849/cm

Swr = 0.08

K = 170 m/d

(c)

0 0.2 0.4 0.6 0.8 1

50

60

70

80

90

100

110

120

dx = 2 cm

dx = 5 cm

(d)

WT

water saturation

Computed steady−state moisture profile

Figure 2. Analysis of appropriate grid spacing to compute capillary rise for different soil-types. Parts (a) and (b) are for a relatively fine

sand, and parts (c) and (d) are for a relatively coarse sand.
7

0 0.2 0.4 0.6 0.8 1

50

60

70

80

90

100

110

120

total liquid saturation

depth (cm)

T = 0, dx = 1cm

T = 0, dx = 2.5 cm

T = 200, dx = 2.5 cm

T = 200, dx = 1cm

WT

(a)

0 0.2 0.4 0.6 0.8 1

40

60

80

100

120

saturation

depth (cm)

WT

Sw, dx = 1cm

Sw, dx = 2.5 cm

Sn, dx = 2.5 cm

Sn, dx = 1cm

time = 5000 (b)

0 100 200 300 400

0

0.5

1

1.5

2

(c)

EARLY TIME

dx= 1cm

dt= 2 s

T = 100

time

mass balance ratio

1000 2000 3000 4000 5000

0

0.5

1

1.5

2

(d)

LATE TIME

dx= 1cm

dt= 2 s

time

mass balance ratio

Figure 3. Results of a one-dimensional, three-phase, DNAPL injection and redistribution simulation, highlighting spatial convergence

and mass balance.
8

Figure 4. Computational analysis of the dissolution model. Parts (a) and (b) illustrate the effect that the rate constant has on the solution.

As the dissolution front sharpens, oscillation appears indicating that a finer grid spacing is required. Parts (c) and (d) illustrate

spatial convergence for ex=24/d. For the parameters chosen a grid spacing of approximately 5 cm is appropriate.

0 0.2 0.4 0.6 0.8 1

0

20

40

60

80

100

120

depth (cm)

NAPL saturation / 0.15

T = 0

T = 3.5 days

dx = 2.5cm

ex = 10/d

ex = 24/d

ex = 100/d

(a)

0 0.2 0.4 0.6 0.8 1

0

20

40

60

80

100

120

depth (cm)

Dissolved concentration / 0.001

T = 0

T = 3.5 days

dx = 2.5cm

ex = 10/d

ex = 24/d

ex = 100/d

(b)

0 0.2 0.4 0.6 0.8 1

0

20

40

60

80

100

120

depth (cm)

NAPL saturation / 0.15

T = 0

T = 3.5 days

ex = 24/d

dx = 12.5cm

dx = 5cm

dx = 2.5cm

(c)

0 0.2 0.4 0.6 0.8 1

0

20

40

60

80

100

120

depth (cm)

Dissolved concentration / 0.001

T = 0

T = 3.5 days

ex = 24/d

dx = 12.5cm

dx = 5cm

dx = 2.5cm

(d)
9

Joseph Guarnaccia and George Pinder are with Research Center for Groundwater

Remediation Design, University of Vermont, Burlington, VT 05401, and Mikhail

Fishman is with Robert S. Kerr Environmental Research Center, Ada, OK 74820.

Thomas Short is the EPA Project Officer (see below).

The complete report, entitled "NAPL: Simulator Documentation," ( Order No. PB95-X;

Cost: $X.00, subject to change) will be available only from:

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22161

Telephone: 703-487-4650

The EPA Project Officer can be contacted at:

U.S. Environmental Protection Agency

National Risk Management Research Laboratory

Subsurface Protection and Remediation Division

P.O. Box 1198

Ada, OK 74820-1198

BULK RATE

POSTAGE & FEES PAID

EPA

PERMIT No. G-35

United States

Environmental Protection Agency

Center for Environmental Research Information

Cincinnati, OH 45268

Official Business

Penalty for Private Use

$300

EPA/600/SR-97/102