xvEPA
United States
Environmental Protection
Agency
Office of Solid Waste
and Emergency Response
Washington DC 20460
July 1986
Solid Watte
Criteria for Identifying
Areas of Vulnerable
Hydrogeology Under the
Resource Conservation
and Recovery Act
Appendix C
Technical Methods for
Calculating
Time of Travel in the
Unsaturated Zone
Interim Final
-------
OSWER POLICY DIRECTIVE NO.
9472-00-2A
GUIDANCE CRITERIA FOR IDENTIFYING
AREAS OF VULNERABLE HYDROGEOLOGY
APPENDIX C
TECHNICAL GUIDANCE MANUAL FOR CALCULATING
TIME OF TRAVEL (TOT) IN THE UNSATURATED ZONE
Office of Solid Waste
Waste Management Division
U.S. Environmental Protection Agency
401 M Street, S.W.
Washington, D.C. 20460
Prepared by
Batelle Project Management Division
Office of Hazardous Waste Management
601 Williams Boulevard
Richland, WA 99352
Under Subcontract to
GCA CORPORATION
GCA/TECHNOLOGY DIVISION
Bedford, Massachusetts 01730
July 1986
-------
OSWER POLICY DIRECTIVE NO.
947200-2A *
DISCLAIMER
This Final Report was furnished to the Environmental Protection Agency
by the GCA Corporation, GCA/Technology Division, Bedford, Massachusetts 01730,
in fulfillment of Contract No. 68-01-6871, Work Assignment No. 20. The
opinions, findings, and conclusions expressed are those of the authors and not
necessarily those of the Environmental Protection Agency or the cooperating
agencies. Mention of company or product names is not to be considered as an
endorsement by the Environmental Protection Agency.
ii
-------
USWtK KULiuY Ul
9472- 00" 2A
EXECUTIVE SUMNBMQT
This appendix to the Guidance Criteria for Identifying Areas of Vulnerable
Hydrogeology describes methods for calculating ground-water time of travel (TOT)
in the unsaturated zone. The methods described in this appendix are intended for
use by hazardous waste facility permit applicants and writers in evaluating the
vulnerability of ground water to contamination.
The appendix presents a review of the general theory of ground-water flow in
the unsaturated zone and describes the processes that control flow. Equations are
presented which describe these processes and illustrate the relationships between
important parameters.
Two general approaches are presented for calculating unsaturated zone TOT. The
first approach involves the use of analytical solutions. These solutions are simp-
lified approaches, appropriate for simple systems, and allow the analyst to directly
y
solve for TOT. Two solutions are described, both for determination of steady state
TOT. The first solution assumes a constant moisture content in the soil profile and
is appropriate for conditions dominated by gravity drainage. The second solution
allows for variable moisture contents and is appropriate for conditions where
factors other than gravity drainage (e.g., capillary forces) are important. The use
of these solutions is described and data requirements and sources of data are
identified.
The second approach involves the use of unsaturated flow models. Two general
classes of models are described, numerical models and water balance models. The
relative complexity of each type of model is described, as are data requirements,
output, and limitations. Methods are presented for determining TOT from model
output for those models where TOT is not directly calculated by the model.
iii
-------
me approaches to calculating TOT are summarized and a decision tree is
presented to aid in selection of the most appropriate approach for specific
applications. Three case histories are presented using data from actual
hazardous waste facilities to illustrate calculation of TOT.
IV
-------
CONTENTS
Executive Summary . . . . . . . iii
Figures . . . . . . vi
Tables ............. viii
Acknowledgements ........... ix
1. Introduction 1.1
2. Technical Background on Unsaturated Flow . . . .2.1
3. Technical Approaches for Determining TOT .... 3.1
Analytical Solution of TOT 3.1
Modeling Solutions of TOT 3.3
Selection of the Appropriate Method to
Determine TOT 3.11
4. Analytical Solutions of TOT 4.1
Data Requirements and Sources . . . . .. .4.1
Description of Analytical Solutions . . . . 4.4
5. Unsaturated^'low Models 5.1
Introduction ......... 5.1
Example Numerical Code - UNSAT1D 5.7
Example Water Balance Code - HELP . . . . .5.14
6. Examples of TOT Determination ....... 6.1
Example 1 .......... 6.1
Example 2 .......... 6.5
Example 3 .......... 6.11
7. References . . . . . . . . . .7.1
Glossary ............. G.I
-------
FIGURES
Number Page
2.0-1 Graph of Moisture Content Versus Suction Head . . .2.4
2.0-2 Graph of Hydraulic Conductivity Versus
Moisture Content 2.5
3.2-1 Cross-Section View of the Node and Element
Discretization of an Unsaturated Flow System
Beneath a Waste Disposal. Facility 3.5
3.2-2 Cross-Section View of the Layer Representation of
an Unsaturated Flow system for Use in a Typical
Water Balance Model 3.6
3.2-3 Example Model Results Showing Migration
of Wetting Front with Time . . . . . . . 3.9
3.3-1 Summary of Selection of Approach for TOT Determination . . 3.12
4.2-1 Schematic of Example Single-Layered System .... 4.7
4.2-2 Schematic of Example Multi-Layered System . 4.10
v
4.2.2-1 Discretization Between Grid Points 4.13
5.2-la Moisture Profile for a Three Layer
Unsaturated Flow System 5.12
5.2-lb Cumulative Water Flux Versus Time Past Several
Elevations in an Unsaturated Flow System . . . .5.12
6.2-1 Plot of Suction Head Versus Moisture
Content for Case Study D . . . . . . .6.7
6.2.2-1 Plot of Hydraulic Conductivity Versus Suction
Head for Case Study D 6.9
6.2.2-2 Steady State Moisture Profile from Iterative
Analytical Solution for Case Study D . . . . . 6.10
6.3.1-1 Soil Water Characteristic Curve for Soil
at Proposed Hanford Hazardous Waste Site . . .6.13
vi
-------
FIGURES (Cont'd.) 7^'OQ
Number . Page
6.3.1-2 Hydraulic Conductivity Versus Water Content
Curve for Soil at Proposed Hanford
Hazardous Waste Site 6.14
6.3.1-3 Simulated Pressure Head Versus Depth with
Time at Proposed Hanford Hazardous Waste Site . . . 6.15
6.3.2-1 Simulated Advance of Wetting Front at Proposed
Hanford Hazardous Waste Site .6.16
6.3.2-2 Simulated Rate of Leachate Discharge at Proposed
Hanford Hazardous Waste Site 6.18
6.3.2-3 Simulated Cumulative Leachate Discharge at
Proposed Hanford Hazardous Waste Site ..... 6.19
vn
-------
TABLES
Number Page
3.0-1 Relative Characteristics of Analytical
Solutions and Unsaturated Flow Modeling . . . .3.1
4.1-1 Representative Values for Saturated _
Hydraulic Conductivity . . . ..... 4.3
4.1-2 Representative Values for Porosity ...... 4.4
5.1-1 Summary of Unsaturated Zone Codes ...... 5.2
5.1-2 Sumnary of Unsaturated Code Documentation
and Availability . . . . . . . 5.6
5.2-1 Important Characteristics and Capabilities
of UNSAT1D .......... 5.8
5.2-2 Summary of UNSAT1D Input Data and Sources . . . . 5.10
5.3-1 Summary of Characteristics and
Capabilities of HELP ........ 5.15
5.3.1-1 Summary of Irvut Data and Sources for HELP .... 5.17
5.3.1-2 Summary of Output Data from HELP ...... 5.19
6.1.1-1 Typical Values for Slope of Soil Moisture
Retention Curve ......... 6.4
viii
-------
ACKNOWLEDGEMENTS
This appendix to the Guidance Criteria for Identifying Areas of Vulnerable
Hydrogeology was prepared by Batelle's Project Management Division, Office of
Hazardous Waste Management under subcontract to GCA Corporation/ GCA/Technology
Division.
The appendix was written by Frederick W. Bond, James M. Doesburg, C. Joseph
English, and Deborah J. Stallings. The authors wish to thank Alfred Leonard,
Charles Young, and Bob Farrell of GCA, and Glen Galen of EPA's Office of Solid
Waste, for their assistance in the preparation of this appendix. Special thanks
are also extended to Peggy Monter, Nancy Painter, Beth Eddy, and Dick Parkhurst, of
'Battelle, for their typing and graphics assistance.
IX
-------
SECTION 1.0
INTRODUCTION
Subtitle C of the Resource Conservation and Recovery Act of 1976 (RCRA)
addresses the management and disposal of hazardous waste. Regulations
developed by the Environmental Protection Agency (EPA) from RCRA legislation
(Regulations for Owners and Operators of Permitted Hazardous Waste
Facilities - 40 CFR 264 and Regulations for Federally Administered Hazardous
Waste Facilities - 40 CFR 270) require owners and operators of hazardous
waste land treatment, storage, and disposal (TSD) facilities to provide
information concerning the design, construction, operation,' and maintenance
of these facilities. This information is provided in the form of a RCRA Part
B permit application. Permit writers in the EPA Regions and authorized
states responsible for writing permits must review this information to
determine if the facility will meet the environmental protection goals
established in the RCRA regulations.
v
A major emphasis of the above environmental protection goals is the
protection of ground-water resources that may be vulnerable to contamination
originating at TSD facilities located in certain hydrogeologic settings.
Therefore, much of the permit application and review process addresses the
adequacy of facility design, construction, operation, maintenance, and
location with respect to ground-water protection. A considerable amount of
guidance has been developed to aid permit writers in evaluating potential
threats to ground water posed by TSD facilities. This appendix presents
methods available for predicting ground-water time of travel (TOT) in the
unsaturated zone.
The intent of RCRA guidance and regulations is to ensure that facilities
are designed, operated, and located such that there will be negligible
migration of contamination beyond the barriers of the facility. Assuming no
release, the velocity of contaminant migration to the ground water beneath
1.1
-------
the site is of little or no importance. The regulations, however, also
recognize the possiblity of contaminant release due to failure of facility
barriers. The contaminant time of travel to the ground water and the
thickness of the unsaturated zone then become important issues in assessing
the potential consequences of site failure.
Consideration of unsaturated TOT may be important for several reasons.
With the exception of certain hydrogeologic settings, all components of a
facility permitted by RCRA will be located above the water table. Therefore,
any release from a facility will migrate through the unsaturated zone before
reaching the water table. The TOT through the unsatuated zone is important
in determining how rapidly, and to what extent, ground-water resources will
be impacted by contaminant release. Consideration of unsaturated TOT may be
necessary to determine appropriate monitoring strategies for detecting
failures as well as for developing appropriate corrective actions.
In some locations, particularly those characterized by arid climates and
deep unsaturated zones, flow through the unsaturated zone may be minimal. In
such cases, the unsaturated zone may effectively form a buffer zone to delay
contaminants when migrating to the water table. Recognition of such
conditions is important in assessing the adequacy of a facility design.
Calculation of TOT in the unsaturated zone is not a trivial task.
»*
Simple graphical methods useful for estimating TOT in the saturated zone
(e.g., construction of flow nets) are not applicable due to nonlinearities
associated with unsaturated flow. Solution of unsaturated TOT must instead
rely on analytical and numerical methods. Available techniques vary
significantly with respect to mathematical complexity and data requirements.
Selection of an appropriate technique must consider the objectives of the
application and the availability of time, resources, and data.
The purpose of this appendix is to acquaint permit writers and
applicants with the techniques available for determination of TOT in the
unsaturated zone. The appendix is divided into six sections. This section
provides an introduction to the need for determination of TOT. Section 2.0
presents technical background material regarding unsaturated flow. A general
discussion of the two technical approaches for TOT determination, analytical
methods, and unsaturated flow models is provided in Section 3.0.
1.2
-------
Sections 4.0 and 5.0 discuss these approaches in more detail, and example
determinations of TOT are presented in Section 6.0.
1.3
-------
SECTION 2
TECHNICAL BACKGROUND ON UNSATURATED FLOW
The unsaturated zone is the transition region between the atmosphere and-
the saturated ground-water system. Passage of water through this zone is
very dynamic and depends on detailed variations in the hydraulic properties
of the water in the soil.
Rainfall, irrigation, and ponded water are the primary sources of water
to the unsaturated zone. Redistribution or downward movement of this water
through the soil occurs under the influence of gravity as long as there is a
sufficient quantity present to overcome the restraining forces of capillary
hydraulic potential.
Water is removed from the surface of the unsaturated zone by the
processes of evaporation and/or transpiration. The rates of both processes
depend directly on available solar energy and surface winds. Water also
moves out the bottom «/ the unsaturated zone as drainage if the soil-water
holding capacity is exceeded. Drained water may possibly enter the water
table depending on its depth.
Water can also move within and be stored in the unsaturated zone. Water
storage is characterized by a water content distribution. Water moves
through and within a soil via two physical mechanisms: capillary Darcian
flow (liquid phase) and vapor diffusion. Darcian flow is described by
hydraulic conductivity and matric potential gradients, both of which are
highly dependent on moisture content. Vapor diffusion controls actual
surface evaporation and results from thermal gradients.
A soil is saturated when all void space (space not occupied by solid
particles) is filled with water. An unsaturated soil contains air-filled
void space as well as water. A measure of the quantity of water contained by
a soil is called the water content which can be defined either on a
volumetric basis (volume of water/total volume of soil, water, and voids) or
a mass basis (mass of water/mass of soil solids).
2.1
-------
Movement of water in the unsaturated zone is always directed from areas
of higher to those of lower potential energy (assuming isothermal
conditions). Total soil water potential ($) is expressed as (Feedes et al..
1978):
where
^ = the pneumatic potential arising from changes in external pressure
^s = osmotic potential arising from the attraction forces of water to a
higher solute concentration
'J' = the matric potential arising from capillary and adsorptive forces
of the soil matrix
^ = the gravitational potential expressing the potential energy of
changes in relative elevation changes.
The negative of the gradient of total potential is the force causing
water movement in a soil.
The gravitation potential is an important component of the driving force
of water downward through the unsaturated zone below a TSD facility. The
gravitational potential of soil water at each point is determined by the
elevation of the poi'v- relative to some arbitrary reference level. The
matric (or capillary) potential is a negative pressure potential resulting
from the adsorptive forces of the soil matrix. The matric potential can be
an important factor, particularly in dry soils. The influence of the
pneumatic and osmotic potential is almost always quite small and, therefore,
they can be disregarded.
The relation between matric potential and soil wetness (water content)
is not generally a unique one due to a phenomena known as hysteresis.
Hysteresis is the phenomenon where the water content of a soil with a given
matric potential can be different depending on whether the soil is wetting
(sorption) or drying (desorption) . The equilibrium water content at a given
suction is greater in desorption than in sorption. The hysteresis effect can
be attributed to several causes which include: 1) the "geometric
nonuniformity of individual pores; 2) entrapped air; and 3) swelling or
2.2
-------
shrinking phenomena. Typically, the hysteresis effect is small and is,
therefore, disregarded in the determination of TOT (Hillel, 1971).
Movement of water through a porous medium is proportional to both the
hydraulic conductivity of the medium and the hydraulic potential gradient
across the medium. The hydraulic conductivity, K, represents the ability of
a soil to transmit water from locations of high hydraulic potential to
locations of low hydraulic potential. In cases of unsaturated flows,
hydraulic conductivity is a function of moisture content (e), such that K can
be represented as K(e). The functional relationship of K(e) will vary from
soil to soil. Typically, K(e) will decrease rapidly by several orders of
magnitude from its maximum saturated value as water content decreases.
As stated above, the moisture content can vary within the unsaturated
flow system. As the moisture content changes, the matric potential (suction
head or negative pressure head) and the hydraulic conductivity also change.
In order to simulate water movement in the unsaturated zone, these
relationships between moisture content and matric potential, and moisture
content (or matric potential) and hydraulic conductivity (soil characteristic
curves) must be known. The only exception to this is for cases where the
moisture content, and therefore the matric head and hydraulic conductivity,
remain constant throuc^out the unsaturated soil profile (unit gradient case).
Example graphs of these two relationships are shown in Figures 2.0-1 and
2.0-2.
Most often the moisture content, matric potential, and hydraulic
conductivity relationships are not available for the soils of interest at a
particular site. If these data are not available, the best means of
obtaining site-specific values is through laboratory measurements. Field
measurements can be made, but unlike those for saturated soils, they are
typically not as accurate and reliable as laboratory tests.
The moisture content versus pressure head relationship can be measured
relatively easily, whereas methods for direct determination of hydraulic
conductivity (K) as a function of moisture content (9) or matric potential
(*m) over the unsaturated range of interest are experimentally difficult. As
a result, the K versus e or «/m relationship is often calculated using
analytical methods such as those presented by Mualem (1976a), Burdine (1953),
2.3
-------
0.4000
0.3000
c
OJ
4-1
c
o
CJ
0.2000
0.1000
0.0000
io5 io6 io7
Suction (cm)
Figure 2.0-1. Graph of Moisture Content Versus Suction
Head
2.4
-------
C
'i
u
>.
4-»
>
4->
U
3
T>
C
O
(J
10°
io:2
io:3
io:4
io:5
io:6
io:7
io:8
io:10
io:11
io:12
io:13
io:14
io:15
io:17-
0.0000
0.1000
0.2000
0.3000
0.4000
Moisture Content
Figure 2.0-2. Graph of Hydraulic Conductivity Versus Moisture Content
-------
and Millington and Quirk (1961). The details of these analytical methods are
discussed in Mualem (1976a).
If site-specific data are not available and cannot be measured in the
field or laboratory, a last resort is to obtain values for representative
soils from the literature. One excellent reference is Maul em (19765) where
example soil characteristic curves for 45 soils for which actual conductivity
measurements as well as moisture content as a function of matric potential
were made.
2.6
-------
SECTION 3.0
TECHNICAL APPROACHES FOR DETERMINING TOT
There are two basic technical approaches for determining TOT in the
unsaturated zone, analytical solutions of TOT, and unsaturated flow modeling.
Both approaches are based on the same fundamental equations, but differ in
the number of simplifying assumptions made in order to solve these equations.
As a result of the simplifying assumptions, the approaches differ
significantly in the time necessary to obtain a solution, in computational
difficulty, and in data requirements. Relative characteristics of the two
approaches are summarized in Table 3.0-1.
Table 3.0-1. Relative Characteristics of Analytical
Solutions and Unsaturated Flow Modeling
Computational Time
- ^
Data Requirements
Complexity of Solution
Time Dependency
Analytical
Solutions
Short
Low to Medium
Simple
Steady State
Unsaturated
Flow Modeling
Medium to Long
Medium to High
Complex
Steady State or Transient
3.1. ANALYTICAL SOLUTION OF TOT
The analytical solution of unsaturated flow TOT is based upon Darcy's
equation for one-dimensional flow
3.1
-------
where
q = flux in the vertical direction
^ = matric potential (suction head or negative pressure head)
K(^) = hydraulic conductivity as a function of matric potential
3i|/3z = hydraulic gradient in the vertical direction
The above relation is identical.to the relation for saturated flow except
that the hydraulic conductivity is not constant.
In unsaturated flow, both hydraulic conductivity and moisture content
are nonlinear functions of pressure head. Pressure head, hydraulic
conductivity, and moisture content need not be constant throughout the soil
column. When these variables are not constant, a direct analytical solution
of Darcy's equation is not possible for unsaturated flow. In order to obtain
an approximate solution, simplifying assumptions must be made. Common
assumptions are:
one-dimensional flow in the vertical direction;
water flow is steady state;
water table conditions exist at the lower boundary;
t the upper boundary condition is constant flux;
soil characteristics (moisture content versus matric potential and
hydraulic conductivity versus matric potential) are constant with depth;
and
y
the hydraulic gradient is vertically down and equals unity (drainage is
3(k_
due strictly to gravity, or - = o).
3^
For nonhomogeneous soils, the constant property assumption can be
approximated by dividing the soil profile into a series of layers, each layer
comprised of soils having approximately the same characteristics.
The unit gradient assumption greatly simplifies the analysis. This
assumption means that the matric potential and, therefore, moisture content
and hydraulic conductivity are constant with depth. Using this assumption,
it is possible to directly solve for moisture content in terms of the flux
through the system and saturated soil properties. Knowing the moisture
content and flux it is possible to calculate the pore water velocity and TOT.
The unit gradient assumption is generally valid if gravitational forces
dominate other forces (e.g., capillary forces).
3.2
-------
If the unit gradient assumption is not made, the analytical solution to
unsaturated flow becomes more complex. In this case, it is necessary to
employ an iterative solution for pressure head and moisture content. This
iterative solution is time-consuming, but can be simplified through the use
of a computer.
The above solutions for TOT are one-dimensional solutions. When
applying these solutions to specific sites, it is important to consider the
horizontal variability of soil characteristics. If soil characteristics vary
spatially, the solution should be applied to the soil profile having the
highest hydraulic conductivity. The solution will then yield the highest
velocity and shortest TOT (e.g., worst case) for the unsaturated flow system.
In summary, analytical solutions provide a means of quickly estimating
TOT. Several assumptions are required to perform the solutions. . These
assumptions and, therefore, the methods themselves, are not appropriate for
all applications. A detailed discussion of two analytical approaches of
calculating TOT are presented in Section 4.0.
3.2 MODELING SOLUTIONS OF TOT
3.2.1 Description of Model Types
Unsaturated flow models provide another means for determining TOT in the
unsaturated zone. ThVs section discusses two general types of unsaturated
flow models and how they can be applied to solve for TOT.
Two general types of unsaturated flow models are available. These are
numerical models and water balance models. Numerical models solve
differential equations describing water movement within the unsaturated zone.
These equations are derived by combining mathematical statements for the
conservation of mass and energy with equations which have been developed to
relate these statements to measurable quantities such as pressure,
temperature, and moisture content. Numerical models are able to account for
the nonlinearities in soil properties and the variations of properties in
space. Water balance models, on the other hand, simulate unsaturated flow
systems such that the flow into the system is equal to the flow out of the
system plus or minus storage within the system for a specific area and for a
specific period of time. These models generally require the direct or
3.3
-------
indirect measurement of soil moisture and other properties which affect water
movement within the unsaturated zone. Simplifying assumptions similar to
those used with analytical solutions are often used with water balance
models.
The two types of models represent very distinct levels of complexity in
their methods of solution and their required input data. Numerical models
typically require more data than water balance models. Numerical solutions
use either a finite difference, finite element, or integrated finite
difference technique, all of which require that the unsaturated flow system
be represented as a series of nodes and elements (Figure 3.2-1). A complex
system may be represented in the models with several hundred nodes, and the
solution complex numerical problems may require several hours of computer
time.
Water balance models, on the other hand, represent the unsaturated flow
system as a series of layers of geologic materials (Figure 3.2-2). In
typical applications, such models seldom have greater than ten layers, and
the model input data defines the properties of each layer. Water balance
models use a "book-keeping" approach to keep track of water entering and
exiting the system, as well as water entering and exiting each layer within
the system. Water balance models can be solved quite rapidly by computers
with solution times usu? ly on the order of minutes.
3.2.2 Determination of TOT from Model Results
The equations used to develop the solutions used in the unsaturated flow
codes do not have velocity as a variable, nor is velocity a primary output of
the models. Therefore, additional analysis of model results is necessary to
derive TOT. The purpose of the following discussion is to present methods of
determining TOT from unsaturated model results.
Four methods of determining TOT from model results are presented below.
One method deals with determining the travel time for a particle along a
travel path in the unsaturated zone. Another method deals with determining
the time required for an instantaneous loading to migrate through the
unsaturated zone to the water table. The other two methods use steady state
solutions for moisture conditions or contaminant concentrations.
3.4
-------
Ground Surface
00
en
Waste
Disposal Facility
Node
Element
Water Table
Figure 3.2-1.
Cross-Section View of the Node and Element Discretization of an
Unsaturated Flow System Beneath a Waste Disposal Facility
-------
Ground Surface
CA>
Water Table
Waste
Clay
"V"
Silty Sand
Sand
Figure 3.2-2.
Cross-Section View of the Layer Representation of an Unsaturated
Flow System for Use in a Typical Water Balance Model
-------
Determination of Particle Travel Time--
This approach to TOT is appropriate for determining TOT associated with
any type of loading to the system (steady state or transient). The approach
is based upon Darcy's Equation for unsaturated flow (Equation 3-1) and
requires that the code be capable of determining pressure head at each node
or layer boundary for each time step. As described by Darcy's Equation, the
velocity between two nodes (or through a layer) will be equal to the
hydraulic conductivity between the two nodes times the hyraulic gradient
between the two nodes. These quantities can be directly determined from
matric potential. The average velocity for the element can be converted to
pore water velocity for a time step, by dividing by moisture content.
Knowing the pore water velocity it is possible to determine how far a
particle travels during that time step.
This approach to TOT determination involves tracking the position of a
particle present at the surface (or water source) at the start of the first
time step. Using the average velocity of the particle's position for the
first time step, the displacement of the particle for that time step is
calculated. After the second time step, the velocity at the particle's new
position is calculated and the particle displaced a distance for 'that time
step. This process is repeated until the particle reaches the water table.
The time required to reach the water table is the unsaturated TOT.
It should also be noted that some unsaturated flow codes have ancillary
programs that will extract appropriate data and perform the above analysis.
If desired, it is possible to modify any numerical code to perform the
analysis.
Solution of Steady-State Travel Time--
An approach similar to the previous approach can be used to determine
steady state travel times through the unsaturated zone. A model is used to
solve for steady state values of hydraulic potential, moisture content, and
hydrauilc conductivity at each node. These nodal values are used to
determine the steady state pure water velocity between each pair of nodes
(i.e., across each element in the model). Elemental velocities are then used
to determine travel times for each element. The sum of all the elemental
travel times along a flow path is the unsaturated TOT.
3.7
-------
The above approach would be appropriate for determining TOT in cases
where steady state conditions will be encountered (e.g., constant surface
flux conditions). A limitation of the approach is that the TOT is valid only
for steady state conditions and does not yield the TOT for the period leading
up to steady state conditions. For example, if a column of soil at one set
of steady state moisture conditions suddenly receives a new or additional
constant flux of water at the surface, the above approach can be used to
determine the TOT after new steady state conditions have been reached. The
method will not determine the TOT associated with water which passed through
the soil column during the time when new steady state conditions were being
established. Therefore, the method will not identify the TOT associated with
the first particle of water to reach the water table.
TOT Associated With Instantaneous Loadings--
The advantage of the above approach is the ability to determine TOT for
any type of loading. The disadvantage is that the solution algorithm may not
be included in the code being used. The following approach is appropriate
for virtually any code, but is limited to cases involving large transient
fluxes at the surface. Such fluxes represent extreme events; including
natural events such as extreme precipitation, or artificial events such as
sudden failure of a landfill or impoundment liner.
Introduction of a^iarge amount of infiltration to the unsaturated zone
will cause the propagation of a wetting front downward through the soil
column. Transient unsaturated flow models may be used to track the progress
of this wetting front through the unsaturated zone. The time required for
the wetting front to reach the saturated zone may be taken as the TOT. This
approach is illustrated in Figure 3.2-3. This figure shows the variation of
moisture content with depth and time. The wetting front appears as an area
where there is a rapid change in moisture content with depth. In the example
shown in Figure 3.2-3, the time for the wetting front to reach a water table
located at 40 m depth would be approximately 90 hours. If graphical model
outputs are not available, the moisture content of each node (layer) should
be examined at the end of each time step. A large increase in moisture
content over a time step signifies the passing of the wetting front.
3.8
-------
10
20
Q.
O)
O
30
1 = 90 min 1.5 hr)
1 = 900 min (15 hr)
I 1 = 1,800 min (30 hr)
I
1 = 3,600 min (60 hr)
1 - 5,400 min (90 hr)
40
0 .05 .10 .15 .20 .25 .30
Moisture Content, 0
Figure 3.2-3. Example Model Results Showing Migration
of Wetting Front with Time
3.9
.35
-------
An important consideration in the above approach is determination of the
input source. The input must be large enough to cause the wetting front to
migrate entirely through the soil column to the saturated zone. For example
in areas with deep unsaturated zones having low moisture contents, even very
large inputs may never reach the saturated zone regardless of how long the
model is run. Unfortunately, there is no general guidance for determining
how large an input is required.
It is very important to note that use of this approach does not imply
that variations of moisture content with time are necessary for unsaturated
flow to occur. This approach is applicable only for transient flow
conditions and only when the water input is large enough to cause a
significant perturbation to existing conditions.
Determination of Contaminant Travel Time
Some unsaturated codes have the capability of modeling the transport of
contaminants through the unsaturated zone. A known concentration of
contaminants can be input to the top of the soil column and the code used to
determine the concentration of contaminants leaving the bottom of the
unsaturated soil column. If the simulation is run for a long period of time,
steady state conditions will be reached where the concentration of
contaminants leaving t^e soil column is equal to the concentration entering
the soil column. IfHhere were no dispersion of contaminants (e.g., plug
flow), the time required to reach steady state conditions would be equal to
the time of travel for the contaminant. Because of dispersion, however, the
average contaminant TOT will be somewhat less than the time to reach steady-
state conditions. An average contaminant TOT can be estimated as some
fraction of the time to reach steady state conditions (e.g., the time to
reach output concentrations equal to one-half the steady state value). The
ground-water TOT can be related to the contaminant TOT by use of a
retardation factor. The retardation factor is equal to the ground-water
velocity divided by the contaminant velocity. Dividing the average
contaminant TOT by the retardation factor will yield the unsaturated ground-
water TOT.
3.10
-------
3.3 SELECTION OF THE APPROPRIATE METHOD TO DETERMINE TOT
A decision tree for selection of the appropriate method for determining
TOT is shown in Figure 3.3-1. The decisions to be made in this figure
consider both the characteristics of the site and the availability of the
data. The representation of the site should be as realistic as possible
depending on the availability of the data. For example, if soil properties
at a site exhibit complex spatial variability but data are not available to
describe this variability, the assumption of simple variability should be
made. Under these circumstances, one has to realize that the assumption(s)
made can significantly impact the results.
The rule of thumb for selecting the appropriate procedures is to choose
the simplest one which can be applied to your specific problem. If the
unsaturated flow system can be represented as a one-dimensional steady state
problem with a single material type, and the flow is controlled by gravity
drainage (i.e., unit gradient), the simplest analytical approach can be used
to obtain a direct solution. If, however, the unit gradient assumption does
not apply but the other assumptions are applicable, an analytical method
using an iterative solution scheme can be applied. All transient problems,
and problems in more than one dimension require the use of simulation
models. Water balvice models are appropriate for simple quasi-two-
dimensional problems where flow is controlled by gravity drainage. Numerical
models should be used for all higher dimension problems having complex
geometry and boundary conditions.
The direct analytical solution technique is quick and easy to apply;
however, the assumptions of the method limit its applicability. The
iterative analytical solution offers more flexibility, but limitations on the
dimensionality and complexity of the problem restrict its applicability. Due
to the limitations of the analytical approaches, simulation modeling can be
the only means of obtaining reliable results for complex problems.
Application of modeling is typically limited by the availability of data
and the availability of time. The more complex a model, the greater its data
requirements (both in time and space). Acquisition of these data often
requires that field and laboratory programs accompany model development.
3.11
-------
Unsaturated
Flow Problem
One-Dimensional
Problem
Yes
Steady State
Yes
Water Table
Boundary
Yes
Spatial
Variability
Iterative
Analytical
Solution
No
Simple
Unit Gradient
Yes
Direct
Analytical
Solution
Spatial
Variability
Simple
Gravity
Drainage
Controlled
No
Water Balance
Model
Figure 3.3-1. Summary of Selection of Approach for TOT Determination
3.12
-------
While complex models may be set up and run with limited data using
simplifying assumptions, such an approach does not take full advantage of the
capabilities of these models and is, therefore, of limited value.
Development and application of models can also be very time consuming
(weeks to months). Calibration and verification of models may take
considerably more time. The time required is usually proportional to the
complexity of the model.
The above limitations relate to modeling in general. Specific
limitations related to the use of models to support permit writing activities
are:
1) Time requirements, especially for development and calibration of
numerical models, may exceed the time available for preparation and
review of Part B applications.
2) Data supplied with Part B applications may be inadequate to develop and
calibrated models, and will almost certainly be inadequate to verify
models.
The above limitations do not imply that the use of unsaturated flow
modeling is inappropriate to support permit writing. Rather, these
limitations are presented to aid the reader in determining the suitability of
modeling for Part B application.
3.13
-------
SECTION 4.0
ANALYTICAL SOLUTIONS OF TOT
This section provides a more detailed discussion of the two analytical
approaches for calculating TOT presented in Section 3.0. The data
requirements and the sources of the data are presented, the methods are
explained, and example calculations are provided. It will be shown that
these methods provide a means of calculating TOT that is easier, and less time
consuming than using unsaturated flow models. However, the application of
these methods is limited due to the assumptions used in their development.
4.1 DATA REQUIREMENTS AND SOURCES
The data required by the analtyic.al solutions for calculating TOT are
listed below:
stratigraphy of the site;
thickness of geol
-------
conductivity (K) and matric potential. These charteristics of the soil(s)
are required for any determination of flow or TOT in the unsaturated zone:
either analytically or through the use of models.
Ideally, these relationships should be measured in the laboratory using
soil samples obtained from the site. If laboratory measurements are not
possible, the following simple analytical relationships between pressure head
and water content, and between conductivity and matric potential (Campbell,
1974) can be used:
K *
where
^g = air entry matric potential;
e$ = saturated water content;
e = field water content;
sat - saturated h-Jraulic conductivity;
b = negative one times the slope of the log-log plot of t versus e; and
n = 2 + 3/b.
Using the above relationships it is necessary to know only the slope of
the log-log plot of 4>m versus e, the saturated hydraulic conductivity, and
the saturated moisture content. The saturated hydraulic conductivity can be
determined in the field or measured in the laboratory. Field methods
are preferred to laboratory methods, and are detailed in Appendix A, Section
1.0. Default values, such as those listed in Table A.1.1, should be used only
as screening factors in choosing a proper field method.since they may
underestimate hydraulic conductivity by several orders of magnitude.
The saturated moisture content (e$) can also be obtained from laboratory
measurements. If measurements are not possible, e can be assumed to be
4.2
-------
Table 4.1-1. Representative Values for Saturated H>draulic Conductivity
(Source: Freeze and Cherry, 1979)
1
ROCKS
Unconsolifloted k k K K K
deposits (dorcy) (crn2, (Cm/s)(m/s) (QQi/aav/f,2)
if
l{_
5
- 1
3 1
~ £ Oj?
£!ss
> (
f 0
i
1
DO
. a!."
" - c
T> O 0
o p fl
O)
^
^ 2 r *E *
^\ !
1 ^ J
gf
3
t
(
c
ft
o
u
T
1
C
e
c
?c
3
*
h
^
'
J
0
t.
3
5
ft
3
;
3
-g
o<3
c
6
E
-10s rtO'3 rl02
to4
IO3
to2
10
1
10"
-1C'2
-10"3
-io-4
- IO"5
-to'6
-to-7
-to'8
-ID'4
-io-s
-io-6
-to-7
-tO'8
-ID'9
-10
- 1
-to"
-ID'2
-to'3
-io-4
-tO'10 - 1C'5
-to"1
- to"1'
- to"3
IO"4
.o"5
10"6
-ID'6
-io-7
-ID'8
-ID'9
- 10"°
-IO"1
r- 1
.1C6
-IO"
-!0'2
-to-3
to-4
-1C'5
ID'6
-to'7
-,o-8
1C'9
10"°
10"'
10"2
IO"3
IO5
1 W
104
-1C3
-IO2
- tn
1 W
t
10-'
ir-2
lu
10'3
to'4
10'5
io-6
io-7
4.3
-------
equal to the total (actual) porosity. Representative values of total
porosity are given in Table 4.1-2.
Table 4.1-2. Representative Values for Porosity
Material
Coarse Gravel
Medium Gravel
Fine Gravel
Coarse Sand
Medium Sand
Fine Sand
Silt
Clay
Porosity
28%
32%
34%
39%
39%
43%
46%
42%
Analytical solutions of travel time assume.steady state flow of moisture
through the unsaturated zone. A simple approximation of steady state flux is
to assume that it is equal to the net infiltration at the site. Net
infiltration is equal to the net precipitation minus actual
evapotranspiration. If such information is not included in the Part B
application, it should 'e obtainable from weather stations or agricultural
research stations.
4.2 DESCRIPTION OF ANALYTICAL SOLUTIONS
Two analytical solutions to unsaturated travel time are presented below.
Both solutions require the steady state assumptions described in Section 3.0.
The first solution assumes that the hydraulic gradient is equal to 1. The
second solution allows for variable moisture content in the soil column
(i.e., hydraulic gradient may not be equal to 1).
4.2.1 Solution for Unit Hydraulic Gradient
The following solution assumes steady state flow and a unit hydraulic
gradient and employs the analytical soil moisture, pressure, and conductivity
relationships described earlier (Campbell, 1974). Utilizing Darcy's equation
and the soil characteristic relationships described by Campbell (1974), it is
4.4
-------
possible to derive the following expression for moisture content as a
function of steady state flux (Heller, Gee, and Myers, 1985)
9 =
\at
where
q = steady state flux;
i/ = saturated hydraulic conductivity;
Nsat
QS = saturated moisture content; and
m = l/(2b + 3), where b is negative one times the slope of the log-log
plot of
-------
Example Problem 1
This first problem provides an example of the procedure for calculating
TOT through an unsaturated zone consisting of a single material type. A
schematic of this single layered system is shown in Figure 4.2-1.
The following parameters must be known (either measured in the field or
laboratory, or obtained from the literature) in order to apply this solution.
Example values of each parameter are provided.
_ Parameter _ Example Value
Flux (q) 0.5 cm/yr
Saturated water content (e ) 0.31 m^/m^
for soil profile
The slope (-b) of a log-log -3.162
plot of *m versus e
Saturated- hydraulic 5.4 x 10 cm/yr
conductivity (Ksat) of the soil
Length of unsaturated 3760 cm
column (L)
Step 1; Calculate m v
m
1
57374 + 3
= 0.107
Step 2: Calculate Steady-State Moisture Content e
e =
_ / 0.5 cm/yr ,0.107 ,Q 31 m3
5.4 x 104 cm/yr
4.6
-------
Figure 4.2-1. v-hematic of Example Single-Layered System
4.7
-------
= 0.09 m3/m3
Step 3: Calculate Travel Time (T)
q
= (3760 cm)(0.09 m3m3)
0.5 cm/yr
= 677 years
4.8
-------
Example Problem 2
This second problem is similar to the first except that it illustrates
the procedure for calculating TOT through a multi-layered unsaturated flow
system. A schematic of this multi-layered system is shown in Figure 4.?-?.
The following parameters must be known in order to apply this solution.
Example values of each parameter are provided.
Parameter
Flux (q)
Value
0.5 cm/yr
Notes
Constant throughout section
Saturated Water Content
e
si
»s2
Js3
Slope of log-log plot
of j, versus e(-b)
l
0.31
0.40
0.42
-3.162
-3.475
-3.610
Saturated hydraulic
conductivity (K .)
satj
K]
K;
5.4 x 10 cm/yr
1.0 x 10 cm/yr
0.5 x 104 cm/yr
Length of unsaturated
column (L)
1
1,000 cm
1,000 cm
1,760 cm
4.9
-------
Figure 4.2-2. Schematic of Example Multi-Layered System
4.10
-------
Step 1: Calculate m for Each Layer
Zb + 3
ml = 0.107
m2 = O-101
m3 = 0.098
Step 2: Calculate Steady-State Moisture Content for Each Layer
e = (r^-f ee
Nsat s
ej = 0.09
e2 = 0.15
93 = 0.17
The total travel time through the section in the sum of the times through
each of the layers.
Step 3; Calculate Travel Time (T)
T = (1.000 cm)(0.09) = ,80
'1 0.5 cm/yr A y
T 9 (1,000 cm)(0.15) _ 300
T2 0.5 cm/yr JUU
T = (1,760 cm)(0.17) = 60Q
T3 0.5 cm/yr buu yr
T = Tl * T2 + T3 = 1'080 yr
4.11
-------
4.2.2 Solution for Variable Moisture Content
Solution of the variable moisture content case is more complex and
requires discretization of the soil profile into a number of node or grid
points as shown in Figure 4.2.2-1. The analytical solution for the case is
(Jacobson, Freshley, and Dove, 1985)
^ = 4»M + A2i (q/K* - 1)
where
^ = pressure head at the upper grid point;
^i-1 = Pressure head at the lower grid point;
q = flux through the soil column;
Az- - elevation difference between grid points; and
K* = harmonic mean hydraulic conductivity between grid points
AZ<
K*
i + 6ZM/KM
Ki'Ki-l ~ hydraulic.,onductivity at the upper and lower grid points,
respectively.
The solution begins with the grid point located at the lower boundary
(water table), where i^_j is known to be 0, and K^_j is known from \J>j_j and
the soil characteristic curve. The solution proceeds iteratively by assuming
a value of \b., determining K*, and then solving for ty-. A new value is
assumed for ty. and the process repeated until there is convergence on a
solution. The calculated value of ^ is then used as \\>^_^ for the next pair
of grid points and the process is repeated.
Once the solution has determined the pressure head at every grid point,
the moisture content and hydraulic conductivity at every grid point can be
obtained from soil characteristic curves. Examples of soil characteristic
curves, how to use them, and where to get them is provided in Section 2.0.
Knowing the moisture content and hydraulic conductivity at two grid points,
the travel time between the grid points is given by
4.12
-------
Grid 1 + 1
Grid i
K. = Hydraulic
1 Conductivity
of Grid
Grid i - 1
Grid i - 2
I
V
'i +1, i + 1
^^"
z ^
* » 1
Z1-l.
1 - 1
Zi - 2,
i-2
A
Node or
_ Grid Point
J« ~
If 1
T AZ. A*.
h ^ 1
.-.L'LljL.
Figure 4.2.2-1. Discretization Between Grid Points
4.13
-------
11 = K* dh.
The above equation is used to determine the travel time between every pair of
grid points. These travel time segments are then summed to obtain the total
travel time through the soil column.
It is possible to perform the above solutions manually for very simple
systems. However, as with all iterative solutions, the process can be very
time consuming. Therefore, the use of a computer is recommended. A computer
code to perform the above solutions for pressure head and travel times has
been developed by Jacobson, Freshley, and Dove (1985).
4.14
-------
SECTION 5.0
UNSATURATED FLOW MODELS
5.1 INTRODUCTION
The two types of unsaturated flow models identified in Section 2.0 are
examined in greater detail in this section. It should be noted that the
intent of this appendix is not to recommend specific computer codes for use
by permit writers or to provide detailed instructions in the use of any
particular codes.
The purpose is to demonstrate the use of two codes which are considered
to be representative of codes in these two categories.
A large number of unsaturated flow codes are presently available.
General characteristics of many of these codes are presented in Table 5.1-1.
A partial list of available unsaturated flow models is also contained in
EPA (1984).
Selection of a cod; should be made by the perspective user considering
such factors as:
1) familiarity with the operation of an appropriate code;
2) availability of data required by the code;
3) applicability of the code to the specific problem (e.g., dimensionality,
complexity of the system);
4) acceptability and documentation of the code; and
5) hardware availability.
Familiarity with a code is perhaps the most important consideration. If
an ana-lyst is already familiar with a particular code, that code should be
used provided it meets the requirements of the application. Availability of
data is the second most important consideration. The available data for an
application should be compared with the data requirements of the codes being
considered. Applicability of the code to the problem is perhaps equally as
5.1
-------
Table 5.1-1. Summary of Unsaturated Zone Codes
Spatial Characteristics
in
ro
Code
Han
AMOCO
ALPURS
BETA- II
BRUTSAERT1
8RUTSAERT2
CMC
COOK
OELAAT
FEMHATER
FLUMP
GRANOAIF
GPS IN
MOMOLS
PORES
REEVES-DUGUIO
SHELL
SSC
STGtfT/MOGHT
SLH2
SUPERMOCK
TRACR30
TRIPH
TRUST
TSIE
UNFLOU
UNSAT1
UNSATIO
UNSAT?
VERGE
VS20
UAFE
Dimensions Discretization Method
Ho. 1
37
21
.40
33
39
38
29
48
24
5
30
31
52 X
49
22
42
41
27 X
46
47
SS
51
4
32
7
28 X
10 X
26
8
45
SO
2
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
3 FOM IFDM FEM Other Special Features
X
X Three- phase oil. water.
gas
X
X X
Two-phase oil £ gas
I Roots, evapotranspiration
X
X. V
X
X X
X
X X
X Radioactive decay
X X
X
X
X Roots, evapotranspiration
X Roots
X X 1- or 2 -phase flow with
tracer in either phase
(air or water)
Freundlich. Langmuir
sorption, radioactive
decay, capillary effects
X 3 -member decay chain
X I
X
X
X
X Roots, evapotranspiration
« Roots
X X
X Roots, evapotranspiration
X Coupled heat » 2-phase
Past Applications
Oil reservoir
Oil reservoir
Experimental
Experimental /Laboratory
Oil reservoir
Oil reservoir
Ground-water extraction
crop production
Underground nuclear
explosions
Experimental
Oil reservoir
Oil reservoir
Oil reservoir
Ground-water e>traction
Ground-water extraction
Tracer flow in unsat.
conditions, Radionu-
clude transport. Tracer
flow in fractured
system
Radioactive waste
disposal
Oil reservoir
Radioactive waste
disposal
Crop studies
Engineering design
Radioactive waste
storage
Confined underground
Principal Contact
AMOCO
Mobil Corp.
Intercomp.
Brutsaert
Brutsaert
CMC
Cook
Oe Laat
Yeh
Narasimhan
Morrison
Exxon
Rojstoc/er
UKAEA
Reeves
Shell Oil Co.
SSC
de Smedt
Oe Laat
Reed
Travis
Gureghian
Narasimhan
Tech. Soft. I Eng.
Pickens
van Genuchten
Bond
Neuman
Verge
Lappala
Travis
Comments
Proprietary Code
Proprietary Code
Proprietary Code
Proprietary Code
Proprietary Code
European Code
Proprietary Code
European Code
Proprietary Code
Proprietary Code
European Code
European Code
Can operate in
1. 2 or
3 dimensions
Proprietary Code
Can operate in
mass transport (air vapor
1 liquid) Accurate treat-
ment of H^O Separate
velocity field phase
radioactive waste dis-
posal . In-situ fossil
energy recovery studies,
2-phase flow and tracer
studies
I or
2 dimensions
KEY: FOM - finite difference method.
IFDH integrated finite difference method.
FEM finite element Method.
-------
Table 5.1-1. Cont'd.
en
CO
Tiame
Spatial Characteristics _
Dimensions Discretization Method
1 7 I FDll'TTBfnTM Oflier
BACHMAT
DUGUIO-REEVES
FECTRA
FEMUASTE
ML IRAN
MMT-DPRU
SCAT1D
SCAT20
TRNHDL
6
19
17
?5
64
44
35
36
2
X
X
X
X
X
X
X
X
X
Otner Special features Pas_t_*J)Plii:aJt'.ons_ ££.'1c_iPJ J Contact Comments
KEY:
FOH = finite difference method.
IFPM = inteqrated finite difference method.
FEM = finite element method.
Surface/ground water
Absorption & decay
li\J order decay, sorption
Ground-water studies
Stochastic velocity field
Stochastic velocity field
Backmat
Ouguid
Baca
Yeh
Reisenauer
Simmons
Oster
Oster
Av-Ron
Middle-east Code
Compatible wi th
Code No. 22
Compatible with
Code No. 22
Compatible with
Code No. 44
Discrete Parcel
Random Walk
Stochastic Code
Stochastic Code
Middle-east Code
-------
Table 5.1-1. Cont'd.
Spatial Characteristics
in
Code
Name
HANKS
MARINO
MCCANN
MOBIOIC
NMOOEL
SEGOL
SHAMTU
SUMATRA- 1
TARGET
TRANS
TRANSONE
TRANSTUO
UNFLW
MAT SOL
UHC
KEY: FOM -
1FON
FEM
Dimensions
No. I 2 3
1 X
S3 X
20 X
14 X
X
9 X
43 X
11 X
18 X
23 X
12 X
13 X
3 X
16 X
34 X
Discretization Method
FDM IFDM FEN Other
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
finite difference method.
integrated finite difference method.
finite element method.
Special Features
Roots
Ist-order reactions
Heat transfer
0- & 1st-order decay
Heat transfer, elegant
numerical solution. 0- &
Ist-order decay. Vari-
able saturation. Radio-
active decay products
Past Applications
Crop production studies
Water loss by
evaporation
Tailings and chemical
waste disposal, radio-
active waste disposal
Salinity studies
Principal Contact
Hanks
Marino
McCann
Couchat
Selim
Segol
Vauclin
van Genuchten
Dames 4 Moore
Walker
van Genuchten
Shapiro
Kapuler
Gaudet
Crooks
Comments
European Code
European Code
Proprietary Code
European Code
Integrated Com-
partment Method
-------
TABLE 5.1.1. Cont'd.
Code
Spatial Characteristics
Dimensions Discretization Method
Name
No. 1 2 3
SESOIL
(Seasonal
Soil Model)
PRZM
(Pesticide
Root Zone
Model)
FDM IFDM FEM Other Special Features
Single constituent
migration through
unsaturated zone;
user-friend ly
Cal^~lates soil
moisture charac-
teristics, crop
root growth,
pesticide
application and
soil transport
Past
Applications
Hydrologic,
Sediment and
pollutant fate
simulation
Pesticide
migration
through
unsaturated
root zone
Principal
Contact
Bona Zountas, M.
A.D. Little
(617) 864-5770
Carsel, R.F.
U.S. EPA
Environmental
Research Lab
Athens, GA
EPA Pub. No.
600/3-84-109
Comments
Analytical
Model
Numerical
Model
-------
important. Applicability involves consideration of such factors as
dimensionality (e.g., one-dimensional flow versus two-dimensional flow) and
complexity (e.g., number of layers, degree of inhomogeneity). Acceptability
may also be an important consideration, particularly within a regulatory
framework. In all cases, an effort should be made to select codes which have
been fully documented and verified against standard solutions. Lastly,
hardware requirements may be important. Codes which require large computer
systems are inappropriate if such systems are not available.
The reader interested in evaluating and selecting an unsaturated flow
code for a particular application is referred to the review of unsaturated
codes prepared by Oster (1982).
Once a code has been selected, detailed instructions on operation and
use should be obtained from the user's manual for the code. Availability of
code documentation is summarized in Table 5.1-2. References to user's
manuals for unsaturated flow codes are provided in Oster (1982).
Models are developed and applied to understand and predict the behavior
of complex physical systems and processes. Physical systems, such as the
unsaturated zone, display characteristic behavior in response to physical
laws. This behavior is described (either exactly or approximately) in terms
of mathematical expressions (e.g., differential equation describing flow).
Many of these expresses are not amenable to analytical solution and are
transformed into approximate solutions in the form of computer codes. These
codes form the framework upon which models are developed. Models are
developed through assignment of representative data to the computer code.
These data are assigned to represent and describe the physical properties of
the system being modeled (e.g., spatial distribution of hydraulic
conductivities).
A model, therefore, consists of two components, the computer code and
the input data for the code. The accuracy of a model is dependent on both of
these components. Codes must adequately describe the processes of importance
for the particular system being modeled. Input data must be provided which
are representative of the properties of the system.
Model results are nothing more than solutions to complex mathematical
expressions. The mere ability of a model to produce results says nothing
5.5
-------
Table 5.1-2. Summary of Unsaturated Code Documentation and Availability
(Source: Oster, 1982)
Documentation
No.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
IS
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
4B
49
50
51
52
S3
54
55
Code
Name
HANKS
TRNHDL
UNFLU
TRUST
FLUW
BACHMAT
UNFLOH
VERGE
SEGOL
UNSAT1D
SUMATRA- i
TRANSONE
TRANS TWO
MOB1DIC
NHOOEL
WATSOl
FECTRA
TARGET
OUGUIO-REEVES
HCCANN
ALPURS
REEVES-OUGUID
TRANS
FEMVATER
FENMSTE
UNSAT2
STGMT/MOGyT
UNSAT1
COOK
GANOALF
GPSIM
TSAE
BRUTSAERT1
UHC
SCAT10
SCAT2D
AMOCO
CMC
BRUTSAERT2
BETA II
ssc
SHELL
SHAMTU
MMT-OPRW
VS20
SUM- 2
SUPERMOCK
OELAAT
PORES
UAFE
TRIPM
NOMOLS
MARINO
MLTRAN
TRACR30
Model
Description
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
V
X
X
X
X
X
X
X
X
X
X
X
X
X
Published
Applications
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
User's
Manual
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
X
Applicat ions
User's Laboratory Field ___^
Data Data Available Proprieta~
Code Availabi1ity_
5.6
-------
947 2 00-2A
with respect to the accuracy or validity of those results. The accuracy of
model results and, hence, the accuracy of the model itself are typically
assessed through the process of calibration. During calibration, the model
input data are adjusted until the model accurately predicts conditions which
are known to exist. Calibration itself, however, is still not proof of model
accuracy. Accuracy can be further tested through the process of
verification. During verification, the calibrated model is used to predict
known conditions at a different time than that tested in the calibration
process.
The processes of calibration and verification are often time-consuming
and expensive. Therefore, many times they are not performed. Lack of
calibration does not necessarily mean that the model is inaccurate, only that
the accuracy of the model has not been established.
Because of the large number of unsaturated flow codes currently
available, and for simplicity of discussion, only one code representative of
each type of model will be examined: UNSAT1D (UNSATurated ^Dimensional), an
example of a numerical code; and HELP Ojydrologic Evaluation of Landfill
Performance), an example of a water balance code. Characteristics of the
codes, as well as the required input data, input data sources, and output
data, will be discussed to provide permit writers and applicants with a
description of the appvjach to unsaturated flow modeling. The limitations
and applicability of each code for the determination of TOT will also be
discussed.
5.2 EXAMPLE NUMERICAL CODE - UNSAT1D
5.2.1 General Characteristics
The UNSAT1D code was originally developed to describe water movement
under typical agricultural conditions (Gupta, et al, 1978). The code and
its auxilliary programs were later revised and incorporated into an
unsaturated flow modeling system (Bond, Cole, and Gutknecht, 1982).
UNSAT1D is a one-dimensional, finite difference code which solves the
differential equation for ground-water flow under saturated and unsaturated
conditions. It simulates infiltration, vertical seepage, and plant root
uptake as a function of the hydraulic properties of a soil, soil layering,
5.7
-------
root growth characteristics, evapotranspiration rates, and frequency, rate,
and amount of precipitation and/or irrigation. UNSAT1D can be used to
estimate ground-water recharge, irrigation and consumptive use of water,
irrigation return flows, and other processes associated with unsaturated and
saturated soils which can be represented as one-dimensional (Bond, Freshley,
and Gee, 1982).
The UNSAT1D modeling system consists of one computer code which solves
the flow equation and several supporting codes which are used for the
preparation of input data and for evaluation and display of model results.
Use of UNSAT1D requires site-specific input data and an understanding of
ground-water flow theory. Input data requirements depend on the problem
being solved, but in most cases include the soil profile description of the
site; the hydraulic properites of each layer of the profile; characteristics
of vegetation at the site; the means by which water is applied to the site;
and climatic data for the site. The model output includes soil-water
potential (suction), moisture content, and water flux at each node (depth
increment) for each time interval considered by the model.
Table 5.2-1 provides a summary of . the important characteristics and
capabilities of UNSAT1D.
TABLE 5.2-1. Important Characteristics and Capabilities of UNSAT1D
1) Simulates partially-saturated ground-water flow.
2) Simulates infiltration, vertical seepage, and plant root uptake.
3) Derives solution using a finite difference, fully implicit method.
4) Describes one-dimensional flow in a vertical or horizontal direction.
5) Accomodates homogeneous, heterogeneous, or layered soil profiles.
6) Simulates up to ten soil layers.
7) Simulates rain, sprinkler or flood irrigation, or constant head
condition for the upper boundary.
8) Simulates lower boundary conditions as water table, dynamic, quasi-
dynamic, or unit gradient.
5.8
-------
OSWER POLICY DIRECTIVE NO.
5.2.2 General Approach to Application Q 4 T Q 0 0 " 2A ^
The first step in application of a numerical model is development of a
conceptual model of the site being considered. A conceptual model must
identify the important features and characteristics of the site and describe,
in a qualitative way, the relationships between the various components and
processes. The conceptual model of a site must be developed before an
analyst can develop a mathematical approximation of the system. For
unsaturated flow, important data required to develop a conceptual model
include stratigraphic data describing soil layering at the site and climatic
data describing net precipitation and evaporation. It should be apparent
that substantial knowledge of geohydrology is required to develop a
conceptual model.
Once a conceptual model has been developed, the analyst must translate
the conceptual model to a mathematical model by supplying appropriate input
data to the computer code. The first step in developing a numerical model is
development of a finite difference or finite element grid network. UNSAT1D
is a one dimensional finite difference code, so grid development involves
specifying a vertical or horizontal array of nodes (depending on the
application).
Once the grid has teen established, input data must be supplied to the
model. Data describing soil properties and characteristics must be supplied
at each node in the grid. Input data may be obtained in the field, measured
in the laboratory, or obtained theoretically. The necessary input data for
UNSAT1D and their sources or methods of estimation are summarized in
Table 5.2-2.
As indicated in Table 5.2-2 a considerable amount of data are required
to operate the UNSAT1D Model Sequence. One of the difficulties of
unsaturated flow modeling is that data requirements often exceed the amount
of available measured data. Various theoretical and .laboratory techniques
may be used to estimate some of these data, and some data may be generated or
estimated using the supporting programs contained within the UNSAT1D Model
Sequence. These programs must be run before UNSAT1D can be used to simulate
a particular unsaturated flow problem.
5.9
-------
Table 5.2-2. Surmrary of UNSAT1D Input Data and Sources
Input Parameter
Depth of soil layers and lower
boundary condition
Soil hydraulic properties
1) soil-water retention relation-
ship, saturated volumetric mois-
ture content (es), and saturated
hydraulic conductivity.
2) hydraulic conductivity vs.
water content.
3) initial moisture content
4) field density
Precipitation and irrigation
with hourly distribution
Potential evapotranspiration with
diurnal variation
Plant growth behavioV
1) leaf-area index
2) root growth and density
3) growing season
Data Source or Estimation3
Must be known or measured via field
dri1 ling
1) laboratory measurements of moisture
contents at various suction heads.
fls may be assumed to be porosity.
2) calculated from soil-water reten-
tion relationship and saturated
hydraulic conductivity.
3) measured from samples or estimated
from water balance history.
4) measured from samples.
Obtain from nearest weather station
and agricultural sources.
Obtain from weather/experimental
station or calculate with detailed
climatic data.
1) published for some plants.
2) published for some plants. May
be assumed over growing season.
3) available from weather service or
agricultural organizations
aThe degree of estimation acceptable depends on accuracy required .in model.
5.10
-------
There are four data preparation programs within the model sequence.
Brief descriptions of these programs and their functions are given below.
HYDRAK - This program estimates the hydraulic conductivity versus water
content or matric potential relationships for each soil.
EXTEND - This program extrapolates additional data points from the high
suction head/low water content end of the soil moisture characteristic curve.
POLYFIT - This program enters the soil moisture characteristic data into the
UNSAT1D code in the form of polynomial expressions, the preferred form for
these data.
FAOPET - This program estimates daily potential evapotranspiration (PET) for
the site, as required by the model when these data are not specifically
available.
Output from UNSAT1D includes the soil water potential (suction),
moisture content, and soil water flux rates at each node for'each time step.
Examples of graphical output from UNSAT1D showing moisture content versus
depth and cumulative drainage or flux versus time past certain elevations
within an unsaturated soil profile are shown in Figures 5.2-la and 5.2-lb.
5.2.3 Determination of TOT from Model Results
The output of UNSAT1D does not include TOT for the unsaturated zone.
Therefore, the model results must be analyzed using one of the techniques
described in Section 372. to determine TOT. For transient simulations, TOT
can be estimated by following the migration of the wetting front. For steady
state simulations, the steady state model solutions for nodal values of
matric potential and moisture content can be used to calculate velocities and
travel times across each model element.
5.2.4 Limitations
UNSAT1D was developed to predict the amount and rate of water entering
and moving through a partially saturated flow system. The code accomplishes
this task by simulating one-dimensional flow through the system using the
differential equation
C(e) 4* - ' (K (e) I?} - S
5.11
-------
Q.
Ol
a
10 -
Volumetric Moisture Content
Figure 5.2-la. Moisture Profile for a Three Layer
Unsaturated Flow System
o>
o
OJ
Figure 5.2-lb.
Time (days)
Cumulative Water Flux Versus Time Past Several
Elevations in an Unsaturated Fl'ow System
5.12
-------
where
C(e) = soil water differential capacity;
e = volumetric water content;
4» = pressure head;
t = time;
z = vertical coordinate;
K(e) = hydraulic conductivity;
H = hydraulic head; and
S = source/sink term.
The above equation is a general solution for one-dimensional flow in the
unsaturated zone. The only limitations to this solution are that it is a
one-dimensional solution and that it does not account for migration of water
in the vapor phase.
The use of a one-dimensional solution may or may not pose limitations,
depending on the characteristics of the site. If materials present at the
site are highly inhomogenous, a two-dimensional model may be more appropriate
(assuming, of coursev that there are adequate data to define the
inhomogeneity). A two-dimensional model may also be more appropriate for
cases where there may be significant lateral flow. For example, leaks from a
surface impoundment located above dry soil would be expected to undergo
significant lateral migration due to capillarity.
Numerical methods of solution require finer resolution of data than
other methods. This increases the amount of data that must be obtained to
set up a model. In addition, the numerical solution of the nonlinear
unsaturated flow problem is quite sensitive to input data. Values of input
data must be reasonably close to actual values in order to obtain a solution.
Therefore, numerical models are difficult to use with "default" data values
since such values may not yield a solution.
The numerical method of solution also requires a great deal of
computational time. The methods are most appropriate for large mini-
5.13
-------
computers or main-frame computers. Even with such computers, solutions may
require several hours of central processing unit (CPU) time.
Lastly, numerical models are very complex and require a good deal of
understanding on the part of the analyst. The analyst must have significant
experience and understanding in order to develop the conceptual model of a
site and set up the finite difference or finite element grid. In addition,
several runs of the model are often required before input data sets are
adequately adjusted to produce a solution. The analyst must have enough
understanding of the workings of the model to be able to adjust and calibrate
the model.
With the exception of the limitation due to one-dimensional flow, the
above limitations apply to all numerical models.
5.3 EXAMPLE WATER BALANCE CODE - HELP
The HELP code was developed and adapted from the EPA's Hydrologic
Simulation Model for Estimating Percolation at Solid Waste Disposal Sites
(HSSWDS) and from the U. S. Department of Agriculture's Chemical Runoff and
Erosion from Agricultural Management Systems (CREAMS) code (Walski et al.,
1983). HELP is a quasi-two-dimensional hydrologic model which rapidly and
economically estimates the amount of runoff, drainage, and leachate that may
be expected to result from operation of landfills. HELP' performs a
sequential daily analyses of water inflow and outflow that takes into account
the effects of runoff, evapotranspiration, percolation, and lateral drainage
on the water balance for a particular site. The code was not developed to
account for lateral inflow and surface runon.
HELP produces daily, monthly, and annual water budgets which describe
both vertical flow through a landfill profile and horizontal flow through its
drainage layers. The model requires climatological data and soil and
landfill design data. Site-specific data should be used for the analysis;
however, if these data are not available a substantial amount of climatic and
soil data are maintained within the model, as well as default options for
vegetative covers. The model's output includes summary data describing the
water balance for each layer in the model. These data can be provided on a
daily, monthly, or yearly basis depending on the needs of the user.
5.14
-------
Table 5.3-1 provides a summary of the important characteristics and
capabilities of HELP.
5.3.1 General Approach to Application
As with numerical models, aoplication of water balance models requires
formulation of a conceptual model of the site and representation of this
conceptual model within the framework of the mathematical model. The HELP
model represents the system as a series of layers. Four types of layers are
allowed:
1) those which only allow vertical percolation;
?.) those which inhibit vertical percolation (barrier layers);
3) those which allow lateral drainage; and
4) waste layers.
Application of the model requires the analyst to review the stratigraphy of
the site and translate the stratigraphy into a series of layers of the
appropriate type.
Table 5.3-1. Summary of Characteristics and Capabilities of HELP
1) Simulates partially-saturated ground-water flow.
2) Performs a sequential daily analysis to determine runoff,
evapotranspiratiorjf percolation, and lateral drainage.
3) Uses a quasi-two-dimensional water budget approach.
4) Describes two-dimensional flow for both vertical flow through the
profile and horizontal flow through drainage layers.
5) Applies to a wide variety of landfill designs.
6) Simulates up to nine layers.
7) Maintains default climatic and soil data, as well as default options for
site vegetation.
8) Assumes gravitational forces to be the most important force for fluid
movement (capillary Forces are ignored) thereby greatly simplifying the
unsaturated flow solution.
5.15
-------
HELP requires input data similar to those for UNSAT1D, but with less
spatial resolution. One advantage of the less stringent data requirements
and less complex numerics of the water balance models is that it makes it
easier to use default values. The HELP model maintains an internal data base
of default data values. This data base includes five years of climatic data
for 102 cities in the United States and default characteristics for 21 soil
types. The model also makes available seven default options for site
vegetation. The necessary input data and their sources and methods of
estimation are summarized in Table 5.3.1-1.
HELP model output consists of a summary of all default or user-provided
input information (except daily precipitation) used for the simulation and a
summary of the analysis computed by the model. The analysis summary includes
a table of annual totals for each year of simulation; a table of average
monthly totals for all years simulated by the model; a table of average
annual totals for all years of simulation; and a table of peak daily values
for all years of simulation. A summary of the information contained in each
table is provided in Table 5.3.1-2. If the user is interested in monthly
output, HELP can produce tables which report monthly totals for all years of
simulation. These tables include:
precipitation;
runoff;
evapotranspiration;
percolation from base of landfill cover;
percolation from base of landfill;
t lateral drainage from base of landfill cover; and
lateral drainage from base of landfill.
If daily output is desired, the model provides daily values for each
Julian date of each year of simulation. In addition to the above data, these
tables include:
head at base of landfill cover;
head at base of landfill; and
soil moisture content of the evaporative zone.
5.16
-------
Table 5.3.1-1. Summary of Input Data and Sources for HELP
Climatologic Data
Daily precipitation values
All can be obtained from data base or measured for each year of
interest (2 - 20 year) -- libraries, universities, agricultural and
climatologic research facilities, and the National Climatic Center
are possible sources.
Mean monthly temperature
Measured for each year or single set of data for all years --
libraries, universities, agricultural and climatologic research
facilities, and the National Climatic Center are possible sources.
Mean monthly solar radiation factors
Measured for each year or single set of data for all years --
agricultural publications, solar heating hand books, and general
reference works are possible -sources.
Winter cover factors
Measured for each year or single set of data for all years --
libraries, universities, agricultural and climatologic research
facilities, and the National Climatic Center are possible sources.
Leaf area indices (LAI)
Measured for e^ch year or single set of data for all years
various refere,Xes, including USDA's publication, "Climate and Man,
Year Book of Agriculture," are possible sources.
Vegetative Cover Data
Root zones or evaporative zone depth
Choose some of seven vegetative cover options from data base, or
must be known or measured/observed on site.
Design and Soil Data
Landfill profile
Modeled from data base or observed/measured on site.
5.17
-------
Table 5.3.1-1. Cont'd.
Soil data
From data base (21 default soil types' "aVailabie";or
observed/measured on site. Observations/measurements must include:
porosity, field capacity, wilting point, hydraulic conductivity,
and evaporation coeffecients for each soil layer of profile.
Soil compaction
From data base or use soil data representative or compacted soil
(from observation of site).
Design data Number of layers and their descriptions type,
thickness, slope, and maximum lateral distance to a drain (if
applicable).
Observed/measured on site.
Design data Whether or not synthetic membranes used in the landfill
cover and/or liner
Observed/measured on site.
5.18
-------
Table 5.3.1-2. Summary of Output Data from HELP
Analysis Summary Table No. 1
Table of annual totals for each year of operation simulated by the
model:
a) precipitation;
b) runoff;
c) evapotranspiration (total of surface and soil evaporation and plant
transpiration);
d) percolation from base of cover;
e) drainage from base of cover;
f) soil water at beginning of year;
g) soil water at end of year;
h) snow water at beginning of year; and
i) snow water at end of year.
Analysis Summary Table No. 2
Table of average monthly totals for all years of operation simulated by
the model:
a) precipitation;
b) runoff;
c) evapotranspiration;
d) percolation from base of cover; and
e) drainage from base of cover.
Analysis Summary Table fl',* 3
^' '
Table of average annual totals for all years of operation simulated by
the model:
a) Precipitation;
b) runoff;
c) evapotranspiration;
d) percolation from base of cover; and
e) drainage from base of cover.
Analysis Summary Table No. 4
Table of peak daily values for all years of operation simulated by the
model:
a) precipitation;
b) runoff;
c) percolation from base of cover;
d) drainage from base of cover;
e) maximum head on base of cover;
f) snow water;
g) maximum soil moisture for vegetative layer; and
h) minimum soil mositure for vegetative layer.
- 5.19
-------
5.3.2 Determination of TOT from Model Results
As discussed earlier, travel times can be determined as particle travel
timer., travel times associated with instantaneous loadings, or steady state
conditions. Application of the first approach requires that model outputs
include the moisture content within a layer and the flux out of the layer.
Since water balance models assume average conditions throughout a layer, the
instantaneous pore water velocity through a layer can be approximated as the
flux out of the layer divided by the average moisture content of the layer.
Determination of the travel times associated with instantaneous (pulse)
loadings may not be appropriate with water balance models. Because of the
averaging that occurs within a layer, there is a loss of resolution.
Therefore, it is difficult to detect the migration of wetting fronts.
Water balance models can be used to solve for steady state moisture
content for each layer. These steady state conditions can be used to
determine steady state travel times through the landfill soil column. As
discussed in Section 3.2, however, this TOT is the steady state TOT and not
the TOT of the first particle of water leaving the site.
Water balance models are not suited' for contaminant transport problems
so that TOT cannot be determined from contaminant TOT.
5.3.3 Limitations
The HELP code was designed to develop long-term water balance models for
landfills to predict generation of leachate from landfills. The soil profile
and landfill are divided into layers and water is budgeted to each layer
based on a mass balance between water flowing into each layer and water
flowing out of each layer. The code allows for lateral drainage from some
layers, giving it quasi-two-dimensional capabilities.
Unlike numerical models, HELP is not based on a general solution to
unsaturated flow and, therefore, has several limitations with respect to its
analytical capabilities. For example, HELP simplifies unsaturated flow by
assuming that gravity is the driving force for all fluid movement; capillary
forces and the effects of vegetation are ignored by the model. Therefore,
the solution obtained by HELP is no more rigorous than those obtained by
analytical methods.
5.20
-------
The water balance approach also lacks the resolution possible with
numerical models. In the water balance solution, soil properties are
assigned by layer and the solution yields the average moisture content for
the entire layer. Using such an approach, subtle effects such as the
migration of a wetting front may not be seen.
The HELP code was developed to simulate the migration of leachate from
landfills. Application of the code to include migration through the
unsaturated zone beneath a landfill will require addition of several soil
layers beneath the landfill. There are internal limits to the number of
layers that can be simulated and internal requirements for certain types of
layers. Therefore, it is not possible to simulate a large number of soil
layers beneath the landfill. Because of limited spatial resolution, the code
is probably best applied to sites in humid areas having thin unsaturated
zones. This limitation was confirmed by a recent comparison between HELP and
7
the UNSAT1D numerical code (Thompson and Tyler, 1983).
5.21
-------
SECTION 6.0
EXAMPLES OF TOT DETERMINATION
This section presents examples of determining unsaturated zone TOT for
several proposed hazardous waste facilities. Examples 1 and 2 demonstrate
the use of analytical solutions for determining steady state TOT at proposed
hazardous waste disposal facilities in the Gulf Coastal Plain and Basin and
Range physiographic regions. Example 3 demonstrates the use of the UNSAT1D
numerical model for determining TOT associated with an accidental spill at a
proposed hazardous waste site in the Columbia-Snake River plateau
physiographic region.
6.1 EXAMPLE 1 - Case Study G
Case study G is a land disposal facility located near the northern edge
of the Gulf Coastal Plain. A review of the data provided in the Part B
permit application for this facility is presented in the Case Study Appendix
to the Phase II Location Guidance (see Appendix C).
The case study G facility is underlain by approximately 21 to 34 m of
fine to medium grained quartz sand, with limited occurrences of silt, clay,
and lignite beds. This sand layer beneath the facility forms the water-table
aquifer. Depths to ground water at the facility range from 6 to 15 m. This
shallow aquifer is recharged by rainfall which averages 119 cm/yr at the
facility.
6.1.1 Description of Method and Data
Unsaturated TOT was calculated using the one-dimensional steady state
analytical solution described by Heller, Gee, and Myers (1985). This
solution assumes that the hydraulic gradient in the unsaturated zone is equal
to one. The unit gradient assumption implies that flow in the unsaturated
zone is dominated by gravity (i.e., capillary forces are negligible).
6.1
-------
Because the site is located in a humid region having a moderately high
rainfall (119 cm/yr), soils in the unsaturated zone should be fairly moist
(i.e., at or above field capacity). Therefore, the unit gradient assumption
is probably valid for this site.
As discussed in Section 4.0, this analytical solution requires
relatively few input data compared to other methods. These data are:
soil profile and depth;
soil characteristics (saturated hydraulic conductivity, saturated water
content, moisture content versus pressure head); and
steady state moisture flux.
The assumptions used to develop input data and the limitations resulting from
these assumptions are described below.
Soil Profile and Depth
The analytical method may be applied to single- or multi-layered soil
profiles. Geologic cross sections of the site identify lenses of clay and
silt within the sand. However, these lenses are not continuous over the
site. Therefore, a uniform profile of sand was assumed. In view of the
higher permeability of sand compared to silt and clay, this is a conservative
assumption (i.e., yields lower TOT).
The thickness of the unsaturated zone was estimated from the geologic
cross sections and repoted ground-water surface .elevations. The minimum
distance from the bottom of the facility to groundwater was estimated to be
6 m. This minimum distance was selected for the analysis to yield a worst
case.
Soil Characteristics--
Saturated hydraulic conductivities for soils at the facility were
determined by aquifer tests. The geometric mean conductivity from tests of
shallow wells was 0.079 cm/sec. For lack of other data, this value was used
for saturated conductivity in TOT calculations.
It should be noted that laboratory permeameter results are the preferred
source of saturated conductivity data. Field measurements from aquifer tests
are easier to obtain, however, and are expected to be the major source of
such data presented in Part B applications. The following limitations to the
use of these data should be recognized:
6.2
-------
t Aquifer test results indicate the hydraulic conductivity of material in
the saturated zone. Unsaturated TOT calculations require the saturated
conductivity of material in the unsaturated zone. In this example,
materials in the two zones are very similar and the use of saturated
conductivities is not expected to be a major source of error.
Aquifer test results represent horizontal saturated conductivity.
Unsaturated TOT calculations require vertical hydraulic conductivity.
There can be significant differences (e.g., order of magnitude) between
vertical and horizontal conductivity for some materials. The
relationship between vertical and horizontal conductivity for the
materials at the site is not known. Use of aquifer test results should
be recognized as a potential source of error.
No saturated moisture content data were presented in the Part B
application. As described in Section 4.0, total porosity is a good
approximation to saturated moisture content. Because no porosity data were
presented in the permit application, the default values presented in Section
4.0 were used. A value of 0.41 was used to represent the average saturated
moisture content based on the default porosities for fine sand and medium
sand.
No data describing the moisture retention characteristics of soils at
the site were provided _n the Part 3 permit application. Typical values of
the slope of the moisture retention curves ("b" values) for different soil
textures are presented by Hall et al. (1977). These values are shown in
Table 6.1.1-1. A value of 4.0 was selected as representative of the sandy
soil at the site.
Moisture Flux
No information was presented in the Part 3 application describing
moisture flow through the unsaturated zone. The yearly average rainfall for
the site was reported to be 119 cm. A conservative assumption would be to
ignore runoff, evaporation, and transpiration and assume that all
precipitation is available for recharge. This assumption would tend to
maximize the unsaturated TOT.
6.3
-------
Table 5.1.1-1. Typical Values for Slope of Soil Moisture
Retention Curve (b) Source: Hall- et al., 1977
Soil Texture
Clay 11.7
Silty Clay 9.9
Silty Clay Loam 7.5
Clay Loam 8.5
Sandy Clay Loam 7.5
Sandy Silt Loam 5.4
Silt,Loam 4.8
Sandy Loam 6.3
Loamy Sand 5.6
Sand 4.0
Summary
The following summarizes the input parameters for the analytical
solution:
depth to ground water, L = 6.0 m;
saturated hydraulic conductivity, K . = 0.079 cm/sec;
saturated moisture content, 9f_. = 0.41;
sac
negative one times the slope of log-log characteristic curve, b = 4.0;
and
moisture flux, q =H19 cm/yr.
6.1.2 Solution of TOT
The following solution follows the same steps as those presented in
Section 4.0.
Step 1: Calculate m
°'091
J 2)(4.0) + 3
Step 2: Calculate Moisture Content e
T-- e
sat sat
6.4
-------
0.091
(119 cm/yr)
(U.O/y cm/sec M-Jit Wo,UUU sec/yr) (0.41)
= 0.16
Step 3: Calculate Travel Time (T)
(6.0 m) (0.16) _ n pi u_ _ onn J,..-
(119 cm/yr)(0.01 m/cm) ' °'81 yr ' 29° days
6.2 EXAMPLE 2 - Case Study D
This Case Study D Facility is a landfill located in the Amargosa Desert,
in the Basin and Range physiographic province. A review of the data provided
in the Part B permit application for this facility is presented in the Case
Study Appendix and the Phase II Location Guidance Manual (see Appendix C).
This review .constitutes the only source of site-specific data for the
unsaturated TOT calculation presented below.
The site is underlain by at least 170 m of alluvial and valley-fill
deposits, primarily sands, gravels, and cobbles of local origin. These
alluvial and valley-fill materials form the water-table aquifer at the site.
The depth to ground-water at the site is approximately 90 m.
The climate at the site is characterized by very low rainfall and high
evaporation. The average rainfall at the site is 11.4 cm/yr, with
evaporation and potential evapotranspiration estimated at 254 cm/yr and
91 cm/yr, respectively.
Analytical methods of unsaturated travel time are based on steady state
flow through the unsaturated zone. If the steady state flux is not known, it
can be estimated as the net recharge at the site. A value of net recharge
for the site of 0.064 cm/yr is reported in Appendix C. Because of this very
low flux, the assumption of a unit hydraulic gradient may not be valid.
Therefore, unsaturated travel time was calculated using both analytical
solutions presented in Section 4.0
6.5
-------
6.2.1 Unit Gradient Analytical Solution
The data requirements for this analytical method were described in
Example 1. Soil characteristic data for the site (Appendix C) were used to
construct soil characteristic curves for a typical soil at the site. A plot
of suction head versus moisture content is shown in Figure 6.2-1. The slope
of the linear portion of this curve was measured to obtain a "b" value of
3.3. The saturated moisture content and saturated hydraulic conductivity of
this soil are 0.40 and 265 cm/day, respectively. From the cross-section data
for the site, the average depth from the bottom of the landfill to the water
table is 76 m. The above data were used to calculate the travel time for the
steady state flux of 0.064 cm/yr.
Step 1: Calculate m
m = -- = = °'10
Step 2: Calculate Moisture Content 6
9 = (^)m9sat
- 0.064 cm/:>r °-10 (0.40)
(265 cm/day) (365 cm/yr)
= 0.10
Step 3: Calculate Travel Time (T)
(76 m) (0.10)
(0.064 cm/yr)(0.01 m/cm)
= 12,000 yrs
6.6
-------
100-
Ld
Z
o
c5
oo
10-
1-
0.1-
Figure 6.2-1.
0.1 1
THETA
Plot of Suction Head Versus Moisture
Content for Case Study D
6.7
-------
The above equation can also be used to solve for a travel time for the
first 100 ft below the facility. This 100 ft travel time is 4,800 years/
which is well above the 100-yr location guidance criterion.
It should be noted that this analytical solution is not strictly
applicable at the site. The steady state moisture content is so low that it
does not fall within the linear, central portion of the curve where the
solution technique is applicable. At this low moisture content (moisture
contents measured at the site were all at or below the wilting point),
capillary forces would be significant and the unit gradient assumption is not
appropriate.
6.2.2 Iterative Analytical Solution
Because of the low moisture contents at the site, the iterative
analytical solution described by Jacobson, Freshley, and Dove (1985) is
probably more appropriate. This solution allows for variable moisture
contents within the soil profile. The data requirement for this method, in
addition to the steady state flux, are soil characteristics curves for
moisture content and hydraulic conductivity. The curve for suction head
versus moisture content was shown previously in Figure 6.2.1-1. A plot of
hydraulic conductivity versus suction head for a typical soil at the site is
shown in Figure 6.2.2-1. The data from these curves were used to generate
tables of suction head cV!d moisture content, and suction head and hydraulic
conductivity for use in the solution.
The above data were used with the iterative solution to calculate travel
time for a steady state flux of 0.064 cm/yr. To employ the iterative
solution, a grid system was constructed to represent the site. The grid
system was constructed to represent the site. The grid consisted, of 251
nodes, uniformly spaced at 1 ft, for a total depth of 250 ft. A boundary
condition of 0 suction head (saturation) was set for the bottom node to
represent the water table. The iterative solution was then applied to solve
for the steady state moisture profile in the soil column. This profile is
shown in Figure 6.2.2-2. Knowing the steady state moisture content and
suction head at each node, the travel time between each pair of nodes was
calculated, and these nodal travel times summed to give the total travel time
6.8
-------
0.01J
0.001J
Q
g O.OOOId
0.00001J
O.OOOOOlHr-rrrr
III) I I I I Illlj F I rriTIT) TTTTTTTTl
0.1 1 10 100 1000
SUCTION HEAD, FT.
Figure 6.2.2-1.
Plot of Hydraulic Conductivity Versus
Suction Head for Case Study D
6.9
-------
250-i
200-
tf
*
O
§
UJ
50-
0.1 0.2 0.3 0.4
MOISTURE CONTENT
Figure 6.2.2-2.
Steady State Moisture Profile from Iterative
Analytical Solution for Case Study D
6.10
-------
through the soil column. The steady state travel time was 16,000 yrs. The
travel times for the first 101 nodes were summed to give the 100 ft travel
time. This travel time is 5,100 yrs, which is much greater than the location
guidance criterion of 100 yrs.
It should be noted that the travel times obtained from this iterative
solution are greater than those obtained from the unit gradient assumption.
The reason for this is the nonuniform distribution of moisture contents in
the soil column. Because of this moisture distribution, the suction head is
not uniform through the soil column and the hydraulic gradient is less than
1, giving longer travel times. Use of the unit hydraulic gradient solution,
therefore, yields a conservative answer.
6.3 EXAMPLE 3 - Pulse Loading to the Unsaturated Zone
The following example uses an actual case study to illustrate how the
UNSAT1D code can be used to estimate TOT through the unsaturated zone
associated with a pulse loading (i.e., transient flow). The example selected
was performed for the Washington Department of Ecology (1979) to determine
the environmental impacts of establishing a hazardous waste disposal facility
at a site in south central Washington.
The study addressed the potential adverse impacts that might occur in a
number of different areac, to include: earth, air, water, flora,'fauna, and
elements of the human Environment. The portion of the study dealing with
water looked at the potential for migration from a surface spill, through the
unsaturated zone, into the saturated zone, and eventually to a discharge
point; in this case the Columbia River. The purpose of the unsaturated flow
modeling was specifically to estimate the TOT from a hypothetical accidental
liquid spill at the surface vertically downward through the unsaturated zone
to the water table.
6.3.1 Simulation Details
The scenario simulated with the UNSAT1D code was a hypothetical
accidental liquid spill of 424,000 liters spread over an area of 930 square
meters. For a one-dimensional simulation, this volume of spill is equivalent
to an initial ponding of 37.5 cm. Based on percolation rates for the soil at
the proposed facility, the depth of the ponding was linearly reduced over a
6.11
-------
10 day period until it had all infiltrated. After infiltration, a condition
of no surface evaporation was assumed in the model, which gave a conservative
drainage prediction.
One-dimensional vertical flow beneath the spill was assumed. Neglecting
the lateral movement of the infiltrate due to capillary forces further
contributed to a conservative estimate of TOT.
The vertical distance between the spill and the constant water table was
52 m. A homogeneous soil profile was assumed. The soil moisture retention
characteristics (Figure 6.3.1-1) and the hydraulic conductivity versus water
content relationship (Figure 6.3.1-2) for the soil were obtained from
laboratory measurements of soil samples taken from a well constructed near
the proposed site. The saturated hydraulic conductivity equaled 1.7 x 10
cm/sec throughout the entire soil profile. For use in the model, the soil
properties (Figures 6.3.1-1 and 6.3.1-2) were fit with logarithmic
polynomials.
The vertical column was defined with 53 nodes in the UNSAT1D model. The
node spacing was uniform at 1 m. The initial pressure conditions (pressure
head) in the model were set to equilibrium (i.e., pressure head equal to
negative one times the elevation above the water table), as shown in
Figure 6.3.1-3.
During the time th; t the spill was infiltrating, changes in pressure
head near the surface were rapid and, therefore, a small time step of 0.02
hours was used. This time step was used for the first 20 days of the
simulation. After 20 days the time step was doubled after every iteration
until a maximum time step of 24 hours was reached. The simulation was run
fo- a total time period of 300 years.
6.3.2 Model Predictions
The UNSAT1D model results, in terms of pressure head versus depth in the
soil profile at various points in time, are shown in Figure 6.3.1-3. This
figure illustrates the advance of the wetting front at 0 and 10 days, and at
10, 100, and 300 years.
Figure 6.3.2-1 illustrates the advance of the wetting front with time.
This figure shows the time required for 5 cm of leachate to seep past a given
6.12
-------
CO
100,000
!
o
o
l-
3
t/>
t-
o.
10,000
1,000
100
10
0.00
0.10
0.20
Moisture Content
0.30
Figure 6.3.1-1. Soil Water Characteristic Curve for Soil at
Proposed Hanford Hazardous Waste Site
(Source: Washington Department of Ecology, 1979)
J
j
J
j
0.40
-------
U
3
o
c
o
CJ
10
_-,
10
-1
10
-2
10
-4
10
-5
10
-7
10
-8
0.00
Figure 6.3.1-2.
0.20
Moisture Content
0.30
... J
0.40
Hydraulic Conductivity Versus Water Content Curve for
Soil at Proposed Hanford Hazardous Uaste Site
(Source: Washington Department of Ecology, 1979)
-------
E
u
0.
O)
o
j
6,000
0
' 0 Days
r-t-r- T T -i r-
1,000
2,000 3,000 4,000
Suction Head (cm)
.-.,. i.. ,
5,000
Figure 6.3.1-3. Simulated Pressure Head Versus Depth with Time
at Proposed Hanford Hazardous Waste Site
(Source: J^hington Department of Ecology, 1979)
. J
6,000
-------
1,000
2,000
E
o
O.
O
o
3,000
4,000
5,000 -
TT
6,000 I............ ._,
0 10
Figure 6.3.2-1
. j.
20
30
r-r-VI I I l-» I-T-J-T TfTTTn T l-| I r I T ri r I :
40
Time (yr)
50
60
1..1-1 - . . I.
70
Simulated Advance of Uetting Front at
Proposed Hanford Hazardous Waste Site
(Source: Washington Department of Ecology, 1979)
i
i
80
-------
depth (depth of leachate in a one-dimensional case is equal to volume of
leachate in a three-dimensional case). Based on the data shown in
Figure 6.3.2-1 the model predicted that approximately 73 years would be
required for 5 cm of leachate from the spill to reach the water table.
The seepage rate of leachate into the water table (depth of leachate per
10 year interval) is shown in Figure 6.3.2-2. Although the first arrival of
leachate occurs in 40 years, the maximum rate occurs at 100 years. The
? o
maximum leakage rate over the entire 930 m area was approximately 3.5 m /yr.
The time required for percentages of the total leachate from the spill
to arrive at the water table is illustrated in Figure 6.3.2-3. The results
indicate that after 300 years (the total simulation period), about 80% of the
leachate has reached the water table.
6.3.3 Summary
The results of the TOT estimates with the UNSAT1D model show that the
first arrival of infiltrate from the spill was approximately 40 years after
the spill occurred, the maximum seepage rate occurred approximately 100 years
after the spill, and more than 300 years are required for all of the
infiltrate to reach the water table. This case study provides an excellent
example of how a numerical model can be used to estimate the travel time for
a pulse loading.
6.17
-------
4.0
3.0
00
>»
o
e
o
0>
+-)
to
2.0 l-
(O
-tr
o
1.0 I
0
100
200
Time (yr)
Figure 6.3.2-2.
Simulated Rate of Leachate Discharge at
Proposed Hanford Hazardous Waste Site
(Source: Washington Department of Ecology, 1979)
j
300
-------
UD
t.
Q)
o
c
o
CD
(O
x:
u
(0
o>
01
D)
i-
KO
.C
u
o-
>
o
O
90
80
70
60
50
40
30
20
10
V
100
200
.. . . J
300
Time (yr)
Figure 6-3.2-3.
Simulated Cumulative Leachate Discharge at
Proposed Hanford Hazardous Waste Site
(Source: uAingtori Department of Ecology, 1979)
-------
SECTION 7.0
REFERENCES
Bates, R. I., and J. A. Jackson. 1980. Glossary of Geology. American
Geological Institute, Falls Church, VA.
Bond, F. W., C. R. Cole and P. J. Gutknecht. 1982. Unsaturated Groundwater
Flow Model (UNSAT1D) Computer Code Manual. EPRI CS-2434-CCM, Electric Power
Research Institute, Palo Alto, CA.
Bond, F. W., M. D. Freshley and G. W. Gee. 1982. Unsaturated Flow Modeling
on a Retorted Oil Shale Pile. PNL-4284, Pacific Northwest Laboratory,
Richland, WA.
Burdine, N.T. 1953. "Relative Permeability Calculations from Size
Distribution Data." Transactions AIME, Vol. 198, pp. 71-78.
Campbell, G. S. 1974. "A Simple Method for Determining Unsaturated
Conductivity from Moisture Retention Data." Soil Science, Vol. 117, pp. 311-
314.
EPA. 1984. Procedures for Modeling Flow Through Clay Liners to Determine
Required Liner ThicknessEPA/530-SW-84-001, U.S.EnvironmentalProtection
Agency, Office of Solidv-Jaste, Washington, DC.
Feedes, R. A., P. J. Kowalik and H. Zaradny. 1978. Simulation of Field
Water Use and Crop Yield. John Wiley and Sons, New York, NY.
Freeze, R. A., and J. A. Cherry. 1979. Groundwater. Prentice-Hall,
Englewood Cliffs, NJ.
Gupta, S. K., K. Tanju, D. Nielsen, J. Biggar, C. Simmons and J. Maclntyre.
1978. Field Simulation of Soil-Water Movement with Crop-Water Extraction.
Water Science and Engineering Paper No, 4013, Department of Land, Air, and
Water Resources, University of California, Davis, CA.
Hall, -D. G. M., A. J. Reeve, A. J. Thomasson and V. F. Wright. 1977. Water
Retention, Porosity, and Density of Field Soils. Soil Survey Technical
Monograph 9, Rothamsted Experimental Station, Harpenden, England.
Heller, P. R., G. W. Gee and D. A. Myers. 1985. Moisture and Texttiral
Variations in Unsaturated Soils/Sediments Near the Hanford Wye Barricade.
PNL-5377,'Pacific Northwest Laboratory, Richland, WA.
7.1
-------
Jacobsen, E. A., M. D. Freshley and F. H. Dove. 1985. Investigations of
Sensitivity and Uncertainty in Some Hydrologic Models of Yucca Mountain and
Vicinity. DRAFT.PNL-5306,SAND84-7212, PacificNorthwest Laboratory,
Rich land, WA.
Millington, R. J., and J. P. Quirk. 1961. "Permeability of Porous Solids."
Transactions Faraday Society, Vol. 57, pp. 1200-1207.
Mualem, Y. 1976a. "A New Model for Predicting the Hydraulic Conductivity of
Unsaturated Porous Media." Water Resources Research, Vol. 12, pp. 513-522.
Mualem, Y. 1976b. A Catalogue of the Hydraulic Properties of Unsaturated
Soils. Haifa, Israel, Israel Institute of Technology.
Oster, C. A. 1982. Review of Ground-Water Flow and Transport Models in the
Unsaturated Zone. NUREG/CR-2917, PNL-4427, U. S. Nuclear Regulatory
Commission, Washington, DC.
Thompson, F. L., and S. W. Tyler. 1984. Comparison of Two Groundwater Flow
Models - UNSAT1D and HELP. EPRI CS-3695, Electric Power Research Institute,
Palo Alto, CA.
Walski, T. M., J. M. Morgan, A. C. Gibson and P. R. Schroeder. 1983. User
Guide for the Hydroloqic Evaluation of Landfill Performance (HELP) Model,
DRAFT. U. S. Environmental Protection Agency, Office of Reasearch and
Development, Municipal Environmental Research Laboratory, Cincinnati, OH.
Washington Department of Ecology. 1979. Environmental Impact Statement,
Proposed Hazardous I
Ecology, Olympia, WA
Proposed Hazardous Waste Site, Supplement Final. Washington Department of
_^_
7.2
-------
GLOSSARY*
Capillarity -- The action by which a fluid, such as water, is drawn up (or
depressed) in small interstices or tubes as a result of surface tension.
Evapotranspiration -- Loss of water from a land area through transpiration of
plants and evaporation from the soil.
Flux -- The specific discharge of water, volumetric flow divide'd by cross-
sectional area (Freeze and Cherry, 1979).
Heterogeneous Spatial variation of physical properties (Freeze and Cherry,
1979).
Homogeneous -- No spatial variation of physical properties (Freeze and
Cherry, 1979).
Hydraulic-Conductivity Boundary The boundary between two materials having
different hydraulic conductivities (Freeze and Cherry, 1979).
Hydraulic Gradient -- The rate of change of total head per unit of distance
of flow at a given point.
Hydraulic Head -- The sum of the elevation at a point and the pressure head
at the point (Freeze and .Vierry, 1979).
Infiltration The flow of a fluid into a solid substance through pores or
small openings.
Inhomogeniety -- Heterogeniety.
Moisture Content The amount of moisture in a given soil mass, expressed as
weight of water divided by weight of oven-dried soil.
Percolation -- See infiltration.
Pore Water Velocity -- See seepage velocity.
Pressure Head -- The portion of the total hydraulic head due to fluid
pressure (Freeze and Cherry, 1979).
Seepage The act or process involving the slow movement of water or other
fluid through a porous material such as soil.
G.I
-------
Seepage Velocity -- The rate at which seepage water is discharged through a
porous medium per unit area of pore space perpendicular to the direction of
flow.
Steady State Flow where the magnitude and direction of the flow velocity
at any point are constant with time (Freeze and Cherry, 1979).
Transient Flow where the magnitude and direction of the flow velocity at
any point varies with time (Freeze and Cherry, 1979).
Unsaturated Zone -- A subsurface zone containing water under pressure less
than that of the atmosphere, including water head by capillarity.
Vadose Zone -- See unsaturated zone.
Water Table -- The surface of a body of unconfined ground water at which the
pressure is equal to that of the atmosphere.
*A11 definitions are after Bates and Jackson, 1980, unless otherwise noted.
G.2
------- |